Estimating Submicron Aerosol Mixing State at the Global Scale With Machine Learning and Earth System Modeling

This study integrates machine learning and particle‐resolved aerosol simulations to develop emulators that predict submicron aerosol mixing state indices from the Earth system model (ESM) simulations. The emulators predict aerosol mixing state using only quantities that are predicted by the ESM, including bulk aerosol species concentrations, which do not by themselves carry mixing state information. We used PartMC‐MOSAIC as the particle‐resolved model and NCAR's CESM as the ESM. We trained emulators for three different mixing state indices for submicron aerosol in terms of chemical species abundance (χa), the mixing of optically absorbing and nonabsorbing species (χo), and the mixing of hygroscopic and nonhygroscopic species (χh). Our global mixing state maps show considerable spatial and seasonal variability unique to each mixing state index. Seasonal averages varied spatially between 13% and 94% for χa, between 38% and 94% for χo, and between 20% and 87% for χh with global annual averages of 67%, 68%, and 56%, respectively. High values in one index can be consistent with low values in another index depending on the grouping of species and their relative abundance, meaning that each mixing state index captures different aspects of the population mixing state. Although a direct validation with observational data has not been possible yet, our results are consistent with mixing state index values derived from ambient observations. This work is a prototypical example of using machine learning emulators to add information to ESM simulations.

To quantify mixing state, Riemer and West (2013) introduced a metric, the mixing state index χ, based on diversity measures derived from the information-theoretic (Shannon) entropy of the chemical species distribution among particles, with per-particle mass fractions as the fundamental quantities. The mixing state index χ varies between 0% (for completely external mixtures) and 100% (for completely internal mixtures) for any given aerosol. The metric has been applied to field observations in different environments, for example, Paris during the MEGAPOLI campaign (Healy et al., 2014), in Northern California during CARES (O'Brien et al., 2015), in central London (Giorio et al., 2015), and in Pittsburgh, PA (Ye et al., 2018). It has provided useful insights into the processes that govern diurnal changes in mixing state, and mixing state changes related to air mass origin. It can also be used for error quantification as shown in Ching et al. (2017), , and Ching and Kajino (2018). Ching et al. (2017) quantified errors in CCN concentration prediction when assuming a fully internal mixture compared to using the realistic mixing state predicted by a particle-resolved model. In a follow-up study,  applied this concept to errors in the prediction of black carbon (BC) mass fraction that is incorporated into cloud droplets by nucleation-scavenging. Considering the deposition of soot particles in the human respiratory system, Ching and Kajino (2018) demonstrated that the mixing state index is also an important metric for health impact evaluation.
The most direct way of predicting aerosol mixing state and its index χ is to employ a particle-resolved model. This modeling approach allows for tracking the composition of individual particles and therefore χ can be directly calculated. However, particle-resolved modeling is extremely computationally expensive and thus not practical for use in large-scale models. Alternatively, Hughes et al. (2018) developed a method that uses the output of a large ensemble of particle-resolved box model simulations combined with machine learning techniques to train a model of the mixing state metric χ. This lower-order model for χ uses as inputs only variables known to the global climate model of interest (in Hughes et al. [2018] it was GEOS-Chem-TO-MAS, which uses a sectional aerosol modeling approach assuming an internal mixture within each size bin). The outcomes of this procedure were global maps of χ based on GEOS-Chem-TOMAS data.
In this work, we revisit the approach by Hughes et al. (2018), but extend the analysis by defining the mixing state index in three different, but complementary, ways (see Section 3): in terms of chemical species abundance, in terms of the mixing of optically absorbing and nonabsorbing species, and in terms of the mixing of hygroscopic and nonhygroscopic species. We discuss their relationship with each other and their interpretation in Section 5.
We used the Community Earth System Model (CESM; Danabasoglu et al., 2020;Hurrell et al., 2013) in conjunction with the Modal Aerosol Module (MAM4; Liu et al., 2016) as a large-scale model that provides the features for the machine learning procedure. We compared and contrasted the spatial distribution and seasonal variation of aerosol mixing state indices at a global scale.

Ensemble of Particle-Resolved Model Scenarios
To generate the training data, we ran an ensemble of particle-resolved model scenarios using the Part-MC-MOSAIC. This section gives an overview of the model and the design of the scenario ensemble.

PartMC-MOSAIC Model Description
PartMC is a particle-resolved tool to model the aerosol mixing state and its impacts under various meteorological and environmental conditions (Riemer et al., 2009). In brief, this box model simulates individual aerosol particles within a representative volume of air, including stochastic coagulation, gas-and particle-phase chemistries, particle-phase thermodynamics, and dynamic gas-particle mass transfer. The Part-MC algorithm has storage cost proportional to the number of particles, computational cost for evaporation/ condensation proportional to the number of particles, and computational cost for coagulation proportional to the number of coagulation events. PartMC has been coupled with the state-of-the-art aerosol chemistry model MOSAIC (Model for Simulating Aerosol Interactions and Chemistry; Zaveri et al., 2008) to provide the chemistry modules. MOSAIC consists of four computationally efficient modules: (1) the gas-phase photochemical mechanism CBM-Z (Zaveri & Peters, 1999); (2) the multicomponent Taylor expansion method for estimating activity coefficients of electrolytes and ions in aqueous solutions (Zaveri, Easter, & Wexler, 2005); (3) the multicomponent equilibrium solver for aerosols (MESA) for intraparticle solid-liquid partitioning (Zaveri, Easter, & Peters, 2005); and (4) the adaptive step time-split Euler method for dynamic gas-particle partitioning over size-and composition-resolved aerosol (Zaveri et al., 2008). Secondary organic aerosol (SOA) from anthropogenic and biogenic precursors is included using the Secondary Organic Aerosol Model (Schell et al., 2001). Overall, MOSAIC treats all locally and globally important gas and aerosol species, including sulfate, nitrate, chloride, carbonate, ammonium, sodium, calcium, primary organic aerosol, SOA, BC, and inert inorganic mass (a surrogate species for mineral dust). The coupled model system, PartMC-MOSAIC, predicts number, mass, and full composition distributions without a priori assumptions about their evolution, which is suitable for use as a numerical benchmark of mixing state for more approximate models (Fierce et al., 2016(Fierce et al., , 2017Kaiser et al., 2014). Here, we used Part-MC-MOSAIC simulations to provide data for training and testing emulators that predict the mixing state indices on the global scale.

Design of Training and Testing Scenarios
The PartMC-MOSAIC training and testing scenarios provide a large number of particle populations with different mixing states aiming at covering the range of conditions found in different environments around the globe. Our strategy to generate this data was to vary the input parameters (primary emissions of different aerosol types, and primary emissions of gas phase species, which served as precursors for secondary aerosol, as well as meteorological parameters; 45 parameters in total), as listed in Table 1. To sample this parameter space efficiently, we employed a Latin Hypercube sampling approach to assemble input parameter ZHENG ET AL.    Hughes et al. (2018). All scenarios were run with 10,000 computational particles. They started at 6:00 a.m. local time and used a simulation time of 24 h, saving the hourly output.
We obtained the daily near-surface temperature from the CESM2 (0.9° by 1.25° resolution) fully coupled historical simulations (Danabasoglu, 2019). For each latitude and day of year (DOY) over the years from 1970 to 2014, we retrieved the maximum and minimum temperatures. For a given pair of latitude and DOY in PartMC-MOSAIC, the temperature was then sampled between the maximum and minimum temperatures at the nearest latitude in CESM2 with the same DOY. Gas initial conditions and emission rates were based on Riemer et al. (2009). We sampled the gas phase emissions using a scaling factor (0%-200%) for the gas phase species listed in Table 1.
We set up five blocks of scenarios (blocks A-E) that differed in the way the aerosol types were combined. Table 2 specifies these simulations blocks using abbreviations as defined in Table 1. Block A consisted of 800 scenarios that included carbonaceous aerosol, dust, and sea salt. We uniformly sampled the relative humidity for these scenarios from a range of 40%-100%, and latitude from a range of 90°S to 90°N. The latitude governs the solar zenith angle and the length of the day, which are important for the photochemical production of secondary aerosol. We included initial concentration (Omori et al., 2017), and emissions (Xu et al., 2016) of the biogenic trace gas dimethyl sulfide (DMS) for the scenarios with sea salt emission, as it is a precursor of sulfate aerosol in marine environments (Lana et al., 2011).
Block B had 400 scenarios and the same setup as Block A, but did not contain dust emissions. Block C had 400 scenarios that included only carbonaceous aerosol and dust (no sea salt). For the Block C scenarios, we uniformly sampled the relative humidity from 10% to 100%, and latitude from 70°S to 70°N. Block D comprised 200 scenarios and only contained carbonaceous aerosol. Otherwise, the setup was the same as for Block C. Block E contained carbonaceous aerosol and sea salt similar to Block B, but had additional sea salt and sulfate emissions in the Aitken mode size range. The purpose of Block E was to capture conditions where particles in the Aitken mode size range were present, including those resulting from new particle formation. The latitude range for the Block E scenarios was limited to high latitudes from 70°N-90°N and 70°S-90°S. This was motivated by the fact that, in these areas, CESM simulations showed the Aitken mode particles accounted for an appreciable fraction of submicron aerosol mass. Aitken mode sulfate particles were introduced into the PartMC-MOSAIC simulation by emission rather than by simulating the process of new particle formation and growth explicitly. While PartMC-MOSAIC includes the process of new particle formation (Tian et al., 2014), 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.
To create aerosol initial conditions with realistic mixing states we adopted the following approach: We performed a first set of runs, starting with the aerosol initial concentrations set to zero for all the simulations (the "initial runs"). We then repeated the same set of runs, but replaced the aerosol initial condition with a randomly sampled population from the initial runs (the "restart runs"). For the training and testing datasets, we only used the results from the restart runs. The size distribution parameters for the aerosol emissions were prescribed as log-normal, with geometric mean diameter (D g ) and geometric standard deviation (σ g ) depending on the particle type. Note. The variables D g , σ g , and E a refer to geometric mean diameter, geometric standard deviation, and particle emission flux, respectively.

Aerosol Mixing State Indices Definition
The mixing state index χ (Riemer & West, 2013) quantifies where an aerosol population lies on the continuum of external to internal mixing, that is, how "spread out" the chemical species are over an aerosol population, see Figure 2 for an illustration. We will focus here on the aerosol mixing state of submicron particles (PM 1.0 ) due to their impact on light scattering and absorption , and contribution to CCN formation (Asmi et al., 2011;Pierce et al., 2015;Yu & Luo, 2009).
The mixing state index χ is given by the affine ratio of the average particle species diversity (D α ) and bulk population species diversity (D γ ) as The diversities D α and D γ are calculated as follows. First, the per-particle mixing entropies H i are determined for each particle using the per-particle species mass fractions: Here, A is the number of distinct aerosol species and a i p is the mass fraction of species a in particle i. These values are then averaged (mass-weighted) over the entire population to find the average particle entropy, H α , and the average particle species diversity, D α , by ZHENG ET AL. , 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 entropy, H γ , and bulk diversity, D γ , is calculated as: .
Importantly, the definition of the "species" to calculate D α and D γ depends on the application or may be dictated by the instrumentation used to estimate per-particle mass fractions. For example, in some previous studies, elemental species (e.g., N, O, and C) have been used (Bondy et al., 2018;Fraund et al., 2017;O'Brien et al., 2015), while others used molecular species (e.g., SO 4 , NH 4 , and NO 3 ) (Healy et al., 2014;Lee et al., 2019;Ye et al., 2018). Yet other studies employed surrogate species where several species with similar properties were grouped together, such as hygroscopic/nonhygroscopic species (Ching et al., 2017;Hughes et al., 2018), or volatile/nonvolatile species (Dickau et al., 2016). The choice of which species to group as one surrogate species is then motivated by the science question under investigation.
Here, we compare and contrast the aerosol mixing state indices defined in three different ways, namely based on model chemical species abundance (χ a ), based on the mixing of optically absorbing and nonabsorbing species (χ o ), and based on the mixing of hygroscopic and nonhygroscopic species (χ h ). Table 3 shows the specific definitions of these aerosol mixing state indices. The index χ a was defined based on all six model aerosol species. For χ o , we considered two surrogate species, BC (absorbing) and all other aerosol species grouped together (assumed to be nonabsorbing). Thus, a lower value in χ o refers to the case where the absorbing species BC and the sum of all other (nonabsorbing) species are more externally mixed. Similarly, χ h was also calculated from two surrogate species. We combined BC, dust, and primary organic matter as one surrogate species, given their comparatively low hygroscopicities, and NaCl, SOA, and sulfate as the other surrogate species. Here, a lower value in χ h represents the case where hygroscopic and nonhygroscopic species tend to be present in separate particles. Figure 2 illustrates the concept of different aerosol mixing state indices corresponding to a particle population which consists of six particles and six chemical species that are present in equal amounts in the bulk.
Moving from left to right, the bulk composition for each population remains the same, but the population becomes more internally mixed. That is, the particles become more similar to each other. This is reflected in an increasing mixing state index, and this applies to all three mixing state indices individually. Another lesson is that the mixing state index depends on how we define the surrogate species. In our example, the population with a mixing state index of χ a = 20% has χ o = 63% and χ h = 67%. For this example, the population is quite externally mixed with respect to all chemical model species, but appears more internally mixed with respect to hygroscopicity or optical properties. A consequence of the definition of the surrogate species as supersets of the chemical model species is that, for the extreme cases of χ a = 100% and χ a = 0%, χ o and χ h are also 100% or 0%, respectively. The reverse, however, is not true. That is, χ o = 100% does not imply χ a = 100%. In general, the exact mapping between χ a , χ h , and χ o is not straightforward and depends on the relative abundance of the different (surrogate) species.
ZHENG ET AL. Note. Six aerosol species (see Table 4) are used in calculating the aerosol mixing state indices based on different definitions. We calculate χ a based on all six aerosol model species. The mixing state indices χ o and χ h are based on two surrogate species.

Feature Definitions
The features for the emulator were selected according to the following two criteria. First, the features should be available within both the PartMC-MOSAIC and CESM simulations, so that we have consistent features for training, validation, and application processes. In other words, quantities that are only available in Part-MC-MOSAIC but not in CESM simulations cannot serve as features for the emulators. Second, the features should be physically or chemically related to the aerosol properties. A subset of the aerosol and gas phase species and environmental variables were selected as features in this study (Table 4). The specific aerosol species were defined by the modal aerosol module (MAM4; Liu et al., 2012Liu et al., , 2016 and PartMC-MOSAIC.

Emulator Development
We Note. Type: A (aerosol), G (gas), E (environmental). The words "SRF," "a1," "a2," and "a4" refer to surface, Accumulation mode, Aitken mode, and Primary carbon mode.  (Friedman, 2001) framework, XGBoost is highly efficient, flexible, and portable. It can not only handle complex nonlinear interactions and collinearity among predictors, but also provides a parallel implementation that solves data science problems (particularly with structured data) efficiently. It has been successfully used in various fields such as physics (Mott et al., 2017), medical research (Abelson et al., 2018), and ornithology (Rosenberg et al., 2019). Here, we developed three emulators, corresponding to the three aerosol mixing state indices outlined in Section 3.
For each aerosol mixing state index, the emulator was trained from the entire 1800 training scenarios with hyperparameters (number of gradient boosted trees, maximum tree depth for base learners, and boosting learning rate) previously determined by using grid search with 10-fold cross-validation. Figure 3 shows the performance of emulators using the testing scenarios, which have never been seen by the training process. Table 5 shows the predictive accuracy of the XGBoost-based emulators compared to the benchmark using ordinary least squares. These metrics are important to keep in mind when interpreting the results in Section 5. The error for any predicted mixing state index can be roughly thought of as ±10%, and somewhat higher for χ h . Figure 4 shows the overall feature importance based on the "gain" metric, which quantifies the improvement in accuracy brought by a feature to the branches it is on. A higher value of this metric when compared to another feature implies it is more important for generating a prediction. In other words, the predictions are more sensitive to the features with higher gain values. This shows that the aerosol species are most important, although the ranking is different for the different mixing state indices. Environmental conditions (relative humidity and temperature) are of medium importance, and the gas phase concentrations are of lower importance.

Emulator Application
Applying the emulator is fast, compared to embedding PartMC-MOSAIC into regional models or ESMs for each grid cell. Note. The metrics are root-mean-square error (RMSE), mean absolute error (MAE), median absolute deviation (MAD), index of agreement (d; Willmott, 1981), Pearson correlation coefficient (PCC), and coefficient of determination (r 2 ). The emulators from linear regression with ordinary least squares (OLS) are examined as a benchmark. The p value (<0.001) applies for both methods using the F-test.

Table 5
Evaluation of Mixing State Emulators Using the Testing Data Set a testbed because CESM is a state-of-the-art and open-source ESM with many substantial science and infrastructure improvements (e.g., improved treatment of dust [Albani et al., 2014;Bogenschutz et al., 2018] and primary organic matter ), and improved historical simulations compared to available observations (Danabasoglu et al., 2020). The aerosol configuration in CESM/CAM has been calibrated by many observations Tilmes et al., 2015). We ran the model for the year 2011 with 6 years (2005-2010) spin-up. See Gantt et al. (2014) and He and Zhang (2014) for details of the spin-up procedure. The simulation was conducted at a resolution of 0.9° latitude × 1.25° longitude with the emission inventories from CMIP6 emissions . We stored the instantaneous outputs every 3 h in the simulation, which yielded 2920 timestamps for each surface-layer grid cell for the entire year of simulation time.
For reference, Figures S1 and S2 illustrate the seasonal variation and spatial distribution of the submicron aerosol species at the global scale. As we discussed in Text S1 (supporting information), the distributions were in line with previous studies (Bond et al., 2004(Bond et al., , 2013Hommel et al., 2011;Lack et al., 2004;Mahowald et al., 2005;Ramanathan & Carmichael, 2008;Sofiev et al., 2011;Tanaka & Chiba, 2006). It is worth mentioning that our simulation does not focus on optimizing the parameterizations of CESM. We used the default setting to conduct the simulations as a basis for estimating the different mixing state indices with the corresponding emulator in each grid cell. However, it should be noted that limitations and systematic biases in CESM (Hodzic et al., 2020;Schwantes et al., 2020) may be translated into the mixing state estimates.

Emulator Validation With Observational Data
To validate our χ emulators requires an observational data set that has both measured χ values as well as measured values for the emulator features (Table 4). Measuring χ requires quantitative per-particle mass fractions, which are difficult to estimate from observations (see Section 5.2). The data set that comes closest to having the required information was produced during the MEGAPOLI winter campaign (Healy et al., 2014) for a site in Paris, France from January 15 to February 11, 2010. We used hourly measurements from this observational data.
Healy et al. (2014) used single particle mass fraction estimates for BC, organic aerosol, ammonium, nitrate and sulfate to calculate single-particle diversities and the mixing state index χ obs . The mass fraction estimates were derived from a combination of single particle mass spectrometer, aerosol mass spectrometer, and multiangle absorption photometer measurements.
Using this data set for validating our emulator poses some challenges. The species list used to calculate χ obs was somewhat different from that used for χ a in this paper (Table 3). Further, not all emulator features (Table 4) were measured during the campaign. The missing features (dust, sea salt, DMS, H2O2, H2SO4, SO2, and SOAG) were set as "NaN" when applying the emulator, which is how XGBoost indicates missing data. XGBoost has explicit support for missing data and it learns default tree branch directions for missing values during training. To partition the measured organic aerosol and create an estimate of the features "pom" and "soa," we assumed a ratio of 1:1, based on CESM results for the same grid cell. We should not expect a perfect agreement between the emulator values, χ a , and observational values, χ obs , due to the difference in species lists used to compute χ and the fact that 7 out of 15 features were missing, including the most important feature "dst" (left panel of Figure 4). Nevertheless, over all MEGOPOLI campaign data the emulator had a mean value of χ a = 63%, which was close to the measured mean value of χ obs = 59%. In terms of the distribution, the 25th percentile, median, and 75th percentile values for the emulator were χ a = 56%, 63%, and 70%, which were also similar to the measured values of χ obs = 54%, 61%, and 65%, respectively. The minimum and maximum values during the campaign were χ a = 42% and 88%, compared to χ obs = 37% and 72%. Although the emulator produced a good match of the distribution compared to the observations, the predictions did not capture the diurnal cycle, which may be due to the differences in species lists and the missing features. Despite the difficulties in validating the emulator with observational data, this distributional agreement provides some support for the accuracy of the emulator.

Global Aerosol Mixing State Indices
In this section, we will first present the seasonal variation and spatial distribution of the three aerosol mixing state indices as obtained from applying the emulators. We will then put these findings in context with available observations, and finish this section with quantifying the relationship between the three indices as they were predicted for the entire model domain.

Seasonal Variation and Spatial Distribution of Aerosol Mixing State Indices
We averaged the 3-h instantaneous values of aerosol mixing state index for the two seasons (December-January-February, DJF, and June-July-August, JJA) to determine their seasonal variations and spatial distributions. We masked the areas where the mass fraction of any one species is higher than 99% (for χ o ) or 97.5% (for χ a and χ h ). The rationale for applying these masks is that only when at least two species are present in a given location does it make sense to quantify mixing state. Figure 5 shows the global maps of the aerosol mixing state index χ a based on chemical species abundance. In DJF, χ a ranged from 19% to 93% over the globe, with a mean of 64%. Of note is the region across the Atlantic from West and Central Africa to the northeast coast of South America, where χ a was smaller than 40%. In this region, dust and primary organic matter were prevalent, with contributions of sea salt ( Figure S2), and so this result can be interpreted to mean that these species are comparatively, although not completely, externally mixed.
In contrast, aerosol species were rather internally mixed (χ a > 70%) over North America, the North Atlantic Ocean, the Arctic Ocean, East Asia, and the North Pacific Ocean. The aerosol composition varied from region to region. Briefly, in DJF, aerosols over the continents of those regions were dominated by primary organic matter, SOA, and sulfate, while sea salt and sulfate governed the aerosol composition over the oceans.
The high values of χ a suggest that the species in these locations are rather internally mixed.
ZHENG ET AL.
10.1029/2020EA001500 11 of 19 In JJA, χ a varied within a range of 13%-94% and a mean of 69%, which is similar to DJF, but exhibited a different spatial pattern compared to DJF. The extent of regions with predominantly externally mixed aerosols (χ a < 40%) decreased, while those with internally mixed aerosols (χ a > 80%) increased. The lowest values of χ a shifted from the North Atlantic Ocean to Western and Northern Africa. Significant differences between seasons existed in the Gulf of Guinea and the Arabian Sea. Compared to DJF, the mixing state index increased greatly over the Gulf of Guinea (12°S-12°N, 12°W-12°E) for JJA, from χ a = 40% to χ a = 54%. The bulk composition of the aerosol also changed in this region between seasons. Combining this information, in DJF the submicron aerosol was mainly primary organic matter and dust, and it was comparatively externally mixed, while in JJA it was mainly primary organic matter, secondary organic matter, and sulfate, and was comparatively internally mixed. In contrast, χ a over the Arabian Sea (0°N-25°N, 50°E-75°E) experienced a drastic decrease of χ a from DJF to JJA, where a more internal mixture of sulfate and primary organic matter was replaced by a more external mixture dominated by sea salt, sulfate, and dust. Over the Southern Ocean (30°S-60°S), χ a values increased during JJA to values of nearly 80%. The bulk composition of the aerosol, however, was mainly sea salt (over 80% mass fraction), with only small contributions of the other aerosol species. Figures 6a and 6b show the global distribution of the aerosol mixing state index χ o (based on the mixing of optically absorbing and nonabsorbing species) for DJF and JJA, and, for reference, Figures 6c and 6d depict the corresponding bulk mass fractions of BC. We masked out the areas where the bulk mass fraction of BC is lower than 1% to only evaluate χ o where BC is present in nonnegligible amounts.
In DJF, χ o ranged from 40% to 93%. This range is smaller than the range for χ a although the global mean (67%) was close to that of χ a . The high values of χ o suggest that in most of the regions around the globe, BC tended to be internally mixed, and even more so in the continental regions, although a complete internal mixture was not reached anywhere. This means that BC is not evenly distributed over all particles, but it is also not exclusively confined to some particles. We can refer to the middle row in Figure 2 for an illustration of what the particle population might look like going from χ o = 63% to 91%. Note that BC-free particles may coexist with BC-containing particles under these conditions. ZHENG ET AL.
10.1029/2020EA001500 12 of 19  (d), respectively, where the areas with mass fraction of black carbon lower than 1% or higher than 99% are masked.
A high degree of internal mixture (χ o larger than 70%, meaning BC is rather spread out over the aerosol population) was predicted over the Arabian Sea, Bay of Bengal, India, and south-west of China. These are areas where the mass mixing ratios of BC were relatively high, and χ a was also relatively high. It suggests that BC was internally mixed with other species, with the most abundant species being primary organic matter and sulfate. However, high values of χ o coincide with low values of χ a over the Gulf of Guinea. This can be attributed to the external mixture of the nonblack-carbon species (dust and SOA) at this location.
In JJA, χ o had a range of 38%-94% and a global mean of 69%, overall similar to the DJF season. BC was more internally mixed in Southern Africa in JJA. In this area, it was mainly primary organic matter that formed an internal mixture with BC. Another pronounced change compared to DJF occurred in Eastern China, where the BC was internally mixed with an elevated level of primary organic matter, sulfate, and SOA. In general, areas with high mass fractions of BC were associated with higher values of χ o , meaning that in polluted urban regions, BC existed in an internal mixture with other aerosol species. In contrast, areas with low BC mass fractions, for example over the oceans, were associated with lower values of χ o (between an internal and external mixture). Figure 7 shows the distributions of the aerosol mixing state index χ h (based on the mixing of hygroscopic and nonhygroscopic species), and the corresponding mass fraction of nonhygroscopic species (BC, dust, and primary organic matter). Similar to Figure 5, we focus on areas where both surrogate species are present with a mass fraction threshold of 2.5%. The mixing state index χ h in DJF ranged from 20% to 82%, with a global mean of 54%. The Bay of Bengal, south-west of China, and Mongolia were areas where hygroscopic and nonhygroscopic species existed in a more internal mixture (χ h larger than 55%). Sulfate and primary organic matter were the most abundant species over these regions ( Figure S2). In contrast, hygroscopic and nonhygroscopic species were more externally mixed over the North Atlantic Ocean (near the equator), Southern Africa, Australia, Eastern Europe, and East China. Here, the dominant aerosol species were mineral dust and sea salt for the ocean, and sulfate, primary organic carbon, as well SOA for the land. In JJA, χ h varied between 26% and 87%, with a mean of 58%. Over Central Africa and Southern Africa, χ h increased to higher values (greater than 55%) than the DJF season. Here, primary organic carbon was the dominant ZHENG ET AL.  (d), respectively, where the areas with mass fraction of black carbon, dust, and primary organic matter lower than 2.5% or higher than 97.5% are masked.
nonhygroscopic species, while SOA and sulfate were both present as hygroscopic species. In contrast, χ h over the Arabian Sea decreased from DJF to JJA, similar to χ a . In this region, a predominantly sea salt/dust mixture that was more externally mixed replaced the more internally mixed primary organic carbon/sulfate mixture present in DJF.
The χ h distribution maps can be related to the framework of CCN error quantification by Ching et al. (2017). Their Figure 6 established a relationship between χ h and the error in CCN concentration prediction when a fully internal mixture was assumed (χ h = 100%). In our case, three-quarters of the grid cells have seasonally averaged χ h values between 20% and 63%, which corresponds to an error in CCN concentration between 50% and 100%. For 11% of the grid cells, χ h is over 75% and the error is expected to be less than 10%. We do not have any grid cells with seasonally averaged χ h values less than 20%, which would be associated with errors larger than 100%. However, instantaneous χ h values smaller than 20% do occur and therefore the seasonally averaged error based on instantaneous χ h values would be higher since the error in CCN concentration is a convex function of χ h .
The mixing state indices varied considerably over the oceans. The outflow of biomass burning or dust emission regions in Africa were areas of comparatively external mixtures in terms of species abundance (quantified by χ a ), owing to the coexistence of mineral dust or biomass burning aerosol and sea salt aerosol. Judging from the maps for χ o , plumes of BC originating from major biomass burning regions in Africa were initially comparatively internally mixed and then became more externally mixed as the plume was transported over the ocean. This makes sense considering that an increasingly external mixture should form as the carbonaceous plumes mix with sea salt aerosol as they are transported away from the coasts. Another finding is that BC is fairly (albeit not completely) internally mixed in highly polluted areas, which is consistent with prior work from modeling (Fierce et al., 2016) and measurements (Wang et al., 2010). These studies found that short aging time scales apply to BC aerosol in polluted regions. However, this may also be due to the crude resolution of horizontal grids in the global model. Very close to the source an external mixture may apply (Ching et al., 2019), but capturing this would require a high spatial resolution down to the kilometer scale.

Comparison to Observational Data
The comparison to observational data is still challenging at this point, as per-particle mass fractions are needed to determine the mixing state indices and they are difficult to estimate from observational data. While many data sets exist where single-particle data was collected in different environments (Ault et al., 2009;Gunsch et al., 2018;Murphy, Thomson, et al., 1998;Qin et al., 2012), the data has rarely been analyzed in a way that yields per-particle mass fractions based on chemical species that are compatible with the species used in chemical transport models. Notable exceptions are the studies by Healy et al. (2014), Ye et al. (2018), andChing et al. (2019) where single-particle mass spectrometers were used that produce estimates of per-particle mass fractions of sulfate, ammonium, nitrate, organics BC, which are species that chemical transport models track.
Here, we attempted to compare at least qualitatively to the work by Healy et al. (2014) A direct comparison of our results and observations like these is generally challenging, since point measurements in an urban environment are not directly comparable to a global modeling grid cell. However, it is encouraging that our results are relatively close to the observations. Clearly, it would be beneficial to have more observations available for model-observation comparisons, especially from regions that are not dominated by local sources. In addition, difficulties in a comparison may arise when the species list used in the observations is not easily mapped onto model species, as is for example the case with methods using electron microscopy or X-ray spectroscopy where χ was calculated but based on elemental composition (Bondy et al., 2018;Fraund et al., 2017).
Our χ o index can be related to observational studies that use SP2 measurements. Raatikainen et al. (2015) conducted measurements in the Finish Arctic during winter 2011-2012 and found that the number fraction of BC-containing particles was on average 24%, and that those BC-containing particles were thickly coated. While this study did not provide quantitative mixing state index calculations, it is an important finding that BC-containing particles (with various amounts of coatings) coexist with BC-free particles, and our result of χ o ≈ 72% for this area is consistent with this. As we see from Figure 2, a χ o value of 72% could correspond to a population where only some particles contain BC, but where BC-containing particles have a substantial nonBC fraction. Representing this mixing state is a challenge for many aerosol models used on the global scale that assume internal mixtures, but important for estimating aerosol optical properties accurately (Oshima et al., 2009). Figure 8 shows the correlation (with p value < 0.001) among the aerosol mixing state indices χ a , χ o , and χ h . The datapoints included in this figure are seasonally averaged values from each surface-layer CESM gridpoint for both DJF and JJA seasons. There is no strong relationship between χ a and χ o , which shows that each mixing state index captures different aspects of the population mixing state. As mentioned above, a population can appear as completely internally mixed when we only consider "black carbon" and "other species" as the two surrogate species. At the same time, the distribution of model species can vary within the "other species" class, resulting in a comparatively externally mixed population (<50%) in terms of χ a . In contrast, a high value for χ a can correspond to a lower χ o in cases where the average particle species diversity ZHENG ET AL. D α (in terms of optical properties) is low. This situation typically occurs when the particles are dominated by the "other species." A small fraction of "black carbon" exerts more influence on χ o than on χ a . Similar principles apply to the relationship between χ h and χ a , and χ h and χ o .

Relationship Amongst the Three Mixing State Indices
Overall, these results illustrate that, when determining mixing state from observed data using a certain set of species (determined by the instrumentation), it is difficult to infer the mixing state based on a different set of species. Mixing state indices for different purposes capture different information about the aerosol. Therefore, care has to be taken regarding the comparison of mixing states in different environments when the measurement techniques to estimate particle mass fractions are different. Figure 9 summarizes the seasonally averaged aerosol mixing indices for grid points with high anthropogenic pollution in the Northern hemisphere. We defined these grid points as those where the mass mixing ration of the sum of BC, POM, SOA, and SO4 is higher than the 90% percentile for both DJF and JJA. Regardless of which mixing state index we consider, extremely external mixtures (χ lower than 20%) were not predicted for the seasonal averages. However, individual 3-h values reached as low as 3%, 11%, and 9% for χ a , χ o , and χ h , respectively. Areas of high anthropogenic pollution in the northern hemisphere, such as central Europe and south-east Asia had larger χ h -values during JJA which is consistent with faster chemical aging during northern hemisphere summer.

Conclusions
In this paper, we presented a framework for estimating submicron aerosol mixing state indices at a global scale. A machine learning model (an emulator) was trained on high-detail simulations, using the particle-resolved PartMC-MOSAIC, and then applied to the output from a coarser model (MAM4 within CESM) to enhance its information content.
We developed three emulators based on the XGBoost algorithm to determine the aerosol mixing state indices for submicron aerosol in terms of chemical species abundance (χ a ), the mixing of optically absorbing and nonabsorbing species (χ o ), and the mixing of hygroscopic and nonhygroscopic species (χ h ) using variables available in CESM. Based on error metrics for the testing data set, the predictions of mixing state metrics have an error of ∼±10% associated with them, and somewhat higher for χ h . The emulators were applied to CESM simulations to produce global maps of aerosol mixing state indices for every 3 h timestamp. For this work, we focused on seasonal variation and spatial distribution of these aerosol mixing state indices.
The seasonally averaged indices varied spatially between 13% and 94% for χ a , between 38% and 94% for χ o , and 20% and 87% for χ h , with global annual averages of 67%, 68%, and 56%, respectively. Hygroscopic and nonhygroscopic species were more externally mixed over the North Atlantic Ocean near the equator, off the coasts of Southern Africa, and Australia. An internal mixing state assumption for those areas could result in errors in CCN concentration of 50%-100%.
High values in one mixing state index can be consistent with low values in another index, depending on the grouping of species and their relative abundance. This indicates that the different indices capture different aspects about the mixing state of an aerosol. When comparing mixing state indices for different environments from observations, care needs to be taken that they use the same choice of species as the basis.
Comparing our results with available ambient observations of mixing state indices in the literature is still challenging, as estimates of per-particle mass fractions needed to determine the mixing state indices are difficult to obtain and long-term records are not yet available. Qualitative comparisons with studies by Healy et al. (2014), Ye et al. (2018), and Ching et al. (2019) show that our results are similar in magnitude to the observations.