Effects of volatility, viscosity, and non-ideality on particle–particle mixing timescales of secondary organic aerosols

Abstract

Different populations of aerosol are constantly mixed throughout the atmosphere. Large-scale models often assume no particle–particle mixing or fast mixing among aerosol populations, so that they stay externally mixed or instantaneously form internal mixtures. We apply the kinetic multilayer model of gas–particle interactions (KM-GAP) to simulate the evaporation of semi-volatile species from one particle population and partitioning into another population with various phase states and nonideal mixing conditions. We find that the particle–particle mixing timescale (τ_mix) is prolonged when the semi-volatile species transport to a population in which it is miscible, as more mass must be transported. Extremes of volatility prolong the τ_mix, as low-volatility species evaporate slowly, while high-volatility species condense slowly. When the bulk diffusivities of the two populations are greater than 10⁻¹⁵ cm² s⁻¹, semi-volatile species mix rapidly; otherwise, the τ_mix can be prolonged beyond 1 h. We apply KM-GAP to particle–particle mixing experiments of H-toluene SOA into D-toluene SOA and limonene SOA, showing that τ_mix is prolonged when toluene SOA is highly viscous, while initial partitioning of gas phase semi-volatile species from toluene SOA into limonene SOA is rapid because of the low viscosity of limonene SOA. Simulations of mixing toluene SOA and β-caryophyllene SOA indicate that the apparent discrepancy of limited mixing under conditions where both are predicted to have low viscosity are explained by limited miscibility of the semi-volatile components. Our study demonstrates that particle–particle mixing timescales are affected by a complex interplay among volatility, diffusion limitations, and non-ideal miscibility.

Editor:

• Cari Dutcher

Introduction

Organic aerosols constitute a major fraction of aerosol mass in the atmosphere, affecting climate, air quality, and public health (Jimenez et al. 2009; Pöschl and Shiraiwa 2015). Organic aerosols contain a huge number of different species from a variety of primary and secondary sources (Goldstein and Galbally 2007; Zhang et al. 2007; Prather, Hatch, and Grassian 2008). Secondary organic aerosol (SOA) is formed via oxidation of volatile organic compounds to form semi-volatile and low-volatility products and subsequent condensation of these oxidation products onto preexisting particles. SOA is a complex mixture of species with varying volatilities and degrees of oxidation. The physicochemical properties of SOA depend on chemical composition as determined by the precursor and the formation conditions, as well as chemical evolution and mixing state with different aerosol populations in the atmosphere.

Measurements have shown that SOA can exist in a highly viscous semisolid or glassy state under low relative humidity and temperature (Liu et al. 2017; Reid et al. 2018; Kasparoglu et al. 2021; Song et al. 2021, 2022; Jeong et al. 2022; Kaur Kohli et al. 2023; Sheldon et al. 2023). Modeling studies have also shown that SOA should be highly viscous under dry and/or cold conditions in the troposphere, influencing multiphase chemical processes and cloud activation (Mu et al. 2018; Shiraiwa et al. 2017; Shrivastava et al. 2017; Schmedding et al. 2020; Rasool et al. 2021). This high particle viscosity corresponds to small bulk diffusion coefficients according to the Stokes–Einstein relation, inducing diffusion limitations within the particle (Ullmann et al. 2019; Shiraiwa et al. 2011, 2013). While SOA is often modeled with instantaneous equilibrium with its surroundings, bulk diffusion limitations can induce long equilibration timescales for partitioning into viscous particles (Shiraiwa and Seinfeld 2012; Booth et al. 2014; Schervish and Shiraiwa 2023). The phase state of SOA also regulates the evolution of size distributions, growth rates, heterogeneous chemistry, ice nucleation, and long-range transport of pollutants (Wang et al. 2012; Berkemeier et al. 2014; Zhang et al. 2018; Zaveri et al. 2020; Bertozzi et al. 2021; Gkatzelis et al. 2021; He et al. 2021; Octaviani et al. 2021; Lei et al. 2022; Shen et al. 2022; Shrivastava et al. 2022; Song et al. 2022; Rasool et al. 2023). For example, it has been shown that viscous organic coatings can facilitate the long-range transport of toxic polycyclic aromatic hydrocarbons (Shrivastava et al. 2017; Mu et al. 2018).

Atmospheric aerosols consist of a complex mixture of organics, inorganics, and water. Mounting evidence indicates that aerosol particles can phase separate into an aqueous inorganic core and a hydrophobic organic shell, indicating a low degree of miscibility between these two phases (Freedman 2020; Huang et al. 2021; Li et al. 2021). The organic aerosol may also phase separate into a high O:C, hydrophilic core and a low O:C, hydrophobic shell even in the absence of inorganics (Song et al. 2012, 2018; Freedman 2017; Huang et al. 2021; Mahrt et al. 2022). In contrast, different species may be favorably miscible, leading to an enhancement in mixing compared to ideality (Kim, Sartelet, and Couvidat 2019; Rusumdar et al. 2020; Roston 2021). According to calculations, these non-ideal interactions are represented with an activity coefficient, which can be orders of magnitude away from 1 (Zuend et al. 2008, 2011; Hyttinen et al. 2020).

Different aerosol populations are constantly being mixed throughout the atmosphere, e.g. indoor-to-outdoor mixing, concentrated pollution plumes mixing with ambient aerosol, or transport and mixing of ambient aerosols from different sources (Riemer et al. 2019; Cummings et al. 2021; Sedlacek et al. 2022). The introduction of two different aerosol populations to each other will affect each population’s equilibrium as semi-volatile components of each population will evaporate and mix with the other population. The timescale for each population to achieve a new equilibrium will depend on the properties of each aerosol population as well as those of the semi-volatile components present (Li and Shiraiwa 2019; Schervish and Shiraiwa 2023). Experiments aimed at probing these mixing timescales have shown that mixing between different SOA components can be rapid at high RH but can be hours at low RH (Robinson, Saleh, and Donahue 2013; Saleh, Donahue, and Robinson 2013; Liu et al. 2016; Ye et al. 2016, 2018; Preston and Zuend 2022). In contrast, Habib and Donahue (2022) observed little kinetic limitations for mixing erythritol and various sugars even at low RH and cold temperatures. Additionally, while Ye et al. (2016, 2018) observed rapid mixing of toluene SOA with deuterated toluene SOA, isoprene SOA, and limonene SOA at high RH, limited mixing was observed between toluene SOA and β-caryophyllene SOA at high RH. These apparent discrepancies suggest complex effects of different volatility distributions of different components, diffusion limitations due to higher particle viscosity at low RH, and possibly miscibility limitations (Wilson et al. 2015). Marcolli et al. (2004) investigated the mixing of two aerosol populations and concluded that equilibrium is rapid for a semi-volatile species present in different amounts in the two populations but only considered diffusion limitations in the gas phase and did not consider bulk diffusion limitations or effects of non-ideality (Marcolli et al. 2004).

Large-scale models often represent no particle–particle mixing or fast mixing between different aerosol populations relative to other atmospheric processes to rapidly form an internal mixture, where aerosols of a given size will contain the same composition because of computational limitations. While the latter assumption is likely valid when particle–particle mixing timescales are less than a typical model timestep between 20 min and 1 h, previous work has shown that the equilibration timescales can be hours or even days for viscous particles (Li and Shiraiwa 2019; Schervish and Shiraiwa 2023). Therefore, understanding under what conditions and what properties of aerosols lead to long mixing timescales provides crucial insight for an improved representation of aerosol dynamics in large-scale models (Mai et al. 2015; O’Meara, Topping, and McFiggans 2016). Here, we investigate the mixing of two particle populations with varying viscosity and probe their mixing timescale. We will simulate a semi-volatile species initially in one population with varying volatility and varying miscibility in the other population.

Methods

The kinetic multilayer model of gas–particle interactions in aerosols and clouds (KM-GAP) is used to simulate the mixing of two different particle populations. The details of KM-GAP, including rate equations and kinetic parameters, can be found in previous publications (Shiraiwa et al. 2012, 2013), so here we will give a brief description and the specific implementation in this study. KM-GAP explicitly treats the processes of gas-phase diffusion, adsorption and desorption, surface-bulk exchange, and bulk diffusion within the particle. KM-GAP consists of a gas phase, a near-surface gas phase, and, historically, one particle population with a sorption layer, a near-surface bulk layer, and a number of bulk layers. KM-GAP does not require simplifying assumptions about steady-state conditions and homogeneous radial mixing, which are often applied in other aerosol models. Mass transport processes are treated as fluxes, and rates of chemical reactions are treated in differential equations to describe the mass or number balance of each molecule in each model layer of gas and particle phases. The size change of particles and the temporal evolution and spatial profile of the concentration of each species can be modeled by solving ordinary differential equations in Matlab software.

In this implementation, we simulate two particle populations, each with these layers and one shared gas phase. As shown in Figure 1, we simulate two particle populations that are initially externally mixed and introduced to each other at time t = 0. The two populations are made up of different nonvolatile organic species, with one population (here referred to as population 1) containing a semi-volatile species. The semi-volatile species, present initially only in population 1, mixes ideally with the nonvolatile component of population 1 with an activity coefficient (γ₁) of 1. It has a diffusivity in population 1 of D_b1 = 10⁻⁸ cm² s⁻¹. The semi-volatile species has a total concentration of 4 µg m⁻³ and is initially partitioned to the gas phase based on its volatility so that the semi-volatile species is in equilibrium between the gas phase and the particle phase of population 1. We assume a surface accommodation coefficient of 1 for the semi-volatile species for each particle population (Julin et al. 2014; Domaros et al. 2020). In this work, we vary the pure-compound saturation mass concentration (C⁰) or volatility of the semi-volatile species in the range C⁰ =10⁻⁴–10⁴ µg m⁻³, which corresponds to the vapor pressure of ∼10⁻⁹–10⁻¹Pa. This range encompasses the volatility of many atmospherically relevant organic species, which are distributed into volatility bins separated by powers of 10, including LVOCs (C⁰ = 10⁻⁴–0.1 µg m⁻³), SVOCs (C⁰ = 1–100 µg m⁻³), and IVOCs (C⁰ = 10³–10⁶ µg m⁻³). Lower volatility species can likely be considered essentially nonvolatile on atmospherically relevant timescales, and higher volatility species are unlikely to exist in large quantities in the particle phase (Donahue et al. 2011).

Figure 1. Schematic of the two populations implementation in KM-GAP. Both populations consist of a near-surface bulk layer (nsb) and the number of bulk layers (b_1-L; L = 10 in this study). Additionally, each particle has a sorption layer containing adsorbed and desorbing species and a near-surface gas phase to account for potential gas phase diffusion, which is not pictured. There is one shared gas phase between the two populations. Each population consists of different nonvolatile species, while population 1 initially contains semi-volatile species in equilibrium with the gas phase.

Both populations are monodisperse with initial diameters of 200 nm, and a total particle number concentration of 10⁴ cm⁻³. This size is chosen as it is roughly the mean diameter of ambient SOA particles in the volume size distributions (Martin et al. 2010; Pöschl et al. 2010; Riipinen et al. 2011). Both nonvolatile species, making up the bulk of each population, and the semi-volatile species have a molecular weight of 200 g mol⁻¹. Therefore, initially, there is exactly the same amount of mass (∼40 µg m⁻³) and number of molecules in each population.

To investigate how particle–particle mixing timescales are affected by different properties of the aerosol and semi-volatile species, the volatility of the semi-volatile species, its miscibility in population 2, and the viscosity of population 2 are all varied. The latter two properties are adjusted via the activity coefficient (γ₂) and the bulk diffusivity (D_b2) of the semi-volatile species in population 2, respectively. γ₂ is varied between 10⁻³ and 10³. While the miscibility of atmospherically relevant organic species leads to activity coefficients often within an order of magnitude of unity, activity coefficients can deviate from 1 by orders of magnitude, for example, for partitioning of highly hydrophilic compounds such as glyoxal and nonpolar compounds like alkanes in the aqueous phase (Erdakos et al. 2006; Zuend et al. 2008, 2011; Shiraiwa et al. 2013; Hyttinen et al. 2020; Wallace et al. 2021; Hilditch et al. 2023). Hilditch et al. (2023) measured activity coefficients for dicarboxylic acids in organic matrix droplets and found activity coefficients between 0.1 and 10. For polycyclic aromatic hydrocarbon benzo(a)pyrene, AIOMFAC predicts activity coefficients greater than 10³ at high water activity (Zhou et al. 2019). For the viscosity of population 2, we simulate a liquid particle (D_b2 = 10⁻⁸ cm² s⁻¹), a semi-solid particle (D_b2 = 10⁻¹⁵ cm² s⁻¹), and an ultraviscous semi-solid particle (D_b2 = 10⁻¹⁸ cm² s⁻¹). These diffusivities are chosen to represent the range of viscosities and diffusivities seen in atmospherically relevant organic or inorganic aerosol (Kiland et al. 2019).

In these simulations, species with high activity coefficients condense into a particle population (population 2) in which they are not miscible, and we assume that they still form one phase in this work. In reality, these species may form their own separate phase within population 2. For the formation of core-shell morphology with the condensing species making up the shell, the condensing species would not need to diffuse to the particle bulk but just need to be accommodated to the particle surface. In this case, the equilibration timescale would be controlled by gas-phase diffusion and surface accommodation. The situation where a condensing species that is immiscible in the shell but mixes ideally with a core phase has been investigated in detail in Schervish and Shiraiwa (2023) with various shell viscosities.

It is important to define timescales carefully. Depending on the quantity of interest, the relevant timescale can differ dramatically. For example, if the quantity of interest is the timescale for the total organic aerosol mass of a single internally mixed population to relax to an equilibrium value with respect to semi-volatile gas-phase constituents, the timescale is mostly determined by the aerosol condensation sink (Saleh, Donahue, and Robinson 2013). However, the timescale for individual molecules to move from one population to another can be very different (as evaporation timescales can be very long even, without diffusion limitations) (Donahue, Chuang, and Schervish 2019). That is the timescale we consider here; the particle–particle mixing timescale (τ_mix) discussed in this work is the e-folding timescale for the concentration of the semi-volatile species to be within 1/e of its equilibrium value relative to its initial value. While each population should reach full equilibrium with the gas phase at the same time, the e-folding timescale can be different for each population. As population 1 is initially in equilibrium with the gas phase, less material needs to evaporate from population 1 compared to the amount that needs to condense to population 2. Therefore, when the bulk diffusivity of semi-volatile species in population 1 is greater than or equal to that in population 2, the equilibration timescale in population 2 will be longer. Hence, the equilibration timescale between population 2 and the gas phase represents the timescale for the whole system to achieve equilibrium and will be referred to as the particle–particle mixing timescale, τ_mix.

As noted in Schervish and Shiraiwa (2023), equilibrium between the particle phase and the gas phase can be achieved long before the particle becomes homogeneous with no concentration gradients in the particle bulk. This occurs when the particle bulk is viscous, and mass can build up in the outer layers of the particle. While the whole system achieves full equilibrium at a single time, this buildup allows the particle to be within 1/e of its equilibrium value relative to its initial value (the condition for τ_mix) long before internal concentration gradients have dissipated. In this work, the $\tau$_mix represents the timescale for the two populations to achieve equilibrium with the gas phase, not necessarily for internal particle-phase concentration gradients to dissipate. The chosen τ_mix is most representative of gas–particle equilibrium in the context of large-scale aerosol modeling and laboratory-probed timescales for the exchange of semi-volatile species between two populations, as both models and experiments often do not resolve or measure concentration gradients within a particle.

Results and discussion

Liquid population 2

We calculate τ_mix when population 2 is liquid with D_b2 = 10⁻⁸ cm² s⁻¹. The pure-compound saturation mass concentration and activity coefficient of semi-volatile species in population 2 are varied from C⁰ = 10⁻⁴–10⁴ µg m⁻³ and γ₂ = 10⁻³–10³. Figure 2a shows an increase in τ_mix with decreasing C⁰ for all cases of γ₂. This occurs for two reasons: (1) as C⁰ decreases, more material must be transported to population 2 for equilibrium to be achieved, and (2) at lower C⁰, evaporation from population 1 becomes slower. Li and Shiraiwa (2019) showed that in an open system, evaporation timescales for low volatility species (C⁰ < 1 µg m⁻³) are consistently longer than 1 h, even in liquid particles. While the simulation condition does not represent a truly open system, population 2 acting as an absorbing medium for the semi-volatile species means that much more needs to evaporate than if population 1 was alone in a closed system, especially when γ₂ is low.

Figure 2. Particle–particle mixing timescale (τ_mix) of two populations where population 1 is liquid (D_b1 = 10⁻⁸ cm² s⁻¹) and population 2 is (a) liquid (D_b2 = 10⁻⁸ cm² s⁻¹), (b) semi-solid (D_b2 = 10⁻¹⁵ cm² s⁻¹), and (c) ultraviscous semi-solid (D_b2 = 10⁻¹⁸ cm² s⁻¹).

The equilibrium concentration of semi-volatile species in population 2 decreases with increasing γ₂, leading to a decrease of τ_mix because less mass needs to be transported to population 2. When γ₂ $\le$ 1 and C⁰ < 10³ µg m⁻³, τ_mix is roughly independent of γ₂ because the equilibrium concentration in population 2 is roughly independent of γ₂. Figure 3 shows the equilibrium concentration in population 2, which is a function of both C⁰ as well as γ₂, but independent of the viscosity of population 2. The equilibrium concentration increases with decreasing C⁰ and decreasing γ₂. When γ₂ = 0.1, there will be approximately an order of magnitude more semi-volatile species in population 2 than in population 1 at equilibrium. When C⁰ = 10 µg m⁻³ and γ₂ = 0.1, approximately 90% of the available semi-volatile species in the system will reside in population 2 at equilibrium. While a further decrease of γ₂ will still lead to an increase in the equilibrium concentration in population 2, the increase is small because of the limited amount of semi-volatile species in the system.

Figure 3. The final equilibrium number concentration of the semi-volatile species in population 2 with different volatilities (C⁰) and activity coefficients of semi-volatile species in population 2 (γ₂).

Overall, when both populations are liquid, equilibrium will be achieved in less than an hour for all cases of miscibility for the semi-volatile species with C⁰ > 0.1 µg m⁻³. Specifically, within the SVOC range, τ_mix is less than or equal to approximately 30 min for all cases of γ₂, and thus, an equilibrium assumption is appropriate for SVOC mixing between liquid populations. For LVOCs, τ_mix may range from a couple of hours for C⁰ = 0.1 µg m⁻³ to several months for C⁰ = 10⁻⁴µg m⁻³ and γ₂ < 1.

Semi-solid population 2

When population 2 is semi-solid with D_b2 = 10⁻¹⁵ cm² s⁻¹, τ_mix is higher than the liquid case. There is very little difference in the trends of τ_mix for LVOCs and SVOCs compared to the liquid case. With only moderate diffusion limitations in population 2, the mixing for the LVOCs is driven mainly by slow evaporation from population 1 and the equilibrium concentration in population 2 (Figure 3).

The remarkable difference from the liquid case is in τ_mix for IVOCs, where τ_mix is controlled by slow condensation and the equilibrium concentration in population 2, as seen in Figure 2b. The uptake of a species into the particle bulk is controlled by competition between bulk diffusion and desorption. As species remain near the surface for longer when bulk diffusion is slow, desorption and re-evaporation are more significant. Therefore, when there is a diffusion limitation in population 2, τ_mix increases with decreasing γ₂ in the IVOC region compared to the liquid case. The increase is more pronounced with low γ₂ because of the higher equilibrium concentration of the species in population 2. Equilibration will take longer, as in addition to the other effects, IVOCs need to move through a more viscous particle bulk in the presence of bulk diffusion limitations.

It is notable that τ_mix is still less than an hour for SVOCs for all cases of miscibility. In the IVOC range, mixing is rapid if the semi-volatile species is immiscible with population 2 with high γ₂, but τ_mix increases with decreasing γ₂ up to ∼5 days when γ₂ = 10⁻³ because of the increase in the equilibrium concentration in population 2 with decreasing γ₂; even when γ₂ = 1, τ_mix is approximately 1.3 h with C⁰ = 10⁴ µg m⁻³. At lower C⁰ (LVOCs and ELVOCs), τ_mix is similar to those of the liquid population 2 case and again can reach several months, as τ_mix is controlled by very slow evaporation from population 1.

Ultraviscous semisolid population 2

When population 2 adopts an ultraviscous semi-solid state with D_b2 = 10⁻¹⁸ cm² s⁻¹, particle–particle mixing is slower in all cases for C⁰ and γ₂. Slower diffusion leads to a very long τ_mix for IVOCs due to slow condensation and for low γ₂ due to the higher equilibrium concentration. An interesting feature is observed for C⁰ between 10⁻² and 10 µg m⁻³ and γ₂ < 0.1 in Figure 2c: while we see the expected trend of decreasing τ_mix when the total equilibrium concentration decreases with increasing γ₂ for γ₂ > 0.1, τ_mix slightly decreases for low γ₂ < 0.1. This occurs due to the accumulation of semi-volatile species in the near-surface bulk in the presence of strong kinetic limitations of bulk diffusion, achieving its’ total equilibrium concentration rapidly. More mass of the semi-volatile species can build up for lower γ₂, and thus, lower γ₂ leads to slightly lower τ_mix in cases where the final equilibrium concentrations are similar. A similar build-up effect was observed in Schervish and Shiraiwa (2023).

This effect does not persist to the IVOC range as eventually it is overtaken by the effect of the difference in the amount of mass needed to be transported, which is larger for the IVOC range as when γ₂ < 1, the equilibrium concentration in population 2 increases substantially with decreasing γ₂ (Figure 3). The equilibrium concentration in population 2 is similar for γ₂ = 10⁻¹ up to C⁰ = 10³ µg m⁻³ and γ₂ < 10⁻² up to C⁰ = 10⁴ µg m⁻³. These C⁰ values correspond to when this reversal of the expected trend based on equilibrium concentration ends in Figure 2c. This effect is also not observed in the LVOC range, as slow evaporation from population 1 limits the mixing timescale in this range. The buildup effect in the region for C⁰ between 10⁻² and 10 µg m⁻³ and γ₂ < 0.1 leads to these miscible species having a τ_mix that is at most 1/3 lower than that when γ₂ = 1. While not an insignificant effect, the τ_mix for these favorably miscible species still exceeds 5 days; thus, this effect does not lead these species to mix rapidly. When population 2 is ultraviscous, τ_mix is always longer than 1 h under investigated conditions of C⁰ and γ₂, indicating that if two aerosol populations are mixing and even just one is ultraviscous semisolid, glassy, or amorphous solid, then the assumption of rapid equilibrium is unlikely to hold.

Viscosity changes

Figure 4a shows the τ_mix of semi-volatile species (C⁰ = 10 µg m⁻³) with various D_b2 and γ₂ in population 2. Population 1 is kept as liquid with D_b1 = 10⁻⁸ cm² s⁻¹. As expected, τ_mix increases with a decrease in D_b2 because of the kinetic limitations of bulk diffusion. Overall, τ_mix is shorter for larger γ₂, as less material needs to be transported. In the presence of diffusion limitations with D_b2 < ∼10⁻¹⁴ cm² s⁻¹, τ_mix can become slightly shorter with decreasing γ₂ when γ₂ < 1 due to the effect of accumulation in near-surface bulk, as explained earlier.

Figure 4. (a) Particle–particle mixing timescale (τ_mix) with various diffusivities and activity coefficients of semi-volatile species (C⁰ = 10 µg m⁻³) in population 2 (γ₂). (b) τ_mix when the bulk diffusivity (D_b) of a semi-volatile species with C⁰ = 10 µg m⁻³ is varied in both populations. In these simulations, the semi-volatile species is ideally miscible in both populations.

All of the above simulation results address only changes in the viscosity of population 2 by keeping population 1 liquid. To investigate the effect of the phase state of population 1, we conducted simulations by varying bulk diffusivity in both populations while assuming that the semi-volatile species (C⁰ = 10 µg m⁻³, still initially present in population 1 and in equilibrium with the gas phase) mixes ideally with both populations (γ₂ = γ₁ = 1). Figure 4b shows τ_mix when D_b in both populations is varied in the range of 10⁻¹⁰–10⁻²⁰ cm² s⁻¹. Note that the equilibration timescale of a population with smaller D_b is systematically longer than a population with higher D_b, while their difference remains within a factor of 2. Hence, we continue to use the equilibration timescale of population 2 to represent the particle–particle mixing timescale (τ_mix), but it is worth noting that the system would take up to twice as long to reach full equilibrium if population 1 is more viscous than population 2. Overall, the τ_mix is very similar if the D_b of populations 1 and 2 are swapped. Particle–particle mixing is driven by evaporation from population 1 and condensation to population 2. For the same D_b, the timescales of evaporation and condensation should be similar, given that condensation is an inverse process of evaporation, as also shown by a previous study (Shiraiwa et al. 2013). Hence, a symmetric distribution of τ_mix is reasonable as the overall particle–particle mixing timescale is limited by a slower process with lower D_b.

A typical timestep of many large-scale aerosol models is 20 min to 1 h (Philip, Martin, and Keller 2016). τ_mix is less than 1 h when both particle populations have D_b > 10⁻¹⁵ cm² s⁻¹, while τ_mix can be several hours to days when both particle populations have lower D_b. Hence, the quick establishment of internal mixing of semi-volatile species would be a valid assumption for low-viscosity particles, but this assumption would be invalid for highly viscous particles. In this case, semi-volatile mixing should be treated kinetically, or an assumption of external mixing would be more appropriate if particles adopt ultraviscous (D_b < 10⁻¹⁷ cm² s⁻¹) or glassy phase states.

Effect of organic aerosol mass concentration

The total organic aerosol mass loading has been shown to affect equilibration timescales for condensation and evaporation, with low mass loadings consistently leading to longer equilibration timescales (Shiraiwa and Seinfeld 2012; Saleh, Donahue, and Robinson 2013; Li and Shiraiwa 2019). While typical chamber experiments are conducted with total organic aerosol concentrations (C_OA) in the range of 10–100 µg m⁻³, ambient organic aerosol mass concentrations in less polluted urban areas are often <10 µg m⁻³ (Hodzic et al. 2020). Therefore, it is valuable to consider lower C_OA concentrations. In the above simulations, both populations have C_OA1 = C_OA2 = 40 µg m⁻³. Here, we maintain all the same properties for the species making up each population as well as the 200 nm particle diameter for each population but change the total particle number concentration to achieve C_OA2 in the range of 0.1–100 µg m⁻³. With this setup, we thus do not capture any potential effect of particle size changes that may occur when particles are formed in the lab or atmosphere to generate different total organic aerosol concentrations but are able to narrow in on the effect of changing only the total aerosol concentration. C_OA1 is fixed at 40 µg m⁻³ as in the previous results, so the total amount of semi-volatile species is the same.

Figure 5 shows the τ_mix for a semi-volatile species with C⁰ from 10⁻³ to 10³ µg m⁻³ when C_OA2 is varied from 0.1 to 100 µg m⁻³. τ_mix for SVOCs decreases with increasing C_OA2 because the condensation sink is larger with a higher surface area at higher C_OA2 (Donahue, Chuang, and Schervish 2019). This decrease is relatively small (less than a factor of 5) within this range of C_OA2. τ_mix for LVOCs are much less affected by the change in mass loading because the increased condensation sink is less effective at decreasing the equilibration timescale for species that condense rapidly, even at low organic aerosol mass loadings. τ_mix for IVOCs is also less affected by C_OA2, as higher volatility compounds mostly reside in the gas phase and their partitioning is less affected by the absorbing mass. These results indicate that the species volatility is a more critical parameter for determining τ_mix. While higher total organic aerosol concentration can shorten the equilibration timescale for SVOCs, the effect of C_OA is minimal for IVOCs and LVOCs.

Figure 5. Particle–particle mixing timescale (τ_mix) when mass loading of population 2 (C_OA2) is varied in the range of 0.1–100 µg m⁻³, while mass loading of population 2 (C_OA1) is held constant at 40 µg m⁻³.

Comparison to laboratory mixing experiments

With this setup of KM-GAP, we investigate experimental results on the particle–particle mixing of two aerosol populations. We simulate the results from Ye et al. (2016, 2018), where two SOA populations were generated separately, allowed to achieve equilibrium, and subsequently mixed. Once mixed, any semi-volatile components of the two populations could evaporate and mix into the other population. In the experiments investigated here, only H-toluene SOA and deuterated toluene SOA contain a significant semi-volatile fraction. The evaporated semi-volatile species from toluene SOA partitions into a different toluene SOA population, limonene SOA, or $\beta$-caryophyllene SOA, representing population 2. In Ye et al. (2016, 2018), SOA particles had mobility diameters between 200 and 700 nm. In the model simulations, the particle diameter was assumed to be 400 nm, with particle number concentrations chosen to achieve the total organic aerosol concentrations of each population in the experimental work.

H-toluene SOA + D-toluene SOA

As H-toluene SOA and D-toluene SOA should mix ideally, these experiments can be used to derive the fraction of SVOC in toluene SOA. To achieve a mass fraction of ∼0.3 of population 1 in population 2 at equilibrium, we set the initial total mass of SVOC in the system to be ∼30% of the organic aerosol concentration of population 1. It is assumed that both H-toluene SOA and D-toluene SOA contain two semi-volatile species with C⁰ of 1 and 10 µg m⁻³. These species are assumed to be present in equal quantities and are initially present in equilibrium with the gas phase. As evaporation from population 2 will affect the experimentally observed property of the mass fraction of population 1 in population 2, we consider this process in simulating H-toluene SOA and D-toluene SOA mixing.

The viscosity of toluene SOA has been measured across wide temperature and RH ranges because of its use as a proxy for anthropogenic SOA. Song et al. (2016) measured toluene SOA viscosity to be >10⁸ Pa s for RH < 20%, and it decreased to 7.8 × 10³ Pa s at ∼50% RH (Song et al. 2016), which is consistent with particle bounce measurements conducted by Bateman, Bertram, and Martin (2015) and Li et al. (2015). Using the Stokes–Einstein equation, we estimated the bulk diffusivity to be 3 × 10⁻¹⁷ cm² s⁻¹ at 10% RH, 1 × 10⁻¹⁵ cm² s⁻¹ at 20% RH, 1 × 10⁻¹³ cm² s⁻¹ at 40% RH, and 2 × 10⁻⁸ cm² s⁻¹ at 90% RH, respectively, which are used as input parameters in our model simulations.

Figure 6a shows the observed mass fraction of H-toluene SOA (population 1) in D-toluene SOA (population 2) by Ye et al. (2016) and the results of the model simulations, showing a very good agreement between measurements and simulations. Ye et al. (2016) conducted experiments under four different relative humidity (RH) conditions. All scenarios were initiated with ∼10% RH. In the first experiment, the RH was maintained at ∼10% throughout the experiment, and there was little to no mixing observed between the two populations. In the second experiment, the RH was increased to 20% at t = 1 h, leading to a small enhancement in mixing. In this experiment, the mass fraction of SVOC from population 1 in population 2 reached ∼0.19 after 4 h. Full equilibrium was not achieved in 4 h under 10% and 20% RH, which is consistent with the highly viscous semi-solid phase state as measured by Song et al. (2016). In this simulation, as well as for all the following simulations, we assumed that the change in bulk diffusivity is instantaneous upon RH increase in the experiment, as previous studies have shown that water diffusion is rapid in the order of seconds or minutes under room temperature depending on relative humidity and particle size (Shiraiwa et al. 2013; Price et al. 2015). Realistically, there may still be a delay in the particle viscosity change as water diffuses through each population; hence, the simulation results show more abrupt changes in mass fraction compared to the observations (Figure 6). In experiment 3, the RH was increased to ∼40% after 2 h, leading to a mass fraction of 0.3 after 4 h. In the final experiment, the RH was increased over the course of the experiment to 90% RH at 1.5 h. The fastest uptake was observed in this experiment, with a mass fraction of 0.3 being achieved in ∼2 h (Figure 6a). The overall agreement between the experimental observations and the simulation results indicates that the limited mixing of toluene SOA can be attributed to its high viscosity at low RH and the increased mixing at higher RH can be attributed to the lower viscosity of toluene SOA as RH increases.

Figure 6. Simulations of toluene SOA mixing experiments and the data from Ye et al. (2016, 2018). Plotted is the mass fraction of SVOCs evaporated from population 1 in population 2. (a) Mass fraction of toluene SOA (population 1) in D-toluene SOA (population 2). RH legend indicates the RH at the end of the experiment. For each experiment, the RH begins at 10%. For RH 20%, the RH is increased at 1 h; for RH 40%, it is increased at 2 h; and for RH 90%, the experimental RH was increased continually, reaching 90% at 1 h, while in the model this is achieved by a stepwise increase in RH. (b) Mass fraction of toluene SOA (population 1) in limonene SOA (population 2) initially at 10% RH with an increase to 30% at 2 h. (c) Mass fraction of toluene SOA (population 1) in $\beta$-caryophyllene SOA (population 2) initially at 10% RH with an increase to 30% at 2 h.

Toluene SOA + limonene SOA

Ye et al. (2018) investigated the mixing of populations of D-toluene SOA and limonene SOA (Ye et al. 2018). As shown in Figure 6b, they observed an initial rapid increase in the mass fraction of toluene SOA present in the limonene SOA followed by a slow increase over the next 2 h at RH < 10%. The RH was then increased to ∼30% after 2 h. The mass fraction of toluene SOA in limonene SOA continued to increase slowly, leveling out at 0.3 after 3 h. They hypothesized that these two SOA populations mix ideally as this is the same mass fraction achieved from mixing toluene SOA with deuterated toluene SOA in Ye et al. (2016) as discussed earlier.

Ullmann et al. (2019) measured bulk diffusion coefficients in limonene SOA using a fluorescent technique to be 5.5 × 10⁻⁹ at 90% RH and 7.1 × 10⁻¹³ cm² s⁻¹ at 33% RH (Ullmann et al. 2019). Petters et al. (2019) derived a viscosity of dry limonene SOA to be ∼10⁶ Pa s, corresponding to D_b of ∼10⁻¹⁴ cm² s⁻¹ for a 0.4 nm diffusing molecule based on the Stokes–Einstein equation. Hinks et al. (2016) measured a viscosity approximately two orders of magnitude lower for dry limonene SOA, corresponding to a D_b of ∼10⁻¹² cm² s⁻¹ (Hinks et al. 2016). As these values are two orders of magnitude apart, we conducted simulations using both 10⁻¹⁴ and 10⁻¹² cm² s⁻¹. We assumed that limonene SOA is entirely nonvolatile, as Ye et al. (2018) observed very little uptake of limonene SOA into toluene SOA. The diffusivity of the toluene SOA is set to be D_b,TOL,RH10 = 3 × 10⁻¹⁷ cm² s⁻¹ and D_b,TOL,RH30 = 10⁻¹³ cm² s⁻¹ as discussed earlier.

A mass fraction of ∼0.3 is achieved in ∼2 h, as observed in Ye et al. (2016) for both D_b cases, but we are only able to reproduce the initial rapid uptake of toluene SOA with the D_b of 10⁻¹² cm² s⁻¹ (Figure 6b). The initial rapid uptake of toluene SOA compounds into limonene SOA is from the vapors present in equilibrium with toluene SOA in the main chamber where the toluene SOA was generated. Therefore, this initial uptake only depends on the D_b of the limonene SOA and indicates that there are likely little diffusion limitations in limonene SOA, even at low RH. Uptake slows down after these vapors are depleted, and gradual uptake continues over the next 1–2 h. While the simulation reproduces the initial rapid uptake of toluene SOA into limonene SOA in the first 5 min, the observations show slower uptake than the model predicts at later times. This may be due to an overestimation of the gas-phase concentration of toluene semi-volatile oxidation products, leading to a larger effect of the initial vapor uptake on the total uptake. This discrepancy may also be caused by the model simplification by treating only two volatility bins with an assumption of an equal fraction of each semi-volatile species. The uptake of the lower volatility species is slower because a larger fraction of its total concentration starts in the viscous toluene SOA. Therefore, a larger fraction of lower volatility species would slow down the overall uptake, potentially leading to better agreement at later times. Additionally, we have assumed both populations are monodispersed as no information on specific size distributions is available. As diffusion limitations would lead to equilibrium being achieved slower in larger particles, faster uptake of toluene SOA vapors into smaller particles and slower uptake into larger particles of the limonene SOA would be expected. Overall, these simulations indicate that the enhanced mixing of toluene SOA and limonene SOA, relative to the H-toluene SOA and D-toluene SOA mixing at low RH, is due to the lower viscosity of the limonene SOA. The rapid initial uptake is due to vapors present in equilibrium with the toluene SOA, which slows down once these vapors are depleted.

Toluene SOA +$\beta$-caryophyllene SOA

Ye et al. (2018) also investigated particle–particle mixing of toluene SOA and β-caryophyllene SOA. As shown in blue markers in Figure 6c, they observed only ∼10% of the toluene SOA is taken up by the β-caryophyllene SOA after 4 h of contact with an increase in RH from 7 to 50%. This extent of mixing is lower than the ∼30% uptake observed for mixing with D-toluene SOA and limonene SOA as above, implying that non-ideality plays a role in suppressing mixing. A recent investigation has indeed observed that mixtures of β-caryophyllene SOA and toluene SOA formed two separate phases, indicating a lack of miscibility between these two populations (Mahrt et al. 2022).

Maclean et al. (2021) determined that D_b of organic molecules in β-caryophyllene SOA to be ∼4.1 × 10⁻¹² cm² s⁻¹ at 48% RH and ∼1.8 × 10⁻¹⁴ cm² s⁻¹ at 0% RH based on viscosity measurements and subsequent application of the fractional Stokes–Einstein equation, which is consistent with another measurement using a fluorescent molecule by Evoy et al. (2021). These lead to characteristic bulk diffusion timescales of $\le$ 1 h for particles smaller than 400 nm, even at 0% RH (Seinfeld and Pandis 2016). Hence, bulk diffusion limitations would be expected to be weak in the β-caryophyllene SOA at 50% RH (Grayson et al. 2016).

In Figure 6c, we simulate the mixing of toluene SOA and β-caryophyllene SOA using D_b of β-caryophyllene SOA to be 10⁻¹⁴ cm² s⁻¹ for the first 2 h of the experiment at 7% RH and increase D_b to 10⁻¹² cm² s⁻¹ at 1.5 h when the RH is increased to ∼20% (Maclean et al. 2021). We consider non-ideality by applying the activity coefficient of SVOCs in β-caryophyllene SOA (γ) of 0.5, 1, 2, 3, and 4. We are able to reproduce the measurements with γ = 4. If we assume ideal mixing with γ = 1, a mass fraction of toluene SOA in β-caryophyllene SOA exceeds 0.1 within 2 h and reaches ∼0.3 within 4 h. To reproduce the measurements with ideal mixing by maintaining a mass fraction of 0.1 for over 4 h, it would require a very low D_b of 10⁻¹⁶ cm² s⁻¹ for the entire experiment, which is unrealistically too low based on measurements. Therefore, our results suggest that the low uptake of toluene SOA compounds into β-caryophyllene SOA is due to miscibility limitations with a high activity coefficient without significant bulk diffusion limitations as speculated in Ye et al. (2018).

Summary and implications

In this work, we investigate the particle–particle mixing timescale of a semi-volatile species in equilibrium with an aerosol population (population 1) when another population of aerosol (population 2) is introduced to that system. We find that particle–particle mixing can be prolonged when even one population is highly viscous and at extremes of volatility due to the slow evaporation of low-volatility species from population 1 and slow condensation of high-volatility species to population 2. Additionally, the mixing timescale is prolonged as the miscibility between the semi-volatile species and population 2 is increased with low activity coefficients. A notable exception to this is when population 2 is ultraviscous, and the semi-volatile species is favorably miscible in population 2 (γ₂ < 1). In this case, the semi-volatile species can build up near the particle surface, allowing the equilibrium concentration in population 2 to be achieved faster.

Our model simulations suggest that mixing of SVOCs from a liquid population to a liquid or semi-solid population would proceed within 1 h for any miscibility investigated in this study. We also find that total organic aerosol concentration, which is known to affect equilibration timescales, does not have a large effect on the particle–particle mixing timescale. Thus, a typical assumption of equilibrium partitioning of SVOCs in aerosol models is likely to hold within 1-h timestep when the mixing populations have D_b < 10⁻¹⁵ cm² s⁻¹. Note that under some conditions including low C⁰ and low γ_2, equilibrium can be prolonged from hours to months and in this case, particle–particle mixing would need to be treated kinetically.

We have applied this model to reconcile apparent discrepancies between experimental measurements of SOA viscosity and direct measurements of particle–particle mixing. We show that mixing between H-toluene SOA and D-toluene SOA is prolonged at low RH because of slow bulk diffusivity, but mixing is achieved faster at moderate and high RH with little kinetic limitations. Additionally, the initial rapid mixing of toluene SOA into limonene SOA at low RH can be explained by gas–particle partitioning, while continuous slow mixing can be explained by slow evaporation from toluene SOA because of diffusion limitations. For mixing toluene SOA and $\beta$-caryophyllene SOA, limited mixing is explained by limited miscibility between these two populations in the absence of kinetic limitations. Through these simulations, we show that these direct experimental observations of particle–particle mixing are consistent with measurements of SOA viscosity. The apparent discrepancy of rapid mixing between toluene SOA and limonene SOA but slow mixing between H-toluene and D-toluene SOA at low RH can be explained by the lower viscosity and little bulk diffusion limitations in limonene SOA. The apparent discrepancy of slow mixing between toluene and $\beta$-caryophyllene SOA but rapid mixing between H-toluene and D-toluene SOA at high RH can be explained by a high activity coefficient of toluene SOA in $\beta$-caryophyllene SOA.

Mixing of different types of organic aerosol from different sources occurs frequently in the atmosphere. For example, biomass burning and wildfires can generate large amounts of organic aerosol (BBOA), diluted by mixing with background aerosols with different chemical compositions and phase states (Liu et al. 2017). Measurements of BBOA viscosity indicate that the viscosity can vary over orders of magnitude with water content. Models that have assumed completely internal mixtures between biomass burning and non-biomass burning components of aerosol have dramatically overestimated the warming effect of black carbon from biomass burning, and assuming external mixtures brings model predictions of single-scattering albedo into better agreement with observations (Cappa et al. 2012; Brown et al. 2021; Jahl et al. 2021). Thus, to accurately represent the climatic effect of biomass burning aerosol, it is critical to understand the extent to which it mixes with background aerosol. Outdoor-to-indoor aerosol mixing can also have implications for the evolution of particle size distributions and multiphase processes occurring between species exchanged between the two populations. Indoor and outdoor aerosols are formed and evolved under very different environmental conditions and can, therefore, exist in different phase states when they come into contact with each other. Therefore, it is also critical to evaluate the timescales for mixing indoor and outdoor aerosol populations, which we intend to investigate in future studies.

Acknowledgments

The authors would like to acknowledge Dr. Qing Ye for helpful discussions and comments on the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was funded by U.S. Department of Energy under Grant DE-SC0022139 and Alfred P. Sloan Foundation under grant G-2020-13912. NMD acknowledges the funding by the National Science Foundation (AGS-2132089).