Quantifying the structural uncertainty of the aerosol mixing state representation in a modal model

Aerosol mixing state is an important emergent property that affects aerosol radiative forcing and aerosol–cloud interactions, but it has not been easy to constrain this property globally. This study aims to verify the global distribution of aerosol mixing state represented by modal models. To quantify the aerosol mixing state, we used the aerosol mixing state indices for submicron aerosol based on the mixing of optically absorbing and non-absorbing species (χ_o), the mixing of primary carbonaceous and non-primary carbonaceous species (χ_c), and the mixing of hygroscopic and non-hygroscopic species (χ_h). To achieve a spatiotemporal comparison, we calculated the mixing state indices using output from the Community Earth System Model with the four-mode version of the Modal Aerosol Module (MAM4) and compared the results with the mixing state indices from a benchmark machine-learned model trained on high-detail particle-resolved simulations from the particle-resolved stochastic aerosol model PartMC-MOSAIC. The two methods yielded very different spatial patterns of the mixing state indices. In some regions, the yearly averaged χ value computed by the MAM4 model differed by up to 70 percentage points from the benchmark values. These errors tended to be zonally structured, with the MAM4 model predicting a more internally mixed aerosol at low latitudes and a more externally mixed aerosol at high latitudes compared to the benchmark. Our study quantifies potential model bias in simulating mixing state in different regions and provides insights into potential improvements to model process representation for a more realistic simulation of aerosols towards better quantification of radiative forcing and aerosol–cloud interactions.

1 Introduction

The direct and indirect climate effects of atmospheric aerosols greatly depend on the particles' spatial distribution in the atmosphere and their climate-relevant properties, including their hygroscopicity, optical properties, and their ability to act as cloud condensation nuclei (CCN) and ice nuclei (Boucher et al., 2013). These properties, in turn, are closely related to the aerosol mixing state (Ching et al., 2012; Cziczo et al., 2009; Fierce et al., 2016, 2017). Aerosol mixing state refers to the way in which different aerosol chemical species are distributed among and within the aerosol particles (Riemer et al., 2019). As shown in many observational field studies, atmospheric aerosols have complex mixing states (Bondy et al., 2018; Healy et al., 2014; Lee et al., 2019; Ye et al., 2018; Yu et al., 2020), ranging between the two extremes of an “internal mixture”, where the composition of all particles within the population is identical (and equal to the bulk composition of the aerosol), and an “external mixture”, where each particle in a population consists of only a single species (which may be different for each particle).

This poses a unique challenge for the modeling of aerosols in Earth system models, which, for the sake of computational efficiency, represent aerosols by simplifying the true aerosol mixing state using various mixing-state-related assumptions. For example, bulk aerosol models predict the abundance of individual aerosol chemical species by tracking the species' mass concentrations, inherently treating the aerosol as external mixtures of, e.g., sulfate, black carbon, organic carbon, sea salt, and dust (Ghan et al., 2012). Univariate sectional models are able to represent size-resolved composition but cannot resolve the diversity of the aerosol within a certain size range. For modal models, the ability to resolve mixing state depends on the definition and the placement of the modes. Different approaches for modal models have been developed, ranging from a small number of internally mixed, non-overlapping modes (e.g., three modes in MAM3, Liu et al., 2012; or CMAQv5.2, United States Environmental Protection Agency, 2017) to a larger number of modes that may overlap in a given size range and separate out different aerosol mixtures (e.g., nine modes in MADE3, Kaiser et al., 2014; or 16 modes in MATRIX, Bauer et al., 2008). For these multi-modal models, the aerosol processes of gas–aerosol partitioning and coagulation make it necessary to define rules for how the modes interact (Wilson et al., 2001). Condensation of secondary aerosol on a mode reserved for a pure species (e.g., black carbon or dust) requires moving mass over to a mixed mode when a critical mass fraction of secondary aerosol is exceeded. Transfer terms due to coagulation of particles in different modes can be calculated analytically (Binkowski and Shankar, 1995), and rules need to be defined regarding the destination mode after coagulation. Generally, the transfer of aerosol mass from smaller modes to larger modes during growth can lead to inaccuracies. The removal of particles due to scavenging by cloud activation is another issue that is difficult to reconcile. Hence, the choice of the number of modes, their compositions, and the criteria for transfer between modes are user-defined, which introduces structural uncertainty in aerosol simulations that still needs to be quantified.

Given that modal models are to some extent mixing-state-aware, the following question arises: how well do modal models represent mixing state? Due to the scarcity of relevant observational data, we are not yet at the point where we can comprehensively validate model output of aerosol mixing state as is done for other aerosol-related quantities, such as bulk mass concentrations or aerosol optical depth. However, higher-detail models can serve as benchmarks to perform a verification of simulated aerosol mixing state. This paper aims to verify the global distribution of aerosol mixing state represented by a modal model by using benchmark simulations from the particle-resolved stochastic aerosol model PartMC-MOSAIC (Riemer et al., 2009; Zaveri et al., 2008). Our usage of the term “aerosol representation” in this paper encompasses the representation of processes that go along with the aerosol representation itself, since the two are in practice tightly coupled.

We used the aerosol mixing state index χ (Riemer and West, 2013) as a metric to quantify aerosol mixing state. The mixing state index χ can be interpreted as a label for particle populations to rigorously characterize where the population lies on the spectrum from external (χ=0 %) to internal (χ=100 %) mixture. This concept has been successfully applied to observational data (Healy et al., 2014; Ye et al., 2018) and for error quantification studies (Ching et al., 2017, 2018). Particularly relevant for this work is the study by Ching et al. (2017), which showed that assuming an internal mixture when the aerosol is actually not completely internally mixed can result in errors of up to 150 % in CCN predictions.

PartMC-MOSAIC tracks the composition of individual particles and therefore resolves aerosol mixing state explicitly (Riemer et al., 2009; Zaveri et al., 2008). However, this modeling approach is computationally very expensive and therefore not practical for large-scale simulations of several months or years of simulation time. To estimate the global spatial distribution of mixing state, we recently developed a machine-learned (ML) model based on high-detail particle-resolved simulations (Zheng et al., 2021) that uses inputs that are known from global model simulations to predict χ. In this paper, we use this ML model to predict the spatial distribution of the mixing index χ and then compare the results with χ values that are derived from the Community Earth System Model version 2 (CESM2 version 2.1.0; Danabasoglu et al., 2020) using the four-mode version of the Modal Aerosol Module (MAM4; Liu et al., 2016).

This paper is organized as follows. In Sect. 2 we introduce the setup of the Earth system model simulations. The definition of mixing state indices and the derivation of aerosol mixing state indices for modal models are given in Sect. 3. Section 4 briefly describes the ML model generated with machine learning and particle-resolved modeling for estimating the benchmark aerosol mixing state indices. Section 5 focuses on the comparison of mixing state indices from the particle-resolved and modal models, and Sect. 6 summarizes our findings.

2 Global model simulations

Here we employed CESM2 to provide the global model simulation data. Specifically, we used the component set FHIST to set up the global simulations with aerosols. This component set represents a typical historical simulation in the Community Atmospheric Model (CAM6; Bogenschutz et al., 2018) using an active atmosphere and land with prescribed sea-surface temperatures and sea-ice extent, as well as a 1^∘ finite-volume dycore with the forcing data available from 1979 to 2015.

MAM4 is the default aerosol module of this component set, which represents the aerosol size distribution with four lognormal modes (Aitken, accumulation, coarse, and primary carbon modes; Liu et al., 2016). MAM4 tracks six aerosol species, and these are distributed over the four modes as follows. The Aitken mode consists of dust, sulfate, secondary organic aerosol (SOA), and sea salt. The accumulation mode includes sulfate, SOA, sea salt, primary organic matter (POM), black carbon (BC), and dust. The coarse mode contains sulfate, dust, and sea salt. The primary carbon mode contains only BC and POM, which are supplied by primary aerosol emissions.

The choice of modes in MAM4 is motivated by the desire to treat the microphysical aging of the primary carbonaceous aerosols in the atmosphere (Liu et al., 2016) similar to other modal models used in regional or global models (Riemer et al., 2003; Vogel et al., 2009). In MAM4, mass and number concentrations of BC and POM in the primary carbon mode are transferred to the accumulation mode by the processes of intermodal coagulation and condensation of SOA and sulfuric acid onto the primary carbon mode. The accumulation mode then represents aged BC and POM, as these species are internally mixed with other aerosol species. The MAM4 treatment of aging is critical for improving the long-range transport of carbonaceous aerosols to remote regions such as the polar region, which suffered from a low bias in a prior version of the model when only three internally mixed modes were used (Liu et al., 2016).

We ran the model for the year 2011 with 6 years (2005–2010) of spinup. The simulation was conducted at a resolution of 0.9^∘ latitude by 1.25^∘ longitude along with emission inventories from CMIP6 emissions (Emmons et al., 2020). We stored the instantaneous outputs every 3 h during the simulation, which yields 2920 timestamps for each surface-layer grid cell for the entire year of simulation time. The surface layer was chosen to be in line with the PartMC-MOSAIC model scenarios that were used as training data for the ML models of mixing state indices (see Sect. 3) and which were designed to represent conditions in the planetary boundary layer.

3 Aerosol mixing state indices: definition and calculation

3.1 Particle-based aerosol mixing state index

The mixing state index χ (Riemer and West, 2013) quantifies where an aerosol population lies on the continuum from external to internal mixing – that is, how spread out the chemical species are over an aerosol population. We will focus here on the mixing state of submicron aerosols (PM_1.0) due to their relevance for light scattering and absorption (Wang et al., 2015) and their contribution to CCN formation (Asmi et al., 2011; Pierce et al., 2015; Yu and Luo, 2009).

To summarize, the mixing state index χ is given by the affine ratio of the average particle species diversity, D_α, and bulk population species diversity, D_γ, as

$$\begin{array}{}\text{(1)}& \mathit{\chi }={\displaystyle \frac{D_{\mathit{\alpha }}-\mathrm{1}}{D_{\mathit{\gamma }}-\mathrm{1}}}.\end{array}$$

The diversities D_α and D_γ are calculated as follows. First, the per-particle mixing entropies H_i are determined for each particle by

$$\begin{array}{}\text{(2)}& H_i=\sum _{a=\mathrm{1}}^A-p_i^a\ln p_i^a.\end{array}$$

Here, A is the number of distinct aerosol species and $p_i^a$ is the mass fraction of species a in particle i. These values are then averaged (mass-weighted) over the entire population to obtain the average particle species diversity D_α by

$$\begin{array}{}\text{(3)}& {\displaystyle }{H}_{\mathit{\alpha }}=\sum _{i=\mathrm{1}}^{N_{\mathrm{p}}}{p}_i{H}_i,\text{(4)}& {\displaystyle }{D}_{\mathit{\alpha }}=e^{H_{\mathit{\alpha }}},\end{array}$$

where N_p is the total number of particles in the population and p_i is the mass fraction of particle i in the population. Finally, the bulk diversity D_γ is calculated as

$$\begin{array}{}\text{(5)}& {\displaystyle }{H}_{\mathit{\gamma }}=\sum _{a=\mathrm{1}}^A-p^a\ln p^a,\text{(6)}& {\displaystyle }{D}_{\mathit{\gamma }}=e^{H_{\mathit{\gamma }}},\end{array}$$

where p^a is the bulk mass fraction of species a in the population.

Note that the definition of “species” for calculating χ is based on application needs. It can be based on operationally defined chemical species (Healy et al., 2014; Ye et al., 2018), elemental composition (Fraund et al., 2017; O'Brien et al., 2015), or species groups such as volatile and nonvolatile species (Dickau et al., 2016) or hygroscopic and non-hygroscopic species (Ching et al., 2017; Hughes et al., 2018). Other possibilities include the propensity for aerosols to undergo heterogeneous reactions, quantified by the heterogeneous reaction rate coefficient for a specific reaction. In this paper we consider three different definitions of χ, which we explain in more detail in Sect. 3.2.

3.2 Mode-based aerosol mixing state index

The framework laid out in Sect. 3.1 can be easily generalized to a modal modeling framework (see Fig. 1). The bulk mixing entropy, H_γ, and the bulk diversity, D_γ, can be calculated using the bulk mass fractions, p^a, of species a from the MAM4 simulation and Eqs. (5) and (6). To calculate the average particle mixing entropy, H_α, and the average particle species diversity, D_α, we use

$$\begin{array}{}\text{(7)}& {\displaystyle }{H}_m=\sum _{a=\mathrm{1}}^A-p_m^a\ln p_m^a,\text{(8)}& {\displaystyle }{H}_{\mathit{\alpha }}=\sum _{m=\mathrm{1}}^M{p}_m{H}_m,\text{(9)}& {\displaystyle }{D}_{\mathit{\alpha }}=e^{H_{\mathit{\alpha }}},\end{array}$$

where $p_m^a$ is the mass fraction of species a in mode m, p_m is the mass fraction of mode m in the population, and H_m represents the per-mode mixing entropies. Finally, the mixing state index, χ, can be calculated using Eq. (1). Note that Eqs. (7) and (8) are analogous to Eqs. (2) and (3). A detailed derivation of these equations is provided in the Appendix A.

Figure 1 Illustration of the mode-based calculation of the aerosol mixing state index. The coarse mode is removed because only modes dominated by submicron particles are used for calculations. Note that the Aitken mode mass fraction is very low compared to the other modes and the caption does not obscure any data.

In this study, we consider the mixing states of submicron aerosols including the Aitken, accumulation, and primary carbon modes, and we do not include the coarse mode because the coarse particles are above 1 µm. Since the mixing entropies are mass-weighted (rather than number-weighted), the mixing state index is more representative of the modes with the larger particles, i.e., the accumulation and primary carbon modes.

3.3 Grouped surrogate species

Here we compare and contrast the aerosol mixing state indices defined in three different ways, namely based on the mixing of optically absorbing and non-absorbing species (χ_o), based on the mixing of primary carbonaceous and non-primary carbonaceous species (χ_c), and based on the mixing of hygroscopic and non-hygroscopic species (χ_h). Table 1 shows the definitions of these aerosol mixing state indices.

Table 1 Aerosol mixing state index definitions. Six aerosol species (bc: black carbon, dst: dust, ncl: sea salt, pom: primary organic matter, soa: secondary organic aerosol, so4: sulfate) are used in calculating the aerosol mixing state indices based on different species groupings. The mixing state indices χ_o, χ_c, and χ_h are based on two grouped surrogate species.

For χ_o, we considered two surrogate species: black carbon (strongly absorbing, assigned a mass absorption coefficient in CESM2 at 533 nm and 0 % RH of 8.144 m² g⁻¹) and the five other aerosol species grouped together (less absorbing or non-absorbing, with mass absorption coefficients in CESM2 at 533 nm and 0 % RH of 0.1442, $\mathrm{9.975}\times {\mathrm{10}}^{-\mathrm{2}}$, $\mathrm{4.703}\times {\mathrm{10}}^{-\mathrm{2}}$, $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$, and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ m² g⁻¹ for POM, SOA, dust, sea salt, and sulfate, respectively). Thus, a lower value in χ_o refers to the case where the strongly absorbing species black carbon (Yang et al., 2009) and the sum of the other species (termed “non-absorbing” here for convenience) are more externally mixed.

The index χ_c is motivated by the primary carbon treatment of MAM4, where the primary particulate organic matter and black carbon are assigned to a separate primary carbon mode (Liu et al., 2016). A lower value in χ_c refers to the situation where the primary carbonaceous species and all other species exist separately in different particles.

Similarly, χ_h was also calculated from two surrogate species. We combined black carbon, primary organic matter, and dust as one surrogate species, given their comparatively lower hygroscopicities (kappa values of ∼0, ∼0, and 0.068, respectively). Accordingly, NaCl (1.16), SOA (0.14), and sulfate (0.507) were grouped as the other surrogate species. Here, a lower value in χ_h represents the case where hygroscopic and non-hygroscopic species tend to be present in separate particles.

4 Machine-learned models of mixing state indices

Aerosol mixing state indices can be calculated directly using particle-resolved modeling, but this comes with large computational costs. Alternatively, Zheng et al. (2021) developed ML models, which integrate machine learning and particle-resolved aerosol simulations to estimate aerosol mixing state indices. To generate the training and testing data sets for developing such ML models, an ensemble of particle-resolved model scenarios was created using the particle-resolved model PartMC-MOSAIC (Riemer et al., 2009; Zaveri et al., 2008). In brief, PartMC-MOSAIC simulates individual aerosol particles within a representative volume of air, including stochastic coagulation, particle-phase thermodynamics, gas- and particle-phase chemistries, and dynamic gas–particle mass transfer. Thus, the composition of the individual particles within a population evolves dynamically, and assumptions about mixing state are not necessary.

The strategy to generate the data was to vary the input parameters (45 in total) for the PartMC-MOSAIC model, including primary emissions of different aerosol types (e.g., carbonaceous aerosol and dust emissions, including contribution from Aitken mode, accumulation mode, and coarse mode size ranges), primary emissions of gas phase species (e.g., SO₂, NO₂, and various volatile organic compounds), and meteorological parameters (see Table 1 in Zheng et al., 2021, for more information). For instance, to vary the gas emissions, scaling factors were sampled from 0 % to 200 % for different gas species, based on the emission rates in Riemer et al. (2009). A Latin hypercube sampling approach was employed to sample the parameter space efficiently for the training and testing data sets. We note that new particle formation and growth was not simulated explicitly, but Aitken mode sulfate particles were introduced into the simulation by emission for a subset of scenarios (Zheng et al., 2021) as a proxy for having particles present that originate from new particle formation. While PartMC-MOSAIC includes the process of new particle formation (Tian et al., 2014), the reason for this simplification was that considerable uncertainty exists regarding the subsequent growth of the freshly nucleated particles (Kulmala et al., 2014), which poses a challenge for a highly detailed aerosol model such as PartMC-MOSAIC. Errors in representing this particle type adequately may result in underestimating the abundance of BC-free particles in some regions (Zhang et al., 2017) and thereby overestimating the degree of internal mixture. This would imply that the error in the MAM4 simulations is even larger than currently indicated. Other processes that are not explicitly included in generating training data are aerosol removal by nucleation-scavenging and other cloud processes. However, for the purpose of this study, the emphasis is on the aerosol state, i.e., having a sufficiently comprehensive set of aerosol populations that can serve as training data, not necessarily that all the processes are included.

The ML models were derived by the machine learning algorithm eXtreme Gradient Boosting (XGBoost; Chen and Guestrin, 2016) from 45 000 particle populations. Each ML model was a tree-based ensemble model that could handle complex nonlinear interactions and collinearity among features. The hyperparameters were determined by grid search with 10-fold cross-validation. The ML models can be expressed as

$$\begin{array}{}\text{(10)}& {\mathit{\chi }}_S(x,y,t)=f_S\left(A\right(x,y,t),G(x,y,t),E(x,y,t\left)\right),\end{array}$$

where ${\mathit{\chi }}_S(x,y,t)$ is the mixing state index (χ_o, χ_c, or χ_h) at location (x,y) in the model layer nearest the surface at time t, and f_S denotes the function for calculating the corresponding mixing state index χ_S. The set names A (aerosol), G (gas), and E (environmental) represent the predictors (features) used for predicting the mixing state index. The choice of features is determined by the overlap of variables that are present in both PartMC-MOSAIC and CESM2. Aerosol species include black carbon, mineral dust, sea salt, primary organic aerosol, secondary organic aerosol, and sulfate. Of note is that we used the bulk (not the per-mode) concentrations of submicron aerosol species as the features. The gas species include dimethyl sulfide, hydrogen peroxide, sulfuric acid, ozone, semi-volatile organic gas, and sulfur dioxide. The environmental variables are air temperature, relative humidity, and solar zenith angle. Table 2 shows the performance of the ML models when predicting the mixing state indices. The mixing state calculation in this study was purely based on the above six aerosol species (excluding other aerosol species) for a fair comparison with the mode-based aerosol mixing state index, which resulted in slightly different performance of the ML model compared to Zheng et al. (2021). The average error of the ML model (using the hold-out testing samples) is about 5 % for χ_o and 8 % for χ_c and χ_h (measured by mean absolute error).

Table 2 Predictive performance of the ML models using the testing data set. Metrics include the mean absolute error (MAE), root-mean-square error (RMSE), median absolute deviation (MAD), index of agreement (d; Willmott, 1981), Pearson correlation coefficient (PCC), and coefficient of determination (r²).

We would like to emphasize that this ML modeling framework cannot compensate for any biases that the global model (here CESM2) might have in simulating the quantities that serve as the features. Instead, what we can expect from this approach is that it provides the most likely mixing state associated with the species concentrations that CESM2 simulates.

5 Results

5.1 Quantitative comparison of mode-based and particle-based mixing state indices

Let ${\mathit{\chi }}_{S,t}^{\mathrm{ML}}$ and ${\mathit{\chi }}_{S,t}^{\mathrm{MAM}\mathrm{4}}$ denote the mixing state indices computed by the ML model and by the MAM4 model for each grid cell at timestamp t, respectively. The corresponding time-averaged values for a certain time interval and for each grid cell are $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_S^{\mathrm{ML}}}$ and $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_S^{\mathrm{MAM}\mathrm{4}}}$. Here we consider the full year as the time-averaging interval. An analysis of the seasonal variation of mixing state indices can be found in Zheng et al. (2021).

To compare the annual mean values, we calculated the mean difference ($\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\chi }}_S}$) and the mean absolute difference ($\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_S\right|}$) for each grid cell of the layer closest to the surface:

$$\begin{array}{}\text{(11)}& {\displaystyle }\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\chi }}_S}={\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T({\mathit{\chi }}_{S,t}^{\mathrm{MAM}\mathrm{4}}-{\mathit{\chi }}_{S,t}^{\mathrm{ML}})=\stackrel{\mathrm{‾}}{{\mathit{\chi }}_S^{\mathrm{MAM}\mathrm{4}}}-\stackrel{\mathrm{‾}}{{\mathit{\chi }}_S^{\mathrm{ML}}},\text{(12)}& {\displaystyle }\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_S\right|}={\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T\left(\right|{\mathit{\chi }}_{S,t}^{\mathrm{MAM}\mathrm{4}}-{\mathit{\chi }}_{S,t}^{\mathrm{ML}}\left|\right),\end{array}$$

where the subscript S refers to the mixing state index (o, c, or h), and the total number of timestamps is T=2920. Since it only makes sense to quantify mixing state when at least two species are present in a given location, areas where the mass fraction of any one surrogate species was higher than 99 % for χ_o (due to the low mass fraction of black carbon) and 97.5 % for χ_c and χ_h were ignored for the calculation and appear as hatched areas in Fig. 3. We will first discuss the overall probability density functions of these quantities (Fig. 2) and then their spatial distributions (Fig. 3).

Figure 2 Probability density functions of annual averaged mixing state indices using the MAM4 model and ML model. The thin black lines refer to their mean values.

Figure 3 Global distribution of annually averaged mixing state indices (χ_o, χ_c, and χ_h) using the ML model, MAM4 model, their mean difference ($\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathit{\chi }}$), and mean absolute difference ($\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }\mathit{\chi }\right|}$). Areas are hatched where the mass fraction of any one surrogate species was higher than 99 % for χ_o (due to the low mass fraction of black carbon) and 97.5 % for χ_c and χ_h.

Figure 2 shows the probability density functions of the annual averaged mixing state indices computed by the ML model ($\stackrel{\mathrm{‾}}{{\mathit{\chi }}_S^{\mathrm{ML}}}$), by MAM4 ($\stackrel{\mathrm{‾}}{{\mathit{\chi }}_S^{\mathrm{MAM}\mathrm{4}}}$), their average difference ($\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\chi }}_S}$), and their average absolute difference ($\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_S\right|}$) for each surface-layer grid cell. The results show large discrepancies in mixing state indices between the ML model and the MAM4 model, without a clear relationship between them (see Fig. 2d–f).

The annual average of the mixing state index χ_o estimated by the ML model, $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}}$, ranged between 55 % and 96 %, with a mean of 73 %. Calculated by the MAM4 model, $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}}$ varied spatially from 46 % to 99.76 %, with a higher mean of 86 %. The similar mean values of $\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{o}}}$ (14 %) and $\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{o}}\right|}$ (18 %) were caused by higher values in ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$ compared to ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}$, which is confirmed below with Fig. 3. The averaged mixing state index $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{ML}}}$ ranged between 31 % and 84 % with a mean of 54 %, while $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{MAM}\mathrm{4}}}$ had a wider range (from 9 % to 99.81 %) with a mean (of 58 %). Similarly, $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{h}}^{\mathrm{ML}}}$ ranged from 21 % to 81 % with a mean of 58 %, while $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{h}}^{\mathrm{MAM}\mathrm{4}}}$ varied between 10 % and 99.85 % with a mean of 63 %. The large discrepancy between the mean difference (4.8 % for $\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{c}}}$ and 4.7 % for $\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{h}}}$) and mean absolute difference (30 % for $\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{c}}\right|}$ and 38 % for $\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{h}}\right|}$) indicates that the errors in χ_c and χ_h were symmetric (positive and negative) but large. The maximal errors in $\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{c}}\right|}$ and $\stackrel{\mathrm{‾}}{\left|\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{h}}\right|}$ between the two methods were up to 59 and 76 percentage points, respectively.

The implications of these discrepancies are more easily discussed with Fig. 3, which illustrates the global spatial distribution of annually averaged mixing state indices predicted by the ML model (first column), MAM4 (second column), their mean difference (third column), and their mean absolute difference (fourth column). The differences in mixing state indices between the ML model and MAM4 varied strongly across the globe.

High values of $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}}$ occurred in the continental regions (77 %) compared to oceans (69 %). Specifically, the ML model predicted high values for χ_o in central Africa (20^∘ S–15^∘ N, 12–30^∘ E), the Arctic (66.5–90^∘ N), and southern Asia (5–38^∘ N, 60–90^∘ E). These are also the regions with relatively larger mass fractions of black carbon (∼5 %; see Fig. A2). The mixing state index ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$ showed a higher degree of internal mixing over the globe (with a median of 90 %) compared to the ML model. The only exceptions were oceans in the Northern Hemisphere at the mid-latitudes (45–60^∘ N, dominated by sea salt, sulfate, and secondary organic aerosol in the accumulation mode) and Antarctica (66.5–90^∘ S, dominated by sea salt and sulfate in the accumulation mode as well as sulfate in Aitken mode), where $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}}$ was 75 %. Qualitatively, the MAM4 model captured the trend that areas with high black carbon concentration (defined here as concentrations above the 95 % percentile) tended to have higher χ_o values.

The ML model estimate $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{ML}}}$ suggested a rather homogeneous spatial distribution of the annually averaged mixing state, with values of approximately 50 %. Compared to $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{ML}}}$, $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{MAM}\mathrm{4}}}$ values were lower (primary carbonaceous aerosol more externally mixed) at high latitudes and higher at low and mid-latitudes (primary carbonaceous aerosol more internally mixed). Note that, while $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{MAM}\mathrm{4}}}$ values were similar in the Arctic and Antarctic, the abundance of primary carbonaceous species was predicted to be higher in the Arctic compared to the Antarctic (see Fig. A2).

The spatial distributions of $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{h}}^{\mathrm{MAM}\mathrm{4}}}$ were similar to $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{MAM}\mathrm{4}}}$. That is, the MAM4 model predicted that the hygroscopic species and non-hygroscopic species were more externally mixed at high latitudes and more internally mixed at low latitudes. In contrast, the spatial distribution of $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{h}}^{\mathrm{ML}}}$ shows qualitative differences compared to $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{ML}}}$ in two aspects. First, $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{h}}^{\mathrm{ML}}}$ was higher than $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{ML}}}$ at high latitudes, meaning that hygroscopic species and non-hygroscopic appeared more internally mixed than primary carbonaceous and non-carbonaceous species in this region. Second, areas over the North Atlantic Ocean (0–20^∘ N, 20–45^∘ W), southern Africa (5–32^∘ S, 5–20^∘ E), and Australia (10–30^∘ S, 100–140^∘ E) appeared rather externally mixed. These are areas where mineral dust is the dominant aerosol species (see Fig. A2).

These two facts lead to the overall finding that χ_h exhibits the largest differences between the two methods. This applies especially to regions where mineral dust was the dominant aerosol species, which points to an important structural issue of the four-mode setup used in MAM4. While the ML model predicted a more external mixture in these regions (dust externally mixed from sea salt and other species), the MAM4 model could not represent this because the accumulation mode included all six aerosol species in an internal mixture. Figure 4 illustrates the relationship of the mean absolute difference of χ_h and the mass fraction of dust for all model grid points. It confirms that grid points with large dust mass fractions were associated with larger mean absolute differences in χ_h. These results confirm the tradeoff discussed in Liu et al. (2012): MAM3 (and MAM4 in Liu et al., 2016) intentionally combines dust and sea salt in the same mode to reduce the computational burden; however, this simplification does not always realistically reflect the aerosol mixing state in the ambient atmosphere.

It is interesting to note that the areas where sea salt is present, but not dust, are not associated with large errors, even though sea salt – just like mineral dust – is a primary aerosol type. The reason for this lies in our surrogate species definitions (Table 1) for computing the mixing state index. Based on our mixing state definitions, sea salt, secondary organic aerosol, and sulfate are always grouped together. Therefore, none of the mixing state indices as defined here tell us how externally mixed sea salt is when it is considered as a single aerosol type.

Figure 4 Dependence of mean absolute difference of χ_h on dust mass fractions for all model grid points.

Figure 5 further demonstrates the zonal mean annual aerosol mixing state indices, highlighting that differences between χ_c and χ_h tended to be zonally structured, where the MAM4 model overestimated at low latitudes, while it underestimated at high latitudes relative to the ML model. In contrast, the MAM4 model overestimated χ_o at all latitudes north of 60^∘ S.

Figure 5 Zonal mean annual aerosol mixing state indices (a) χ_o, (b) χ_c, and (c) χ_h using the MAM4 model and ML model. The bands refer to the standard deviation.

5.2 Interpretation of findings

From Sect. 5.1, the following picture emerges: MAM4 overestimates the mixing state index χ_o except in regions at high latitudes in the Southern Hemisphere. At the same time, χ_c and χ_h are overestimated at low latitudes to mid-latitudes and underestimated at high latitudes. These findings point towards too rapid a transfer from the carbonaceous mode to the accumulation mode at low latitudes to mid-latitudes and too slow a transfer at high latitudes.

To conceptually illustrate these relationships, here we use χ_o and χ_c as examples and contrast the conditions for high and low latitudes. Figure 6a–f show conditions representative of high latitudes. A grid cell sampled from the CESM2/MAM4 simulation (73^∘ N, 151^∘ W) contains 15 % BC and 37 % POM, distributed over the accumulation and primary carbon mode as shown in Fig. 6a and d. The corresponding value for χ_o is 80 %. Figure 6b depicts particle population that was sampled from the MAM4 population in Fig. 6a. All particles, except for the smallest ones (corresponding to Aitken mode particles), contain BC, which results in the relatively high mixing state index value for χ_o. Note that in MAM4 BC is not included in the Aitken mode by definition. Considering the same particle population, but now evaluating the mixing state metric χ_c, which quantifies the degree of mixing of primary carbon and other species, yields the following observation. The entire primary carbon mode, by definition, consists of POM and BC, which results in an appreciable number of particles that contain only primary carbon (BC + POM), giving a mixing state index χ_c of only 27 %.

We now compare the MAM4-sampled particle populations above to particle populations that were sampled from our PartMC scenario library. We searched for populations with similar mass fractions of BC and POM as in the MAM4 populations and that were simulated at a similar latitude as the grid point location of the CESM2/MAM4 model output. Figure 6c shows that the PartMC results have comparatively more BC-free particles, and Fig. 6f shows that comparatively more particles are mixtures of primary carbon and other species. Overall, this means that in MAM4 BC appears too internally mixed (because irrespective of whether BC is placed in the primary carbon or accumulation mode, it is by design mixed with other species) and that at high latitudes the primary carbon mode is not transferring mass to the accumulation mode as quickly as is the case in PartMC simulations.

The reason why MAM4 behaves in this way can be explained by the aging process treatment in MAM4. Aging in MAM4 is formulated using a threshold criteria. That is, BC and POA mass is transferred from the primary carbon mode to the accumulation mode when a certain threshold of sulfate and SOA has condensed. In MAM4 this threshold is set to a relatively large value. This is done to prevent BC from being removed too quickly by wet deposition – because the primary carbon mode has a lower hygroscopicity than the accumulation mode and thus a lower wet scavenging efficiency – thereby counteracting a low bias in BC concentrations in the Arctic regions. From Wang et al. (2018) we already know that using such a high threshold may not be appropriate. However, the global model also has biases in other processes that contribute to the low BC bias in the Arctic, and setting the threshold to a high value compensates for these errors. Our results are a reflection of this fact. While adjusting the threshold criteria in MAM4 to a lower value may improve the agreement with the ML simulations in some regions, it may deteriorate the overall results in other areas. This is a good example how structural uncertainty manifests itself, namely by the fact that adjusting a parameter does not fundamentally fix the issue.

Figure 6g–l show conditions representative of low latitudes. A grid cell sampled from the CESM2/MAM4 simulation (20^∘ N, 120^∘ E) contains 11 % BC and 24 % POM, distributed over the accumulation and primary carbon mode as shown in Fig. 6g and j, with most of the mass in the accumulation mode. The corresponding value for χ_o is therefore 99 %, an almost complete internal mixture. For the same reason, χ_c is also very high. Similarly to the high-latitude case, Fig. 6i and l show that the comparable PartMC population has comparatively more BC-free particles and more particles that contain very low amounts of primary carbonaceous material, leading to lower values of both χ_o and χ_c compared to the MAM4 results.

Figure 6 Illustration to explain the differences in mixing state representation between MAM4 and the ML model at high and low latitudes.

5.3 Comparison to observational data

The question that arises from Sect. 5.1 and 5.2 is of course the following: which spatial distribution of aerosol mixing state reflects reality more closely? The validation of simulated mixing state indices with observational data is still challenging since per-particle mass fractions of species are required for calculating the mixing state indices (see Sect. 3.1). These are in principle obtainable from in situ deployments of single-particle mass spectrometers or by using electron microscopy techniques, but their quantitative derivation comes with challenges and is not routinely done, so that only very few data sets exist that allow for a meaningful comparison (Riemer et al., 2019). Keeping these limitations in mind, Zheng et al. (2021) reported a qualitative comparison of available measurements of mixing state metrics in locations in developed countries (Paris, France; Pittsburgh, USA; various locations in Japan) (Healy et al., 2014; Ye et al., 2018; Ching et al., 2019) with seasonally averaged results from the ML model based on particle-resolved simulations. This showed that the ML model was able to capture the range of values that is consistent with the observations.

We further compared the ML model estimates using recent observations from China. Specifically, we compared ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}$ and ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$ with χ values from Taizhou (Zhao et al., 2021) and Beijing (Yu et al., 2020) derived from Single Particle Soot Photometer (SP2) measurements. For both locations, ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$ overestimated the observed χ values, while ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}$ was in the range of the observations. Specifically, the χ measured at a suburban site Taizhou from 26 May to 18 June 2017 ranged from 62 % to 82 %. During the same time period (but in the year 2011), the values of ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}$ were between 63 % and 84 %, while ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$ was between 84 % and 96 %. The χ values at the urban site of Beijing ranged between 55 % and 70 % in winter (from 10 November to 10 December 2016) and varied between 60 % and 75 % in summer (from 18 May to 25 June 2017). Using our simulations of the year 2011, ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}$ varied from 60 % to 88 % in winter and from 59 % to 83 % in summer. As a comparison, ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$ ranged from 92 % to 97 % in winter and from 87 % to 95 % in summer. A caveat when comparing ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}$ and ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$, respectively, with the observations reported in Zhao et al. (2021) and Yu et al. (2020) is that the definition of ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{ML}}$ and ${\mathit{\chi }}_{\mathrm{o}}^{\mathrm{MAM}\mathrm{4}}$ included BC-free particles, while the χ values in the measurements by Zhao et al. (2021) and Yu et al. (2020) were calculated only considering the subpopulation of BC-containing particles. This might introduce a bias in the mixing state index between the χ_o index used in this paper and the observations (depending on the fraction of the BC-free particles present at any given location).

We can also relate our χ_o index qualitatively to the SP2 measurements in the Finnish Arctic during winter 2011–2012 (Raatikainen et al., 2015). Although this study did not provide quantitative mixing state index calculations, it is an important finding that BC-containing particles (with various amounts of coatings) co-existed with BC-free particles. As we saw in Sect. 5.2, this condition can easily be represented with a particle-resolved approach. However, the modal model with modes configured as in MAM4 puts black carbon in all accumulation-sized particles (Fig. 6), which is not consistent with the observations.

6 Conclusions

In this paper we present a framework for evaluating the error in submicron aerosol mixing state induced by aerosol representation assumptions, which is one of the important contributors to structural uncertainty in aerosol models. We quantitatively compared mixing state indices for submicron aerosol predicted by the modal model MAM4 within the global model CESM to a machine-learned model based on high-detail particle-resolved simulations. We focused on the mixing of optically absorbing and non-absorbing species (χ_o), the mixing of primary carbonaceous with other aerosol species (χ_c), and the mixing of hygroscopic and non-hygroscopic species (χ_h).

For χ_o, the MAM4 modal representation generally overestimated the degree of mixing of BC with other aerosol species. This overestimation is due to the fact that MAM4's choice of modes does not allow for representing BC-free particles in the accumulation and primary carbon modes. This is in contrast to field observations by Brown et al. (2021), which showed that BC and POM may be externally mixed near sources. The implication of this is that, if optical properties are calculated based on the aerosol composition, absorption will be overestimated.

For χ_c and χ_h, the error tended to be zonally structured, where the MAM4 model overestimated the mixing state indices at low latitudes and underestimated them at high latitudes compared to the ML model. This behavior could be explained by modeling choices in MAM4, in particular that (1) BC is always emitted with POM, (2) no BC-free particles exist in the submicron modes, and (3) dust is always internally mixed with other aerosol species.

Mixing state is an important emergent property that affects the aerosol radiative forcing and aerosol–cloud interactions, but it is not easy to constrain this property globally. To the best of our knowledge, this is the first study that evaluated the spatial distribution of aerosol mixing state as predicted by a global model. Since errors in mixing state predictions propagate into errors in aerosol climate impacts, our findings provide a framework and reference for Earth system model developers and users regarding simulation reliability. For example, this framework can be used to (1) quantify model bias in simulating mixing state in different regions, identifying model structural deficiencies, and (2) provide insights into potential improvements of model process representations for a more realistic simulation of aerosols.

Appendix A: Derivation of mode-based aerosol mixing state index

Table A1 details the notation for aerosol mass and mass fractions to calculate H_α using modal information.

To explain how to obtain Eqs. (7) and (8) from Eqs. (2) and (3), let us assume that each mode m contains N_m particles and the number of species in the population is A. The mixing entropy of particle i in mode m, H_m,i, is given by

$$\begin{array}{}\text{(A1)}& H_{m,i}=\sum _{a=\mathrm{1}}^A-p_{m,i}^a\ln p_{m,i}^a.\end{array}$$

The average particle mixing entropy of the entire population (summed over all modes), H_α, is

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{H}_{\mathit{\alpha }}& =\sum _{m=\mathrm{1}}^M\sum _{i=\mathrm{1}}^{N_m}{p}_{m,i}{H}_{m,i}=\underset{m=\mathrm{1}}{\underbrace{p_{\mathrm{1},\mathrm{1}}{H}_{\mathrm{1},\mathrm{1}}+p_{\mathrm{1},\mathrm{2}}{H}_{\mathrm{1},\mathrm{2}}\mathrm{\cdots }+p_{\mathrm{1},N_{\mathrm{1}}}{H}_{\mathrm{1},N_{\mathrm{1}}}}}\\ & +\mathrm{\cdots }+\underset{m=M}{\underbrace{p_{M,\mathrm{1}}{H}_{M,\mathrm{1}}+p_{M,\mathrm{2}}{H}_{M,\mathrm{2}}\mathrm{\cdots }+p_{M,N_M}{H}_{M,N_M}}}.\end{array}\end{array}$$

Given that each mode is assumed to be internally mixed, particles within the same mode have the same composition, and we have

$$\begin{array}{}\text{(A3)}& p_{m,i}^a={\displaystyle \frac{{\mathit{\mu }}_{m,i}^a}{{\mathit{\mu }}_{m,i}}}={\displaystyle \frac{{\mathit{\mu }}_m^a}{{\mathit{\mu }}_m}}=p_m^a.\end{array}$$

This results in

$$\begin{array}{}\text{(A4)}& H_{m,i}=\sum _{a=\mathrm{1}}^A-p_{m,i}^a\ln p_{m,i}^a=\sum _{a=\mathrm{1}}^A-p_m^a\ln p_m^a=H_m.\end{array}$$

Therefore, based on Eq. (A4) and the fact that $p_m={\sum }_{i=\mathrm{1}}^{N_m}{p}_{m,i}$, Eq. (A2) can be rewritten as

$$\begin{array}{}\text{(A5)}& H_{\mathit{\alpha }}=\underset{m=\mathrm{1}}{\underbrace{p_{\mathrm{1}}{H}_{\mathrm{1}}}}+\mathrm{\cdots }+\underset{m=M}{\underbrace{p_M{H}_M}}=\sum _{m=\mathrm{1}}^M{p}_m{H}_m.\end{array}$$

With the mode-based H_α, the other mixing state quantities can be computed as described in Sect. 3.2.

Table A1 Aerosol mass and mass fraction definition and notation. The number of modes is M (M=3 for MAM4 without the coarse mode), the number of particles in mode m is N_m, and the number of species is A.

Figure A1 Aerosol species mixing ratio (µg kg⁻¹). Accumulation mode: a1; Aitken mode: a2; primary carbon mode: a4. The coarse mode (a3) is not used in this study and therefore omitted in this figure. Black carbon: bc; dust: dst; sea salt: ncl; primary organic matter: pom; secondary organic aerosol: soa; sulfate: so4.

Figure A2 Fraction of aerosol species mixing ratio (%). Accumulation mode: a1; Aitken mode: a2; primary carbon mode: a4. The coarse mode (a3) is not used in this study and therefore omitted in this figure. Black carbon: bc; dust: dst; sea salt: ncl; primary organic matter: pom; secondary organic aerosol: soa; sulfate: so4.