\section{Introduction and Motivation}
The concept of aerosol mixing state describes how different aerosol
chemical species are distributed among the particles in a population
\citep{RIEMERETAL19}. A completely ``externally mixed'' population
contains only one species per particle, while a completely
``internally mixed'' population contains particles with the same
composition as the bulk population. Many different intermediate mixing
states are possible between those two extremes, and these commonly
occur in the ambient atmosphere. Aerosol mixing state influences the
particles' properties such as hygroscopicity
\citep{FIERCEETAL17,HOLMGRENETAL14}, optical properties
\citep{FIERCEETAL16,JACOBSON01}, cloud condensation nuclei activity
\citep{CHINGETAL12,WANGETAL10}, ice nucleation potential
\citep{KNOPFETAL18}, the aerosols' lifetime \citep{KOCHETAL09} in the
atmosphere, and their deposition in the human respiratory system
\citep{CHINGKAJINO18}. However, Earth System models (ESMs) or regional
climate models usually hold simplified mixing state assumptions. These
can influence how accurately physicochemical properties of aerosols
are predicted \citep{CHINGKAJINO18, FIERCEETAL17}, and thereby limit
the accuracy of estimating aerosol impacts on climate and health.

\citet{RIEMERWEST13} proposed an entropy-based diversity metric to
quantify aerosol mixing state. This metric, the mixing state index
($\chi$), has a range from 0\% (a fully externally mixed population)
to 100\% (a fully internally mixed population). The index $\chi$ 
has been applied to both modeling \citep{ZHUETAL16} and experimental
\citep{HEALYETAL14,OBRIENETAL15,FRAUNDETAL17,BONDYETAL18,YEETAL18}
work. For example, $\chi$ has been used for
quantifying errors in CCN concentration \citep{CHINGETAL17} prediction,
as well as the prediction of soot particles depositing in the human
respiratory system \citep{CHINGKAJINO18}.

Knowing the global distribution of aerosol mixing state index is
desirable, as it can be used for error quantification at the large
scale. But deriving this distribution directly by using a benchmarking
particle-resolved model \citep{RIEMERETAL09} at the global scale would
be computationally prohibitive over the time scales of interest. To
overcome this limitation, supervised-learning emulators were developed
for predicting the spatial distribution aerosol mixing state
across the globe \citep{HUGHESETAL18,ZHENGETAL20}. However, the
spatial delimitation based on aerosol mixing state indices remains
unclear, which hinders further spatiotemporal analysis.

The research questions for this paper are: (1) Is it possible to
regionalize the global mixing state indices, i.e., define aerosol
mixing state zones, similar to distinct climate zones
\citep{FOVELLFOVELL93,MITCHELL76}? (2) What are spatiotemporal
patterns of the regionalized mixing state indices? This paper
describes an effort to regionalize the global mixing state indices
using a simple but effective unsupervised learning approach.

\section{Methods}
\subsection{Mixing state indices definition}

The mixing state index $\chi$ \citep{RIEMERWEST13} is given by the affine
ratio of the average particle species diversity ($D_{\alpha}$) and bulk
population species diversity ($D_{\gamma}$) as
${\chi} = \frac{D_{\alpha}-1}{D_{\gamma}-1}$. 
The diversities $D_{\alpha}$ and $D_{\gamma}$ are calculated based on
per-particle species mass fractions and mass fraction of particles as
described in detail in \citet{RIEMERWEST13}.

This study focused on submicron aerosols because they are the dominant
category of particles for radiation interactions and the provision of
cloud condensation nuclei. We defined the mixing state indices in
three different ways,
namely based on chemical species abundance ($\chi_\textrm{a}$), based
on the mixing of absorbing and non-absorbing species
($\chi_\textrm{o}$), and based on the mixing of hygroscopic and
non-hygroscopic species ($\chi_\textrm{h}$). The index
$\chi_\textrm{a}$ was defined based on all the model aerosol
species. A lower value for $\chi_\textrm{a}$ indicates that the
species tend to be present in separate particles. For
$\chi_\textrm{o}$, we considered two surrogate species, black carbon
(absorbing) and all other aerosol species grouped together (assumed to
be non-absorbing). Thus, a lower value for $\chi_\textrm{o}$ indicates
the absorbing species black carbon and the sum of all other
(non-absorbing) species are more externally mixed. Similarly,
$\chi_\textrm{h}$ was also calculated from two surrogate species. We
combined black carbon, dust, and primary organic matter as one
surrogate species, given their comparatively lower hygroscopicities,
and salt, secondary organics, and sulfate as the other surrogate
species. Correspondingly, a lower value for $\chi_\textrm{h}$
represents the case where hygroscopic and non-hygroscopic species
tended to be present in separate particles. Note that mixing state
indices defined for different purposes capture different aspects of
the overall mixing state, and they are uncorrelated
\citep{ZHENGETAL20}.

\subsection{Data}
We used the supervised-learning surrogate models developed by
\citet{ZHENGETAL20} to predict the mixing state indices from ESM
(Earth System Model) output. The surrogate models were trained on an
ensemble of particle-resolved model simulations generated by the
particle-resolved model PartMC-MOSAIC
\citep{RIEMERETAL09,ZAVERIETAL08}, with the training labels for $\chi$
calculated directly from the particle-resolved data. The strategy to
generate the training and testing data was to vary 45 input parameters
for the PartMC-MOSAIC model scenarios, including primary emissions of
different aerosol types, primary emissions of gas phase species, and
meteorological parameters. A Latin Hypercube sampling approach was
employed to sample the parameter space efficiently. The surrogate
models were trained by using XGBoost (eXtreme Gradient Boosting,
\citet{CHENGUESTRIN16}), and can be expressed as:
\begin{equation}
\chi_S = f_S(A,G,E).
\end{equation}
Where $\chi_S$ is the mixing state index ($\chi_\textrm{a}$,
$\chi_\textrm{o}$, or $\chi_\textrm{h}$) and $f_S$ denotes
the surrogate function for the corresponding mixing state index.
The variants $A$ (aerosol), $G$ (gas), and $E$
(environment) represent the features (variables) from the PartMC
simulations that are used for predicting the mixing state index. These variables are
also available from the ESM.  A detailed definition of the features for the
surrogate models is given in \citet{ZHENGETAL20}.

We used the simulations from the Community Earth System Model Version
2 \citep[CESM2 version 2.1.0;][]{DANABASOGLUETAL20} with MAM4
\citep{LIUETAL16} to provide variables at the global
scale. Specifically, we used the component set ``FHIST'' for the
global simulation configuration. This component set represents a
typical historical simulation using an active atmosphere and land with
prescribed sea surface temperatures and sea-ice extent, and a 1~degree
finite volume dycore with the forcing data available from 1979 to
2015. We ran the model for the year 2011 with 6~years (2005-2010)
spinup. The simulation was conducted at a resolution of 0.9$^\circ$
latitude by 1.25$^\circ$ longitude along with emission inventories
from CMIP6 emissions \citep{EMMONSETAL20}. We stored the instantaneous
outputs every three hours in the simulation, which yielded
2920~timestamps for each surface-layer grid cell for the entire year
of simulation time.

With the surrogate models and the ESM simulation, we predicted the
mixing state indices for each grid cell at each timestamp (i.e. every
3~hours). For each grid cell, the mixing state indices were averaged by
month.  Therefore, each grid cell has a vector of stacked mixing state
indices with a length of 36 (12 months $\times$ 3 mixing state types),
represented as $\boldsymbol{x} =
(\chi_{\textrm{a},1},\dots,\chi_{\textrm{a},12},\chi_{\textrm{o},1},\dots,\chi_{\textrm{o},12},\chi_{\textrm{h},1},\dots,\chi_{\textrm{h},12})$.

\subsection{Regionalization strategy}
Here we interchangeably use the terms ``region'' and ``cluster'' to
refer to a group of grid cells. We used the $k$-means unsupervised
clustering algorithm to partition $n$ grid cells into $k$ clusters,
minimizing the within-cluster variances. The $k$-means method has
been applied to environmental sciences for ecoregion delineation
\citep{KUMARETAL11}, environmental risk zoning of the chemical
industrial area \citep{SHIZENG14}, and clustering haze trajectory of
peatland fires \citep{KHAIRATETAL17}, among other applications.

In our case, given $n=55\,296$ (192 latitude $\times$ 288 longitude) grid cells
$(\boldsymbol{x}_1, \boldsymbol{x}_2, \dots, \boldsymbol{x}_{55296})$, where
each grid cell contains a 36-dimensional vector (three mixing state indices over
12 months), $k$-means clustering aims to partition these grid cells into $k$
($k \le n$) sets $\boldsymbol{S} = \{S_1,S_2,\dots,S_k\}$ so as to
minimize the within-cluster variance. The objective is to find:
\begin{equation}
\mathop{\arg\min}_{\boldsymbol{S}} \ \ \ \sum_{i=1}^{k}\sum_{\boldsymbol{x}\in S_i}=\Vert{\boldsymbol{x}}-\boldsymbol{\mu}_i\Vert^2
\end{equation}
where $\boldsymbol{\mu}_i$ is the mean of the points in $S_i$.

The $k$-means clustering method requires the a priori specification of
the number of clusters ($k$). As the value of $k$ increases, the
variance $V_k$ asymptotically approaches 0 until $k$ equals $n$. Here
we used the relative difference in variance between $k-1$ and $k$ to
determine the number of clusters. Starting from $k=2$, we computed the
criteria of relative difference $V_{k-1} - V_{k} \le 0.05V_{k-1}$ to
determine the cluster number $k$. In this study, since
$V_{10} - V_{11} \le 0.05V_{10}$, we selected $k=11$
clusters.

\section{Results}

The regional clusters of mixing state are shown in
Figure~\ref{fig:map_clusters}. Unlike with supervised learning, the
performance of unsupervised learning are challenging to evaluate since
there is no ground truth.  However, the clustering process of this
study was merely based the information of mixing state indices,
without the direct guidance such as aerosol species distribution or
other explicit spatial information. Therefore a success indicator of
the spatiotemporal clustering is whether spatially-contiguous regions
emerge from this procedure. The clusters in
Figure~\ref{fig:map_clusters} exhibit notable regional patterns among
the eleven clusters (mixing state zones), suggesting that this
approach has indeed identified meaningful clusters. Each cluster has a
unique spatiotemporal distribution of three mixing state indices, and
captures the spatiotemporal variations that could not be detected by
the global overall averages. For example, ocean regions (e.g.,
clusters 8, 9, 10 in the Southern Hemisphere) tend to display a more
pronounced seasonal cycle than land regions (e.g., cluster 7) and
overall averages.

Clusters~1 (the Arctic) and~11 (the Antarctic) were characterized by
high mixing state indices (with monthly means of at least 67\%)
throughout the year. The high values of all three mixing state indices
means that the different aerosol species are rather internally mixed
in these regions. Clusters~5 (tropical oceans), 8 (the South China
Sea, Indian Ocean and the Pacific Ocean in the northern hemisphere
near the equator, and the northern part of the Southern Ocean) and~9
(the southern part of the Southern Ocean) were ocean areas. Similar
bulk aerosol composition was predicted in these regions, however the
annual cycle of mixing state varies more for clusters~8 and~9 than
for~5. Clusters~2, 3, and~4 covered oceans at mid and high latitudes
in the Northern Hemisphere, and lands including Northern and Southern
America, Europe, Middle East, Southeast Asia, as well as
Australia. Clusters~6 (the region
across the Atlantic from Western and Northern Africa to the northeast
coast of South Africa, as well as part of the Indian Ocean) and~7
(Middle, Eastern, and Southern Africa, as well as Asia) were
characterized by their higher values in $\chi_\textrm{o}$ compared to
$\chi_\textrm{a}$ and $\chi_\textrm{h}$, meaning that black carbon and
non-absorbing species are comparatively internally mixed.

\begin{figure}
  \centering
  \includegraphics[width=\textwidth,trim=50 0 0 0]{./mixing_state_zones.pdf}
  \caption{Aerosol mixing state zones based on unsupervised
    regionalization. Also shown are monthly average mixing state
    indices and annual averages of aerosol species mass fractions
    (bulk composition).}\label{fig:map_clusters}
\end{figure}

\section{Conclusions}
In this paper, we developed a simple but effective unsupervised
learning approach to regionalize global aerosol mixing state
indices. We used the monthly averages of the spatially-varying aerosol
mixing state indices as input data. Each grid cell was then assigned
to a region using $k$-means clustering without explicit spatial
information as input. This approach resulted in eleven distinct
regions over the world with specific spatial aggregation
patterns. Each region exhibited a unique distribution of mixing state
indices, suggesting that the unsupervised regionalization approach had
identified meaningful regions. To the best of our knowledge, this is
the first study to define aerosol mixing state zones, which we
anticipate will be helpful for future studies of the global aerosol
burden.

\begin{ack}
We would like to acknowledge high-performance computing support from
Cheyenne (\url{doi:10.5065/D6RX99HX}) provided by NCAR's Computational and
Information Systems Laboratory, sponsored by the National Science
Foundation. The CESM project is supported primarily by the National
Science Foundation.
This research used resources of the Oak Ridge Leadership Computing
Facility, which is a DOE Office of Science User Facility supported
under Contract DE-AC05-00OR22725. This research was supported in part
by an appointment to the Oak Ridge National Laboratory ASTRO Program,
sponsored by the U.S. Department of Energy and administered by the Oak
Ridge Institute for Science and Education. We also acknowledge funding
from DOE grant DE-SC0019192 and NSF grant AGS-1254428.
This research is part of the Blue Waters sustained-petascale computing
project, which is supported by the National Science Foundation (awards
OCI-0725070 and ACI-1238993) the State of Illinois, and as of
December, 2019, the National Geospatial-Intelligence Agency. Blue
Waters is a joint effort of the University of Illinois at
Urbana-Champaign and its National Center for Supercomputing
Applications. Louisa Emmons is thanked for thoughful comments on the
CESM2 simulations and the manuscript. We thank AWS for providing AWS
Cloud Credits for Research.
\end{ack}

\clearpage
