A single-column particle-resolved model for simulating the vertical distribution of aerosol mixing state: WRF-PartMC-MOSAIC-SCM v1.0

The PartMC-MOSAIC particle-resolved aerosol model was previously developed to predict the aerosol mixing state as it evolves in the atmosphere. However, the modeling framework was limited to a zero-dimensional box model approach without resolving spatial gradients in aerosol concentrations. This paper presents the development of stochastic particle methods to simulate turbulent diffusion and dry deposition of aerosol particles in a vertical column within the planetary boundary layer. The new model, WRF-PartMCMOSAIC-SCM, resolves the vertical distribution of aerosol mixing state. We verified the new algorithms with analytical solutions for idealized test cases and illustrate the capabilities with results from a 2-day urban scenario that shows the evolution of black carbon mixing state in a vertical column.


Introduction
Aerosol particles impact the Earth's radiative budget directly by scattering and absorbing shortwave radiation (Mc-Cormick and Ludwig, 1967;Charlson and Pilat, 1969;Charlson et al., 1992), and indirectly by modifying cloud microphysical properties (Twomey, 1977;Albrecht, 1989;Rosenfeld, 2000). The magnitude of these impacts on climate depends not only on the bulk amount of aerosol material in the atmospheric column, but also on its vertical distribution within the column (Haywood and Shine, 1997;Schulz et al., 2006;Zarzycki and Bond, 2010;Samset and Myhre, 2011;Ban-Weiss et al., 2012), and on the microphysical characteristics of the aerosol population, such as the size distribution of the particles and the aerosol composition on a per-particle level (McFiggans et al., 2006;Moffet and Prather, 2009;Zelenyuk and Imre, 2009;Zelenyuk et al., 2010). For the purposes of this paper we use the term "aerosol mixing state" to refer to the distribution of chemical species across the aerosol population (Riemer and West, 2013;Winkler, 1973). This is distinct from the use of the term "mixing state" for the arrangement of components within a particle (e.g., homogeneous mixture or core-shell arrangements).
Observational evidence shows that aerosol mixing state varies with altitude. For example, Pratt and Prather (2010) used the aircraft aerosol time-of-flight mass spectrometer to measure vertical profiles of single-particle composition over Wyoming and northern Colorado and found that carbonaceous particles were mixed with ammonium, nitrate, and sulfate at low altitudes, while they were mixed with sulfate and sulfuric acid at higher altitudes. Measurements with the single-particle soot photometer (SP2) showed that the fraction of black carbon particles that are heavily coated increases with altitude (McMeeking et al., 2011).
These experimental findings confirm that the composition of individual aerosol particles constantly changes during the particles' lifetime as a result of aging processes such as coagulation (Fassi-Fihri et al., 1997), condensation (Pósfai et al., 1999), and photochemical processes (Kotzick and Nießner, 1999), and that this is intimately linked to the transport processes in the atmosphere. To better represent these processes in chemical transport models, several twodimensional sectional models have been developed, such as MADRID-BC (Oshima et al., 2009a, b), the MS-Resolved WRF-Chem (Matsui et al., 2013), and WRF-Chem/ATRAS-MOSAIC (Matsui et al., 2014). A common feature of these models is that they use a two-dimensional sectional framework to represent black-carbon-containing particles, with one dimension being dry diameter and the other dimension being black carbon mass fraction. Building on previous two-dimensional sectional frameworks, the MOSAIC-MIX model (Ching et al., 2016) adds an additional dimension to represent hygroscopicity and shows that this optimizes the calculations of CCN concentrations and aerosol optical properties. The SCRAM model (Zhu et al., 2015(Zhu et al., , 2016a, also a two-dimensional sectional model, uses an alternative discretization based on both size and composition where composition is tracked by mass fractions of different chemical groups such as inorganic hydrophilic, organic hydrophilic, organic hydrophobic, black carbon, and dust. From the application of the different types of aerosol models described above within spatially resolved threedimensional chemical transport models, we learn that it is important to track the aerosol mixing state in order to accurately predict particle aging, the associated aerosol optical properties, and the resulting heating rates (Riemer et al., 2003;Matsui et al., 2013;Zhang et al., 2014;Matsui, 2016;Zhu et al., 2016a, b). While some chemical transport models focus on representing black carbon mixing state (Matsui et al., 2013;Matsui, 2016), other models have allowed for more general mixing state representations (Zhang et al., 2014;Zhu et al., 2016a). However, further extension of aerosol bin schemes to include additional dimensions to capture greater mixing state detail eventually becomes computationally prohibitive.
In contrast to the distribution-based models mentioned here, particle-resolved aerosol models simulate a representative group of particles distributed in composition space, treating coagulation, condensation/evaporation, and other important processes on an individual particle level. Relative particle positions within this computational volume are not tracked, but processes such as coagulation are instead simulated stochastically, following the approach pioneered by Gillespie (1975). Particle methods are attractive, because they resolve the full aerosol mixing state without any ad hoc assumptions. The storage cost of these models is proportional to the number of particles, the computational cost for evaporation/condensation is proportional to the number of particles, and the computational cost for coagulation is proportional to the number of coagulation events (Riemer et al., 2009).
For the large number of computational particles needed for atmospheric simulations, we developed efficient algorithms for coagulation (Riemer et al., 2009;Michelotti et al., 2013) and for appropriately weighting computational particles (DeVille et al., 2011). These were implemented in the Particle Monte Carlo (PartMC) model for simulating atmo-spheric aerosol dynamics and coupled with the state-of-theart MOSAIC aerosol chemistry model (Zaveri et al., 2008), which simulates the gas-and particle-phase chemistries, particle-phase thermodynamics, and dynamic gas-particle mass transfer in a deterministic manner. The coupled model system, PartMC-MOSAIC, predicts number, mass, and full composition distributions, and is therefore suited for applications where any or all of these quantities are required. The particle-resolved approach eliminates any errors associated with artificial numerical diffusion in composition space. As a result, its treatment of aerosol mixing state dynamics and chemistry makes PartMC-MOSAIC suitable for use as a numerical benchmark of mixing state for more approximate models (McGraw et al., 2008;Kaiser et al., 2014).
In previous work PartMC-MOSAIC has been used as a box model (Zaveri et al., 2010;Tian et al., 2014;Fierce et al., 2016;Ching et al., 2016), and hence it was not possible to resolve spatial gradients in aerosol mixing state. To overcome this limitation, we have now coupled PartMC-MOSAIC with the Weather Research and Forecast (WRF) model to allow transport of aerosol particle populations and gas species concentrations. In this paper we present the model development that couples PartMC-MOSAIC with the WRF Single Column Model, resulting in a fully coupled one-dimensional atmospheric-dynamics/aerosol-particle model that not only resolves the particle mixing state on a per-particle level, but also resolves the vertical structure of the atmosphere. This paper is structured as follows. In Sect. 2 we write the governing equations for the coupled gas-aerosol onedimensional column model. In Sect. 3 we discuss the specifics of the coupled model including numerical approximations to model the vertical transport of aerosols. In Sect. 4 we present test-case verification of the two new model processes of turbulent transport and particle removal by dry deposition. Section 5 shows the multidimensional particleresolved results for an idealized scenario, focusing on the evolution of black-carbon-containing particles.

Coupled aerosol-gas governing equations
In this section we describe the model equations that govern the evolution of aerosol particles and trace gases in a vertical column. We include gas-phase chemistry, gas-toparticle conversion, coagulation of aerosol particles, emission of aerosol and gases, and the transport of aerosol particles and trace gases in the vertical column. We ignore horizontal diffusion and advection of trace gases and aerosol particles into and out of the column by assuming horizontal homogeneity.
An aerosol particle contains mass µ a ≥ 0 of species a, for a = 1, . . ., A, so that the particle composition is described by the A-dimensional vector µ ∈ R A . The cumulative aerosol number distribution at height z with constituent masses µ at time t is N (z, µ, t) (m −3 ). The aerosol number distribution at height z and time t with constituent masses µ is then defined by The concentration of gas-phase species i at height z and time t is given by g i (z, t), for i = 1, . . ., G, so that gas-phase concentrations form the G-dimensional vector g(z, t) ∈ R G . We assume that the first C aerosol and gas species undergo gas-to-particle conversion and are indexed in the same order so that gas species i partitions with aerosol species i for i = 1, . . ., C. Additionally, species C + 1 is assumed to be water.
The evolution of the multidimensional aerosol number distribution is given by where w(z, t) (m s −1 ) is the vertical velocity, K h (z, t) (m 2 s −1 ) is the diffusion coefficient of heat, K(µ, µ ) (m 3 s −1 ) is the coagulation rate between particles µ and µ , n emit (z, µ, t) (m −3 kg −A s −1 ) is the number distribution rate of aerosol emissions which can be specified at any height, c a (kg mol −1 ) is the conversion factor from moles of gas species a to aerosol species a, I a (µ, g, t) (mol s −1 ) is the condensation flux of gas species a, c w (kg mol −1 ) is the conversion factor for water, and I w (µ, g, t) (mol s −1 ) is the condensation flux for water. The turbulent transport term is written using the gradient of mixing ratio rather than the gradient of concentration to account for the vertical variations in density that are present in the atmosphere (Eq. 6, Venkatram, 1993). Equation (2) does not contain a term for gravitational sedimentation since we focus our test case on submicron particles for which the settling velocities are very small. As an example, over the course of the 48 h simulation period a 1 µm particle would only settle by about 10 m, which is less than the smallest vertical grid size used here. Gravitational settling should be included in scenarios that involve larger particles such as sea salt and dust or simulations with finer vertical resolution. The evolution of trace gas concentrations is given by whereġ emit,i (z, t) (mol m −3 s −1 ) is the emission rate of gas species i and R i (g(z, t)) (mol m −3 s −1 ) is the concentration growth rate of gas species i due to gas-phase chemical reactions. For Eqs. (2) and (3) we use reflective boundary conditions at the top of the domain and partly reflecting, partly absorbing boundary conditions at the surface. For the aerosol distribution in Eq. (2), this is given at the top of the domain (z = h) by and at the surface by Here V d (µ) (m s −1 ) is the dry deposition velocity, which depends on particle size and composition. Dry deposition velocities for aerosols are computed using the size-dependent dry deposition scheme described in Zhang et al. (2001) by their Eqs. (1), (2), (3), (5), (6), (7c), (8), and (9). Instead of using Eq. (4) in Zhang et al. (2001), which describes the calculation of the aerodynamic resistance, the aerodynamic resistance computed by the WRF model is used, as described by McRae et al. (1982). Further, we do not need to use the parameterization for the correction of particle size for high relative humidity conditions given in Eq. (10) in Zhang et al. (2001), since we explicitly compute the water content of the aerosol particles and hence directly calculate the particles' wet diameters.
For the gas-phase concentrations of Eq.
(3) the boundary condition at the top of the domain is given by and at the surface by where V d,i (m s −1 ) is the dry deposition velocity of gas species i. The dry deposition velocity of each gas species is determined by WRF/Chem as described in Grell et al. (2005) with the use of the surface resistance parameterization from Wesely (1989).

Model discretization
We coupled the different model components (WRF, PartMC, and MOSAIC) by using the operator splitting (Press et al., 2007, Sect. 20.3.3) which allows for the use of independent numerical methods to solve each portion. The Advanced Research Weather Research and Forecasting (WRF-ARW) Model ( WRF t ) is used to solve for the meteorological variables, details of which are further described in Skamarock et al. (2008). WRF computes temperature, pressure, eddy diffusivity, aerodynamic resistance, and dry deposition velocity for gases, which are then used in Eqs.
(2) and (3). The WRF model is discretized using a terrain-following hydrostatic-pressure coordinate system η which is constant in time. The aerosols and gas species are evolved in the transport step ( Trans t ) on a geometric height coordinate system z that is computed from the geopotential field and changes over time due to column pressure changes.
The vertical advection terms found in Eqs. (2) and (3) are entirely due to pressure-induced grid changes and do not cause particles to be transported across grid cell edges. At every time step the z grid is moved by WRF, and this accounts for vertical advection.
The Particle-resolved Monte Carlo (PartMC) model ( PartMC t ) is used to treat the coagulation term and the emission term in Eq. (2). Emissions are simulated by stochastically sampling a finite number of particles at each time step, approximating the continuum emission distribution. Coagulation is efficiently simulated using a fixed time step method and a binned acceptance procedure. These approaches are further described in Riemer et al. (2009), DeVille et al. (2011), and Michelotti et al. (2013. The Model for Simulating Aerosol Interactions and Chemistry (MOSAIC) ( MOSAIC t ) is used to solve gas-phase and gas-to-particle chemistry terms in Eqs. (2) and (3). The composition of each particle can change due to evaporation and condensation of chemical species to and from the gas phase. MOSAIC (Zaveri et al., 2008) deterministically treats gasphase chemistry and gas-particle partitioning. MOSAIC consists of gas-phase mechanism Carbon-Bond Mechanism version Z (CBM-Z) (Zaveri and Peters, 1999), the Multicomponent Equilibrium Solver for Aerosols (MESA) for aerosol solid-liquid partitioning (Zaveri et al., 2005a), the multicomponent Taylor expansion method (MTEM) for estimating activity coefficients of electrolytes and ions in aqueous solutions (Zaveri et al., 2005b), and the adaptive step timesplit Euler method (ASTEM) for gas-particle partitioning (Zaveri et al., 2008). MOSAIC uses the Secondary Organic Aerosol Model (SORGAM) scheme (Schell et al., 2001) for the treatment of secondary organic aerosol. The transport model ( Trans t ) solves the turbulent transport terms for gas species and particles deterministically for gases and stochastically for particles as described in Sect. 3.1. The finitevolume method used for particle transport is also applied to gas transport to provide consistency between the transport of gases and aerosols.

Stochastic aerosol transport algorithm for turbulent diffusion
This section details the treatment of the turbulent transport term of Eq.
(2). We use a stochastic sampling approach for moving particles between grid boxes, rather than tracking the exact location of each simulated particle. This is done for computational efficiency. By using a stochastic sampling approach, only a fraction of particles will have their grid cell positions updated in each time step. In many cases, this transported fraction is rather small. Within each grid cell k of the model domain we represent the aerosol state k with N k p particles, where k = {µ 1 , µ 2 , . . ., µ N k p }, and the particle order is not significant. There may be multiple identical particles in a single grid cell, so k is a multiset in the sense of Knuth (1998, p. 473). Figure 1. Schematic of the single-column domain centered on grid cell k with neighboring grid cells k − 1 and k + 1.

Each particle is an
is the mass of species a in particle i, for a = 1, . . ., A and i = 1, . . ., N k p . The aerosol population k can be thought of as populating a computational volume V k which is smaller than the physical volume of the grid cell. The computational volume is representative of the mean properties of that grid cell. The value of the computational volume is the ratio of the number of computational particles contained in the grid cell over the number concentration of the grid cell. The assumption that the modeled aerosol represents the aerosol throughout the grid cell is the same as the one used for the box model version of PartMC-MOSAIC presented in Riemer et al. (2009), which simulated particles within a given computational volume that was representative of the well-mixed boundary layer during the day, and of the residual layer during the night. For simplicity we use here a flat weighting for all computational particles in the sense of DeVille et al. (2011).
Since PartMC-MOSAIC resolves a finite population of particles k in a given volume V k within grid cell k, we must determine the finite number of particles that are transported in and out of grid cell k due to turbulent transport. To determine the set of particles to remove from particle set k as well as to add to k from neighboring grid cells, we discretize the vertical turbulent transport term of Eq. (2) with respect to space, time, and particle number. This is in contrast to conventional models, which transport scalar variables such as mass mixing ratios. Therefore, the discretization process requires an additional step to transport particles.
We first present the discretization in space and time in terms of deterministic number concentrations and particle numbers (Sect. 3.1.1). The resulting equation for turbulent transport of deterministic particle numbers is then discretized for stochastic particle number as shown in Sect. 3.1.2, and time step selection is presented in Sect. 3.1.3.
3.1.1 Discretization in space and time in terms of deterministic particle number Figure 1 introduces our notation for the spatial discretization of the single-column model domain. The vertical grid spacing z is non-uniform and also varies over time. The variation over time is a result of WRF using the pressure-based vertical coordinate η, while the physical height z of grid cells is computed from geopotential height.
To account for the variation in z both with respect to height and to time, we define the distance from the top to the bottom edge of grid cell k at time t s to be and the distance between the center of grid cells k and k + 1 at t s by z s To obtain transported number concentrations and eventually a discrete number of transported particles, we define the total aerosol number concentration N (z, t) at height z and time t as where n(z, µ, t) is the multidimensional aerosol number distribution. We then apply the turbulent transport process in terms of total number concentration as where K h (z, t) is eddy diffusivity of heat, ρ(z, t) is density of dry air, and q(z, t) is the mixing ratio defined by An explicit, second-order accurate discretization scheme was selected for this work. While an explicit method simplifies the parallel implementation because it only requires communication between neighboring grid cells, there exists no theoretical reason why other numerical schemes may not be used, including higher-order, semi-implicit, or implicit methods. Following a finite-volume discretization, we arrange the derivation (see Appendix B for details) to isolate gain and loss terms of the number concentration N s k , giving The four transport terms are illustrated in Fig. 2. The arrows in Fig. 2 show the transported number concentrations from and to grid cell k. The first subscript indicates the origin grid cell, and the second subscript indicates the destination grid cell. The superscripts indicate either gain (G) for the destination grid cell or loss (L) for the origin grid cell at time step s. For example, N G,s k+1,k is the number concentration transported from grid cell k + 1 to grid cell k at time s, representing a gain for grid cell k. N L,s k,k+1 is the number concentration transported from grid cell k to k + 1 at time s, representing a loss for grid cell k. These transport terms can also be expressed in terms of the product of a coefficient β and the number concentration in the origin grid cell.
Eventually, we want to perform turbulent transport of discrete particles; therefore, as the next step we express Eq. (14) in terms of real-valued particle number instead of number concentration. The deterministic real-valued particle number of grid cell k at time s is related to number concentration by where V s k is the computational volume of grid cell k at time t s . Applying Eq. (15) to Eq. (14) and multiplying by V s k yields the following equation for the real-valued particle number where naming conventions are used similarly to Eq. (14), with particle number C replacing number concentration N and coefficient p replacing β.

Stochastic turbulent transport of particles
Equation (16) expresses the gains and losses of grid cell k in terms of deterministic real-valued particle number C k . However, the PartMC model simulates a finite set of particles for each grid cell rather than a number concentration or deterministic real-valued particle number. Therefore, an additional step is required to transform equations from the de-terministic real-valued particle number to an integer number of particles that are lost and gained from the finite particle population k . We use the same subscript and superscript notation as previously. To determine the discrete particle gain and loss sets, we discretize the deterministic real-valued particle number gains and losses of Eq. (16) by applying binomial (or multinomial) sampling of the particle set, where a binomial sample of the finite particle set with a probability p of selecting each particle is denoted by Binom( , p) for 0 ≤ p ≤ 1. For example, the discretization of a real-valued particle number from cell k + 1 to k, C G,s k+1,k , with coefficient p G,s k+1,k is given by which is a stochastic discretization of the corresponding deterministic equation C G,s k+1,k = p G,s k+1,k C s k+1 . From Eq. (16), the particle numbers involved in transport from grid cell k to grid cell k + 1 are Comparing Eqs. (18) and (19), we see that these gain and loss numbers are generally different from one another due to two factors, namely the difference in the vertical grid cell sizes and the difference in the computational volumes containing the particles as illustrated in Fig. 3. When we sample the sets of gain and loss populations, we wish to ensure that the gain and loss populations have as many particles in common as possible. This is done in order to minimize particle duplications and removals. To accomplish this, it is convenient to define the ratio of loss to gain by This ratio contains the quantities that cause the particle number lost from grid cell k to k + 1 to be different from the particle number gained by grid cell k + 1 from k. Note that To construct particle sets that are as similar as possible, we manipulate the gain and loss terms, C G,s k,k+1 and C L,s k,k+1 , as follows: where the maximum transport term C T,s k,k+1 and transport probability p T,s k,k+1 are defined as This ensures that C G,s k,k+1 and C L,s k,k+1 are always less than or equal to C T,s k,k+1 , with at least one of them being equal. Applying binomial sampling to Eq. (26) for each grid cell simultaneously, we obtain the finite transport particle sets. We can sample particles transported in both directions from particle population s k as where each is a finite set of particles, T,s k,k+1 is the set of particles which are candidates for addition to k + 1 and removal from k , T,s k,k−1 is the set of particles which are candidates for addition to k − 1 and removal from k , and Mult ( , p 1 , p 2 , p 3 ) is a multinomial distribution that samples particles into three subpopulations according to the three probabilities, which must sum to 1: p 1 +p 2 +p 3 = 1. In practice, this is evaluated as the equivalent series of binomials 4064 J. H. Curtis et al.: A single-column particle-resolved model for simulating aerosol mixing state Figure 3. Schematic of the transport process for particle number from grid cell k to k + 1. Panel (a) shows the scenario when the particle number removed from grid cell k is equal to the particle number added to grid cell k + 1. Panel (b) shows the scenario when the particle number lost from grid cell k is not equal to the particle number gained by k + 1, which arises when the gain and loss probabilities are unequal due to differences in grid cell sizes and computational volumes. by G,s k,k+1 ∼ Binom T,s k,k+1 , min 1, 1/γ k,k+1 , L,s k,k−1 ∼ Binom T,s k,k−1 , min 1, γ k,k−1 .
Some of the particles initially sampled into the transport sets T,s k,k+1 and T,s k,k−1 are not lost and will remain in cell k: Finally, the sets above can be combined to give the new set of particles s+1 k in cell k by where is the multiset sum.

Selection of sub-cycle time step
To maintain numerical stability with the explicit finitevolume scheme, sub-cycle time steps are taken for vertical transport that differ from the model time step for other processes. Within each sub-cycle time step t T , particles will only transition between immediate grid cell neighbors. However, within a full model time step t, particles may be transported more than one grid cell away.
To determine an appropriate sub-cycle time step, the transfer rates are first computed for all grid cells in the column. The sub-cycle time step must be chosen such that the total particle transition probabilities are always in [0, 1]. To achieve this, the critical time step is taken to be We then determine the number of time steps required to reach the full model time step t by where · is the ceiling function, and defines the sub-cycle time step for the transport as This ensures that the sum of the transport probabilities computed for time step t T is always in [0, 1]. That is, ( t T / t)(p T,s k,k+1 + p T,s k,k−1 ) ≤ 1 for all k. To see this, we observe that Eqs. (44) and (45) imply that t T ≤ t c . From Eq. (43) we have that t c p T,s k,k+1 ≤ t/2, and similarly for p T,s k,k−1 , giving the desired result. Note that scaling the probabilities by the ratio t T / t is the same as computing the probabilities with time step t T because the probabilities are linear in t.
The complete algorithm for stochastic turbulent transport of finite particle sets is given in Algorithms C1 and C2 (Appendix C1). The Binom( , p) binomial samples reflect the fact that each particle has an equal probability of being transported. A Binom( , p) sample can be implemented by first sampling a scalar binomial function n ∼ Binom(N p , p), where N p = | | is the number of particles in population , and then choosing n particles uniformly from .

Rebalancing of computational particle number
During a given simulation, the number of computational particles changes as particles are added due to emission, are transferred from one grid cell to another due to turbulent transport, and are removed by coagulation and dry deposition. When the number of computational particles falls below half the initial prescribed number in a given grid cell, in order to maintain an adequate statistical sample, we duplicate every particle and double the computational volume. When the number of particles is twice the initially prescribed number, in order to alleviate the higher computational cost, half the computational particles are discarded and the computational volume is halved. This strategy has been previously used for particle populations in the zero-dimensional PartMC-MOSAIC box model (Riemer et al., 2009) as well as other Monte Carlo simulations for particle dynamics (Efendiev and Zachariah, 2002;Maisels et al., 2004).

Aerosol dry deposition algorithm
Particles near the surface are subject to removal by the process of dry deposition. This is parameterized by evaluating a dry deposition velocity for each particle in the aerosol population of the lowest grid cell 1 . The parameterization presented in Zhang et al. (2001) is applied here on a per-particle basis. The dry deposition velocity is dependent on the perparticle diameter and density, in addition to the meteorological conditions and surface characteristics. Given a dry deposition velocity for particle i, denoted V d,i , a loss rate probability, d,i , can be expressed as where z 1 is the lowest model layer thickness. Algorithm C3 in Appendix C shows the procedure for the removal of particles from the particle population 1 , where each particle i is tested for removal with the associated dry deposition loss rate probability d,i .

Verification of the stochastic particle algorithms
In the following section, we will present separate numerical verifications of the new algorithms for particle transport by turbulent diffusion (Sect. 4.1) and particle removal by dry deposition (Sect. 4.2). These verification scenarios are motivated by atmospheric applications, but have somewhat artificial grid structures (see Fig. 4), chosen to enable smooth refinement studies. Therefore, the numerical values in this section will differ from those in actual atmospheric case studies.

Particle transport by diffusion
For verification of the stochastic transport method of particles presented in Algorithm C1, the particle transport code was implemented to solve the one-dimensional diffusion equation given by where N (z, t) is total aerosol number concentration as defined by Eq. (11). Here, the eddy diffusion coefficient K h is taken to be constant in both time and space so that the equation can be solved analytically. To isolate diffusion, all other aerosol processes that may contribute to the evolution of the aerosol state were excluded from the simulation. Reflective boundary conditions were imposed at the surface and at the top of the domain. The model was initialized with an instantaneous area source in the x − y plane, with an initial thickness z, centered at altitude z and with uniform perturbation number concentration N 0 and a background number concentration N back . The analytical solution for comparison to model results is where the method of images was applied to impose boundary conditions, with imaginary sources at z a and z b . In this diffusion test case the background particle concentration was N back = 3.2 × 10 3 cm −3 . The instantaneous finite-thickness particle cloud was placed at the altitude z = The variable grid spacing was determined by with C 0 as a domain scaling parameter given by C 0 = 15 000/(3n z ). The grid cell edges are determined by for k = 1, . . ., n z . Figure 4 shows the variation in vertical spacing of grid cells resulting from Eq. (50) where the small grid cells, located at 7500 m, are approximately half as large as the largest grid cells, located at the domain edges. We use a sequence of grids, each with twice as many grid cells as the last, giving n z = 10, 20, 40, 80. Figure 5 shows the aerosol particle number concentration evolution for a single simulation with n z = 20 grid cells, using N p = 10 4 computational particles. The simulated number concentration is compared to the analytical Gaussian solution and shows good agreement. The noise in the number concentration is a result of stochastic sampling and could be further reduced by averaging several independent simulations to form an ensemble mean.
To verify the convergence of the transport algorithm to the analytical solution, we quantified the error in the total number concentration for ensemble member j by the weighted L 2 error: where n z is the total number of grid cells,N t,j k is the stochastic number concentration for grid cell k at time t for ensemble member j ,N t k is the average analytical number concentration over grid cell k, and z k is the grid cell size of cell k. The average analytical number concentration of grid cell k is given bȳ where N (z, t) is the number concentration given by the analytical solution Eq. (48). The total relative error for ensemble member j is given by The ensemble mean errorē t is given bȳ (e t,j ) 2 (54) and the standard deviation in the error σ t e given by (55) Figure 6 shows the error convergence behavior as the number of computational particles, N p , and average grid cell size, z, vary. The average grid cell size is given by the domain height divided by the number of grid cells, so z = (15 000 m) /n z . The ensemble size was n run = 20 and the number of computational particles N p ranged from 20 to 10 7 . Error bars represent the 95 % confidence interval.
We expect that the stochastic particle solutionN t k will converge to the finite-volume solution N t k (see Eq. B7) as N p → ∞. In turn, we expect N t k to converge to the analytical solutionN t k as z → 0. Thus, we anticipate the convergencē To understand the rates of convergence, we decompose the error as O z 2 finite-volume error→0 as z→0 (58) Here we see that the total error is bounded by the stochastic and finite-volume errors. The stochastic error is O(1/ N p ), while the finite-volume error is O( z 2 ) since the spatial discretization is second order accurate. Figure 6a shows convergence of total error for fixed z as N p → ∞, with the expected −1/2 slope until the finitevolume error dominates. Figure 6b shows the convergence as the grid size z decreases for large N p . Each data point in Fig. 6b corresponds to the converged value of a line in Fig. 6a, taken with N p = 10 7 . In this figure we see the expected slope of 2.

Particle dry deposition
To verify Algorithm C3 for dry deposition we developed a test case that only considered the removal of particles by dry deposition. The simulation was initialized with two monodisperse particle populations with different diameters, 1 and 10 µm, and identical densities of 1800 kg m −3 . The number concentration of the 1 µm particle population was 10 3 higher than the number of concentration of the 10 µm population, so that the mass concentrations of the initial particle populations were equal. The dry deposition velocities of the two populations did not change over time due to the absence of coagulation and condensation. Each population was expected to decay at a rate based on their computed dry deposition velocities. As a result, the deposition process can be represented as a first-order decay equation where M D is the aerosol mass concentration of a given population of particles with diameter D. The loss rate of the particles, λ D , for population with diameter D is given by the deposition velocity of particles in that population, V d,D , and the reference height z ref . The analytical solution is where M D,0 is the initial aerosol mass concentration and M D (t) is the aerosol mass concentration at time t. Figure 7 shows the evolution of initial identical mass concentrations for two particle populations, one with particles with diameter of 1 µm and the other with diameter of 10 µm. The results are an average of 10 independent model runs. The average mass concentration and 95 % confidence interval are shown to be in good agreement with the analytical concentrations as given by Eq. (60). The mass concentration associated with particles with diameter of 10 µm decays more quickly because 10 µm particles have a higher settling velocity than 1 µm particles. Particles with diameter of 1 µm experience a very slow loss in mass concentration as a result of ineffective removal by any of the processes represented by dry deposition. These particles are too large to be removed effectively by Brownian diffusion and too small to be removed by gravitational settling.
5 Application of a single-column model with an idealized scenario

Setup of idealized scenario
We constructed an idealized scenario to illustrate the model capabilities of WRF-PartMC-MOSAIC-SCM. The scenario is similar to the box model study presented in Riemer et al. (2009) and Zaveri et al. (2010), but for the first time we now gain insight into aerosol mixing state as it varies spatially with altitude. We simulated a 48 h episode, starting at 06:00 local standard time (LST). Initial gas mixing ratios were based on initial conditions given by Riemer et al. (2009) and decreased linearly with height to 3.5 km. Gas-phase emissions were specified only at the surface and were also based on the urban plume case described in Riemer et al. (2009), adapted from the Southern California Air Quality Study (SCAQS) simulation (26-29 August 1988 period) of Zaveri et al. (2008). Table 1 shows the initial aerosol distributions and aerosol emissions used in this scenario with two aerosol modes, an Aitken mode and an accumulation mode. Both aerosol modes consisted of particles that contained ammonium sulfate and primary organic aerosol. Initial aerosol number concentration was constant with height.
The model allows for the inclusion of aerosol emissions within any grid cell in the column and has flexibility in the choice of the parameters of the size distribution as well as the particle composition of the emitted particles. Carbonaceous aerosols were emitted at the surface from three different sources: diesel vehicles, gasoline vehicles, and meat cooking. Due to the importance of the timing of atmospheric turbulent mixing and emissions, we applied a diurnal cycle to the particle emission rates. This is in contrast to Riemer et al. (2009), where the particle emission rates were held constant with time. The gasoline and diesel emission source strengths were varied over time by redistributing the mean aerosol emissions from Riemer et al. (2009) based on the weekday traffic distribution fractions as described in Marr et al. (2002). The resulting 48 h time series used for this scenario are shown in Fig. 8. We consider this set of simplified surface emissions, as the underlying assumption of the onedimensional column setup is horizontal homogeneity.
WRF-PartMC-MOSAIC-SCM was initialized with 59 vertical levels, logarithmically spaced with 16 levels in the lowest 1 km of the domain. The Mellor-Yamada-Janjic (MYJ) planetary boundary layer scheme (Janjic, 1994) was used to model turbulent transport and parameterize the diffusion coefficient for the particle transport scheme. The presented model formulation requires the use of local boundary layer schemes such as MYJ and Mellor-Yamada-Nakanishi-Niino (MYNN) Niino, 2006, 2009). Nonlocal schemes such as the Asymmetric Convection Model 2 Scheme (ACM2) (Pleim, 2007) may be included in future work. The model was initialized with approximately 25 000 computational particles for each level, resulting in approximately 1.5 million particles in the column. The number of computational particles within each level fluctuated during the simulation due to emission, coagulation, transport, and dry deposition, and was restricted to a range between half and double the initial number of particles (12 500 to 50 000) to maintain accuracy while avoiding higher computational costs as described in Sect. 3.1.4.

Aerosol distribution functions
Given the complexities of the multidimensional aerosol distribution, we must project the distribution for purposes of dis- playing results. We take N (D) to be the cumulative number distribution where N (D) is the number of particles per volume that have a diameter less than D. Given the cumulative number distribution, we define the number distribution n(D) as To characterize the aerosol mixing state, we define the perparticle mass fraction of species a as w a,dry = µ a µ dry , where µ a is the mass of species a in a given particle and µ dry is the total dry mass of the particle. Here species a can be a single aerosol species such as black carbon (BC), sulfate (SO 4 ), or nitrate (NO 3 ), or it consists of a group aerosol species such as secondary organic aerosol (SOA). The one-dimensional cumulative number distribution N(w a,dry ) is the number concentration of particles with dry mass fraction of species a less than w a,dry . The corresponding number distribution n(w a,dry ) is defined as n(w a,dry ) = ∂N (w a,dry ) ∂w a,dry .
We can also construct two-dimensional number distributions with respect to dry diameter, D dry , and dry mass fraction of species a, w a,dry . The two-dimensional number distribution n(D dry , w a,dry ) is defined as n(D dry , w a,dry ) = ∂ 2 N (D dry , w a ) ∂log 10 D dry ∂w a , where N (D dry , w a ) is the number concentration of particles that have dry diameter less than D dry and dry mass fraction of species a less than w a .

Simulation results
In this section we present an illustration of model output and will focus on the evolution of the mixing state of blackcarbon-containing particles within the boundary layer. However, before we discuss the results of the aerosol mixing state in detail, we will provide a brief description of the bulk quantities of the scenario. For this 48 h scenario, the temperature and relative humidity varied over time, as simulated by the WRF model and shown in Fig. 9. Figure 10 shows the evolution of O 3 and NO 2 profiles for the 2-day simulation period. During the daytime we observed production of O 3 , with the highest mixing ratio of 137 ppb occurring in the afternoon of the second simulation day and a peak surface O 3 mixing ratio of 132 ppb at 15:00 LST. NO 2 reached a maximum mixing ratio of 38.8 ppb found at the surface at 07:20 LST on the second day. Figure 11 gives an overview of the evolution of aerosol bulk properties. Figure 11a shows the time evolution of vertical profiles of black carbon mass concentration. Blackcarbon-containing particles were emitted at the surface and vertically mixed in the boundary layer by turbulent diffusion. As the stable boundary layer developed around 18:00 LST on each day, black carbon emissions accumulated within that layer with a depth of ∼ 250 m, resulting in higher surface concentrations. A maximum black carbon mass concentration of 4.63 µg m −3 was found at the surface at 08:20 LST as vehicle emissions began to increase and before the mixing layer began to deepen. Later in the morning the boundary layer height increased, allowing black carbon concentrations to be dispersed vertically and become well mixed. Figure 11b-d show the bulk aerosol mass concentrations of nitrate, sulfate, and SOA. Ammonium nitrate formation was responsible for the majority of the aerosol mass, with a peak mass concentration of 27.7 µg m −3 . Maximum sulfate and SOA concentrations were 3.2 and 7.7 µg m −3 , respectively, at the surface. Figure 11e shows the evolution of the total number concentration. Aerosol number concentrations were impacted by emissions, coagulation, deposition, and turbulent transport. During the afternoon, the boundary layer remained well mixed with respect to number concentration. During the nighttime, the aerosol number concentration decreased with time, most noticeably from the top of the stable boundary layer (∼ 250 m) to the top of the residual layer due to coagulation and lack of transport of surface emissions. The maximum number concentration at the surface was ∼ 14 000 cm −3 at 19:30 LST, when vehicle emissions became trapped near the surface as a result of the development of the nocturnal boundary layer. Number concentrations in the stable boundary layer decreased with time overnight due to coagulation and dry deposition in combination with relatively lower emission rates over that period. Figure 11f shows how the number mean wet diameter varied with height and time. The mean wet diameter was largest in the residual layers of each night due to the particle populations containing high amounts of ammonium nitrate, which resulted in water uptake of particles in the high relative humidity environment as indicated in Fig. 9b.
To understand how the mixing state of black-carboncontaining particles evolves in time and with respect to height, Fig. 12 shows the two-dimensional number distribution n(D dry , w BC,dry ) at 06:00 and 12:00 LST on day 2 at heights of 25, 241, and 552 m. Fresh emissions occurred at the surface and appear as horizontal lines, with diesel emissions prescribed as w BC,dry = 70 % and gasoline emissions as w BC,dry = 20 %.
At 06:00 LST on the second day, the horizontal lines representing fresh emissions were most pronounced at the surface. As a result of the stable boundary layer limiting the vertical extent of turbulent mixing, the fresh emissions were contained to levels near the surface. By 12:00 LST, the height of the boundary layer had grown to ∼ 750 m. As a result, particles with high BC mass fraction were vertically transported to higher levels. However, the greatest number concentrations of the fresh particles were still found near the surface. Diagonal band structures of high number concentrations were a result of condensation of nitrate, which gradually shifted particles to larger diameters and lower BC mass fraction.
While Fig. 12 only shows the BC-diameter distribution of the aerosol, the simulation results contain the full highdimensional distribution over all constituent species, thus permitting the calculation of any desired mixing state measures or visualizations.

Computational costs
The WRF-PartMC-MOSAIC-SCM model has computational cost for evaporation/condensation proportional to the number of computational particles, computational cost for coagulation proportional to the number of coagulation events, and computational cost for transport proportional to the number of sampled particles. As a result, the evaporation/condensation of secondary species is the dominant cost of the model, comprising more than 97 % of the computational cost for the simulation presented here. By contrast, the particle transport is computationally inexpensive, typically representing less than 2 % of the cost.
The WRF-PartMC-MOSAIC-SCM model is both computationally and memory intensive and therefore benefits greatly from parallelization. While the host WRF model utilizes a horizontal decomposition that is not applicable to the single-column model, the WRF-PartMC-MOSAIC-SCM model features a vertical domain decomposition where the gas and aerosol domain is distributed across multiple cores. Since the dominant cost is evaporation/condensation, which is a per-particle process and requires no communication with neighboring grid cells, the model scales efficiently even when the domain is decomposed to a single grid cell per core.

Conclusions
In this paper we presented the development and application of the WRF-PartMC-MOSAIC-SCM model. This model, for the first time, resolves the aerosol composition on a perparticle level in an Eulerian single-column domain and couples the aerosol-and gas-phase chemistry with the meteorology.
We developed and implemented two new algorithms, a stochastic aerosol transport algorithm to treat vertical tur-bulent diffusion, and a stochastic aerosol dry deposition algorithm. Both model processes were verified with test cases and performed as expected when compared to analytical solutions. A stochastic sampling strategy, which does not track particle position, was used to reduce the computational burden of particle transport. This method relies on an explicit discretization of the diffusion equation and the use of nearestneighbor diffusion. Potential future model development includes the implementation of other numerical methods for turbulent diffusion, such as higher-order and/or semi-implicit schemes, and non-local boundary layer schemes such as ACM2.
To illustrate the newly coupled model capabilities, an idealized urban scenario was developed. This 48 h simulation showed the evolution of the black carbon mixing state due to coagulation, secondary aerosol formation, particle emission, dry deposition, and turbulent transport. In the presented scenario, freshly emitted diesel and gasoline particles existed in the highest concentrations near the surface where they were emitted. As particles were vertically mixed due to turbulent transport, emitted particles experienced changes in composition due to coagulation with aged particles as well as due to condensation of secondary aerosol species. While we focused on the composition of black-carbon-containing particles to demonstrate the model capabilities, we do store the full composition for each computational particle, so a similar analysis can be made for other aerosol species. Future applications of the model include quantifying the impact of aerosol mixing state on secondary aerosol formation and on climate-relevant aerosol properties, such as aerosol absorption and CCN concentration, and comparing these findings to existing studies (Matsui et al., 2013;Zhang et al., 2014;Zhu et al., 2016a).
Code availability. The box model version of PartMC is available from http://lagrange.mechse.illinois.edu/mwest/partmc/ under the GNU General Public License (GPL). The version of WRF-PartMC-MOSAIC-SCM presented here, excluding coupling with MOSAIC, is available upon request from Nicole Riemer (nriemer@illinois.edu) and is available under the GNU GPL. To couple chemistry, MOSAIC may be obtained upon request from Rahul Zaveri (rahul.zaveri@pnl.gov).

Appendix A: List of symbols
See Table A1 for a list of symbols used in this paper.
Given that we seek to represent the loss of particles from each grid cell to the neighboring cells and the gain of particles from those neighboring cells, we need a form of the equation in terms of gains and losses rather than a flux formulation. We rearrange Eq. (B8) to isolate gain and loss terms of the number concentration N s k , giving (B9) Algorithm C2 Stochastic aerosol particle transport from a grid cell.
function PARTICLESAMPLE Π, p G , p L , pprev → Π G , Π R , Π U , pnew Input: Π is the starting particle population of the source grid cell p G is the particle gain probability p L is the particle loss probability pprev is the sum of previous transfer probabilities Output: Π G is the particle population to be gained by the destination grid cell Π R is the particle population that is temporarily removed from the source grid cell Π U is the unsampled particle population in the source grid cell pnew is the updated sum of transfer probabilities pnew ← pprev + p T end function 33 Algorithm C3 Stochastic aerosol particle dry deposition for one time step Input: Π 0 1 is the initial particle population in the grid cell closest to surface at the start of the time step { d,i} N p,1 i=1 are the particle loss rate probabilities for each particle in Π 0 1 ; see Eq. (46) Output: Π1 is the updated particle population with particle removals due to dry deposition after the time step