DDA-BNN v1.0: a morphology-aware surrogate model for the optical properties of black carbon-containing particles

Black carbon (BC) is the most strongly absorbing component of atmospheric aerosol and significantly impacts Earth's energy balance. The optical properties of BC-containing particles depend on particle-level variability in size, chemical composition, and internal morphology. Such particle-level details are not easily represented in large-scale atmospheric models. Existing parameterizations typically assume idealized particle geometries (e.g., homogeneous spheres or concentric core–shell spheres) and homogeneous mixing, which can yield biased predictions and provide no quantitative estimate of model-form uncertainty at the single-particle level. In this work, we present a probabilistic framework for predicting the optical properties of individual BC-containing particles using a hybrid Bayesian neural network (BNN) model trained on discrete dipole approximation (DDA) simulations. The hybrid BNN is a flexible combination of deterministic and Bayesian layers allowing for more realistic treatment of particle optical properties and quantification of uncertainty. The hybrid BNN model predicts extinction efficiency, single-scattering albedo, and asymmetry parameter and returns predictive uncertainty that can be decomposed into aleatoric (data-driven variability) and epistemic (uncertainty due to limited training coverage) components. We show that the hybrid BNN outperforms homogeneous-sphere and core–shell Mie approximations for calculating extinction and scattering-sensitive quantities (Q_ext, SSA and g), while maintaining comparable accuracy for absorption-related metrics. We further demonstrate how epistemic uncertainty highlights under-sampled regions of particle parameter space, enabling targeted design of future DDA simulations that most effectively reduce model uncertainty. This uncertainty-aware surrogate provides a practical pathway for incorporating realistically complex particle morphologies into parameterizations of aerosol optical properties, which will ultimately improve the reliability of model-based assessments of BC impacts on the atmosphere.

1 Introduction

Predicting the optical properties of black carbon (BC)-containing particles presents a unique and longstanding challenge in atmospheric models and remote sensing algorithms (Bond and Bergstrom, 2006). BC-containing particles represent a small fraction of atmospheric aerosol mass but contribute significantly to total light absorption (Bond et al., 2013), so modest errors in BC optics can translate into disproportionately large errors in aerosol direct radiative effect and inferred emissions. Despite decades of study, a central problem remains: particle-resolved morphology and mixing state strongly control absorption and scattering, yet large-scale models require computationally cheap optics schemes that reduce this complexity into a small set of bulk parameters.

BC-containing particles are distinct from many other aerosol types because of their strong and highly variable absorption properties (Samset et al., 2018; Laskin et al., 2015; Bond and Bergstrom, 2006; Cappa et al., 2012; Fierce et al., 2020; Liu et al., 2017b, 2015). BC-containing particles are emitted as fractal aggregates of primary spherules, and their morphology evolves after emission through coagulation and condensation of semi-volatile materials (China et al., 2013). Accumulation of coatings (or water during cloud processing) can restructure aggregates toward more compact configurations via surface tension and capillary forces (Corbin et al., 2023; Beeler et al., 2025), while the coating itself can range from weakly absorbing to strongly absorbing (“brown carbon”, BrC). Laboratory and ambient observations have shown that these combined changes in internal structure, coating amount, and coating absorptivity lead to broad variability in single-particle extinction, single-scattering albedo (SSA), and asymmetry parameter (g) (Hu et al., 2021; Liu et al., 2020; Bond et al., 2013).

To achieve computational efficiency, large-scale atmospheric and radiative transfer models typically represent BC-containing particles using simplified shapes (homogeneous spheres or core–shell spheres) and assume uniform composition within size modes (Liu et al., 2016). While these assumptions simplify calculations, they introduce substantial error by neglecting complex particle structures and compositional diversity (Fierce et al., 2016, 2020; Liu et al., 2015). Detailed models of single-particle light scattering, such as those based on the discrete dipole approximation (DDA), represent realistically complex particle morphologies and compositions (Yurkin and Hoekstra, 2011), but are too computationally expensive for routine use in regional or global modeling and are impractical for embedding directly into iterative workflows such as inverse modeling, data assimilation, or uncertainty propagation. Recent work has explored the use of machine learning (ML) to approximate optics calculations for BC-containing particles with much lower computational cost (Qin et al., 2025; Romshoo et al., 2024; Li et al., 2023; Li and May, 2022). However, these studies rely on multiple-sphere T-matrix models, which are more computationally efficient than DDA, but generally require idealized representations (e.g., non-overlap constraints and idealized coatings) (Liu et al., 2017a). Additionally, most existing ML surrogates return only point predictions and do not separate uncertainty due to irreducible variability in the training data from uncertainty due to sparse coverage of particle parameter space. This is a key limitation because the relevant BC parameter space is high-dimensional and observational constraints remain incomplete.

As a path forward, this work introduces a hybrid Bayesian Neural Network (BNN) approach (Bonnet et al., 2023; Olivier et al., 2021; Gal and Ghahramani, 2016) to predict the absorption and scattering properties of individual BC-containing particles, as well as confidence intervals in modeled quantities. The hybrid BNN combines deterministic and Bayesian layers to balance computational efficiency with uncertainty quantification. The model is trained on single-particle optics simulations that have been informed by detailed laboratory measurements of the mixing state-dependent morphology of BC-containing particles (Beeler et al., 2025). The hybrid BNN predicts the extinction efficiency, single-scattering albedo (SSA), and asymmetry parameter (g) and associated uncertainties for individual particles, allowing for the decomposition of total model uncertainty into aleatoric and epistemic components. We compare the performance of the hybrid BNN against homogeneous-sphere and core–shell Mie approximations commonly used in atmospheric models. Finally, we demonstrate how epistemic uncertainty maps can be used to design targeted future DDA experiments that reduce surrogate uncertainty most efficiently. This hybrid BNN lays the groundwork for an uncertainty-aware, computationally efficient surrogate to enable more complete representation of BC optical properties in larger-scale models.

2 Training Set Generation

The training dataset comprises two ensembles of simulated BC-containing particles designed to span contrasting representations of the coupling between particle coatings and aggregate morphology. Ensemble I (“non-compacting”) is taken from Beeler and Chakrabarty (2022) and intentionally decouples coating amount from aggregate compaction, representing scenarios where coating does not induce restructuring. Ensemble II (“compacting”) includes particles whose compaction depends on coating amount following laboratory-informed scaling relationships (Beeler et al., 2025). These particles represent scenarios where liquid coatings induce restructuring via surface tension and capillary forces (Corbin et al., 2023). Together, the two ensembles allow the surrogate to learn both morphology-driven and coating-driven variability while retaining physical plausibility for compacting particles. The following sections provide details of the particles in each ensemble and calculation of their optical properties.

2.1 Ensemble I: Non-Compacting Black Carbon-Containing Particles

Non-compacting BC-containing particles are taken from Beeler and Chakrabarty (2022), and full details on their generation can be found therein. Briefly, these particles have three morphologies (fresh, partially collapsed, fully collapsed), with each morphology ranging from uncoated to thickly coated (Fig. 1). These particles are representative of BC-containing particles with decoupled morphologies and coating thicknesses. Physically, decoupled morphologies and coating thicknesses can arise if coatings evaporate from BC aggregates, coating is added via coagulation, or coating accumulates under high supersaturation conditions (Corbin et al., 2023). External coatings are quantified using the ratio of total particle volume to BC volume, also known as the volume growth ratio ($V/V_{\mathrm{0}}$). The amount of BC in each particle is quantified by the number of primary particles in the aggregate (N_pp). All particles in this group are made up of primary particles with radii of 20 nm. The degree of compaction is quantified using the core fractal dimension (D_f,c), with fresh aggregates having low D_f,c (1.8) and compact aggregates having high D_f,c (>2.5) (Koeylue et al., 1995). The D_f,c of each aggregate is calculated from the structure factor. The structure factor, which governs the angular dependence of scattered light, is the Fourier transform of the density autocorrelation function. The structure factor exhibits an inverse power-law dependence as a function of the scattering vector in the fractal regime, with an exponent equal to D_f,c (Sorensen, 2001).

Figure 1 Examples of simulated particles used to train the hybrid Bayesian neural network model. The training dataset is made up of two ensembles of particles (I and II). Particles in Ensemble I (top row) are taken from a previous publication (Beeler and Chakrabarty, 2022) and have statistically independent ($⟂⟂$) total-to-black carbon volume ratio ($V/V_{\mathrm{0}}$) and core fractal dimension (D_f,c). The particles in Ensemble II (bottom row) represent compacting particles and have $V/V_{\mathrm{0}}$, where the coating amount is selected to reflect the relationship between D_fc and $V/V_{\mathrm{0}}$ described in Beeler et al. (2025). Therefore, $V/V_{\mathrm{0}}$ and D_f,c are not statistically independent for these particles (D_f,c $\cancel{⟂}⟂V/V_{\mathrm{0}}$).

It is important to note that D_f,c and $V/V_{\mathrm{0}}$ are statistically independent in Ensemble I. We simulated 900 particles in Ensemble I, with N_pp ranging from 13–1189, D_f,c ranging from 1.8–3.0, and $V/V_{\mathrm{0}}$ ranging from 1–75 (see Table 1). The optical properties are calculated at wavelengths (λ) of 405, 532, 880, 1200, and 2500 nm. The refractive index of BC is assumed constant at 1.95+0.79i, and the coating refractive index is $\mathrm{1.55}+[\mathrm{0.00},\mathrm{0.01},\mathrm{0.05}]i$ (Bond and Bergstrom, 2006).

Table 1 Overview of discrete dipole approximation (DDA) simulations used to train the hybrid Bayesian neural network. The training dataset is made up of two ensembles. Particles in Ensemble I are taken from a previous publication and have independent total-to-black carbon volume ratio ($V/V_{\mathrm{0}}$) and core fractal dimension (D_f,c), while these are coupled for particles in Ensemble II. Other parameters used as inputs to the DDA model include the number of primary particles in each aggregate (N_pp), the radius of primary particles (r_pp), the wavelength of incident light (λ), the refractive index of BC (RI_BC), the coating refractive index at 550 nm (RI_{coat,550 nm}), and the absorption Ångstrom exponent of the coating (AAE_coat).

2.2 Ensemble II: Compacting Black Carbon-Containing Particles

Previous studies have developed scaling laws for D_f,c as a function of $V/V_{\mathrm{0}}$ using detailed laboratory measurements and single-particle simulations (Beeler et al., 2025). Therefore, the D_f,c of particles in this ensemble is dependent on $V/V_{\mathrm{0}}$ (following the formulation of Beeler et al., 2025). The N_pp of compacting particles is evenly sampled in log-space from 20 to 1000, and the $V/V_{\mathrm{0}}$ of each particle is randomly selected between 1 and 100 (see Table 1). The D_f,c of a particle for a given $V/V_{\mathrm{0}}$ is given by

$$\begin{array}{}\text{(1)}& D_{\mathrm{f},\mathrm{c}}-D_{\mathrm{f},\mathrm{0}}=c_{\mathrm{1}}\mathit{\{}\exp \left[c_{\mathrm{2}}\left(x_c-x_s\right)\right]-\exp \left[c_{\mathrm{2}}\left(x_{\mathrm{0}}-x_s\right)\right]\mathit{\}},\end{array}$$

where x_c is given by

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{x_c-x_{\mathrm{f}}}{x_{\mathrm{0}}-x_{\mathrm{f}}}}=\exp \left[-\mathrm{0.799}\left({\displaystyle \frac{V}{V_{\mathrm{0}}}}-\mathrm{1}\right)\right],\end{array}$$

and x_f is given by

$$\begin{array}{}\text{(3)}& x_{\mathrm{f}}=\mathrm{1.07}{x}_{\mathrm{0}}{N}_{\mathrm{pp}}^{-\mathrm{0.026}}.\end{array}$$

Here, $D_{\mathrm{f},\mathrm{0}}=\mathrm{1.8}$, x₀=0.46, $x_s=\mathrm{1}/\mathrm{3}$, c₁=11.66, and $c_{\mathrm{2}}=-\mathrm{41.39}$. BC-containing particles with given N_pp are generated using diffusion-limited cluster-cluster aggregation; compaction to the desired D_f,c (Eq. 1) is simulated using a spring-van der Waals compaction model; and coating to the desired $V/V_{\mathrm{0}}$ is simulated using a capillary coating model. These methods are described in detail in Beeler et al. (2025). The D_f,c of each particle in Ensemble II is then calculated using its structure factor.

Ensemble II included 990 particles that are representative of atmospheric BC-containing particles that have undergone coating-induced compaction. Half of the particles in Ensemble II are coated with non-absorbing coatings and have a coating refractive index of 1.6+0.0i. The remaining half of Ensemble II particles are coated with brown carbon (BrC). The imaginary part of the BrC refractive index (κ) at the wavelength of interest is given by

$$\begin{array}{}\text{(4)}& \mathit{\kappa }\left(\mathit{\lambda }\right)={\mathit{\kappa }}_{\mathrm{550}\mathrm{nm}}{\left({\displaystyle \frac{\mathit{\lambda }}{\mathrm{550}\mathrm{nm}}}\right)}^{-{\mathrm{AAE}}_{\mathrm{coat}}}.\end{array}$$

Here, κ_550 nm is the imaginary part of the refractive index at 550 nm and AAE_coat is the coating absorption Ångstrom exponent. Previous studies have shown that there is a relationship between κ_550 nm and AAE_coat, given by

$$\begin{array}{}\text{(5)}& {\mathrm{AAE}}_{\mathrm{coat}}=\mathrm{7.6}\exp \left[-\mathrm{12.4}{\mathit{\kappa }}_{\mathrm{550}\mathrm{nm}}\right].\end{array}$$

This relationship represents the measured properties of BrC across absorptivity (Saleh et al., 2018). The optical properties of Ensemble II particles are calculated at wavelengths randomly selected in 50 nm-spaced bins between 300 and 1200 nm. The refractive index of BC is assumed constant at 1.95+0.79i, the radius of primary particles is 20 nm, κ_550 nm is randomly sampled between 0.0 and 0.25, and the coating refractive index at the chosen wavelength is given by Eqs. (4) and (5).

2.3 Optical Properties of Simulated Black Carbon-Containing Particles

We use the ADDA code to calculate the optical properties of each particle in the training set (Yurkin and Hoekstra, 2011). This model operates by breaking each aggregate into sub-volumes that are much smaller than the wavelength of incident light (dipoles). Each dipole is then treated as a point scatterer that interacts with surrounding dipoles. This method has been shown to accurately calculate the optical properties of arbitrarily shaped objects and is highly accurate in capturing near-field interactions among the dipoles present in touching BC monomers. A full description of the ADDA configuration can be found in Appendix A. The optical properties of particles in Ensemble I are calculated at 405, 532, 880, 1200, or 2500 nm. We calculate the optical properties of Ensemble II particles at wavelengths randomly distributed among 50 nm-spaced bins ranging from 300 to 1200 nm.

3 Hybrid Bayesian Neural Network Description

Accurate prediction of the optical properties of BC-containing particles is of the utmost importance for determining their behavior in and effects on the atmosphere. However, BC-containing particles have variable compositions and morphologies, both of which have a large but uncertain effect on their optical properties. Here, we train a hybrid BNN model to infer the optical properties of BC-containing particles coated by absorbing and non-absorbing materials, at different stages of compaction. This model is unique from other similar models (such as traditional neural networks) as it enables uncertainty quantification of the inferred optical properties. Uncertainty quantified by the hybrid BNN can further be separated into aleatoric (irreducible noise inherent to the training dataset) and epistemic (error due to limited data). Uncertainty quantification in single-particle optical properties can then be propagated to particle populations, allowing for more complete representation of BC optical properties in atmospheric models. The hybrid BNN framework developed here is a sequence of deterministic and Bayesian layers trained on the dataset detailed in Sect. 2. The following sections outline data pre-processing routines, model architecture, uncertainty decomposition, and model evaluation.

3.1 Data Pre-Processing

Because several optical properties of BC-containing particles are correlated with one another (such as absorption and extinction efficiencies), we have chosen the extinction efficiency (Q_ext), single-scattering albedo (SSA), and asymmetry parameter (g) as the target variables of the hybrid BNN model. These optical properties are also commonly supplied to radiative transfer algorithms. Similarly, several properties of BC-containing particles are also correlated with one another, such as volume-equivalent diameter and particle mass. To avoid potential correlations, we have selected five microphysical properties of BC-containing particles as input variables to the hybrid BNN model: N_pp, $V/V_{\mathrm{0}}$, D_f,c, the volume-equivalent size parameter (χ_ve), and the imaginary part of the coating refractive index (κ_coat). These inputs are chosen as they are either directly calculated or easily derived by large-scale atmospheric and radiative transfer models. Additionally, the collection of these parameters fully describes the amount of BC (N_pp), the amount of coating ($V/V_{\mathrm{0}}$), the shape of particles (D_f,c), the coating material (κ_coat), and the wavelength of light (χ_ve). The optical properties of a volume-equivalent homogeneous sphere are also used as input variables in the hybrid BNN (Q_ext,hs, SSA_hs, and g_hs). The optical properties of a homogeneous volume-equivalent sphere are calculated using Mie theory, with refractive index given by the volume-weighted average of the BC and coating refractive indices (Sumlin et al., 2018). These were used as input variables because they are strong predictors of the target variables and are currently calculated in large-scale models.

Prior to training, each component of the target vector ${\mathit{y}}_{\mathrm{t}}=\left(Q_{\mathrm{ext}},\mathrm{SSA},g\right)$ is mapped to an unbounded real domain using a transformation T that preserves physical bounds (i.e., $T\left({\mathit{y}}_{\mathrm{t}}\right)=\left(T\left(Q_{\mathrm{ext}}\right),T\left(\mathrm{SSA}\right),T\left(g\right)\right)$). For strictly positive targets (e.g., optical efficiencies), we use a logarithmic transform. For bounded targets in (0,1) (SSA and g) (Koshelev, 2015), we use the logit transform. While the theoretical limits of g are $(-\mathrm{1},\mathrm{1})$, it is extremely rare for atmospheric particles to have g<0 (Zhao et al., 2018). Therefore, we use the logit transform for g. Specifically, T is given by:

$$\begin{array}{}\text{(6)}& T\left(y\right)=\left\{\begin{array}{ll}\ln \left(y\right),& y\in (\mathrm{0},\mathrm{\infty })\\ \text{logit}\left(y\right)=\ln \left(\frac{y}{\mathrm{1}-y}\right),& y\in (\mathrm{0},\mathrm{1})\end{array}\right.\end{array}$$

These transformations preserve constraints of the underlying physical variables while allowing the model to operate in an unconstrained latent space to ensure appropriately bounded uncertainties. After T transformations have been applied, all input variables (x) defined above and transformed target variables (y_t) are independently standardized using Z-score normalization (Fei et al., 2021) to give inputs and target variables in unconstrained latent space ($\stackrel{\mathrm{̃}}{\mathit{x}}$ and ${\stackrel{\mathrm{̃}}{\mathit{y}}}_{\mathrm{t}}$, respectively). More details on Z-score normalization can be found in Appendix B.

3.2 Model Architecture

The proposed neural network (NN_θ) is made up of an arbitrary sequence of deterministic (D) or Bayesian (B) layers (Kendall and Gal, 2017; Gal and Ghahramani, 2016). For example:

$$\begin{array}{}\text{(7)}& {\mathrm{NN}}_{\mathit{\theta }}=\left[D,B,D\right],\end{array}$$

defines a three-layer network with a deterministic input layer, a Bayesian hidden layer, and a deterministic output layer. Each layer (z) is generally defined as

$$\begin{array}{}\text{(8)}& \mathit{z}\mapsto \mathit{\varphi }\left(\mathbf{W}\mathit{z}+\mathit{b}\right),\end{array}$$

where W and b are the model weights and biases, and ϕ is the activation function. In a deterministic layer, $\mathit{\theta }\in \left(\mathbf{W},\mathit{b}\right)$ are the fixed parameters learned during training. Conversely, in a Bayesian layer the weights and biases are treated as random variables with an approximate (variational) posterior distribution, e.g., $q_{\mathit{\varphi }}\left(\mathit{\theta }\right)=\mathcal{N}({\mathit{\mu }}_{\mathit{\theta }},{\mathit{\sigma }}_{\mathit{\theta }}^{\mathrm{2}})$, parameterized by a learnable mean μ_θ and scale (standard deviation) σ_θ, from which θ is sampled during training and inference. We model each Bayesian-layer parameter independently, so this approximation does not represent correlations between parameters.

The outputs of the model, for a given set of sampled weights and biases θ and latent-space inputs $\stackrel{\mathrm{̃}}{\mathit{x}}$, are a predicted mean vector ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)$ and a vector of per-target aleatoric standard deviations ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)$ in latent space. These outputs are mapped back to physical units via the inverse Z and T transforms (Sect. 3.1) using Monte Carlo sampling (Sect. 3.3). We assume a diagonal aleatoric covariance (independent targets) and leave prediction of the full aleatoric covariance matrix for future work. The proposed pipeline is visualized in detail in Fig. 2.

Figure 2 Example of the hybrid Bayesian neural network architecture. Inputs are standardized using Z-score normalization, and target variables are mapped to an unbounded latent space using the transform T followed by Z-score normalization. The neural network (NN_θ) consists of a combination of deterministic and Bayesian layers producing a latent-space mean (${\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}$), aleatoric uncertainty (${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\theta }}$), and epistemic uncertainty (${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{e},\mathit{\theta }}$). Predictions are propagated to physical space by sampling in latent space and then applying the inverse transforms (Z⁻¹ and T⁻¹) to obtain optical properties in physical units.

3.3 Uncertainty Decomposition

During inference we seek both a point prediction and a decomposition of predictive uncertainty into aleatoric ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a}}$ and epistemic ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{e}}$ components. Suppose two observers are trying to predict the outcome of the next flip of an unknown coin. The first observer has recorded the results of 1000 previous flips and estimates that the probability of heads is 0.6. Their remaining uncertainty about the next flip is dominated by aleatoric uncertainty; even with perfect knowledge of the coin's bias, each flip is still random. The second observer has not seen any previous flips and does not know whether the coin is fair or biased. Their uncertainty includes epistemic uncertainty arising from their lack of knowledge about the coin, in addition to the aleatoric uncertainty inherent in the flipping process. Therefore, aleatoric uncertainty represents irreducible variability in the data (e.g., unresolved physical variability and numerical noise in the training set), whereas epistemic uncertainty reflects uncertainty in the model parameters due to limited data and/or model capacity. The initial aleatoric and epistemic decompositions are computed in latent space under Gaussian assumptions (Sect. 3.3.1). We then calibrate the uncertainty (Sect. 3.3.2) before propagating to physical space by sampling and applying the inverse transforms, from which we report predictive quantiles in physical units and estimate physical-space aleatoric and epistemic standard deviations via a nested Monte Carlo approximation (Sect. 3.3.3).

3.3.1 Latent-Space Predictive Distribution

At inference time we draw S independent weight and bias samples θ^(s)∼q_ϕ(θ) from the learned approximate posterior. For each draw, the neural network NN_θ directly outputs ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right):={\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{{\mathit{\theta }}^{\left(\mathrm{s}\right)}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a}}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right):={\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},{\mathit{\theta }}^{\left(\mathrm{s}\right)}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)$. The latent-space predictive mean is estimated by Monte Carlo averaging,

$$\begin{array}{}\text{(9)}& {\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)={\displaystyle \frac{\mathrm{1}}{S}}\sum _{s=\mathrm{1}}^{\mathrm{S}}{\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right).\end{array}$$

Aleatoric uncertainty in latent space is estimated as the mean of the predicted variances,

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\theta }}^{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)={\displaystyle \frac{\mathrm{1}}{S}}\sum _{s=\mathrm{1}}^{\mathrm{S}}{\left[{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a}}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)\right]}^{\odot \mathrm{2}},\end{array}$$

and epistemic uncertainty is estimated from the variations of the sampled means across draws,

$$\begin{array}{}\text{(11)}& {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{e},\mathit{\theta }}^{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)={\displaystyle \frac{\mathrm{1}}{S}}\sum _{s=\mathrm{1}}^S{\left[{\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)-{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)\right]}^{\odot \mathrm{2}}.\end{array}$$

Finally, the total latent-space predictive uncertainty is taken as the sum of epistemic and aleatoric variances,

$$\begin{array}{}\text{(12)}& {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{tot},\mathit{\theta }}^{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)={\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\theta }}^{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{e},\mathit{\theta }}^{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right).\end{array}$$

3.3.2 Uncertainty Scaling

Although Eqs. (10)–(11) provide an aleatoric-epistemic decomposition in latent space, the overall predictive uncertainty can be miscalibrated due to non-linear activation functions, leading to prediction intervals that are systematically too wide or narrow. We therefore apply per-target uncertainty scaling in latent space, using multiplicative factors τ fitted on a validation set (Guo et al., 2017; Kuleshov et al., 2018).

We fit τ by minimizing the Gaussian negative log-likelihood in latent space, using the latent predictive mean $\left({\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\right(\stackrel{\mathrm{̃}}{\mathit{x}}\left)\right)$, the total latent predictive standard deviation $\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{tot},\mathit{\theta }}\right(\stackrel{\mathrm{̃}}{\mathit{x}}\left)\right)$, and the validation targets transformed into latent space $({\stackrel{\mathrm{̃}}{\mathit{y}}}_{\mathrm{t}}=ZT({\mathit{y}}_{\mathrm{t}}\left)\right)$. For each target dimension j, we solve

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathit{\tau }}_j& =\arg \underset{\mathit{\tau }>\mathrm{0}}{min}{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\mathrm{NLL}\left({\stackrel{\mathrm{̃}}{y}}_{\mathrm{t},j}\left({\mathit{x}}_i\right);{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta },j}\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_i\right),\right.\\ & \left.\mathit{\tau }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{tot},\mathit{\theta },j}\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_i\right)\right).\end{array}\end{array}$$

At inference time, the calculated τ values are applied in latent space before propagation to physical space, scaling both epistemic and aleatoric uncertainty components. Details of the fitted scaling factors along with their effect on predictive-interval coverage and reliability are reported in Appendix C.

3.3.3 Propagation to Physical Space

Because T⁻¹ is nonlinear, the assumed Gaussian distribution in latent space does not necessarily map to a Gaussian distribution in physical space. We therefore propagate uncertainty by Monte Carlo sampling in latent space and applying inverse transforms. Specifically, starting from the S posterior draws $\mathit{\left\{}{\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{\left(\mathrm{s}\right)}\right(\stackrel{\mathrm{̃}}{\mathit{x}}),{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a}}^{\left(\mathrm{s}\right)}(\stackrel{\mathrm{̃}}{\mathit{x}}){\mathit{\}}}_{s=\mathrm{1}}^S$, we apply the latent-space scaling factors τ and then generate physical-space predictive samples via nested Monte Carlo sampling. To scale the epistemic component (the spread across S), we keep the Monte Carlo mean fixed and scale the deviations around it,

$$\begin{array}{}\text{(14)}& {\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\tau }}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)={\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)+\mathit{\tau }\odot \left({\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)-{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)\right).\end{array}$$

We then directly scale the corresponding per-sample aleatoric standard deviations.

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\tau }}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)=\mathit{\tau }\odot {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a}}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right).\end{array}$$

Next, for each posterior draw s we generate L latent samples by adding Gaussian noise,

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{\mathit{y}}}^{(\mathrm{s},\mathrm{\ell })}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)={\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\tau }}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\tau }}^{\left(\mathrm{s}\right)}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)\odot {\mathit{\eta }}^{(\mathrm{s},\mathrm{\ell })},\\ & {\mathit{\eta }}^{(\mathrm{s},\mathrm{\ell })}\sim \mathcal{N}(\mathbf{0},\mathbf{I}),\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L.\end{array}\end{array}$$

Finally, each latent sample is mapped back to physical space,

$$\begin{array}{}\text{(17)}& {\mathit{y}}^{(\mathrm{s},\mathrm{\ell })}\left(\mathit{x}\right)=T^{-\mathrm{1}}\left(Z^{-\mathrm{1}}\right({\stackrel{\mathrm{̃}}{\mathit{y}}}^{(\mathrm{s},\mathrm{\ell })}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)\left)\right),\end{array}$$

with predictive summaries (e.g., means, percentiles) computed from the resulting physical-space samples {y^(s,ℓ)(x)}.

We then estimate aleatoric and epistemic uncertainties in physical units using a nested Monte Carlo approximation to the law of total variance (Kendall and Gal, 2017). Let 𝔼_ℓ[⋅] and Var_ℓ[⋅] denote the mean and variance over $\mathrm{\ell }=\mathrm{1},\mathrm{\dots },L$ at fixed s, and let 𝔼_s[⋅] and Var_s[⋅] denote the mean and variance over $s=\mathrm{1},\mathrm{\dots },S$. We compute

$$\begin{array}{}\text{(18)}& {\displaystyle }{\widehat{\mathrm{Var}}}_{\mathrm{a}}(\mathit{y}\mid \mathit{x})\approx {\mathbb{E}}_{\mathrm{s}}\left[{\mathrm{Var}}_{\mathrm{\ell }}\left({\mathit{y}}^{(\mathrm{s},\mathrm{\ell })}\left(\mathit{x}\right)\right)\right],\text{(19)}& {\displaystyle }{\widehat{\mathrm{Var}}}_{\mathrm{e}}(\mathit{y}\mid \mathit{x})\approx {\mathrm{Var}}_{\mathrm{s}}\left({\mathbb{E}}_{\mathrm{\ell }}\left[{\mathit{y}}^{(\mathrm{s},\mathrm{\ell })}\left(\mathit{x}\right)\right]\right),\end{array}$$

element-wise per target. The physical-space aleatoric and epistemic standard deviations (σ_a and σ_e) are then given by the square root of ${\widehat{\mathrm{Var}}}_{\mathrm{a}}(\mathit{y}\mid \mathit{x})$ and ${\widehat{\mathrm{Var}}}_{\mathrm{e}}(\mathit{y}\mid \mathit{x})$, respectively. During inference, the computational burdens scale with the product of S and L. Therefore, to balance computational costs and statistical accuracy, we set $S=L=\mathrm{50}$ during inference. A sensitivity analysis shows that this results in approximately 400-fold decrease in computation time and less than 4 % change in uncertainty when compared to S and L of 1000. More details on sensitivity estimates and inference throughput are detailed in Appendices D and E.

3.4 Loss Function

Model training minimizes the negative evidence lower bound (ELBO), which balances data fit and model complexity,

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\mathcal{L}}_{\mathrm{ELBO}}& =-{\mathbb{E}}_{\mathit{\theta }\sim q_{\mathit{\varphi }}\left(\mathit{\theta }\right)}\left[\log p\left(\stackrel{\mathrm{̃}}{\mathit{y}}\mid \stackrel{\mathrm{̃}}{\mathit{x}},\mathit{\theta }\right)\right]\\ & +\mathrm{KL}\left[q_{\mathit{\varphi }}\left(\mathit{\theta }\right)\Vert p\left(\mathit{\theta }\right)\right],\end{array}\end{array}$$

where the expectation is approximated with Monte Carlo samples from the variational posterior q_ϕ(θ) (Blundell et al., 2015). The first term is the (Gaussian) negative log-likelihood in transformed space, and the second is the KL (Kullback–Leibler divergence) term applied to the Bayesian layers, regularizing the approximate posterior toward the prior.

For a given training pair $(\stackrel{\mathrm{̃}}{\mathit{x}},\stackrel{\mathrm{̃}}{\mathit{y}})$ and a weight and bias sample θ∼q_ϕ(θ), we assume a Gaussian likelihood with diagonal covariance where ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)$ are the network outputs in latent space,

$$\begin{array}{}\text{(21)}& p\left(\stackrel{\mathrm{̃}}{\mathit{y}}\mid \stackrel{\mathrm{̃}}{\mathit{x}},\mathit{\theta }\right)=\mathcal{N}\left(\stackrel{\mathrm{̃}}{\mathit{y}};{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathit{\theta }}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right),\mathrm{diag}\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{a},\mathit{\theta }}(\stackrel{\mathrm{̃}}{\mathit{x}}{)}^{\odot \mathrm{2}}\right)\right).\end{array}$$

For training stability, we use a deterministic warm-up of E_warm epochs in which the variational posterior over Bayesian-layer parameters is $q_{\mathit{\varphi }}\left(\mathit{\theta }\right)=\mathcal{N}({\mathit{\mu }}_{\mathit{\theta }},{\mathit{\sigma }}_{\mathit{\theta }}^{\mathrm{2}})$ with σ_θ fixed to a near-zero constant σ_warm≈0. In this regime, sampling yields θ≈μ_θ, so Bayesian layers behave nearly deterministically. After warm-up, we allow σ_θ to be learned and initialize it to a small but nonzero value σ_init, enabling epistemic uncertainty to develop during subsequent training epochs.

3.5 Model Optimization

We optimized the performance of the hybrid BNN model using a grid search over seven hyperparameters. For each hyperparameter set, we trained the model using 80 % of the ADDA data points, 10 % for validation during model selection, and held out the remaining 10 % for testing; performance was evaluated using the negative log-likelihood (NLL) term of the ELBO loss function. To verify that random splitting did not introduce significant information leakage, we computed test-to-train nearest-neighbor Euclidean distances in standardized input-feature space (excluding the homogeneous sphere inputs); the calculated distances had a minimum of 0.029, a 1st percentile of 0.040, and a median of 0.176, indicating that test samples are generally well separated from their nearest training neighbors. Appendix F provides full details of the hyperparameter search ranges, optimal values, and stopping criteria. Additional training diagnostics and test-set performance plots are provided in Appendix G.

4 Hybrid Bayesian Neural Network Outperforms Homogeneous Sphere and Core-Shell Approximations

We evaluate the hybrid BNN in two complementary ways. First, we assess its capability as a surrogate by comparing its predictions with ADDA reference values on the held-out test set, establishing whether the surrogate can reproduce the underlying morphology-resolving physics. Second, we compare the surrogate against homogeneous-sphere and core-shell approximations, which are the baselines used in large-scale atmospheric models and treat BC-containing particles as homogeneous or concentric volume-equivalent spheres. While the accuracy advantage over these baselines is inherited from ADDA physics, the distinct contributions of the hybrid BNN framework are the computationally efficient access to ADDA-quality optics along with per-particle predictive uncertainty with aleatoric-epistemic decomposition.

We further analyzed the accuracy of each model for predicting the optical properties of BC-containing particles in three groups based on their shape. We quantified particle shape using the fractal dimension (D_f), which is distinct from D_f,c because it accounts for the shape of the BC core and coating, as opposed to just the BC core. The structure factor of each coated particle was used to calculate D_f (Sorensen, 2001). We separated particles into three groups based on their D_f. The first group has D_f≤2.0 and represents uncoated and thinly coated particles that have not undergone significant core compaction. The second group includes particles with $\mathrm{2.0}<D_{\mathrm{f}}\le \mathrm{2.5}$, and represents moderately coated uncompacted particles and thinly coated particles that have undergone some degree of core compaction, but their surface morphologies are not spherical. The final group has D_f>2.5, and represents particles with morphologies approaching spherical. These particles include thickly coated particles that have not undergone core compaction and moderately to thickly coated particles which have undergone significant core compaction.

Figure 3 Accuracy of hybrid BNN compared with homogeneous sphere and core-shell approximations for calculating (A) extinction efficiency (Q_ext), (B) single-scattering albedo (SSA), and (C) asymmetry parameter (g). The accuracy of each model is quantified for particles with varying fractal dimension (D_f). Overall, homogeneous sphere and core-shell approximations have improved accuracy when D_f>2.5, while the hybrid BNN is accurate across the range of D_f. Lines represent median values, boxes represent the interquartile range, and error bars represent 5th and 95th percentiles.

The hybrid BNN was more accurate than homogeneous sphere and core-shell approximations for particles with D_f≤2.5 (Fig. 3a). However, the accuracy of homogeneous sphere and core-shell approximations for predicting Q_ext improved as D_f increased, a trend which has been demonstrated in previous studies (Schwarz et al., 2008; Lack and Cappa, 2010; Shamjad et al., 2012). The hybrid BNN predicted Q_ext of particles with D_f≤2.0 with median error of 0.02, while homogeneous sphere and core-shell assumptions had median errors of 0.56 and 0.35, respectively. The hybrid BNN also outperformed homogeneous sphere and core-shell approximations when predicting Q_ext of particles with $\mathrm{2.0}<D_{\mathrm{f}}\le \mathrm{2.5}$ (median errors of 0.05, 0.21, and 0.17, respectively). Conversely, the hybrid BNN was comparable to homogeneous sphere and core-shell approximations when predicting Q_ext of particles with D_f>2.5 (median errors of 0.03, 0.05, and −0.05, respectively). The hybrid BNN replicates the observational trend of increased accuracy of homogeneous sphere and core-shell approximations for particles with D_f>2.5, but is more accurate than these approximations for BC-containing particles with D_f≤2.5.

The hybrid BNN also gave the most accurate predictions of SSA for particles with D_f≤2.0 (Fig. 3b). The hybrid BNN predicted SSA with median error of 0.02 for D_f≤2.0, while homogeneous sphere and core-shell approximations both had median errors > 0.17. The reduced accuracy of Mie theory for BC-containing particles with D_f≤2.0 is rooted in their complex morphologies. Given that scattering is sensitive to the morphology of particles (Bohren and Huffman, 2008), approximating particles with D_f≤2.0 as spheres leads to large errors when calculating their scattering properties. The hybrid BNN accurately captured differences in the scattering properties of thinly coated BC-containing particles because it was trained using particles with realistically complex morphologies. The accuracy of homogeneous sphere and core-shell approximations improved for particles with $\mathrm{2.0}<D_{\mathrm{f}}\le \mathrm{2.5}$ (median errors of 0.08), but the hybrid BNN was still more accurate for these particles (median error of −0.01). Similar to predictions of Q_ext, the hybrid BNN was comparable to homogeneous sphere and core-shell approximations for predicting SSA when D_f>2.5 (median errors all < 0.01).

The hybrid BNN was also the most accurate for predicting g, with median error between −0.02 and 0.01. This is a substantial improvement over homogeneous sphere and core-shell approximations, which both systematically underestimate g (median errors between −0.38 and −0.13) (Fig. 3c). Unlike predictions of SSA, g was also underestimated by homogeneous sphere and core-shell approximations even for BC-containing particles with D_f>2.5. This is because directional scattering properties are highly sensitive to the internal morphology of particles. Non-spherical aggregates introduce complex anisotropic internal fields and coherence patterns that are not present in symmetric spheres (Bohren and Huffman, 2008). Additionally, adjoining monomers introduce near-field coupling that also affects directional scattering. Incomplete representation of the internal structure of BC-containing particles thus leads to inaccurate estimates of g under homogeneous sphere or core-shell approximations. The lower errors produced by the hybrid BNN relative to ADDA-calculated optical properties show that the hybrid BNN is able to effectively learn the results of the morphology-resolving ADDA physics code without the need to run the more computationally expensive model. Therefore, the advantage of the hybrid BNN model is rooted in its ability to provide mean and uncertainty estimates for optical properties that are reflective of morphology-resolving physics without the need to run an expensive morphology-resolving optics model.

The hybrid BNN outputs also allow for calculation of scattering and absorption efficiencies of particles (Q_abs and Q_sca). Because BC-containing particles play an important role in total absorption, their Q_abs is of particular interest for radiative transfer modeling. However, Q_abs is often not supplied directly to radiative transfer models, and is instead given by the product of Q_ext and 1−SSA. Therefore, we have trained the hybrid BNN to predict Q_ext and SSA and focused our analysis more heavily on these variables. Nevertheless, the hybrid BNN has similar accuracy to homogeneous sphere and core-shell approximations when Q_abs is calculated from Q_ext and SSA across particle D_f (see Fig. H1). Previous studies have shown that absorption by BC-containing particles is largely insensitive to particle shape (Beeler and Chakrabarty, 2022; Fierce et al., 2020), a trend which is reflected by hybrid BNN predictions. The hybrid BNN was significantly more accurate for predicting Q_sca than homogeneous sphere and core-shell approximations, regardless of particle shape. A more detailed comparison of the hybrid BNN model and core-shell and homogeneous sphere approximations for Q_abs and Q_sca can be found in Appendix H.

Finally, predictions from the hybrid BNN were also compared with those using the Rayleigh-Debye Gans (RDG) approximation (Sorensen et al., 2018). The RDG approximation accounts for the non-spherical shape of BC-containing particles, but it is not commonly implemented in large-scale atmospheric models. The hybrid BNN was more accurate than the RDG approximation for calculating Q_ext and g, and had similar accuracy to the RDG approximation for calculating SSA. More details on the RDG approximation calculations and comparisons with hybrid BNN predictions can be found in Appendix I.

5 Mean and Uncertainty Predictions using Hybrid Bayesian Neural Network

The hybrid BNN model allowed for calculation of the mean optical cross-sections of BC aggregates (with absorbing and non-absorbing coatings) along with their associated uncertainties (Fig. 4). We calculated the absorption cross-section (C_abs), scattering cross-section (C_sca), and asymmetry parameter (g) of particles with various amounts of BC, which is quantified by the core volume equivalent diameter (core d_ve). For these calculations, the D_f,c that was input to BNN was calculated from $V/V_{\mathrm{0}}$ following Eqs. (1)–(3). It should be noted that all particles included in these calculations have $D_{\mathrm{f},\mathrm{c}}<$ 2.5. The difference between hybrid BNN-calculated optical properties and those calculated using homogeneous sphere or core-shell approximations should be interpreted as the contribution of particle morphology to the optical property of interest. The calculated C_abs of BC-containing particles was larger from the hybrid BNN than from homogeneous sphere and core-shell Mie theory for uncoated and thinly coated particles (Fig. 4a–c). This trend agrees well with previous work which showed that the absorption efficiency of compact particles is smaller when compared to equivalent fractal particles (Beeler and Chakrabarty, 2022).

Figure 4 Absorption and scattering predictions at 532 nm using hybrid BNN compared with homogeneous sphere and core-shell approximations. The absorption cross-sections (C_abs) predicted by the hybrid BNN are slightly larger than those predicted by the homogeneous sphere and core-shell approximations (A–C). The scattering cross-sections (C_sca) predicted by the hybrid BNN are smaller than those predicted by the homogeneous sphere and core-shell approximations (D–F). The asymmetry parameters (g) predicted by the hybrid BNN are consistently larger than those predicted by the homogeneous sphere and core-shell approximations (G–I). Lines show predicted mean values and shaded areas represent the total standard deviation (computed from the sum of aleatoric and epistemic variances).

Conversely, the predicted C_sca of pure BC aggregates ($V/V_{\mathrm{0}}=\mathrm{1}$) was smaller using the hybrid BNN than for homogeneous sphere and core-shell approximations across core sizes (Fig. 4d). The observed lower C_sca for uncoated particles is likely due to the optical effects arising from the porosity of fractal BC particles, which has been demonstrated in previous studies (Sorensen, 2001). The C_sca of a non-spherical particle is proportional to the effective refractive index (m_eff) of a particle, given by

$$\begin{array}{}\text{(22)}& \mathit{\varphi }\left({\displaystyle \frac{m^{\mathrm{2}}-\mathrm{1}}{m^{\mathrm{2}}+\mathrm{2}}}\right)=\left({\displaystyle \frac{m_{\mathrm{eff}}^{\mathrm{2}}-\mathrm{1}}{m_{\mathrm{eff}}^{\mathrm{2}}+\mathrm{2}}}\right),\end{array}$$

where ϕ is the volume packing fraction and m is the refractive index of the material. The value of $|m_{\mathrm{eff}}-\mathrm{1}|$ is therefore inversely proportional to the porosity of a particle, and higher values of $|m_{\mathrm{eff}}-\mathrm{1}|$ lead to higher C_sca. Aggregates have much higher porosity than equivalent spheres, leading to lower $|m_{\mathrm{eff}}-\mathrm{1}|$ and lower C_sca for aggregates when compared to equivalent spheres. During coating, the ϕ and m_eff of particles evolve as coating fills voids within aggregates and compaction occurs. This leads to complex relationships between core d_ve and $V/V_{\mathrm{0}}$ (Fig. 4e). Finally, the m_eff of BC-containing aggregates approaches that of an equivalent sphere as particles become thickly coated, and relationships between core d_ve and C_sca approach homogeneous sphere and core-shell approximations (Fig. 4f).

The calculated g from the hybrid BNN is larger than those using homogeneous sphere and core-shell approximations regardless of size or coating amount (Fig. 4g–i). This trend is also rooted in the sensitivity of directional scattering to the internal morphology of particles, which is not accurately captured by homogeneous sphere and core-shell approximations.

6 Input Feature Importance

We calculated Shapley additive explanations (SHAP) values for the eight input features to examine the influence of each input feature on the predicted optical properties using the hybrid BNN model (Fig. 5). A broad SHAP value distribution indicates that the given feature has a large influence on the given output, while narrow distributions indicate small influence on the output value (Lundberg and Lee, 2017). High feature values clustered with high SHAP values indicate that increasing the value of the given feature tends to increase the value of the output (and vice versa), while randomly distributed feature values and SHAP values indicate that the feature has a small influence on the output. Our analysis showed that χ_ve had a significant effect on Q_ext, and that increasing χ_ve led to increased Q_ext (Fig. 5a). This agrees with well-established trends in aerosol optical properties, which state that extinction efficiency increases as particle size parameter increases (Bohren and Huffman, 2008; Kokhanovsky, 2008). The SHAP value distributions also showed that increasing the amount of coating (increasing $V/V_{\mathrm{0}}$) led to decreased Q_ext. This trend is likely rooted in the subtleties of the refractive index of particles. As $V/V_{\mathrm{0}}$ increases, the real and imaginary parts of the particle refractive index (volume-averaged) decrease. Decreasing the real and imaginary parts of the particle refractive index tends to reduce particle-light interactions and therefore decreases Q_ext. However, the exact relationship is not strictly monotonic and depends on other factors such as particle shape and χ_ve. Similarly, increasing κ_coat while keeping $V/V_{\mathrm{0}}$ constant increases the absorption efficiency of particles, which in turn increases Q_ext.

Figure 5 Shapley additive explanation (SHAP) values for extinction efficiency (Q_ext), single-scattering albedo (SSA), and asymmetry parameter (g) in hybrid BNN. (A) Calculated Q_ext is most sensitive to χ_ve and Q_ext,hs (broad SHAP distributions). (B) Calculated SSA is most sensitive to input $V/V_{\mathrm{0}}$, κ_coat, and SSA_hs. (C) Calculated g is most sensitive to χ_ve and Q_ext,hs.

We also found that $V/V_{\mathrm{0}}$ had a large influence on SSA (Fig. 5b), but that the effect of $V/V_{\mathrm{0}}$ on SSA diminished at high $V/V_{\mathrm{0}}$. This trend agrees with previous laboratory studies which show that the addition of scattering material leads to an increase in SSA (Kalbermatter et al., 2022; Xu et al., 2018). The diminishing effect of $V/V_{\mathrm{0}}$ on SSA at high $V/V_{\mathrm{0}}$ is also in good agreement with our current understanding of particle optics. As $V/V_{\mathrm{0}}$ increases, the extinction efficiency of particles becomes dominated by scattering while absorption effectively saturates (Beeler and Chakrabarty, 2022), both of which lead to plateauing SSA. However, increasing κ_coat increases the contribution of absorption to total extinction, leading to decreased SSA, an effect that was also well represented in the SHAP value distributions.

Finally, SHAP distributions show that χ_ve had a large effect on g (Fig. 5c). This reflects fundamental principles in directional light scattering, which predict that forward scattering is enhanced as the size of particles relative to the wavelength of incident light increases (Bohren and Huffman, 2008; Kokhanovsky, 2008). The coating absorptivity also had a significant effect on g, with increasing κ_coat leading to increased g. This trend also agrees well with fundamental particle optics, which predict that when χ_ve is < 1, increasing κ_coat effectively decreases large-angle scattering while leaving small-angle scattering relatively unchanged. This leads to increased g when κ_coat increases and χ_ve< 1. However, complex coherence and interference patterns when χ_ve> 1 lead to more complex relationships between g and κ_coat.

The input homogeneous sphere optical descriptors (Q_ext,hs, SSA_hs, and g_hs) also had a significant impact on all outputs of the hybrid BNN. For instance, a large input value for Q_ext,hs led to a larger predicted Q_ext value. To evaluate the relative importance of the combined microphysical descriptors with the homogeneous sphere optical descriptors on the target variables, we calculated the total group importance, given by the mean of ${\sum }_{j\in G}\left|{\mathit{\psi }}_j\right|$, where j is the feature in group G and ψ_j are the corresponding SHAP values. The total group importance gives the cumulative influence of the group on the prediction. This analysis showed that the total importance of the microphysical descriptors was 0.38, 0.08, and 0.07 for Q_ext, SSA, and g, respectively. The importance of the homogeneous sphere optical descriptors on the predictions was 0.80, 0.13, and 0.10 for Q_ext, SSA, and g, respectively. These findings collectively show that the homogeneous sphere optical descriptors are very important for hybrid BNN predictions. Including the homogeneous sphere optical properties as an input provides the hybrid BNN with a large fraction of the underlying physics involved in light scattering and absorption, which allowed the hybrid BNN to be trained on a relatively small number of data points. The importance of Q_ext,hs, SSA_hs, and g_hs for hybrid BNN prediction also indicates that the hybrid BNN may be accurate for refractive indices and monomer sizes that were not represented in the training set if the inputs are reflective of the new refractive index or monomer size.

7 Extrapolation of Model Predictions Across Refractive Indices and Monomer Sizes

The hybrid BNN training set included optical properties of BC-containing particles with monomer radii of 20 nm, coatings with real refractive indices (Re(RI_coat)) of 1.55 or 1.6, and BC with refractive index (RI_BC) of 1.95+0.79i. We performed a perturbed-parameter ensemble to evaluate the accuracy of the hybrid BNN for particles with monomer radii, Re(RI_coat), and RI_BC outside of these values. We randomly selected 10 points from the held-out test set and systematically changed the input monomer radii in ADDA simulations. The same process was repeated for simulations in which Re(RI_coat) and RI_BC were altered, resulting in a total of 90 additional ADDA simulations. We then used the BNN to predict the optical properties of the particles and compared the BNN predictions to the ADDA simulations. We benchmarked model error across Re(RI_coat), RI_BC, and monomer radii using the average difference between calculated and true values, where true values are given by ADDA simulations (Fig. 6). The shaded regions in Fig. 6 denote the parameter values represented in training.

Figure 6 Mean error in extinction efficiency (Q_ext), single-scattering albedo (SSA), and asymmetry parameter (g) calculated using the hybrid BNN, homogeneous sphere approximation, and the core-shell approximation for a perturbed parameter ensemble. Here, error is given by the difference between predicted and true values. (A–C) The hybrid BNN shows varying accuracy as the real part of the coating refractive index (Re(RI_coat)), BC refractive index (RI_BC), or monomer radius deviate from training set values. However, it still has comparable accuracy to homogeneous sphere and core-shell approximations. (D–I) The hybrid BNN is more accurate than homogeneous sphere and core-shell approximations when calculating SSA and g across all Re(RI_coat), RI_BC, and monomer radius. Shaded areas in each panel indicate values that were used in the training set and lines show average values.

Deviations in Re(RI_coat), RI_BC, or monomer radius from training values affected the accuracy of the hybrid BNN for predicting Q_ext, but the hybrid BNN still had comparable errors to homogeneous sphere and core-shell approximations (Fig. 6a–c). These results indicate that the hybrid BNN loses predictive accuracy for calculating Q_ext under high relative humidities, where the refractive index of water would contribute significantly to Re(RI_coat). However, it is likely that the hybrid BNN would still be more accurate than homogeneous sphere or core-shell approximations under dry to moderate relative humidities, assuming dry coatings have Re(RI_coat) between 1.5 and 1.6 (Bond and Bergstrom, 2006). The hybrid BNN was more accurate than homogeneous sphere and core-shell approximations for calculating SSA and g across Re(RI_coat), RI_BC, and monomer radii (Fig. 6d–i).

It should be noted that the monomer radius, Re(RI_coat), and RI_BC are not directly included as inputs to the hybrid BNN. However, changing the monomer radius alters χ_ve, and changing the refractive indices of BC or the coating material changes the Mie theory-calculated optical properties (Q_ext,hs, SSA_hs, and g_hs), both of which are used as inputs to the hybrid BNN. Including Q_ext,hs, SSA_hs, and g_hs as hybrid BNN inputs allows for extrapolation of BNN predictions to monomer sizes and refractive indices that were not used to train the model with little loss in accuracy.

8 Future Experiment Design Using Epistemic Uncertainty Predictions

One of the unique features of the hybrid BNN is the ability to separate aleatoric and epistemic uncertainty in predicted optical properties. To demonstrate that epistemic uncertainty can identify sparse regions in input-parameter space, the hybrid BNN was retrained after masking approximately 12 % of the training data in the region $\mathrm{70}\le N_{\mathrm{pp}}\le \mathrm{180}$ and $\mathrm{8}\le V/V_{\mathrm{0}}\le \mathrm{25}$. Test particles within this masked region were retained for evaluation but were not used during re-training. Epistemic uncertainty increased more inside the masked region than outside it for all three target variables (see Appendix J). Because epistemic uncertainty increases in regions with sparse training data, patterns in epistemic uncertainty can be used to identify regions of the parameter space that would most benefit from additional training data. Therefore, uncertainty estimates from the hybrid BNN can be leveraged to guide the design of future simulations that would most effectively improve model performance.

Figure 7 Epistemic uncertainty in predicted optical properties can be used for targeted future experiments aimed at reducing model uncertainty. (A–C) Epistemic uncertainty in predicted Q_ext is highest for large, thickly coated particles (N_pp>800 and $V/V_{\mathrm{0}}>\mathrm{40}$) at 405 nm. Epistemic uncertainty for SSA and g is largest for small, thinly coated particles at 532 and 880 nm. (D–I) For SSA and g, epistemic uncertainty across wavelengths is largest for small, thinly coated particles with N_pp<100 and $V/V_{\mathrm{0}}<\mathrm{10}$.

We used the hybrid BNN to predict aleatoric and epistemic uncertainty in Q_ext, SSA, and g at three wavelengths commonly used in optical instruments (405, 532, and 880 nm). Here, we quantified error in optical properties as a percentage of the predicted mean value. Overall, we found that aleatoric uncertainty is highest for smaller, thinly coated aggregates, which is due to the high variability in the morphology of these particles as they compact due to the accumulation of coatings (see Appendix K). We also found that epistemic Q_ext uncertainty at 405 nm is highest for large, thickly coated particles with N_pp>800 and $V/V_{\mathrm{0}}>\mathrm{40}$ (Fig. 7a). Epistemic Q_ext uncertainty at 532 nm is highest for thinly coated particles ($V/V_{\mathrm{0}}<\mathrm{5}$) and N_pp<20 or N_pp>700 (Fig. 7b). At 880 nm, epistemic Q_ext uncertainty is highest for small, thinly coated particles with N_pp<50 and $V/V_{\mathrm{0}}<\mathrm{10}$ (Fig. 7c). For SSA and g, epistemic uncertainty across wavelengths is largest for small, thinly coated particles with N_pp<100 and $V/V_{\mathrm{0}}<\mathrm{10}$ (Fig. 7d–i).

Given that ADDA simulations are extremely computationally expensive, this analysis can be used to identify regions of the input parameter space that are either sparsely represented or where there is significant variability in the optical properties of particles. Including additional simulations that fall within these areas of the input parameter space would be the most beneficial for improved model performance. By targeting these sparsely represented input parameters, future simulations can be designed more efficiently, maximizing model improvement while reducing computational burdens. These results also show that future laboratory experiments with BC-containing particles with N_pp<100 and $V/V_{\mathrm{0}}<\mathrm{10}$ would also aid in constraining model uncertainty. Previous studies have shown that particles in this region are in the midst of moderate restructuring that is dependent on N_pp (Beeler et al., 2025). Therefore, these particles are expected to have more morphological variability than particles with N_pp>100 and $V/V_{\mathrm{0}}>\mathrm{10}$. Additional laboratory experiments targeting these particles would therefore be particularly valuable for constraining their optical properties. Measurements that jointly constrain aggregate size, coating amount, core and whole-particle morphology, and wavelength-resolved optical properties would provide especially useful constraints for evaluating and informing future versions of the model.

9 Discussion and Conclusions

Efficiently predicting the optical properties of BC-containing particles remains a longstanding challenge in atmospheric models and remote sensing algorithms (Ramanathan and Carmichael, 2008; Bond and Bergstrom, 2006). Traditional atmospheric models employ simplifications to BC morphology and mixing with other aerosol components to reduce computational burdens (Liu et al., 2016). However, these simplifications often do not capture the complexity of light scattering and absorption by irregularly shaped BC-containing particles. Previous studies have attempted to improve predictions of the optical properties of BC-containing particles using machine learning-based approaches (Qin et al., 2025; Romshoo et al., 2024; Luo et al., 2018). However, these studies use single-particle optics models which do not accurately capture near-field interactions between adjacent monomers in BC aggregates and do not quantify uncertainty in model predictions due to variability in particle morphology.

We developed a hybrid BNN to predict the optical properties of BC-containing particles and their associated uncertainties. The model was trained using ADDA simulations of single-particle optical properties, which accurately represent the morphology and near-field interactions between monomers in BC-containing particles (Yurkin and Hoekstra, 2011). The hybrid BNN was able to effectively learn the morphology-resolving physics of ADDA simulations while reducing computation time by 2–6 orders of magnitude. The hybrid BNN takes in five inputs that are easily tracked or derived in large-scale models and fully describe the size, shape, and composition of BC-containing particles. The hybrid BNN also uses the optical properties of a volume-equivalent homogeneous sphere as an input feature. The optical properties of a volume-equivalent homogeneous sphere can be calculated with low computational cost and inform the hybrid BNN of the complex physics involved in particle-light interactions. Our results showed that the hybrid BNN more accurately predicts mean BC optical properties than traditional homogeneous-sphere or core–shell approximations, with the largest improvements for SSA and g. The hybrid BNN also quantifies uncertainty in the predicted optical properties of BC-containing particles.

We found that the size of the particle relative to the wavelength, imaginary part of the coating refractive index, and the optical properties of a homogeneous volume-equivalent sphere were the most influential input features for model predictions. We also showed that the hybrid BNN is able to extrapolate to monomer sizes, BC refractive indices, and coating refractive indices that were not represented in the training set with little loss in accuracy. This is because inclusion of homogeneous sphere optical descriptors encodes much of the physics of these parameters on target variables, and therefore makes hybrid BNN predictions applicable to various assumed BC and coating refractive indices and monomer sizes. Finally, uncertainty decomposition in the hybrid BNN allows for identification of regions in the input parameter space that are underrepresented in the training dataset or have significant variability in the optical properties. This allows for targeted future additions to the training dataset that will be most effective in improving model predictions.

It should be noted that the hybrid BNN introduced here is a single-particle model that is not suitable for direct incorporation into large-scale atmospheric models. Multiple lines of evidence show that large-scale models do not accurately represent the mixing state of BC-containing particles, which is extremely important for predicting their optical properties (Fierce et al., 2020, 2016). The hybrid BNN introduced here does not address this gap, and therefore is not likely to improve BC optics calculations in large-scale models. However, our hybrid BNN can be used as a computationally efficient surrogate model for calculating BC optics in particle-resolved models, which often contain thousands of particles. Bulk optical properties from particle-resolved simulations can then be used to evaluate optics schemes in large-scale models, which will be the subject of future studies. Additionally, the hybrid BNN framework could be extended to predict scattering quantities relevant for remote sensing applications, such as phase function or depolarization. This will also be the subject of future studies.

Overall, our results present the first machine learning-based model for BC optical properties with uncertainty quantification. It is also the first machine learning-based model for BC optical properties that is trained using ADDA simulations, to our knowledge. Our hybrid BNN accurately predicts the optical properties of BC-containing particles across the morphology-composition space. Quantification of uncertainty in single-particle optical properties allows for propagation to uncertainty in bulk optical properties, which are important for comparison against measurements and for probabilistic predictions of BC behavior in large-scale atmospheric models. Improved predictions of BC optical properties will be a major step forward for improved representation of anthropogenic and biomass burning emissions in large-scale atmospheric models.

Appendix A: ADDA Run Configuration

Version 1.5.0-alpha2 of the ADDA code was used to calculate the optical properties of simulated BC-containing particles. The ADDA user manual recommends keeping the number of dipoles per wavelength ($\mathit{\lambda }/d$) above 10 to keep errors in calculated optical properties to within a few percent. Since computation time increases steeply with the number of dipoles (Fig. A1a), particles were discretized using a dipole length of 0.0033 µm to balance simulation accuracy with computational cost ($\mathit{\lambda }/d>\mathrm{90}$ for all simulations). ADDA simulations were performed on the Deception high-performance computing cluster at Pacific Northwest National Laboratory using a single AMD EPYC 7502 CPU core for each simulation. Finally, we calculated the optical properties of particles for 150 randomly generated orientations and used the average optical properties over the 150 orientations for hybrid BNN training and testing. We randomly selected 40 particles from the hybrid BNN test set and calculated their optical properties for 300 random orientations. We found that using the average optical properties over 150 orientations deviated by 0.1 % or less compared to the average value over 300 orientations (Fig. A1b). We therefore expect errors in the calculated optical properties from ADDA due to discretization and orientation-averaging to be negligible.

Figure A1 Timing and errors from ADDA simulations. (A) Computation time per orientation scales approximately as a power law with the number of dipoles. (B) Optical properties converge after approximately 100 orientations.

Appendix B: Z-Score Normalization

Z-score normalization ensures that each variable has mean zero and standard deviation one in the normalized space. The forward Z-score transformation is

$$\begin{array}{}\text{(B1)}& Z\left(x\right)={\displaystyle \frac{x-{\mathit{\mu }}_x}{{\mathit{\sigma }}_x}},\end{array}$$

for input variables, and

$$\begin{array}{}\text{(B2)}& Z\left(y_{\mathrm{tf}}\right)={\displaystyle \frac{y_{\mathrm{tf}}-{\mathit{\mu }}_{y_{\mathrm{tf}}}}{{\mathit{\sigma }}_{y_{\mathrm{tf}}}}},\end{array}$$

for transformed target variables (Fei et al., 2021). Here, μ and σ are the mean and standard deviation of the corresponding variable computed on the training set only; the same training-set statistics are then used to normalize the validation and test sets. The input and transformed target variables are mapped and stored as $\mathit{x}\stackrel{Z}{⟶}\stackrel{\mathrm{̃}}{\mathit{x}}$ and $\mathit{y}\stackrel{T}{⟶}{\mathit{y}}_{\mathrm{tf}}\stackrel{Z}{⟶}\stackrel{\mathrm{̃}}{\mathit{y}}$.

Appendix C: Uncertainty Scaling Calibration

We applied post-hoc uncertainty scaling to calibrate the predictive uncertainty of the hybrid BNN for each target (Q_ext, SSA, and g). Specifically, we fit a per-target multiplicative scaling factor τ on the validation set by minimizing the Gaussian negative log-likelihood in latent space (Sect. 3.3.2). Table C1 reports the fitted τ values and the empirical coverage of nominal 90 % prediction intervals before (Pre) and after (Post) calibration on the validation and test sets. Values of τ<1 indicate that the raw BNN uncertainties were overestimated. Figure C1 shows reliability diagrams comparing expected versus observed coverage; curves closer to the 1:1 line indicate better-calibrated predictive uncertainty.

Figure C1 Reliability diagrams for hybrid BNN predictive uncertainty for Q_ext, SSA, and g. The dashed line indicates ideal calibration (observed coverage equals expected coverage), and the points/solid line show observed coverage. Deviations from the 1:1 line indicate miscalibration.

Table C1 Uncertainty scaling calibration results. All fitted values have τ<1, indicating that the raw BNN overestimated uncertainty.

Appendix D: Uncertainty Convergence

Table D1 Deviation of predictive uncertainty at representative (S,L) combinations with respect to the converged reference at $(S,L)=(\mathrm{1000},\mathrm{1000})$.

Predictive uncertainties reported in the main text are computed by nested Monte Carlo sampling, with S posterior parameter draws and L aleatoric samples each. We characterize the convergence of the two components with respect to S and L using ${\mathit{\sigma }}_{\mathrm{e}}^{\ast }$ and ${\mathit{\sigma }}_{\mathrm{a}}^{\ast }$ to denote the epistemic and aleatoric standard deviations at a converged reference $(S,L)=(\mathrm{1000},\mathrm{1000})$. All deviations in this section are expressed with respect to that reference.

Figure D1 shows the two components as functions of S (left column, with L pinned at 1000) and L (right column, with S pinned at 1000), for each of the three targets. As the paper operates at $S=L=\mathrm{50}$ on both axes simultaneously, we additionally report the joint deviation at that point along with several representative combinations in Table D1. Values are pooled-median deviations across 10 test particles and 3 seeds. At the paper default, all six deviations remain within ±4 %.

The choice of $S=L=\mathrm{50}$ is a cost–accuracy trade-off. Nested-MC sampling has O(S⋅L) cost in both compute and memory, so moving from $(S,L)=(\mathrm{50},\mathrm{50})$ to the converged reference at (1000,1000) increases per-particle wall-clock and memory footprint by a factor of ∼400 while tightening the reported uncertainties by less than 4 % on any single target. The paper default therefore trades a small, quantified loss of precision for two orders of magnitude of throughput. Absolute inference timings and batch scaling are reported in Appendix E.

Figure D1 Monte Carlo convergence of predictive uncertainty. Curves show the pooled median across 10 test particles and 3 seeds of the ratio $\mathit{\sigma }/{\mathit{\sigma }}^{\ast }$, where σ^∗ is the value at the converged reference $(S,L)=(\mathrm{1000},\mathrm{1000})$. Shaded bands show the interquartile range across the pooled sample; the dashed vertical line marks the paper default $S=L=\mathrm{50}$.

Appendix E: Inference Throughput

Table E1 Hybrid BNN inference cost as a function of batch size at $S=L=\mathrm{50}$. Per-particle time and throughput are computed from wall clock and batch size.

We report the wall-clock cost of the hybrid BNN at the paper default parameters of $S=L=\mathrm{50}$. All timings were measured on a single Intel Core i7-11850H workstation (8 physical cores, 16 GB RAM, CPU-only PyTorch 2.4 with Pyro 1.9) using the median wall clock of 20 repeated inference passes after a single warm-up call. Sensitivity of the reported uncertainties to the (S,L) sampling choice is documented separately in Appendix D.

Single-particle inference (batch =1) is dominated by fixed per-call overhead. As batch size grows the overhead is distributed across particles and the per-particle cost drops toward the asymptotic compute cost of the network forward passes and nested-MC sampling. The per-particle time plateaus around ∼1 ms for batches of 10³–10⁴ particles, delivering approximately 10³ particles s⁻¹ on a single CPU core.

We note that the current inference pipeline is CPU-only, as the Bayesian sampling stack (Pyro's variational inference) does not straightforwardly support CUDA in our configuration. A GPU-enabled implementation would be expected to sustain substantially higher throughput at large batch sizes; extending the pipeline to GPU is left as future engineering work.

Appendix F: Hybrid BNN Hyperparameter Selection

Table F1 summarizes the hyperparameters explored during model selection and the optimal configuration. Layers denotes the total number of hidden layers in the network, and Bayesian layers denotes the number of variational (Bayesian) layers used to represent weight uncertainty. Bayesian layer position indicates where the Bayesian layer(s) were placed within the network (middle or last). The learning-rate scale factor multiplies the base learning rate used by the optimizer for any Bayesian layers (post E_warm), and the prior standard deviation sets the scale of the Gaussian prior on Bayesian-layer weights.

The warm-up parameter E_warm denotes the number of initial epochs used for the deterministic warm-up described in the main text. Training used early stopping based on the validation loss: if the validation loss did not decrease by at least ${\mathrm{\Delta }}_{min}=\mathrm{0.001}$ for 100 consecutive epochs, training was terminated.

For the optimal model, the hidden-layer widths (all using ReLU activations) were 128, 128, 128 (Bayesian), 64, and 32 units. Additional training diagnostics and test-set performance plots corresponding to the selected hyperparameter configuration are provided in Appendix G.

Table F1 Hyperparameters explored during hybrid BNN model selection (Sect. 3.5), including sampled values and the optimal configuration.

Appendix G: Hybrid BNN Training and Test-Set Performance

Figure G1 shows the evolution of the negative log-likelihood (NLL) and KL-divergence components of the hybrid BNN objective during training for the training and validation sets. The similar trajectories of training and validation NLL indicate stable optimization and limited overfitting over the training period.

Figure G2 compares test-set predictions from the hybrid BNN (points with 90 % predictive confidence intervals) and the homogeneous-sphere approximation (orange markers) against ADDA reference values (dashed 1:1 line). Points closer to the 1:1 line indicate higher accuracy relative to ADDA, and the BNN error bars summarize predictive uncertainty for individual test samples.

Figure G1 Training history for the hybrid BNN loss components. Left: negative log-likelihood (NLL) for training and validation sets. Right: KL divergence term for training and validation sets.

Figure G2 Test-set comparison of predicted versus ADDA reference values for Q_ext, SSA, and g. Orange markers show the homogeneous-sphere approximation. Blue points show hybrid BNN predictions with 90 % predictive confidence intervals. The dashed line indicates the 1:1 relationship with ADDA.

Appendix H: Absorption and Scattering Efficiencies from Hybrid Bayesian Neural Network

Predictions of Q_ext and SSA were combined to calculate the absorption (Q_abs) and scattering (Q_sca) efficiencies of particles (see Fig. H1). The hybrid BNN calculated Q_abs with median error of 0.01, 0.03, and 0.01 for D_f≤2.0, $\mathrm{2.0}<D_{\mathrm{f}}\le \mathrm{2.5}$, and D_f>2.5, respectively. The BNN was slightly more accurate than homogeneous sphere and core-shell approximations for D_f≤2.0 and D_f>2.5, but was less accurate for $\mathrm{2.0}<D_{\mathrm{f}}\le \mathrm{2.5}$. This shows that homogeneous sphere and core-shell approximations can predict Q_abs with reasonable accuracy, but the hybrid BNN model also predicts Q_abs with similar accuracy. Conversely, the hybrid BNN greatly outperforms homogeneous sphere and core-shell approximations when calculating Q_sca. Median error for Q_sca is −0.02 to 0.05 using the hybrid BNN model, 0.06–0.49 using the homogeneous sphere approximation, and −0.02 to 0.51 using the core-shell approximation. Core-shell and homogeneous sphere approximations again have improved accuracy as D_f increases. Overall, homogeneous sphere and core-shell approximations may be adequate for calculating absorption, but are not accurate for calculating the scattering properties of BC-containing particles with D_f≤2.5.

Figure H1 Accuracy of hybrid BNN, homogeneous sphere approximation, and core-shell approximation for calculating absorption efficiency (Q_abs) and scattering efficiency (Q_sca) as a function of particle shape (quantified by D_f). (A) The three methods have comparable accuracy for predicting Q_abs across particle shapes. (B) The hybrid BNN is more accurate than homogeneous sphere and core-shell approximations when predicting Q_sca, regardless of particle shape. Lines represent median values, boxes represent the interquartile range, and error bars represent the 5th and 95th percentiles.

Appendix I: Comparison of Hybrid BNN Predictions with the RDG Approximation

Figure I1 Comparison of hybrid BNN accuracy with the core-shell RDG approximation. The hybrid BNN was more accurate than the RDG approximation for Q_ext and g, and had similar accuracy to the RDG approximation for SSA.

We compared the accuracy of hybrid BNN predictions with those using the Rayleigh-Debye Gans (RDG) approximation. The RDG approximation accounts for the non-spherical shape of BC-containing aggregates by assuming that the monomers comprising a BC-containing particle are within the Rayleigh optics limit (Sorensen et al., 2018). The optical properties of particles in the test set were calculated using Eqs. (16)–(23) of Sorensen et al. (2018). To account for the effects of coating on RDG approximation calculations, the optical response of each coated monomer is represented using the electrostatic polarizability of a concentric core-shell sphere (β_cs). Defining the relative permittivity of the core, shell, and surrounding medium as ε_c, ε_s, and ε_m, respectively, and the core volume fraction (f) as

$$\begin{array}{}\text{(I1)}& f={\left({\displaystyle \frac{a_{\mathrm{c}}}{a_{\mathrm{o}}}}\right)}^{\mathrm{3}},\end{array}$$

where a_c and a_o are the core and outer radii. The polarizability is then calculated using

$$\begin{array}{}\text{(I2)}& {\mathit{\beta }}_{\mathrm{cs}}={\displaystyle \frac{({\mathit{\epsilon }}_{\mathrm{s}}-{\mathit{\epsilon }}_{\mathrm{m}})({\mathit{\epsilon }}_{\mathrm{c}}+\mathrm{2}{\mathit{\epsilon }}_{\mathrm{s}})+f({\mathit{\epsilon }}_{\mathrm{c}}-{\mathit{\epsilon }}_{\mathrm{s}})({\mathit{\epsilon }}_{\mathrm{m}}+\mathrm{2}{\mathit{\epsilon }}_{\mathrm{s}})}{({\mathit{\epsilon }}_{\mathrm{s}}+\mathrm{2}{\mathit{\epsilon }}_{\mathrm{m}})({\mathit{\epsilon }}_{\mathrm{c}}+\mathrm{2}{\mathit{\epsilon }}_{\mathrm{s}})+\mathrm{2}f({\mathit{\epsilon }}_{\mathrm{s}}-{\mathit{\epsilon }}_{\mathrm{m}})({\mathit{\epsilon }}_{\mathrm{c}}-{\mathit{\epsilon }}_{\mathrm{s}})}},\end{array}$$

with

$$\begin{array}{}\text{(I3)}& {\mathit{\epsilon }}_i={\left({\displaystyle \frac{m_i}{m_{\mathrm{m}}}}\right)}^{\mathrm{2}}.\end{array}$$

Here, m_i denotes the complex refractive index of the core or shell and m_m is the refractive index of the surrounding medium (1.0+0.00i). This expression reduces to the homogeneous-sphere Lorentz–Lorenz factor when the core and shell have identical optical properties. We therefore use β_cs as opposed to the Lorentz-Lorenz factor when calculating F(m) and E(m). It should also be noted that the fractal dimension and fractal prefactor were assumed fixed at 1.78 and 1.35, respectively, in RDG approximation calculations.

The hybrid BNN is more accurate than the RDG approximation for Q_ext, with median errors between (0.02,0.05) and (−1.4, 0.46), respectively. The hybrid BNN had similar accuracy to the RDG approximation for calculating SSA, with median errors between (−0.01, 0.02) and (−0.04, −0.01), respectively. The accuracy of RDG approximation for calculating SSA but not Q_ext indicates that RDG approximation tends to underestimate total absorption and/or scattering cross-sections of coated aggregates, but approximately captures their relative contributions to total extinction. Finally, the hybrid BNN is also more accurate than the RDG approximation when calculating g for coated aggregates, with median errors between (−0.02, 0.01) and (−0.20, −0.15), respectively.

Appendix J: Sensitivity of Predictive Uncertainty to Data Sparsity

A defining feature of the hybrid BNN is that its epistemic uncertainty σ_e represents the model's own uncertainty in underrepresented regions during training. As a direct test of that interpretation, we retrained the hybrid BNN after masking approximately 12 % of the training data in the region $\mathrm{70}\le N_{\mathrm{pp}}\le \mathrm{180}$ and $\mathrm{8}\le V/V_{\mathrm{0}}\le \mathrm{25}$. Test particles within this masked region were retained for evaluation but were not used during retraining. All other hyperparameters and the random seed were identical to the baseline model.

Figure J1 shows the per-particle changes in epistemic uncertainty (σ_e) and aleatoric uncertainty (σ_a) between the baseline model (trained on the full dataset) and the retrained model (trained on the masked dataset). For test particles inside the masked region, epistemic uncertainty increased by 14.4 %, 28.6 %, and 11.9 % for Q_ext, SSA, and g, respectively. For test particles outside the masked region, the corresponding increases were 3.9 %, 5.0 %, and −1.4 %. The increase in epistemic uncertainty was larger inside the masked region than outside for all three target variables, supporting the interpretation that epistemic uncertainty responds to reduced training coverage. We also examined the corresponding changes in aleatoric uncertainty. For test particles inside the masked region, aleatoric uncertainty increased by 6.9 %, 27.9 %, and 3.6 % for Q_ext, SSA, and g, respectively, compared with changes of 0.8 %, 8.7 %, and −0.4 % outside the masked region. Because the aleatoric and epistemic components are learned jointly through the same predictive model, changes in the training distribution can affect both components. We therefore interpret this experiment primarily as a test of whether epistemic uncertainty responds to reduced training coverage.

Figure J1 Per-particle changes in epistemic uncertainty (σ_e) and aleatoric uncertainty (σ_a) when the hybrid BNN is retrained after masking approximately 12 % of the training data in the region $\mathrm{70}\le N_{\mathrm{pp}}\le \mathrm{180}$ and $\mathrm{8}\le V/V_{\mathrm{0}}\le \mathrm{25}$. Points show the relative change in each uncertainty component with respect to the baseline model for each held-out test particle.

Appendix K: Aleatoric Uncertainty Estimates

We used the hybrid BNN to predict aleatoric uncertainty in Q_ext, SSA, and g at 405, 532, and 880 nm. Predicted uncertainties in optical properties are expressed as a percentage of the predicted mean value. We found that aleatoric uncertainty is highest in regions where particles have the highest morphological variability (uncoated or thinly coated particles). This shows that aleatoric uncertainty is able to identify regions of the parameter space where uncertainty due to particle variability is highest.

Figure K1 Aleatoric uncertainty as a function of BC core size (quantified by N_pp) and coating amount (quantified by $V/V_{\mathrm{0}}$) at three wavelengths commonly used in aerosol optical instruments. We found that aleatoric uncertainty is generally higher for thinly coated aggregates, which reflects high morphological variability of fractal BC aggregates before coating-induced compaction.