Microphysical evolution and column loading drive nonlinear regional contrast in black carbon top-of-atmosphere forcing

Black carbon (BC) aerosols remain among the most uncertain contributors to anthropogenic climate forcing, as their radiative impact depends sensitively on microphysical evolution and atmospheric loading. This study presents a physics-informed, machine learning (ML) approach to estimate clear-sky BC top-of-atmosphere direct radiative forcing (BC TOA) at high spatial-temporal resolution while retaining physical interpretability. The study derives necessary optical properties for radiative transfer modeling (RTM), by constraining them with multi-platform, multi-waveband observations and their associated uncertainties. The RTM outputs are then used to train the ML surrogates and applied over two contrasting urban agglomerates-Xuzhou, China, and Dhaka, Bangladesh. The ML framework closely reproduces physics-based regional climatological mean ($-\mathrm{17.6}\pm \mathrm{2.2}$ W m⁻² versus $-\mathrm{17.4}\pm \mathrm{2.6}$ W m⁻² over Xuzhou; $-\mathrm{14.9}\pm \mathrm{1.1}$ W m⁻² versus $-\mathrm{15.0}\pm \mathrm{1.2}$ W m⁻² for Dhaka), while achieving high predictive fidelity R²>0.95; RMSE ∼ 1.5–1.8 W m⁻² and strong cross-regional consistency (r>0.9). SHAP based predictor attribution indicates that BC TOA estimates are strongly associated with BC aerosol optical depth (BCAOD), column number density, and mixing state, with their relative contributions varying non-linearly across cooling-to-warming regimes. Crucially, similar BC loading can yield contrasting absorption-scattering dynamics across region, which are not captured by simplified forcing parameterization. To test transferability, the combined ML model (trained in Xuzhou, China and Dhaka, Bangladesh) was evaluated zero-shot on two additional regions with contrasting aerosol microphysical conditions represented by Delhi, India (urban and agricultural burning sources) and Mongu, Zambia (strong savanna fires). While transference to Delhi is reasonable (Adj. R²=0.91, RMSE = 2.3 W m⁻²), there is a systematic underestimate at Mongu (Adj. R²=0.83; MBE = −4.2 W m⁻²). Feature-space overlap analysis attributes this degradation to a distributional mismatch in key microphysical predictors. Retraining on an expanded dataset including all four regions preserves urban performance while reducing Mongu RMSE by 68 % and bias from −4.2 to −0.8 W m⁻². Together, the physics-informed ML framework and the multi-domain evaluation provide an efficient and transferable tool for constraining BC radiative impacts across real-world heterogeneity. The study also offers new mechanistic insight into how regional properties reshape BC radiative forcing.

1 Introduction

Among the major climate forcers, aerosols represent one of the largest sources of uncertainty in Earth's radiative budget, primarily through their interactions with solar and terrestrial radiation. These radiative perturbations arise through key aerosol–radiation interactions and manifest as the direct aerosol effect (Chung, 2012; Yu et al., 2006); the indirect aerosol effect (Lohmann and Feichter, 2005; Nenes et al., 2002) and the semi-direct effect (Koch and Del Genio, 2010; Randles and Ramaswamy, 2010). Despite an estimated net global cooling of around $-\mathrm{1.3}\pm \mathrm{0.7}$ W m⁻² from anthropogenic aerosols−offsetting ∼33 % of the GHGs warming (Li et al., 2022), their radiative forcing quantification remains highly uncertain and poorly constrained. The diverse optical and microphysical characteristics of aerosol species (Li et al., 2022), adds to the substantial spread in aerosol's direct radiative forcing (DRF), amplifying uncertainties in climate sensitivity.

Black carbon (BC) aerosols, emitted primarily from incomplete combustion of carbon exerts strong solar-radiation integrated absorption, impacting regional radiative forcing and heating, where it locally competes with greenhouse gases like CO₂ and CH₄ (Hansen and Sato, 2001; Jacobson, 2001). Yet, its DRF estimates vary widely, for instance, one study reported a global BC instantaneous radiative forcing of +0.9 W m⁻² (+0.4 to +1.2 W m⁻²) (Ramanathan and Carmichael, 2008). In contrast, AeroCom Phase II yielded a lower estimate of +0.23 W m⁻² (Myhre et al., 2013). AERONET-derived absorption aerosol optical depth (AAOD) analyses suggested +0.51 W m⁻² (+0.06 to +0.91 W m⁻²) (Bond et al., 2013), while the IPCC AR5 synthesized expert judgment estimating +0.4 W m⁻² (+0.05 to +0.8 W m⁻²) (Bond et al., 2013; Myhre et al., 2013). The recent IPCC AR6 reports a ∼50 % reduction in BC warming and an associated top of the atmospheric (TOA) increase by 0.1° (Forster et al., 2023) due to adjusted rapid climate responses (Everett et al., 2022), however, several studies emphasize that emission products used by the CMIP-6 models (used in the IPCC AR6) did not adequately account for BC loading as well as varying absorption incurred from dynamic variability in morphological complexities (Chen et al., 2022; Kelesidis et al., 2022; Ramachandran et al., 2023).

Both ground and satellite-based observational communities are rigorously working to improve quantification of BCDRF regionally. Ground-based BC mass measurements and/or sun-sky photometer's aerosol optical depth (AOD) provide some essential constraints for optical modeling. These are typically used in iterative Mie model frameworks to derive single scatter albedo (SSA), phase function and asymmetry coefficient (ASY) for radiative transfer model (RTM) simulations (Mehrotra et al., 2024; Zeb et al., 2020). However, these approaches are spatially limited, rely on limited size bins and lacks ability to resolve internal mixing. Furthermore, few observational studies have either interpolated or extrapolated absorption or extinction between a small number of distinct wavebands (where observations are taken) and/or use fixed mass absorption coefficients (MAC of 7.5±1.2 and 8.0±0.7 m² g⁻¹ at 550 nm) (Liu et al., 2020a; Randles et al., 2017; Zangmeister et al., 2018) to compute BC DRF. Satellite-based approaches extend the spatial coverage to better capture regional variability of regional aerosol DRF (Bellouin et al., 2005; Kahn et al., 2023; Matus et al., 2019; Wu et al., 2024). The Clouds and the Earth's Radiant Energy System  CERES) and CloudSat (2B-FLXHR-LIDAR) provides a comprehensive dataset of TOA fluxes which are widely used (Gautam et al., 2023; Loeb et al., 2021; Subba et al., 2020, 2022). While these approaches have analyzed total aerosol DRF, however, accurate isolation of BC signals from co-existing aerosol types within the satellite data is challenging. Furthermore, computing no-aerosol fluxes to derive net TOA forcing adds to the existing challenge. Studies often assume either linear AOD-flux relationships (computing TOA_{no-aerosol-flux} at AOD = 0) or using low-aerosol pixels for pristine flux (Biswas et al., 2017; Christopher and Zhang, 2002; Sundström et al., 2015). However, such approaches can introduce significant uncertainties in regions undergoing rapid industrialization or experiencing persistently high aerosol loading. such as South and East Asia, where truly pristine conditions are rare and the aerosol–radiation interactions are highly nonlinear. Finally, many satellite derived DRF approaches rely on single-wavelength absorption measurements (commonly at 500 or 550 nm), which increases DRF uncertainty as it neglects the spectral dependence of aerosol–radiation interactions across the solar spectrum. Complementarily, reanalysis products like MERRA-2 and CAMS assimilate multi-source data and provide different BC specific properties. However, these products also incorporate simplified aerosol representations like using pre-defined morphology with median radii and standard deviations that are climatologically fixed limiting reanalysis accuracy (Chen et al., 2023; Chin et al., 2002; Ding and Liu, 2022). Moreover, emissions datasets underlying these assimilations are often outdated (Li et al., 2024), particularly over many rapidly developing and industrial areas in the Global South (Dasari et al., 2020; Ren et al., 2025), and many employ externally mixed tracers (Fernandes et al., 2023; Xu et al., 2020) rather than internal mixtures.

Overall, these approaches have advanced our understanding of aerosol DRF, but still leave persistent gaps in capturing BC's dynamic radiative behavior at regional scales. The underlaying uncertainties in BC DRF quantifications are further aggravated at regional level, due to BC's short atmospheric lifetime, rapidly evolving microphysics and multispectral absorption interaction. Reducing these discrepancies requires a geospatially resolved approach that uses multi-waveband information, considers both columnar loading and per-particle properties in tandem (Tiwari et al., 2023). Additionally, an improved accounting of emissions also requires accurate knowledge of microphysical evolution (Fierce et al., 2020; Liu et al., 2017), morphological dynamics (Wu et al., 2018), and vertical lofting of BC (Cohen and Wang, 2014). Accurate regional BC DRF quantification is therefore critical, not only to resolve local radiative impacts but also to reduce uncertainties in global forcing estimates.

In light of persistent challenges in regional BCDRF quantification and building on recent advances in observationally constrained BC microphysical retrievals, including the core-shell Mie model optimization (COSMO), this study extends and refines multi-platform framework of an earlier approach developed by Tiwari et al. (2025). That framework combined multi-wavelength observations with a flexible core-shell morphology to derive probabilistic BC microphysical solutions and column densities, jointly constrained by multi-band SSA. Here, we introduce an additional observational constraint through multi-satellite total AOD inputs, further optimizing the estimation of BC size, mixing state, and column loading, along with their associated optical properties. These physically consistent ensembles are then coupled with RTM to generate high resolution clear-sky BC top-of-atmosphere direct radiative forcing (BC TOA; used interchangeably with BC TOA DRF). Finally, we train a set of linear, parametric and interpretable machine learning (ML) surrogate models on these COSMO-RTM outputs to test their predictive skills, provide a fast and transferable tool for estimating BC TOA DRF and diagnosing how BC microphysics modulate radiative forcing. Unlike most previous assessments which have generally not incorporated dynamically varying size-mixing states (e.g. Luo et al., 2022; Romshoo et al., 2024; Yang et al., 2025) or jointly constrained microphysics and column loading in TOA forcing calculations (e.g. Tiwari et al., 2023), this unified approach integrates these aspects to enhance understanding of regional BC TOA DRF while providing a physically interpretable, transferable and computationally efficient modeling tool and explore new aspects of BC's microphysical and column property influence on clear-sky instantaneous BC TOA DRF .

2 Methodology

2.1 Computing BC TOA DRF from observationally constrained microphysics and RTM

Tiwari et al. (2025) used a multi-platform COSMO approach to derive BC aerosol size, mixing state, and column densities over Xuzhou, China and greater Dhaka, Bangladesh. Greater Dhaka is one urban agglomeration representing a rapidly expanding monocentric megacity in a tropical river delta along with its per-urban settlements. They are characterized by intense emissions from traffic, brick kilns, and small-scale industries along the inland river networks. In contrast, Xuzhou, as the largest city in the Huaihai region is a polycentric urban agglomeration on the North China Plain, where it acts as a central node in a network of cities with significant heavy industry, energy production, and agriculture. Comparing these two regions provides unique insights into how different urban forms and regional economic structures influence the physical and optical properties of light-absorbing aerosols.

The original framework combined multi-waveband single-scattering albedo (SSA) from TROPOMI (388 nm) and AERONET/SONET (440–1020 nm) observations for these regions and used a two-component Mie inversion to generate ensembles of physically consistent microphysical solutions for BC core and non-absorbing shell sizes while accounting for both waveband offsets and cross-platform SSA uncertainties at a grid-by-grid and day-by-day level. A brief description of the Mie model framework and assumed morphological set up, including key assumptions and their rationale, is provided in Appendix A. In essence, Tiwari et al. (2025) first constructed multi-waveband SSA ranges at each grid based on their uncertainties and waveband offsets. Since each BC particle interacts with radiation across all wavebands simultaneously, the COSMO was used to generate paired core-shell size solutions that are physically consistent across the spectrum. Each candidate size pair was evaluated by comparing its modeled SSA at multiple wavebands to the uncertainty-bounded SSA ranges. Only those solutions with SSA values falling within the SSA uncertainty range at all wavebands simultaneously were retained. This filtering approach has been previously implemented using both ground and satellite SSA observations, both in tandem as well as separately (Liu et al., 2024a, b, c; Wang et al., 2021) . The multi-waveband SSA constraining ensures that the resulting microphysical ensembles reflect only BC core/sulfate shell pairs capable of reproducing the observed, uncertainty-aware SSA constraints across the full spectral domain. For more details on the geographical context of the urban agglomerates including, differences in emission characteristics of Xuzhou and greater Dhaka region, and the application of the multi-platform COSMO with SSA constraint to derive BC microphysics and derivation of BC column loading, readers are referred to Tiwari et al. (2025).

This study extends and refines the framework, where it uses the probabilistic solutions of core and shell sizes and total column mass density derived from earlier SSA constraint for computing BC aerosol optical depth (BCAOD) across multiple wavelengths. The BCAOD calculation follows a forward Mie-based physical formulation. Specifically, the diameters of the BC core (D_core) and total diameter (D) of the particle are used to determine the particle's geometric cross-sectional area (σ_g), which reflects the amount of light that the particles interact with. Next, the volumes of the BC core (V_core), SUL shell (V_shell), and total particle volume (V_total) are used to establish the effective density of the aerosol ρ_total based on the volume weighted of BC (ρ_BC) and sulfate (ρ_SUL). Subsequently, the mass of a single particle (m_p) is calculated by multiplying its total volume (V_total) by its effective density (ρ_total). The mass extinction efficiency (α_ext) is then computed to assess how effectively the particle contributes to light extinction per unit mass using extinction efficiency (Q_ext). Finally, the modeled BCAOD (τ_BC) was deduced by multiplying the mass extinction efficiency by the total aerosol mass (M_c) in the column. The mathematical formulation of BCAOD computation is summarized in Eq. (1) through Eq. (8).

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{g}}={\displaystyle \frac{\mathit{\pi }{D}^{\mathrm{2}}}{\mathrm{4}}}\text{(2)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{core}}={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\times \mathit{\pi }{\left({\displaystyle \frac{D_{\mathrm{core}}}{\mathrm{2}}}\right)}^{\mathrm{3}}\text{(3)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{total}}={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\times \mathit{\pi }{\left({\displaystyle \frac{D}{\mathrm{2}}}\right)}^{\mathrm{3}}\text{(4)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{shell}}=V_{\mathrm{total}}-V_{\mathrm{core}}\end{array}$$

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{total}}={\displaystyle \frac{\left({\mathit{\rho }}_{\mathrm{BC}}\times V_{\mathrm{core}}+{\mathit{\rho }}_{\mathrm{SUL}}\times V_{\mathrm{shell}}\right)}{V_{\mathrm{total}}}}\text{(6)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{p}}=V_{\mathrm{total}}\times {\mathit{\rho }}_{\mathrm{total}}\text{(7)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{ext}}={\displaystyle \frac{Q_{\mathrm{ext}}\times {\mathit{\sigma }}_{\mathrm{g}}}{m_{\mathrm{p}}}}\text{(8)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{BC}}=M_{\mathrm{c}}\times {\mathit{\alpha }}_{\mathrm{ext}}\end{array}$$

The SSA-based constraint ensures internal spectral consistency, but it does not by itself rule out microphysical solutions that generate unrealistically strong BC extinction. To address this limitation, we introduce a physically binding requirement that the modeled τ_BC must always be less than total AOD measured in the column. Since total AOD represents the combined extinction from all aerosol species, not only BC, any case where τ_BC exceeds the total AOD indicates a violation and is therefore physically unrealistic. This requirement acts as an additional layer of observational constraint to refine microphysical and column loading solutions.

The Mie model outputs, τ_BC, SSA and ASY, which are computed at three wavebands 388, 470, and 550 nm. To implement this total AOD constraint, we derive specific upper uncertainty bounds based on AOD measurements from Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm-based AOD (470, 550 nm) (Lyapustin et al., 2018; Martins et al., 2017) and TROPOMI (388 nm). The upper limits for AOD are defined in Eqs. (9)–(11). based on Lyapustin and Wang (2022) and de Graaf (2022):

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }\mathrm{Upper}{\mathrm{AOD}}_{\mathrm{388}}={\mathrm{AOD}}_{\mathrm{388}}+\max \left(\mathrm{0.1},\mathrm{0.25}\times {\mathrm{AOD}}_{\mathrm{388}}\right)\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{Upper}{\mathrm{AOD}}_{\mathrm{470}}={\mathrm{AOD}}_{\mathrm{470}}+\mathrm{0.05}+\left(\mathrm{0.1}\times {\mathrm{AOD}}_{\mathrm{470}}\right)\text{(11)}& {\displaystyle }& {\displaystyle }\mathrm{Upper}{\mathrm{AOD}}_{\mathrm{550}}={\mathrm{AOD}}_{\mathrm{550}}+\mathrm{0.05}+\left(\mathrm{0.1}\times {\mathrm{AOD}}_{\mathrm{550}}\right).\end{array}$$

The above violation condition was implemented whereby τ_BC must remain below the upper-bound total AOD observed across all wavebands simultaneously. If τ_BC exceeds the total AOD for even one waveband, that dataset is excluded from further analysis. After implementing this quality control, a total of 178 and 314 d of data were used for analysis for Xuzhou and greater Dhaka respectively from the period of 30 April 2018 to 30 April 2022, consisting of a total of 105 569 and 273 225 individual grid points in space and time.

The τ_BC, SSA and ASY at multiple wavebands, are used within the RTM model to generate BC DRF. Additional necessary datasets, including MODIS Level 2 precipitable water vapor (PWV) (MYD05/MOD05), surface spectral albedo from the MODIS Bidirectional Reflectance Distribution Function product (MCD43A1), and total column ozone from the TROPOMI S5P_L2__O3_TOT_HiR product, were processed and spatio-temporally harmonized to align with the BC data grid resolution. The RTM model used herein is the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) (Ricchiazzi et al., 1998) model which generates the TOA and surface radiative forcing. The description and processing of each ancillary dataset is provided in Appendix B and mathematical formulation for calculating BC TOA DRF from SBDART derived fluxes are provided in Sect. S1 in the Supplement. A schematic overview of the workflow, encompassing BCAOD computation, multi-waveband satellite data filtering, and radiative transfer modeling, is provided in Fig. 1.

Figure 1 Overview of the COSMO-coupled RTM framework for deriving BC radiative forcing.

Overall approach integrates TROPOMI and AERONET/SONET SSAs using waveband offset adjustment and further employs observationally constrained multi-waveband inverse-forward MIE modeling to generate size, mixing state and associated optical and physical properties as detailed in Tiwari et al. (2025). This study further refines the established approach using an additional layer of observational constraint from multi-waveband aerosol optical depth. The resulting set of microphysical and associated optical properties are integrated with PWV, surface spectral albedo, and total column ozone are used to perform SBDART RTM runs and BC radiative forcing at a day-by-day and grid-by-grid resolution over Xuzhou and greater Dhaka are derived.

2.2 Surrogate models for rapid BC TOA DRF estimation

High-fidelity RTMs such as SBDART, when simulated with morphologically dynamic BC aerosol properties and environmental variables can enhance the accuracy of BC DRF estimates. However, computational demands render the approach impractical for large-scale chemical transport or climate modeling applications. As a means of coping with the computational issue, the atmospheric science community often adopts a linear radiative efficiency approach specifically for total aerosol DRF, by assuming the TOA as a linear function of AOD (García et al., 2012; Sena et al., 2013; Sundström et al., 2015). While this method enables rapid regional-to-global scaling of DRF estimates using satellite AOD products, it overlooks non-linear modulations in BC TOA DRF caused by microphysical variability (e.g., particle size, mixing state), column number loading, and their role in both absorption and extinction efficiency. Similarly, BC TOA forcing regimes, characterizing the transition from strong cooling to warming effects were defined by stratifying SBDART-derived BC TOA ensemble for both the regions into percentile-based categories. Low BC TOA (≤ 25th percentile), Mid-Low BC TOA (25th–50th percentile), Mid-High BC TOA (50th–75th percentile), and High BC TOA (≥ 75th percentile).

This study explores this gap by developing three surrogate models for fast BC radiative forcing calculation using the full physical SBDART-derived BC TOA DRF outputs:

• Linear model: mimicking the community's approach of scaling TOA with AOD at 550 nm for total aerosol DRF (Matus et al., 2019; Sena et al., 2013; Sena and Artaxo, 2015), we derive a regression slope (β) between SBDART-generated BC TOA forcing and BCAOD₅₅₀.

  $$\mathrm{BC}\mathrm{TOA}\sim {\mathrm{BC}}_{\mathrm{AOD}\mathrm{550}}\cdot \mathit{\beta }$$

• Parametric multiple linear regression (MLR) model: expanding on the previously proposed MLR equation (Tiwari et al., 2023) we incorporate additional predictors; column mass (M_c) and number concentration (N_c), BC size, and mixing state along with BCAOD₅₅₀ into the MLR framework

  $$\begin{array}{ll}{\displaystyle }\mathrm{BC}{\mathrm{TOA}}_{\mathrm{MLR}}& {\displaystyle }={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{1}}\cdot {\mathrm{BCAOD}}_{\mathrm{550}}+{\mathit{\alpha }}_{\mathrm{2}}\cdot N_{\mathrm{c}}+{\mathit{\alpha }}_{\mathrm{3}}\\ {\displaystyle }& {\displaystyle }\cdot M_{\mathrm{c}}+{\mathit{\alpha }}_{\mathrm{4}}\cdot \mathrm{BC}\mathrm{size}+{\mathit{\alpha }}_{\mathrm{5}}\cdot \mathrm{mixing}\mathrm{state}.\end{array}$$

• Machine learning (ML) model: while MLR improves BC TOA prediction by accounting for additional variables, it assumes fixed, first order or pseudo-linearized, additive contributions of predictors (Cohen and Prinn, 2011). The ensemble learning algorithm herein employs a random forest approach, which helps explore the complex relationships and inter-dependencies between the predictors and physically derived TOA, providing a more robust framework for BC TOA estimation. Specifically, random forest regressor was used to predict COSMO-RTM derived BCTOA from the COSMO-derived predictor variables. The predictors were BCAOD₅₅₀, aerosol column number density (N_c), aerosol column mass density (M_c), BC core size, and mixing state (defined as the ratio of BC core size to total size). These variables are derived from the observation-constrained COSMO retrieval framework described in Sect. 2.1. For each region's RTM output BC TOA is subjected to a 70–30 train-test split, with 70 % of the data used to train the model and 30 % withheld as an independent test set. This approach allowed the model to learn relationships between predictors and BC TOA while ensuring a robust evaluation on unseen data. The Random Forest hyperparameters were set as follows: n_estimators = 100, max_depth = 10, min_samples_split = 5, min_samples_leaf = 3, and max_features = “sqrt”. These regularized settings were used to provide a stable ensemble prediction while limiting overfitting through restricted tree depth, minimum node sizes, and randomized feature selection. The statistical metrics – coefficient of determination (R²), mean bias error (MBE), mean absolute error (MAE) and root mean square error (RMSE) – for both the training and test data are presented, along with the corresponding predicted versus actual plots for each set along with their statistics and error metrics provided in Fig. S1 and Table S1 in the Supplement respectively. The ML-model demonstrated robust predictive capabilities in analyzing the relationship between BC TOA and key predictors in both Xuzhou and Dhaka, with R², RMSE, and MAE of 0.97, 1.7 W m⁻², and 1.1 W m⁻² respectively on the training set and 0.96, 1.8 W m⁻², and 1.2 W m⁻² respectively on the test set, capturing approximately 96 % of the variance in BC TOA in Xuzhou. Similarly, in Dhaka, the model achieved an R², RMSE, and MAE of 0.96, 1.5 W m⁻², and 1.1 W m⁻² respectively for the training and set and 0.96, 1.6 W m⁻², and 1.1 W m⁻² respectively for the test dataset. In all cases, the MBE remained near zero, indicating negligible bias. Cross-validation using a 5-fold split yielded a mean R² of 0.96 with a standard deviation of 0.001, reflecting the model's consistency across different training subsets. Additionally, the learning curve analysis observed converging curve for RMSE as the datasets were added in subsequent runs, with training and validation curve closely converging. These metrics and observations assure confidence to apply the trained model across all the datasets and subset of these datasets under varying conditions.

Currently, to the best of our knowledge, no observational products directly quantify BC TOA DRF while fully accounting for dynamic BC microphysics and using multi-wavelength optical properties in tandem. Although satellite and ground-based platforms provide key aerosol optical parameters (e.g., AOD, AAOD and SSA), they do not provide BC TOA DRF for direct comparison. Therefore, we use physically based COSMO-RTM model outputs, that accounts for BC microphysics and optical properties, as the reference standard to compare our surrogate models. This ensures comparison against a comprehensive, physically consistent benchmark in the absence of direct observational data. As an additional check on the physical credibility of this COSMO-RTM BC TOA DRF target, an AERONET-based forcing-efficiency consistency test was performed over Dhaka using Version 3 Level 2.0 inversion radiative forcing products (Sect. S2). The analysis shows that COSMO-RTM BC TOA forcing efficiency exhibits regime-dependent consistency with AERONET total-aerosol TOA forcing efficiency under absorbing, spectrally BC-like conditions (Fig. S2), and that observed proportional scaling falls within the expected BC-to-total TOA forcing efficiency range (Table S2).

Furthermore, to better understand how different predictors influence BC TOA DRF, each surrogate model is tested across different pollution and morphological regimes. Pollution conditions were categorized based on BCAOD₅₅₀ values, while different morphological regimes were categorized based on BC size and mixing state (ratio of BC core size to total size) conditions. These conditions were generated following percentile-based division as shown in Table S3. A series of performance and error metrics – adj R², root mean square error (RMSE), mean Absolute and bias error, mean absolute percentage error (MAPE), bias-corrected RMSE and confidence intervals etc. – are computed and compared, to assess the performance of different models detailed in.

2.3 SHAP (SHapley Additive exPlanation) analysis

To further explain the feature contributions, Shapley Additive Explanations (SHAP) analysis is used, where SHAP values, derived from cooperative game theory (Lundberg and Lee, 2017), quantify the contribution of each feature by considering all possible combinations. These values capture both individual and joint feature contributions, offering a fair and consistent explanation of the model's predictions. Additionally, violin plots are employed to visualize the distribution of SHAP values across predictions, providing insights into how feature contributions vary.

Figure 2 Conceptual schematic for conducting cross-regional (ML^*) and combined (ML${}^{**}$) analysis.

2.4 Testing ML models transferability: cross-regional and combined approach

Finally, to test if developed surrogate ML models can serve as a fast BC radiative transfer tool in other geographic regions, a cross-regional and combined model testing analysis was carried out. First, the ML model generated in one region (e.g., Dhaka) was applied to predict BC TOA DRF in the other region (e.g., Xuzhou), and vice versa. The specific cross-regional model (ML^*) validation tested the transferability of region-specific models informed by local BC properties into different regions with their own unique local BC properties, providing critical insights into the extent to which these models perform across distinct environmental, microphysical, and total column conditions using the same sets of predictors. Additionally, a combined model (ML${}^{**}$) was developed by implementing stratified regional testing, combining data from both the regions as inputs, hence learning the variability in how BC TOA is physically computed based on relationships using predictors based on environmental conditions across both regions in tandem. Subsequently ML${}^{**}$ is again tested over both Xuzhou and Dhaka and is compared with SBDART BC TOA. The conceptual overview of this framework is summarized in Fig. 2.

Finally high resolution (0.05°×0.05°) long-term average spatial maps of SBDART, ML^* and ML${}^{**}$ BCTOA were generated and spatial correlation are analyzed to assess how well the surrogate models capture spatial variability over a long-term average time frame. Additionally, comparison between the local ML and the ML${}^{**}$ model is performed over these spatial aggregated climatologies, by computing error improvement metrics and spatial difference in absolute prediction errors evaluated against SBDART.

3 Results

3.1 Comparative analysis of black carbon microphysical, optical, and columnar properties in Xuzhou and Dhaka: key predictors of BC TOA DRF

Optical, microphysical and column abundance properties like BCAOD, BC size (radius in nm) and mixing state (ratio of BC size to total size), N_c, M_c, BC core size directly or indirectly influence aerosol–radiation interactions, ultimately impact TOA forcing. Figure 3, below presents the distribution of these parameters via probability density plots over both the regions.

Figure 3 Probability density distributions of key predictors for aerosol radiative forcing analysis in Xuzhou and Dhaka: (a) Black Carbon (BC) core size, (b) mixing state (ratio of BC size to total size), (c) number column density (N_c), (d) total mass column density (M_c), and (e) BCAOD at 550 nm, (f) BC core to Sulfate shell size ratio (BCAOD₅₅₀ > 75th percentile), (g) BC core size for polluted condition (BCAOD₅₅₀ > 75th percentile), (h) BCAOD₅₅₀ distribution (BC core size < 25th percentile or under barely coated conditions).

The regional differences in per-particle BC size and mixing state between Dhaka and Xuzhou revealed contrasting behavior. First, BC core size (Fig. 3a) showed a broader and mostly higher range in Dhaka (from 50 to 320 nm, with mean/25th/50th/75th percentiles distributed as 162.5/130/156/190 nm) than in Xuzhou (from 50 to 260 nm, with corresponding statistics of 140.5/110/134.4/170 nm). The mixing state distribution was also distinct (Fig. 3b), with the value over Dhaka both narrower and higher (from 0.16 to 0.29 with mean/25th/50th/75th observed as 0.23/0.21/0.23/0.25) than over Xuzhou (0.12 to 0.28, with corresponding statistics of 0.2/0.17/0.2/0.22). In general BC aerosols over Xuzhou are emitted as smaller particles which then more rapidly coat, while over Dhaka, larger BC cores are emitted and the coating is lesser. These results are consistent with the wide range of co-emitted sources and higher in-situ secondary aerosol and other atmospheric processing occurring in Xuzhou, and the less efficient combustion leading to more and larger soot monomers combined with less oxidative atmosphere in Dhaka leading to less overall secondary aerosol production. These factors play a crucial role in shaping the microphysical properties of BC, and subsequent influence on the solar radiation stream and overall atmospheric effects.

N_c exhibited distinct distribution patterns (Fig. 3c) with the value in Xuzhou generally higher on average but less extreme (from 6.4×10⁹ to 2.7×10¹² # m⁻², with the mean/25th/median/75th values distributed as 3.0×10¹¹/1.6×10¹¹/2.3×10¹¹/3.5×10¹¹ # m⁻²) than in Dhaka (from 4.5×10⁹ to 3.6×10¹² # m⁻², with corresponding statistics 2.7×10¹¹/1.5×10¹¹/2.2×10¹¹/3.5×10¹¹ # m⁻²). Xuzhou exhibits a consistently higher baseline for N_c, as indicated by its minimum, 25th percentile, and median values, consistent with more sustained yet regulated aerosol sources from known industrial and urban centers, as well as influence from substantial regional transport. Dhaka shows slightly lower mean and baseline values, although its maximum N_c is substantially higher, highlighting the occurrence of extreme aerosol loading events, which are consistent with strong yet episodic sources, including extensive contributions from biomass burning, and brick kiln operations, as well as influence from stagnant atmospheric conditions, which exacerbate aerosol accumulation in the region (Rahman et al., 2019; Salam et al., 2021; Zaman et al., 2023). Interestingly, the 75th percentile values of N_c are comparable across both regions indicating similar levels of moderately high loading.

M_c exhibited a broader range and higher loading across all percentiles in Xuzhou compared to Dhaka (Fig. 3d). Over Xuzhou, M_c ranged from 2.8 to 1668.2 mg m⁻² with mean/25th/median/75th values of 254.2/151.8/222.1/309.6 mg m⁻², while over Dhaka M_c ranged between 1.9 to 1052 mg m⁻² with corresponding values of 211.9/113.3/181.2/284.9 respectively. This is concordant with a previous study which pointed out that the higher M_c in Xuzhou is mostly associated with the larger aerosol coating compared with Dhaka (Tiwari et al., 2025). This discrepancy is consistent with lower combustion efficiency in Dhaka and significant influences of agricultural and biomass burning, leading to larger cores but also slower secondary aerosol growth. Additionally, such conditions lead to reduced local solar radiation required for rapid secondary aerosol production of inorganic aerosols compared to Xuzhou. In contrast, Xuzhou's elevated BC levels are attributed to its higher total power consumption and diverse industrial emissions, which arise from both local and regional transport activities. This finding is consistent with observations from other urban and suburban regions throughout central and western China, where diverse emission sources, effective secondary aerosol formation processes, and significant co-emissions have been identified as contributing factors (Liu et al., 2024c). Consequently, these long-range transport events in Xuzhou may also further contribute to the higher M_c.

Notably, BCAOD at 550 nm showed comparable values between Xuzhou and Dhaka (Fig. 3e), ranging from 0.01 to 2.2 with mean/25th/median/75th values distributed as 0.5/0.3/0.4/0.5 in Xuzhou and 0.4/0.3/0.4/0.6 in Dhaka. Detailed comparison reveals distinct regional differences in BC dynamics, with Xuzhou experiencing more extreme aerosol events contributing to BC loading, as reflected by its higher maximum and elevated mean values. These episodic events likely enhance the BC loading and aging, in turn enhancing net extinction from inorganics, contributing to the longer tail in the distribution. In contrast, Dhaka exhibits a more uniform distribution of moderately high values, as evidenced by its slightly higher 75th percentile but lower overall maximum, suggesting high but consistent aerosol loading, potentially driven by continuous emissions from urban and regional sources. The elevated BCAOD₅₅₀ in Xuzhou is likely influenced by both local emissions and regional loading, as supported by Qin et al. (2022), which observed transported aloft aerosols promote rapid formation of secondary inorganics like secondary sulfate over Xuzhou, which when coated with BC cores can enhance the overall extinction. Additionally, Chen et al. (2019) previously reported high BC to PM ratio in Xuzhou often associated with both southern and northern winds originating from major economic areas in the lower Yangtze River Delta and Beijing-Tianjin-Hebei provinces. Source attribution has pointed to Shandong as well as central provinces, Henan and Anhui, as significant contributors to BC loading over Xuzhou.

Further analysis of high polluted conditions (BCAOD₅₅₀ > 75th percentile) demonstrates distinct BC core size distributions (Fig. 3f and g). The polluted-condition core-to-shell ratio distribution indicates relatively thicker sulfate coatings in Xuzhou, whereas Dhaka exhibits larger core-to-shell ratios, consistent with comparatively larger BC cores and thinner relative coatings. The polluted-condition BC core size PDF further shows that Xuzhou is dominated by a more compact small-to-medium core-size population, while Dhaka exhibits a broader distribution extending toward larger BC core sizes.

The corresponding two-dimensional core-shell KDE analysis further supports this interpretation (Fig. S3). Xuzhou shows a compact distribution with BC cores mostly between 100–120 nm and sulfate shells of 500–600 nm, plus a smaller cluster of 80–90 nm cores with 390–410 nm shells, indicating well-coated small to medium cores. Dhaka's primary cluster features 100–120 nm cores with 390–450 nm shells, and a secondary cluster of larger 140–160 nm cores paired with 550–620 nm shells. The BCAOD₅₅₀ probability density function for small cores (Fig. 3h) associates these sizes with elevated pollution in both areas, while core-to-shell ratios highlight a dominance of thicker shell coatings in Xuzhou, reflecting strong secondary inorganic formation. Conversely, Dhaka's narrower small-core distribution and larger core-to-shell ratios imply frequent pollution events involve larger BC cores with relatively thinner coatings. These contrasts emphasize secondary inorganic growth drives Xuzhou's polluted conditions, while primary emissions drive Dhaka's polluted conditions.

3.2 Evaluating model performance: overall, and across varied pollution condition and morphological regimes

As outlined in Sect. 2.2, the three surrogate models span a deliberate complexity gradient, from AOD-only linear scaling to multivariate parametric regression, and ultimately nonlinear ensemble learning, each model tier is designed to probe a distinct aspect of BC–radiation interactions. Their performance is assessed not only in aggregate but across distinct pollution, coating, and size regimes, where the physical drivers of BC TOA forcing are expected to differ most.

Figure 4 Comparison of BC TOA DRF predictions from various models (Linear, MLR, and ML) with SBDART BC TOA DRF for Xuzhou (a–c) and Dhaka (d–f). The scatter plots (a, d) display the relationship between model predictions and SBDART values, with color maps indicating the density of data points.

The performance of the Linear, MLR, and ML models in predicting SBDART TOA were visualized through 2D density plots for both Xuzhou (Fig. 4a–c) and Dhaka (Fig. 4d–f). The ML-model, MLR, and linear models reproduce the SBDART-derived TOA radiative forcing (mean: Xuzhou = $-\mathrm{17.2}\pm \mathrm{9.2}$ W m⁻²; Dhaka = $-\mathrm{15.1}\pm \mathrm{7.7}$ W m⁻²) with moderate-high fidelity, as evidenced by near-identical means and low percentage errors (ML: 7.8 %–8.5 %; MLR: 13.0 %; Linear: 16.0 %–20.4 %) (Fig. 4). The progressively tighter clustering around the 1:1 line from Fig. 4a–c and d–f confirms that the ML-model exhibits superior accuracy with the lowest RMSE (Xuzhou: 1.7 W m⁻²; Dhaka: 1.5 W m⁻²) and highest R² (0.96–0.97) (Fig. 4c and f), reflecting its improved capacity to capture non-linear BC–radiation interactions. The linear model underestimates extreme BCTOA values (e.g., Xuzhou 1st percentile: Linear = −49.0 W m⁻² versus SBDART = −45.4 W m⁻²), while ML-model closely matches SBDART distributions (Xuzhou 1st percentile: ML = −43.9 W m⁻² versus SBDART = −45.4 W m⁻²), particularly at tails as shown in Fig. S4. In Dhaka, ML's narrower 95 % confidence intervals (−0.002 to 0.009 W m⁻²) and lower MAE (1.1 W m⁻²) highlight its robustness in capturing TOA variability. Both regions show systematic model divergence at warming conditions, where linear and MLR models underestimate warming (e.g., Xuzhou 99th percentile: Linear = −7.2 W m⁻² versus SBDART = −4.1 W m⁻²), whereas ML marginally minimizes bias (ML = −4.5 W m⁻²) (Fig. S4).

The ML, MLR, and linear models exhibit distinct performance patterns across BCTOA radiative forcing regimes, revealing significant variations in accuracy, precision, and robustness. In Xuzhou, the linear model exhibited the highest errors in all BC TOA categories, with a mean error of 1.98 W m⁻² (95 % CI: 1.92–2.03) in Low BC TOA (more cooling scenarios) and −2.9 W m⁻² (95 % CI: -2.89 to -2.83) in High BC TOA (more warming scenarios), alongside persistent bias (RMSE = 3.4 W m⁻² in High BC TOA). The MLR model showed intermediate performance, reducing errors compared to linear (e.g., MAE = 2.8 W m⁻² in Low BC TOA versus Linear's 4.4 W m⁻²) but still underperforming ML (MAE = 1.6 W m⁻²). In High BC TOA, MLR's bias corrected RMSE (1.8 W m⁻²) matched Linear's, indicating unresolved bias, while ML achieved the lowest RMSE and bias corrected RMSE (1.2 and 1.1 W m⁻²) and highest R² (0.73). In Dhaka, similar trends emerged but with smaller errors. Linear model had a mean error of 1.1 W m⁻² in Low TOA and −1.2 W m⁻² (95 % CI: −1.2 to −1.1) in High BC TOA, with RMSE and bias corrected RMSE as 3.5 and 3.3 W m⁻² and 2.0 and 1.6 W m⁻² in Low and High BC TOA scenarios. MLR improved over Linear (e.g., MAE = 2.1 W m⁻²in Low BC TOA versus Linear's 2.8 W m⁻²) but lagged behind ML (MAE = 1.5 W m⁻²). In High BC TOA, MLR's R² (0.58) was higher than Linear's (0.47) but lower than ML's (0.78), with ML achieving the lowest RMSE (1.0 W m⁻²).

An interesting pattern is observed over the range of TOA values transitioning between strong cooling and slight cooling/warming, hereafter considered the transition- BCTOA regimes. For both Mid-Low and Mid-High BC TOA categories, all the models demonstrated lower R² despite similar or lower error metrics than the Low and High BC TOA across some models. In both Xuzhou and Dhaka, the Mid-Low and Mid-High BC TOA categories revealed a striking inconsistency between error metrics and variance explained. For instance, in Xuzhou's Mid-Low category, the ML model achieved a low MAE of 1.1 W m⁻² and RMSE of 1.6 W m⁻², comparable and better to its performance in extreme BC TOA ranges (e.g., Low BC TOA: MAE = 1.6 and RMSE = 2.4 W m⁻² and High BC TOA: MAE = 0.9 W m⁻², RMSE = 1.2 W m⁻²). Similarly, in Dhaka's Mid-Low category, ML had MAE = 1.2 W m⁻² and RMSE = 1.6 W m⁻², better than its performance in the low BC TOA range (MAE = 1.5 W m⁻², RMSE = 2.0 W m⁻²). However, R² values in mid-range categories were significantly lower (e.g., Xuzhou Mid-Low; R²=0.56 versus High BC TOA: R²=0.73 and Low BC TOA: R²=0.89; Dhaka Mid-Low: R²=0.50 versus High BC TOA: R²=0.78 and Low BC TOA: R²=0.87). This pattern persisted across models. In Xuzhou's Mid-High category, MLR exhibited a mean error near zero (−0.4 W m⁻²) but explained only 26 % of variance, while linear model had a 11.9 % MAPE; lower than Low and High BC TOA regimes, but R² of 0.13. In Dhaka's Mid-High category, MLR achieved a low MAPE of 13.2 % when R² was 0.33, compared to its performance at High TOA (MAPE = 21.5 % R²=0.58). These results highlight a bias-variance trade-off, meaning, the models tend to capture the average radiative forcing magnitude fairly well (low RMSE/MAE) but failed to explain data variability in mid-range BC TOA conditions. This discrepancy likely arises from the complex, non-linear interplay of aerosol microphysics and atmospheric dynamics that dominate the transition between cooling to warming BC TOA condition, as well as daily-scale variability which occurs underlying these special events as they crossover from one regime to the other. For instance, transient processes such as mixing-state evolution (e.g., sulfate coating on BC cores; Cappa et al., 2012; Fierce et al., 2020; gradual growth in per particle size and mass; McFiggans et al., 2006, and impacts of changes in precursor species on the size distribution; Seinfeld et al., 2003), can introduce high intra-category variability, or spatial and temporal variation, which linear and MLR framework cannot resolve and ML only partially captures. These regimes also involve competing absorption-scattering balances, where minor shifts in aerosol composition or morphology disproportionately alter per-particle properties (Guan et al., 2026) and even forcing directionality, amplifying unexplained variance (Samset et al., 2018b).

To assess the model performance under different and microphysical conditions, performance comparison across different mixing state and BC core size conditions were analyzed and tabulated in Table 1.

Table 1 Performance comparison of Linear, MLR, and ML surrogate models across BC core size and mixing state regimes (all model differences significant at p<0.05). Reported metrics are Adj. R², errors (RMSE, MAE, MBE (all in W m⁻²), and MAPE (%)).

For heavily coated BC, ML achieves superior accuracy in both Xuzhou (Adj. R²=0.97) and Dhaka (Adj. R²=0.93), reducing RMSE by 59 % and 46 % relative to linear and MLR models, with near-zero mean bias (MBE: 0.1 versus 2.2 in Xuzhou; 0.05 versus 1.4 in Dhaka). Under moderately coated conditions, ML remains dominant, lowering MAE 46 %–50 % in Xuzhou (1.3 W m⁻² versus 2.4 W m⁻²) and MAPE by 19 %–24 % in Dhaka (8.2 % versus 10.7 %). Even for minimally coated BC, ML demonstrates strong predictive skills, reducing RMSE by 35 % compared to MLR over Xuzhou and lowering MAPE by 42 % in Dhaka. Notably, for barely coated BC, linear modeling fails entirely in Xuzhou (Adj. $R^{\mathrm{2}}=-\mathrm{0.20}$), whereas, ML retains robust performance (Adj. R²: 0.81, RMSE: 1.4 W m⁻², MAPE: 11.0 %), underscoring its resilience under low-coating conditions. In Dhaka, ML's MAPE (10.0 %) reducing MAPE by 38 % lower than MLR's (16.1 %), with near-perfect bias neutrality (MBE −0.02 W m⁻²). Linear and MLR models exhibit stronger systematic bias in these regimes including underprediction (MBE: −2.7 W m⁻²) in Xuzhou, and slight overprediction (MBE: 0.5 W m⁻²) in Dhaka. Regarding BC core size, ML achieves high precision for small cores in both cities (Adj. R²=0.98/0.96, RMSE = 1.7/1.4 W m⁻², MAPE = 4.9 %/5.0 % in Xuzhou/Dhaka, reducing RMSE errors by 45 %–60 % over Xuzhou and 30 %–46 % over Dhaka compared to MLR and linear models. Transitioning to very big BC cores, the linear frameworks perform poorly in Xuzhou (Adj. R²=0.29; MAPE = 38.6 %), while ML retains robust accuracy (Adj. R²=0.84, MAPE = 12.3 %). In Dhaka, ML reduces RMSE by 42 % and nearly halves MAPE relative to MLR. Similarly, model performances across different pollution conditions also demonstrates superior accuracy and robustness of ML compared to its contemporaries; detailed evaluation including Taylor diagram analysis in Fig. S5 within Sect. S3.

Analysis of the underlying microphysical properties reveals that ensemble learning models consistently outperform linear and MLR approaches in reproducing BC TOA, across diverse pollution conditions (clean to polluted), mixing states (heavily coated to barely coated), and size categories (small to very large) in both Xuzhou and Dhaka. Statistically, RF achieves superior accuracy (Adj. R² spanning 0.73–0.98), near-neutral bias ($\left|\mathrm{MBE}\right|\le \mathrm{0.1}$), and error reductions of 30 %–67 % (RMSE) compared to linear models, with performance gaps widening under complex conditions-such as heavily polluted scenarios (Xuzhou: RF RMSE 2.26 W m⁻² versus MLR 3.62 W m⁻²) or minimally coated BC (for Dhaka; RF MAPE 9.1 % versus MLR 15.7 %). Physically, ML's dominance stems from its ability to resolve non-linear interactions between BC core size, coating thickness, and atmospheric aging processes, which linear frameworks oversimplify, including the transition states between negative and nearly zero or even positive BC TOA. The observed bias-variance trade-off in mid-range BC TOA regimes further underscores the complex interplay of BC microphysics and atmospheric dynamics and other sources of non-linearities which ML models partially capture but remains a challenge for conventional approaches.

3.3 Interpreting BC forcing drivers using SHAP Analysis of ML

To interpret how the ML surrogate reproduces COSMO-RTM-derived BC TOA, we applied SHAP analysis to decompose the model predictions into feature-level attributions (Fig. 5). This provides a diagnostic of how BCAOD₅₅₀, column number density, mixing state, BC size, and column mass contribute to predicted BCTOA across Xuzhou and Dhaka, while explicitly accounting for non-linear feature interactions.

Figure 5 Overall feature importance and feature contributions for predicting BC TOA in Xuzhou and Dhaka. The bar plots (left panels) display the overall feature importance based on the mean SHAP value, and the violin plots (right panels) show the distribution of SHAP values for each feature, indicating the direction and magnitude of feature contributions to the model's output. Blue corresponds to low feature values, and red indicates high feature values, with the SHAP values reflecting model-attributed shifts in predicted BCTOA toward relatively warmer positive or cooler negative forcing contributions. All SHAP values represent the impact on predicted BC TOA DRF (in W m⁻²).

Across both regions the ML models exhibit broadly consistent feature rankings, with BCAOD₅₅₀ showing the largest mean absolute SHAP attribution. The mean SHAP associated with BCAOD₅₅₀ is larger in Xuzhou (3.4, 62.1 % positive) than in Dhaka (2.9, 57.3 % positive). However, SHAP distributions reveal pronounced non-linearity across TOA forcing regimes. Lower-to-moderate BCAOD₅₅₀ values are more frequently associated with positive SHAP values, whereas higher BCAOD₅₅₀ values are more frequently associated with negative SHAP values. with BCAOD₅₅₀ contributing positively under clean to moderately polluted conditions (BCAOD₅₅₀; 0.01–0.53), but increasingly negative under moderately to heavily polluted conditions (BCAOD₅₅₀; 0.37–1.9). This sign change in SHAP likely reflects a transition in the modeled radiative response, wherein increases in total extinction and scattering enhancement outweigh absorption enhancement at higher aerosol loadings. The overlap region (0.37–0.53), where both positive and negative BCAOD₅₅₀ SHAP values occur, further indicates that column loading and total extinction alone are insufficient to separate weaker-cooling/warming from stronger-cooling regimes. Instead, the ML interpretation indicates that per-particle microphysical properties-represented here through dynamic coating ratios and BC core size, simultaneously modulate absorption, total scattering, and angular redistribution, producing regime-dependent forcing behaviour. This microphysics-driven BC TOA sign variability, previously demonstrated by Tiwari et al. (2023), likely arises from non-linear changes in absorption, scattering and scattering angle that occur simultaneously with small changes in the mixing ratio and core size in tandem. While such interactions and sensitivity are often high AOD contributing to more negative TOA forcing is consistent with the total attenuation under the conditions where scattering enhancement is larger than the absorption enhancement. In this context, the SHAP analysis further shows that per-particle microphysics and column loading do not map to a single TOA response, instead depending on their coupled state, they are associated with both positive or negative contributions to predicted BC TOA. This behaviour highlights that commonly used assumptions of external mixing and reliance on AOD alone maybe insufficient, as it may yield either enhanced scattering or absorption while systematically missing critical microphysical effects that govern the net TOA response In the range of BCAOD values between 0.37 and 0.53, the per particle microphysics and overall column loading can contribute either positively or negatively to TOA forcing, suggesting that commonly used external mixing assumption and AOD in tandem are insufficient (Mahajan et al., 2013; Mehrotra et al., 2024, 2025), either yielding enhanced scattering or absorption, thus assessing and demonstrating BC mostly exerting a positive TOA missing actual microphysical effects.

Followed by BCAOD₅₅₀, N_c emerges as the second most important predictor in both the regions, with comparable mean SHAP of (1.95 over Xuzhou and 1.85 over Dhaka). In Xuzhou, higher N_c (> 65th percentile) coincides with increasingly negative SHAP contributions, implying the model predicts stronger cooling (more negative BC TOA) under elevated N_c. Dhaka observes a similar finding, but with the range of N_c exhibiting positive SHAP (warming tendency) being narrow (< 40th percentile). While SHAP explains model sensitivity rather than emission source characteristics, this behaviour is consistent with BC populations that remain weakly mixed over a broader range of lower particle numbers in Dhaka, consistent with limited microphysical aging and enhanced absorption efficiency. In contrast, the narrower range of N_c associated with positive TOA contributions likely suggests more rapid microphysical evolution, plausibly facilitated by co-emitted non-absorbing aerosols that promote condensation, coagulation, and coating growth.

The SHAP results further reveals a moderate but systematic role of mixing state in modulating the predicted BCTOA forcing. While the mean SHAP magnitude associated with mixing state is similar in both regions (∼1.5) the sign and distribution vary. The distribution of SHAP values indicates that positive contributions occur more frequently in Dhaka (58.6 %) than in Xuzhou (44.8 %), reflecting differences in how mixing-state variability modulates the modeled TOA response across the two urban regions. Higher mixing state values (≥0.21 in Xuzhou; ≥0.23 in Dhaka), corresponding to thinner coatings and/or larger BC core are generally associated with enhanced BC TOA warming, whereas, lower mixing states (<0.17 in Xuzhou, <0.16 in Dhaka), correspond to thicker coatings and/or smaller core fractions, coincide with negative BCTOA (cooling) via increased scattering. These regional contrasts highlight that mixing-state effects on BCTOA forcing are strongly context-dependent and emerge through interactions with column loading and particle number, rather than acting as an independent control.

In contrast, M_c and BC size exhibit comparatively weaker influence on predicted BC TOA DRF. In Xuzhou, M_c shows a higher mean SHAP than Dhaka (1.24 versus 0.43), with only the upper ∼31 % of column mass loading (>450 mg m⁻²) associated with cooling in the model prediction. This pattern is consistent with conditions where scattering by co-emitted aerosols dominates absorption. In Dhaka, BC core size marginally outweighs M_c in importance, though its overall contribution remains minimal, with only 9 % of larger BC cores (>220 nm) coinciding with weakly enhanced BC TOA warming. This weak sensitivity suggests that, within the observed range, BC size variations exert a secondary influence compared to column optical extinction, abundance and mixing state.

When stratified across distinct BC TOA forcing regimes-from cooling to warming, SHAP analysis reveals that BCAOD₅₅₀, consistently exerts the largest influence on the predicted BC TOA DRF, while mixing state and N_c, provide a secondary but systematic contributions with regime dependent importance. This evolving hierarchy of feature contributions highlights the strongly non-linear nature of aerosol impacts on TOA forcing, consistent with transition behaviour arising from coupled changes in column optical loading and per-particle microphysics, rather than by bulk optical properties alone. The detailed SHAP analysis across BC TOA regimes is presented in Fig. S6 within Sect. S4.

3.4 Adaptability of cross-regional and combined ML surrogates under domain shift

Understanding the cross-regional adaptability of ML models is critical for assessing their transferability in terms of predicting BC TOA over different regions from where the ML model was trained. The performance of the cross-regional ML models is evaluated by applying the Xuzhou-trained model to Dhaka data and vice versa. In parallel, model adaptability is evaluated by using combined model (ML${}^{**}$) surrogate trained on samples drawn from both the regions, and performing comparison with region-specific ML surrogates with SBDART DRF as reference. This analysis revealed that the cross-regional application results in slight performance degradation while still retaining non-trivial predictive skill, as visualized in Fig. 6.

Figure 6 Comparison of SBDART-derived BC TOA forcing with cross-regional (ML^*) and combined (ML${}^{**}$) ensemble learning models for Xuzhou (a, b) and Dhaka (c, d). PDF comparisons between SBDART BC TOA and ML-predicted BC TOA, with statistical summaries of BCTOA distributions at different percentiles (e, f).

When the Dhaka-trained ML model is applied over Xuzhou, the predicted BCTOA forcing remains strongly correlated with SBDART estimates (R=0.96, p<0.05), but exhibits reduced variability (predicted standard deviation = 8.7 W m⁻² versus 9.2 W m⁻²), indicating incomplete capture of Xuzhou's TOA extremes. This is accompanied by a systematic negative bias (MBE = −0.46 W m⁻², 95 % CI = −0.47 to −0.44), corresponding to an overestimation of cooling and an elevated RMSE (2.5 W m⁻²), approximately 24 % higher than that of the local Xuzhou ML model (RMSE = 1.9 W m⁻²). The associated percentage error (14 %) and residual standard deviation (2.5 W m⁻² after bias removal) indicate the performance degradation arises from both reduced variance representation and systematic bias. These differences may primarily be linked to shifts in the distributions of BCAOD₅₅₀, N_c and mixing-state, particularly exerting stronger modulation of BCTOA under warming conditions. Conversely, employing the Xuzhou-trained ML model to Dhaka also retains high correlation (R=0.96, p<0.05) but underestimates TOA variability (predicted standard deviation = 7.3 W m⁻² versus observed standard deviation = 7.7 W m⁻²) and systematically over-predicts warming in Dhaka (MBE = +0.82 W m⁻², 95 % CI = 0.81 to 0.83). These inconsistencies likely stem due to the differences in BC microphysical regimes between the two regions, including variations in BC core size distributions and coating evolution, which are not fully captured when models are transferred across regions. Despite these discrepancies, ML^* models outperform local linear models, and performed comparably against the local MLR models. Over Xuzhou, Dhaka-trained model achieves a higher capability to capture variance (R²=0.92) compared with local linear (R²=0.85) and local MLR, with a 30 % lower percentage error than linear (14.0 % versus 20.4 %). Similarly in Dhaka, the Xuzhou-trained ML (R²=0.91) surpasses local linear (R²=0.88) and is similar to local MLR (R²=0.92). The combined ML${}^{**}$ model, exhibits improved generalization across both regions with minimal loss of regional accuracy. In Xuzhou, the combined model attains an Adj. R²=0.96 (RMSE = 1.90 W m⁻², MAE = 1.4 W m⁻²), closely matching the local ML model (R²=0.97, RMSE = 1.70 W m⁻² and MAE = 1.2 W m⁻²) while significantly outperforming local linear (RMSE = 3.6 W m⁻² and MAE = 2.7 W m⁻²) and local MLR (RMSE = 2.5 W m⁻², MAE = 1.8 W m⁻²). Similarly, in Dhaka, the combined model achieves R²=0.96 (RMSE = 1.6 W m⁻², MAE = 1.2 W m⁻²), nearly equaling local ML (R²=0.96, RMSE = 1.5 W m⁻²) and surpassing local linear (R²=0.88, RMSE = 2.6 W m⁻²) and local MLR (R²=0.92, RMSE = 2.2 W m⁻²). Although the combined model does not consistently outperform region-specific ML surrogates, it substantially reduces biases associated with cross-regional transfer.

3.5 Spatial climatological evaluation and robustness of the ML${}^{**}$ model

Swath-by-swath pixel comparisons between the cross-regional (ML^*) and combined (ML${}^{**}$) models provide insight into short-term, local-scale predictability of BC radiative effects, but do not directly characterize the models' ability to reproduce long-term spatial structure and climatological variability. To assess regional long-term average and spatial generalization, we therefore compare the 0.05° weighted BC TOA forcing climatology simulated by ML^* and ML${}^{**}$ against SBDART benchmarks. This spatial analysis serves two complementary purposes: (i) to assess the fidelity with which each surrogate reproduces long-term spatial patterns and climatological variability of BC TOA forcing (Fig. 7), and (ii) to examine whether the improved cross-regional adaptability of the ML${}^{**}$ model translates into spatial performance relative to region-specific ML surrogates (Fig. 8).

Figure 7 Spatial Evaluation of BCTOA DRF Simulations: Comparison of SBDART, Cross-Regional (ML^*), and Combined Data (ML${}^{**}$) surrogate ML models Over Xuzhou and Dhaka. Panels (a) and (b) show weighted climatology maps of BCTOA DRF simulated by SBDART (reference model), cross-regional (ML^*) model, and combined (ML${}^{**}$) model over Xuzhou (a) and Dhaka (b) computed over valid days between April 2018–April 2022. Panels (c) and (d) present grid-by-grid correlation coefficients between SBDART and ML^* and ML${}^{**}$ for Xuzhou (c) and Dhaka (d). Color bars indicate TOA DRF values (W m⁻²) and correlation coefficients. Stars mark key observational sites (SONET/AERONET stations) over Xuzhou and Dhaka respectively.

Figure 8 Spatial comparison of combined and region-specific ML performance. Panels (a) and (c) show maps of the difference in absolute error ($|{\mathrm{ML}}^{**}-\mathrm{SBDART}|-|\mathrm{ML}-\mathrm{SBDART}|$) over Xuzhou and Dhaka, respectively. Positive values indicate grid cells where the combined ML (ML${}^{**}$) yields lower absolute error than the region-specific ML model, while negative values indicate underperformance. Magenta-outlined grid cells mark locations where ML${}^{**}$ achieves strictly lower absolute error. Panels (b) and (d) present radar plots comparing domain-aggregated performance metrics (RMSE, MAPE, and bias) for ML and ML${}^{**}$ over Xuzhou and Dhaka, respectively. Metrics are normalized using min-max scaling, where 0 denotes the best observed performance and 1 the worst. Stars indicate SONET/AERONET sites.

Over Xuzhou (Fig. 7a), both ML models capture the general spatial distribution of BC forcing, including the broad patterns of cooling (negative forcing) over central and northern regions and warming (positive forcing) in the south. The overall mean BC TOA DRF simulated by SBDART across the domain was $-\mathrm{17.4}\pm \mathrm{2.6}$ W m⁻², with ML^* showing a mean of $-\mathrm{17.8}\pm \mathrm{2.1}$ W m⁻² and ML${}^{**}$ at $-\mathrm{17.6}\pm \mathrm{2.2}$ W m⁻². Despite this consistency in weighted climatological mean, localized discrepancies emerge in regions characterized by strong spatial gradients. Specifially, SBDART simulates an intense, spatially confined cooling feature north-northwest of the Xuzhou_CUMT station, which both ML^* and ML${}^{**}$ reproduce only in a spatially diffused form, underestimating its magnitude and sharpness. This smoothing likely reflects a combination of spatial averaging, regularization, and error propagation inherent to ML surrogates, which preferentially preserve coherent regional patterns while attenuating localized extremes. Over Dhaka (Fig. 7b) differences between ML^* and ML${}^{**}$ become more pronounced. SBDART simulates a weighted climatological BCTOA mean of $-\mathrm{15.0}\pm \mathrm{1.2}$ W m⁻² closely matched by ML${}^{**}$ ($-\mathrm{14.9}\pm \mathrm{1.1}$ W m⁻²), whereas ML^* exhibits a weaker cooling ($-\mathrm{14.3}\pm \mathrm{1.1}$ W m⁻²). In terms of spatial coherence, ML${}^{**}$ demonstrates superior agreement with SBDART, particularly north of the Dhaka_University AERONET station, where SBDART resolves localized cooling pockets with sharp gradients. ML${}^{**}$ preserves both the spatial extent and intensity of these features, while ML^* underestimates their magnitude. Similarly, south of the station, SBDART indicates relatively warmer BCTOA forcing, which ML${}^{**}$ reproduces with high fidelity, whereas ML^* overestimates warming in these regions.

A shared limitation of both ML models is evident in their inability to fully reproduce a pronounced cooling pocket centered near 90.5–91° E and 22.5–23° N. This feature, prominent in SBDART simulations, is absent in ML^* and only feebly represented in ML${}^{**}$, suggesting that region-specific non-linear relationships between predictors and TOA forcing are not fully captured by cross-regional or combined training alone. Such discrepancies are consistent with a well-recognized trade-off in ML applications to atmospheric and environmental sciences; wherein spatial regularization improves overall coherence and reduces random noise but can smooth localized extremes and hotspots (Colin et al., 2018; Rasp et al., 2018; Reichstein et al., 2019; Soriano and Maxwell, 2023). Grid-by-grid spatial correlations further demonstrate that ML${}^{**}$ consistently outperforms ML^* across the majority of the domain in both regions (Fig. 7c and d), achieving statistically significant correlations (p<0.05) for more than 90 % of grid cells.

Having established that ML${}^{**}$ outperforms the cross-regional ML^* in reproducing climatological spatial structure, we next evaluate its robustness relative to fully trained, region-specific ML models (Fig. 8). This comparison assesses whether the combined model can generalize across diverse BC microphysical and emission regimes without sacrificing local predictive skill.

In Xuzhou (Fig. 8a), despite complex aerosol environments characterized by high BC column number density, more internally mixed particles, and varying BC core-to-total mass ratios, ML${}^{**}$ achieves lower absolute error than the local ML model in approximately 32 % of grid cells. These improvements are most evident over grids influenced by large industrial sources, where spatial variability in scattering enhancement and coating evolution is pronounced. While ML${}^{**}$ underperforms the local ML model in other parts of the domain, its overall accuracy remains robust (RMSE ∼ 1.1 W m⁻²; MAPE∼ 4.5 %). Over Dhaka (Fig. 8c), ML${}^{**}$ outperforms the local ML model over a larger fraction of the domain, achieving lower absolute error in 53.5 % of valid grid cells. Spatially aggregated metrics indicate comparable performance between ML and ML${}^{**}$, with low bias and RMSE (∼0.6 W m⁻²) and similar MAPE (∼2.8 %), demonstrating that the combined model retains strong predictive capability in a region dominated by distinct BC microphysical regimes.

This grid-level behaviour is further reflected in the domain-aggregated performance metrics (Fig. 8b and d). Over both domain RMSE, MAPE, and bias show that the ML${}^{**}$ performs comparably to the region-specific ML model in both Xuzhou and Dhaka, with no systematic degradation in any metric. In Xuzhou, ML${}^{**}$ exhibits slightly higher bias than the local ML model, consistent with its underperformance in portions of the domain, while near similar RMSE and MAPE. In Dhaka, the two models show nearly identical normalized metrics, confirming that the ML${}^{**}$ retains robust predictive skill while improving performance over a larger fraction of grid cells.

Overall, these results indicate that combined training enables ML${}^{**}$ to capture substantial spatial variability and retain predictive skill across regions with contrasting BC microphysical-optical characteristics. While the combined model does not uniformly outperform region-specific surrogates, it achieves comparable or superior performance over substantial portions of both domains, supporting its utility as a unified framework for transferable BC radiative forcing estimation. At the same time, the results underscore an inherent trade-off between generalization and regional specificity; blending heterogeneous distributions improves robustness to domain shifts but may partially dilute localized interactions that are optimally captured by region-specific models.

3.6 Diagnosing transferability limits: zero-shot evaluation on complex urban and biomass burning regions and improving ML${}^{**}$

While the preceding analysis demonstrates ML${}^{**}$ robust performance across Xuzhou and Dhaka, the question of transferability to regimes with stronger biomass burning influence remains unanswered. To address this, ML${}^{**}$ was in a zero-shot manner (i.e., without any retraining or fine-tuning) on two additional regions with contrasting dominant aerosol types. Delhi, India during the post-monsoon crop residue burning period (October–November 2024), representing an urban-biomass burning mixed regime and Mongu, Zambia during the peak dry season (July–October 2024), a remote AERONET site representative of intense savanna fire emissions.

Figure 9 Zero-shot transferability evaluation of the combined ML${}^{**}$ model (trained exclusively on Xuzhou and Dhaka data) over two contrasting biomass burning-influenced regions.

Figure 9 presents the zero-shot evaluation results. ML${}^{**}$ transfers well to Delhi (Adj. R²=0.91, RMSE = 2.3 W m⁻², MAE = 1.9 W m⁻², MBE = −1.7 W m⁻²), demonstrating effective generalization to a heavily polluted urban-influenced environment during a period when the aerosol loading is significantly influenced by surrounding agricultural burning. Performance at Mongu, however, degrades substantially (Adj. R²=0.83, RMSE = 4.4 W m⁻², MAE = 4.2 W m⁻², MBE = −4.2 W m⁻²), with a systematic underprediction bias more than twice that observed at Delhi. The Adj. R² of 0.83 indicates that ML${}^{**}$ correctly tracks the relative variability in BC TOA DRF. However, the MBE of −4.2 W m⁻² reveals systematic underprediction bias.

To diagnose the cause of this degradation, feature-space overlap analysis was conducted comparing the distribution of key predictors between the training dataset (Xuzhou and Dhaka combined) and each test region. For Delhi, all predictors showed >97 % overlap with training ranges, with Cohen's d values < 0.30 and strong IQR overlap, confirming distributional compatibility. For Mongu, the overlap analysis reveals a systematic, multi-predictor distributional mismatch across the key predictors. Shell thickness shows the most severe divergence (60.7 % overlap, Cohen's $d=-\mathrm{2.83}$), followed by number concentration (80.8 % overlap, Cohen's $d=+\mathrm{1.18}$), BC core size (84.2 % overlap, Cohen's $d=-\mathrm{1.67}$), and total mass (88.3 % overlap, zero IQR overlap). The high M_c overlap is driven solely by distributional tails; the zero IQR overlap reveals that the bulk of Mongu's aerosol mass occupies a range completely disjoint from the training data's central tendency. The mixing state predictor shows a subtler but equally important diagnostic, despite 99.2 % overlap, Cohen's $d=+\mathrm{1.02}$ indicates a large systematic shift. This apparent contradiction is attributable to the fact that the same mixing state ratio can arise from different microphysical configurations. In Mongu, the higher mixing state reflects a small BC core (∼103 nm) with a thin coating (∼203 nm). In the training data, a similar ratio corresponds to a larger core (∼120 nm) with a proportionally thicker coating (∼480 nm). The ratio overlaps numerically, but the underlying particle architecture, and thus the radiative behaviour is physically distinct.

These microphysical distinctions are not merely statistical artefacts; they reflect genuinely different aerosol regimes that can be independently validated against in-situ observations. The COSMO-retrieved BC core median size for Mongu (∼103 nm) is smaller than the training dataset median used in ML${}^{**}$ (Xuzhou+Dhaka). This is consistent with savanna fire BC particles that have undergone limited coagulation and growth compared to relatively more aged urban BC characteristic of Xuzhou and Dhaka. In-situ SP2 measurements of African biomass burning smoke consistently report BC core count median diameters in the range of ∼100–200 nm, with values varying by transport age and fire type (Wu et al., 2021). The Mongu retrieval falls squarely in the middle of this range, as expected for a regional receptor site dominated by fresh to moderately aged savanna smoke. Similar physical consistency extends to the coating thickness. The thinner coatings retrieved at Mongu align with SP2-based measurements of West African savanna fire plumes reporting BC coating thicknesses change by ∼13–50 nm over the first 12 h of transport (Wu et al., 2021), reflecting limited secondary processing in remote savanna environments compared to urban BC. The lifecycle of biomass burning aerosol coatings further supports this, with rapid initial growth followed by slower accumulation or net loss with continued ageing (Sedlacek et al., 2022), keeping coatings thinner relative to urban pollution. The elevated column number concentration and low total aerosol mass retrieved at Mongu are equally physically consistent. Dense particle number concentrations are characteristic of optically thick southern African savanna smoke (Holanda et al., 2020; Wu et al., 2020), while the low total aerosol mass reflects the sulfate-poor nature of savanna fire BC, where sulfate mass fractions remain minimal in the absence of industrial sulfur emissions. Study by Wang et al. (2025) has shown that this area has a very high amount of NO_x emitted during this time, which would provide a way through co-emitted NO_x to rapidly form nitrate aerosols, which are optically identical to sulfate, and hence could fill the gap between the thinner shell sizes observed using COSMO and the even thinner shell sizes from the above observational studies.This indicates that the COSMO retrieved properties represent a realistic biomass burning regime rather than retrieval artifacts. Consequently, the feature-space divergence from the training distribution reflects a genuine aerosol type underrepresented in ML${}^{**}$ training, providing a physical explanation for the systematic underprediction at this site.

To improve transferability to biomass burning regimes, we retrained the model on an expanded dataset including all four regions - Xuzhou, Dhaka, Delhi, and Mongu using a region-stratified 70/30 split. The test set (30 % from each region, 124 691 samples total) was held out entirely during training to ensure rigorous evaluation on unseen data. Internal validation (80/20 split of training data) showed minimal train-validation gaps (R² difference < 0.005, RMSE difference < 0.05 W m⁻²), confirming no overfitting. This improved combined model is denoted by ML${}^{***}$, when applied over the unseen held-out test datasets achieved strong and consistent performance across all four regions (Fig. S7). For the urban sites, Xuzhou (N=31 670) yielded Adj. R²=0.95, RMSE = 2.0 W m⁻², and MBE = −0.2 W m⁻²; Dhaka (N=81 967) yielded Adj. R²=0.95, RMSE = 1.7 W m⁻², and MBE = +0.2 W m⁻². For the biomass burning-influenced site, Delhi (N=6776) achieved Adj. R²=0.94, RMSE = 1.6 W m⁻², and MBE = −1.1 W m⁻²; Mongu (N=4278) achieved Adj. R²=0.89, RMSE = 1.4 W m⁻², and MBE = −0.8 W m⁻². The near-identical performance between urban sites and the modestly lower but still robust performance at Mongu demonstrate that ML${}^{***}$ performs effectively across diverse aerosol regimes without degrading urban accuracy.

Finally, the ML${}^{***}$ was applied to the full datasets across all four regions and the performance against the ML${}^{**}$ model was compared (Table 2). This analysis aids in assessing whether adding biomass burning predictor regimes improve performance at Mongu and Delhi, while simultaneously preserve accuracy at Xuzhou and Dhaka.

Table 2 Performance comparison between the ML surrogates combined model ML${}^{**}$ and the revised four-region combined ML model (ML${}^{**}$) for predicting COSMO-RTM BC TOA DRF.

The results show that, urban performance remains virtually unchanged; Xuzhou and Dhaka show negligible differences, confirming that incorporating extreme biomass burning cases from Mongu and mixed urban-biomass burning cases from Delhi does not compromise the model's ability to capture urban BC–radiation interactions. ML${}^{***}$ substantially improves performance at the biomass burning-influenced sites. For Delhi, RMSE improves by 30 % and bias reduces by 35 %. For Mongu, RMSE decreases by 68 % and bias improves from −4.2 to −0.8 W m⁻². Consistent performance across region by ML${}^{***}$ demonstrates that the revised combined model transfers effectively across diverse aerosol regimes, from polluted urban environments with signatures of biomass burning (Xuzhou, Delhi, Dhaka) to remote savanna fire (Mongu) environments.

4 Conclusions

This study advances the COSMO framework by refining the characterization of BC microphysical and columnar properties derived in Tiwari et al. (2025) and integrating them into a large ensemble of physically consistent radiative transfer simulations. By coupling observation-constrained microphysics with more than 400 000 RTM runs, we generate high-resolution estimates of BC DRF that resolve both spatial and temporal variability beyond what is achievable with conventional mass-based approaches.

A central scientific finding is the emergence of grid-scale BC TOA DRF that spans both positive and negative values within the same regional domains. While positive BC TOA is widely reported, the occurrence of negative BC TOA highlights the importance of resolving microphysical and column-dependent absorption-scattering competition, which is inadequately represented in many current climate-model parameterization. This demonstrates that BC radiative effects are not monotonic with loading and that realistic aerosol–radiation interactions require explicit treatment of particle size, mixing state, and column abundance. The contrast between Xuzhou and Dhaka provides a compelling illustration of this complexity. Similar BCAOD₅₅₀ magnitudes mask fundamentally different absorption-scattering dynamics driven by differences in BC microphysics and column abundance, leading to divergent radiative outcomes and regional climate implications. These results underscore the limitations of mass-based or bulk optical approaches, including methods that infer PM_2.5 or radiative effects from AOD alone (Chen et al., 2020; Kang et al., 2025; Ma et al., 2022; Wei et al., 2023), which inherently overlook the microphysical processes influencing the BC absorption efficiency and mass absorption cross-section variability.

ML surrogates trained on the RTM ensemble outperforms linear and parametric approaches across a wide range of pollution and microphysical regimes. Despite modest bias-variance trade-offs in transitional BC TOA regimes, the ML framework retains robust predictive ability. SHAP analysis of the trained surrogate further indicates that model learns parameterization in which the predicted BC TOA DRF depends strongly on non-linear interactions among BCAOD₅₅₀, mixing state, and N_c, with no universal threshold separating cooling and warming regimes in the learned mapping. Instead, a transitional regime emerges in which scattering increases more rapidly than absorption at high aerosol loading, with region-dependent moderation and overall expression. In industrial Xuzhou, enhanced scattering dominates at high loadings consistent with rapid aging and internal mixing with co-emitted species including sulfate and nitrate, qualitatively illustrated by the pronounced directional haze contrast in Fig. S8. In contrast, in biomass-burning-dominated Dhaka, absorption remains comparatively strong consistent with larger, less-coated BC cores, consistent with observations from other combustion-intensive regions (Guan et al., 2026). These findings demonstrate that no universal parameterization of BC forcing is currently feasible and that BC climate impacts are inseparable from source-specific microphysics and local environmental conditions.

Transferability tests over Delhi and Mongu show that ML surrogate performance depends strongly on how well the target aerosol microphysical and optical regimes are represented within the training feature space. The original Xuzhou-Dhaka-trained combined model (ML${}^{**}$) reproduced Delhi reasonably well, indicating that geographically distinct regions can still be predictable when their BC microphysical and optical predictors remain within the learned domain. In contrast, the systematic underprediction over Mongu reflected a distinct savanna-fire regime characterized by smaller BC cores and thinner coatings outside of the main ranges of the training dataset, as well as higher column number abundance near the top of the range of the training dataset, and shifted predictor distributions relative to the ML${}^{**}$ training data. Incorporating Delhi and Mongu into the revised four-region model (ML${}^{***}$) substantially reduced this bias, particularly over Mongu, while preserving performance over Xuzhou and Dhaka, and slightly improving performance of Delhi. These results indicate that transferability is not universal, but is diagnosable and can be improved by expanding the RTM training ensemble to include underrepresented aerosol regimes. Therefore, surrogate deployment should be preceded by a quantitative feature-space compatibility test. Simple inclusion within the training min-max range is not sufficient, because distributional tails can mask systematic shifts in the central predictor space. Instead, target-region predictors should be compared with the training domain using complementary diagnostics: (i) the fraction of samples within the training 2nd–98th percentile range, (ii) Cohen's d effect size to identify systematic distributional displacement, and (iii) interquartile-range overlap to determine whether the bulk of the target distribution lies within the learned regime. This also demonstrates that the traditional approach of prescribing fixed optical and microphysical properties for urban and biomass burning aerosol types – which global models currently heavily rely upon (Brown et al., 2021; Chen et al., 2023; Das et al., 2024) is likely incompatible with the diversity of microphysical properties revealed within each type by observationally constrained retrievals.

While this framework captures key microphysical controls on BC TOA, several assumptions define its present scope. The use of a core-shell representation with a non-absorbing sulfate proxy reflects a deliberate optical constraint. Although urban aerosols arise from diverse emission sources, extensive field evidence across urban, industrial, and biomass-burning regions in East and South Asia indicates that BC is frequently internally mixed and coated by optically similar inorganic and organic material (Kompalli et al., 2021; Peng et al., 2016; Sun et al., 2024; Zamora et al., 2019), supporting the physical relevance of this assumption. However, particles with irregular shapes, including fractal aggregates, agglomerated structures, are known to occur prominently in freshly emitted BC, are not explicitly addressed in this study. Future extensions could reconcile this limitation through hybrid frameworks that couple fractal or irregular morphologies at emission with spherical coated particles at advanced aging stages. Next, the treatment of organic aerosols follows an optical rather than chemical constraint, with purely scattering organic carbon treated similarly to nitrate and sulfate, while the absorption associated with absorbing organic carbon (aOC) is not explicitly prescribed but is filtered through a strict multi-band SSA consistency requirement. As described in previous studies (Liu et al., 2024b, c), the multi-waveband optical consistency inherently favours particles exhibiting stable cross-band absorption characteristic of BC, while disfavouring wavelength-confined absorbers whose absorption is confined to shorter wavelengths, such as aOC (Flores et al., 2014; Forrister et al., 2015; Lambe et al., 2013; Park et al., 2010). Additional limitations arise from the conservative spatial representativity of merged AERONET-TROPOMI SSA constraints, where pixel exclusion prioritizes physical admissibility over spatial completeness; denser ground-based networks and expanded geostationary coverage offer a pathway to relax these constraints. Furthermore, vertical variability of BC remains unresolved, reflecting both methodological choices and the maturity of current observing systems. A companion CALIPSO-constrained sensitivity analysis over Xuzhou and Dhaka showed that the dominant aerosol layers are largely boundary-layer confined, and that redistributing high-BCAOD, low-SSA cases into CALIPSO-derived elevated layers produced only modest BC TOA DRF changes, with median relative differences below 1 % and 95th-percentile differences below approximately 2 % (Tiwari et al., 2026a). This suggests that the baseline vertical-profile assumption is unlikely to dominate the BC TOA DRF uncertainty over these two regions, although the retrieval remains column-integrated and does not resolve the true vertical profile of BC extinction. Emerging multispectral GRASP-based TROPOMI products (Chen et al., 2024; Litvinov et al., 2024), together with new global BC microphysical datasets (Liu et al., 2026), provide a direct route to extend this framework without substantially altering its core physical foundations.

Finally, this study computes clear-sky BCDRF at the TOA, surface, and within the atmosphere. However, the analysis presented here emphasizes BC TOA, because it represents the net clear-sky TOA perturbation associated with the observation-constrained BC-containing aerosol state. As the primary contributor to the instantaneous radiative energy imbalance caused by BC, BC TOA provides an unambiguous metric for assessing BC-specific climate forcing and feedbacks under clear-sky conditions. Follow-up work will extend this framework to jointly analyse the spatio-temporal variability of TOA, surface, and net atmospheric BC DRF components and provide a more complete assessment of BC radiative impacts over these contrasting urban regions.

Appendix A

Mie model

Mie theory provides a classical solution to Maxwell's equations describing the scattering of electromagnetic waves by homogeneous spherical particles, widely used to model aerosol scattering properties (Bohren and Huffman, 1998; Mie, 1908) such as extinction efficiency (Q_ext), scattering efficiency (Q_scat), single scattering albedo (SSA), and asymmetry factor (g). The Mie model used in Tiwari et al. (2023, 2025) characterizes BC aerosol particles as coated spheres, consisting of a black carbon (BC) core with sizes ranging from 50 to 500 nm and $m=\mathrm{2.0}+\mathrm{1.0}i$ (Schuster, 2005), coated by a non-absorbing sulfate (Sul) shell sized between 10 and 800 nm and $m=\mathrm{1.52}-\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}i$ (Aouizerats et al., 2010). This core-shell configuration allows us to analyze scattering properties at the single-particle level while varying particle size and mixing state, constrained qualitatively by multi-waveband observational data and their uncertainties.

While the Mie model's assumption of spherical, core-shell particles does not capture the full complexity of aerosol morphology – such as agglomerated or irregularly shaped particles – numerous studies have demonstrated that this approach remains representative for BC aerosols in regions dominated by urban, industrial, and biomass burning sources, particularly across South Asia and China (Kompalli et al., 2021; Liu et al., 2020b; Peng et al., 2016; Wang et al., 2021; Zamora et al., 2019). The core-shell assumption simplifies the optical modeling while effectively capturing key interactions between incident radiation and particle morphology relevant to these environments.

Tiwari et al. (2025) derived size, mixing state, column mass and number density over Xuzhou, China and Dhaka, Bangladesh using this morphological setup adopting sulfate as a proxy for the non-absorbing shell material due to its optical similarity to other common scattering aerosol coatings such as nitrate and ammonium. Literature shows that differences in extinction efficiency and SSA among these species are minimal, typically within 2 %, and further reduced when mixed with aerosol water (Dai et al., 2024). Another study demonstrated that ammonium sulfate's refractive index remains nearly constant under varying humidity, with negligible absorption changes (Erlick et al., 2011). Thus, categorizing sulfate, nitrate, ammonium, and similar species collectively as scattering aerosols is justified, based on their low absorption characteristics, supporting the use of sulfate as an effective scattering shell proxy in the used core-shell model. Furthermore, the exclusion of secondary organic aerosol (SOA) particularly absorbing organic carbon (aOC) from analysis is grounded within the methodological screening. Our approach emphasizes identifying core-shell size combinations that exhibit absorption consistent with all observations made across the UV, VIS, and NIR bands in tandem, a trait commonly associated with aerosols like black carbon (Bond and Bergstrom, 2006; Samset et al., 2018a). SOA and aOC exhibit strong UV absorption but diminishing visible and near-infrared absorption, leading to instability in SSA across wavebands and failure to meet the filtering criteria derived from observations. Furthermore, aOC undergoes rapid photochemical bleaching and mass absorption cross-section reduction within hours to a single diurnal cycle after emission, transitioning from absorbing to primarily scattering behavior (Forrister et al., 2015; Lambe et al., 2013; McMeeking et al., 2011; Park et al., 2010; Saleh et al., 2013). Sensitivity tests incorporating aOC within the shell layer of BC cores was conducted by Liu et al. (2024b), revealing lower aerosol absorption retrievals and reduces the number of valid pixels raising question on the physical realism of aOC absorption in aged aerosols, particularly in urban environments where chemical aging accelerates. Consequently, our model prioritizes particles with, cross-band absorption characteristics consistent with BC and sulfate, effectively excluding SOA from the core-shell assumption. This method prioritizes consistency with BC aerosol type and variable size and mixing state, potentially excluding some cases with similar loading but differing optical characteristics. Although conservative, this exclusion enhances the reliability of the inverted results ascertaining the pixel for BC. Future expansion of observational networks and studies will be essential to refine the selection of optimal solutions from the probabilistic solution developed using this approach, improving the accuracy of BC concentration monitoring in varied urban settings and deepening insight into emissions amid dynamic urban-industrial changes.

Appendix B

B1 Satellite and reanalysis dataset processing for SBDART model simulation

B1.1 Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm-based AOD

MODIS aerosol products, widely utilized for regional and global analyses (Tian et al., 2023; Wei et al., 2020), are derived primarily from three algorithms: Dark Target (DT), Deep Blue (DB), and MAIAC (Qin et al., 2021). The DT algorithm is limited by its limited ability to retrieve aerosol optical properties over bright surfaces such as urban areas, deserts, and the Gobi region. The DB algorithm, meanwhile, assumes that short-wave infrared bands, that exhibit minimal reflectance influence from vegetation (Hsu et al., 2013; Sayer et al., 2013), which can cause errors in aerosol optical depth (AOD) retrievals, especially under polluted conditions.

The MAIAC algorithm retrieves AOD over land by applying time series analysis to distinguish aerosol signals from surface reflectance, while explicitly accounting for the bidirectional reflectance distribution function (BRDF) of the terrain (Lyapustin et al., 2011). By integrating temporal and spatial image processing, MAIAC effectively corrects atmospheric effects, enabling reliable AOD retrievals even over bright surfaces and densely vegetated areas (Lyapustin et al., 2011, 2018). Compared to the DT and DB algorithms, MAIAC offers enhanced spatial coverage with a spatial resolution of 1 km × 1 km daily data. This study employs MCD19A2 AOD data at 470 and 550 nm wavelengths, accessed from the https://search.earthdata.nasa.gov/search/granules?p=C2324689816-LPCLOUD&pg[0][v]=f&pg[0][gsk]=-start_date&q=MCD19A2 (last access: 27 June 2026). To ensure data integrity, pixels were rigorously filtered following the MAIAC Collection 6.1 quality assurance protocols, including cloud masking, adjacency effect screening, and AOD quality flags. Only pixels flagged as “Clear”, “Normal Condition/Clear”, and “Best Quality with No glint” were retained for analysis. MAIAC AOD products have been extensively validated at both regional and global scales, including over South Asia, demonstrating superior accuracy compared to DT and DB algorithms (Mhawish et al., 2019). Global study has reported strong agreement with ground-based AERONET observations, with correlation coefficients exceeding 0.8 at over two-thirds of monitoring sites (Qin et al., 2021).

B1.2 Total column precipitable water vapor from MODIS (MYD05 and MOD05)

This study predominantly utilizes the MODIS PWV Level 2 product MYD05_L2, with MOD05_L2 data employed as a substitute when MYD05_L2 is unavailable. The MODIS PWV is derived from the ratio of three water vapor absorbing bands (0.905, 0.936, and 0.930 µm) to two atmospheric window bands (0.865 and 1.240 µm) (King et al., 1992). While both IR and NIR channels provide PWV information, this study opts for the NIR product due to its lower sensitivity to atmospheric temperature compared to atmospheric water vapor content (Bai et al., 2021; Zhou and Cheng, 2025). The NIR product offers a daily data at a spatial resolution of 1×1 km.

Given the sensitivity of MODIS PWV to cloud presence in the field of view (He and Liu, 2019), a series of rigorous quality assurance measures were implemented to retain only the most reliable pixels. The study explicitly checks for Cloud Mask QA flags, retaining only those pixels with “Cloud Mask Status Flag” marked as “determined” and “Cloudiness Flag” as “Confident clear”. Furthermore, each pixel undergoes scrutiny for product quality and retrieval processing flags. Only pixels with PWV (NIR) “Usefulness Flag” designated as “Useful” and “Confidence Flag” as “Good” or “Very Good confidence” are incorporated into the study. For a comprehensive understanding of the physical background of the retrieval method, the retrieval algorithm, adapted inversion technique, uncertainty estimates, and validation, readers are directed to the detailed study conducted by He and Liu (2019). Since the MOD05/MYD05_L2 data only has geolocation scientific data set for PWV(IR) at 5×5 km, this study separately downloaded MYD03_L2/MOD03_L2 1 km Geolocation data and collocated MYD05 and MYD03. Both MOD/MYD03_L2 and MOD/MYD05_L2 datasets utilized in this study are publicly accessible via the NASA's Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center (LAADS DAAC) (URL: https://ladsweb.modaps.eosdis.nasa.gov/search/order/2/MYD03--0,MYD05_L2--0, last access: 27 June 2026).

B1.3 Surface spectral albedo from MODIS BRDF albedo model (MCD43A1)

In this study, we utilize the MCD43A1 dataset, which includes MODIS, Terra and Aqua BRDF/Albedo model parameters on a daily basis at a 500 m resolution, to compute the actual blue-sky albedo. This measured albedo is employed in the radiative transfer model as the surface albedo, replacing the predefined values used by the model.

The MODIS BRDF/albedo algorithm, based on the three-parameter Ross-Thick-Li-Sparse-Reciprocal model, estimates kernel parameters from clear-sky surface reflectance collected by Terra and Aqua MODIS sensors (Lucht et al., 2000; Román et al., 2010; Schaaf et al., 2002). These parameters include the isotropic parameter (f_iso), the kernel for volume scattering (K_vol), and the kernel derived from surface scattering and geometric shadow theory (K_geo). They are available for each MODIS spectral band and for three broad bands: 0.3–0.7, 0.7–5.0, and 0.3–5.0 µm. These parameters can be utilized in the forward model to reconstruct surface anisotropic effects, allowing for the correction of directional reflectance. From the MCD43A1.061 files, layers of BRDF Albedo Parameters for visible, near infrared and shortwave bands along with their respective albedo band mandatory quality flags are used. The weighting parameters described above are used to compute integrated directional hemispherical reflectance also known as Black-sky Albedo (albedo in the absence of a diffuse component and is a function of solar zenith angle) at local solar zenith angle, and bihemispherical reflectance or White-sky albedo (albedo measured without direct sunlight when the diffuse component is isotropic). Together, these two types of albedos represent the extremes of completely direct and completely diffuse illumination. Finally, the actual blue-sky albedo is derived as an interpolated value between these extremes (Jia et al., 2022), depending upon the fraction of diffuse skylight present. This fraction, itself a function of the solar optical depth (Giannaklis et al., 2023). For each short-wave, visible and near-infrared band, we compute the actual blue-sky albedo using the above method and finally spectral average of the blue-sky albedo is taken across these wavebands which serves as the input for radiative transfer model (RTM). This approach aligns with previous methods, where surface albedo values are averaged across the four operational AERONET wavebands. This spectral averaging has been demonstrated to provide a reliable approximation for estimating direct aerosol radiative forcing at the surface (García et al., 2012; Logothetis et al., 2021).

The fundamental equations for calculating black-sky and white-sky albedo using the BRDF kernels, as well as the global assessment/validation of these albedo products, have been extensively covered in earlier literature assuring the quality of the product (Che et al., 2017; Wang et al., 2015). NASA's Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) platform, offers efficient tool for extracting the required layers of data for desired temporal and spatial range. The MCD43A1.061 data can be sourced from https://appeears.earthdatacloud.nasa.gov/task/area (last access: 27 June 2026).

B1.4 Total column ozone from Sentinel-5P TROPOMI (TROPOspheric monitoring instrument)

The total column ozone data used in this study were obtained from the Sentinel-5P TROPOMI (TROPOspheric Monitoring Instrument) S5P_L2__O3_TOT_HiR product, which provides data at a high spatial resolution of 5.5×3.5 km (Inness et al., 2019; Sullivan et al., 2022). This study uses dataset derived from offline and reprocessed products, which employs the GOME-type Direct FITting (GODFIT) algorithm. This advanced retrieval algorithm enhances the accuracy of total ozone measurements by directly fitting simulated radiances in the Huggins bands (325–335 nm) to observational data (Garane et al., 2019). The adjustment process involves varying key atmospheric parameters, including total ozone concentration, effective scene albedo, and effective temperature. To retain only the best pixel, for total column ozone, the study implements a pixel-by-pixel data quality assurance assuring that qa_value is greater than 0.5 following the recommendation of S5P Mission Performance Centre OFFL/RPRO Total ozone Readme Manual.

The S5P_L2__O3_TOT_HiR can be sourced from: https://disc.gsfc.nasa.gov/datasets/S5P_L2__O3_TOT_HiR_2/summary?keywords=Total column ozone (last access: 27 June 2026). This study used the spatio-temporal data subsetting tool provided by Goddard Earth Sciences Data and Information Services Centre (GES DISC) to source the data for desired spatio-temporal range. A comprehensive validation of this product, using both global ground-based and contemporary satellite observations, was conducted by Garane et al. (2019). In a more recent study, various factors – such as solar zenith angle, surface albedo, and measurement parameters like scan angle – and their impact on the uncertainty in TROPOMI's ozone profile retrievals were investigated (Keppens et al., 2024).