Estimation of secondary organic aerosol viscosity from explicit modeling of gas-phase oxidation of isoprene and α-pinene

Secondary organic aerosols (SOA) are major components of atmospheric fine particulate matter, affecting climate and air quality. Mounting evidence exists that SOA can adopt glassy and viscous semisolid states, impacting formation and partitioning of SOA. In this study, we apply the GECKO-A (Generator of Explicit Chemistry and Kinetics of Organics in the Atmosphere) model to conduct explicit chemical modeling of isoprene photooxidation and α-pinene ozonolysis and their subsequent SOA formation. The detailed gas-phase chemical schemes from GECKO-A are implemented into a box model and coupled to our recently developed glass transition temperature parameterizations, allowing us to predict SOA viscosity. The effects of chemical composition, relative humidity, mass loadings and mass accommodation on particle viscosity are investigated in comparison with measurements of SOA viscosity. The simulated viscosity of isoprene SOA agrees well with viscosity measurements as a function of relative humidity, while the model underestimates viscosity of α-pinene SOA by a few orders of magnitude. This difference may be due to missing processes in the model, including autoxidation and particle-phase reactions, leading to the formation of high-molar-mass compounds that would increase particle viscosity. Additional simulations imply that kinetic limitations of bulk diffusion and reduction in mass accommodation coefficient may play a role in enhancing particle viscosity by suppressing condensation of semi-volatile compounds. The developed model is a useful tool for analysis and investigation of the interplay among gas-phase reactions, particle chemical composition and SOA phase state.

1 Introduction

Secondary organic aerosols (SOA) are ubiquitous in the atmosphere and represent a major component of fine particulate matter, affecting air quality, climate and public health (Jimenez et al., 2009; Pöschl and Shiraiwa, 2015). Due to their complexity, SOA represent a large source of uncertainty in current understanding of global climate change and air pollution (Tsigaridis et al., 2014; Ciarelli et al., 2019). Development of SOA models represents one of the most challenging and demanding problems in atmospheric chemistry (Shrivastava et al., 2017). Formation of SOA is initiated by gas-phase oxidation of biogenic and anthropogenic volatile organic compounds (VOCs) (Kroll and Seinfeld, 2008). Typically, multigenerational oxidation of VOCs in the gas phase leads to the formation of a myriad of semi-volatile and low-volatility compounds which can condense on pre-existing particles (Ziemann and Atkinson, 2012; Nozière et al., 2015) or contribute to nucleation and new particle formation (Tröstl et al., 2016). As gas-phase oxidation is a driving step of SOA formation, there is a strong need for a computational tool that can generate exhaustive gas-phase chemical mechanisms. The GECKO-A (Generator of Explicit Chemistry and Kinetics of Organics in the Atmosphere) model is to date the most extensive generator of gas-phase chemical schemes; based on established reaction pathways and structure–activity relationships, it automatically generates detailed gas-phase mechanisms involving thousands to millions of oxidation products from a given VOC precursor (Aumont et al., 2005, 2012; Lee-Taylor et al., 2011). Most aerosol models treat SOA particles as homogeneous and well-mixed liquids with a dynamic viscosity (η) below 10² Pa s. Recent studies provide accumulating evidence that SOA can adopt glassy solid (η≥10¹² Pa s) or amorphous semi-solid phase states (10${}^{\mathrm{2}}\le \mathit{\eta }<{\mathrm{10}}^{\mathrm{12}}$ Pa s) depending on chemical composition, relative humidity (RH) and temperature (Koop et al., 2011; Reid et al., 2018; Virtanen et al., 2010). The particle phase state is impacted by SOA chemical composition, hygroscopicity and water content, as water can act as a plasticizer lowering SOA viscosity (Mikhailov et al., 2009). It has been observed that ambient SOA particles can bounce off from an impactor stage depending on ambient conditions, indicating non-liquid states of organic particles (Virtanen et al., 2011; Bateman et al., 2017; Slade et al., 2019). Measurements of the viscosity of SOA bulk materials derived from the oxidation of limonene (Hinks et al., 2016), toluene (Song et al., 2016), α-pinene (Renbaum-Wolff et al., 2013; Kidd et al., 2014; Bateman et al., 2015; Zhang et al., 2015; Grayson et al., 2016; Hosny et al., 2016) and isoprene (Song et al., 2015) have confirmed that viscosity of SOA particles can vary depending on temperature and RH (Petters et al., 2019). In addition, previous studies have shown that SOA chemical composition can be affected by experimental conditions upon SOA formation such as RH and particle mass concentrations (Kidd et al., 2014; Hinks et al., 2018; Grayson et al., 2016; Jain et al., 2018). These results imply that SOA viscosity and phase state are dynamic properties in response to chemical processing of organic aerosols and variations in RH and temperature in the atmosphere.

In liquid particles, molecules diffuse quickly in the particle bulk, leading to rapid establishment of gas–particle equilibrium. In glassy and viscous particles, kinetic limitations of bulk diffusion can significantly retard gas–particle partitioning, prolonging equilibration timescales (Li and Shiraiwa, 2019; Mai et al., 2015; Vaden et al., 2011), while relatively fast particle–particle mixing was observed at intermediate and high RH (Ye et al., 2016, 2018). Model simulations suggest that an assumption of instantaneous equilibrium partitioning can result in substantial overestimation of particle mass concentration and underestimation of gas-phase mass concentration (Shiraiwa and Seinfeld, 2012). A proper consideration of particle phase state is essential for simulating particle size distribution dynamics in SOA growth (Shiraiwa et al., 2013; Zaveri et al., 2018, 2020). Heterogeneous and multiphase reactions (Kuwata and Martin, 2012; Shiraiwa et al., 2011; Davies and Wilson, 2015; Zhang et al., 2018; Marshall et al., 2018; Zhou et al., 2019) as well as activation to cloud droplets and ice crystals (Slade et al., 2017; Knopf et al., 2018) can also be impacted by glassy and amorphous semisolid states. Thus, it is important to consider the effects of particle phase state on SOA processes for accurate representation of the fate of SOA particles in aerosol models.

Recently, we have developed parameterizations to predict glass transition temperature (T_g) of organic compounds based on molar mass and atomic O:C ratio (Shiraiwa et al., 2017), elemental composition (i.e., number of carbon, hydrogen, oxygen, nitrogen and sulfur atoms) (DeRieux et al., 2018), and volatility (Li et al., 2020). T_g is a characteristic temperature at which a phase transition from a glassy solid state to an amorphous semi-solid state occurs (Koop et al., 2011). Applying the Gordon–Taylor mixing rule for mixtures of SOA multiple components and water, the viscosity of SOA can be estimated by the T_g-scaled Arrhenius plot of viscosity (DeRieux et al., 2018). This method has been successfully applied to estimate particle viscosity based on elemental composition obtained from high-resolution mass spectrometry (Schum et al., 2018; Ditto et al., 2019; Song et al., 2019). Gervasi et al. (2020) simulated viscosity of aqueous SOA surrogate mixtures derived from oxidation of α-pinene, isoprene and toluene using the Aerosol Inorganic–Organic Mixtures Functional groups Activity Coefficients viscosity framework (AIOMFAC-VISC) model combined with the T_g parameterization of DeRieux et al. (2018), demonstrating the capability and flexibility of the group-contribution method in predicting the viscosity for organic mixtures of varying degrees of complexity.

There is a strong need for the development of a modeling tool for further elucidation of the impacts of various aerosol properties and processes on particle viscosity. A better understanding of the interplay among gas-phase reactions, chemical composition and viscosity of SOA is essential for further development of model representation of the evolution of SOA in the atmosphere. In this study, the T_g and viscosity prediction methods are implemented into the GECKO-A and box model. Model simulations with explicit chemical mechanisms are conducted to reproduce chamber experimental conditions for SOA generation and subsequent viscosity measurements. The effects of chemical composition and mass accommodation on SOA partitioning and viscosity are also investigated. The objective of this work is to develop a useful tool for analysis and investigation of SOA chemical composition and phase state and to expand our fundamental understanding of such properties.

2 Methods

2.1 GECKO-A

GECKO-A generates detailed gas-phase chemical oxidation schemes and the associated gas–particle mass transfers for SOA formation. It produces explicit highly detailed chemical mechanisms starting from experimental data and structure–activity relationships (SARs) (Aumont et al., 2005). Implemented into a box model, these explicit chemical mechanisms simulate the oxidation of parent VOC precursors to oxidation products, their subsequent gas-phase chemistry, as well as partitioning into the particle phase based on their vapor pressures (Camredon et al., 2007). GECKO-A is frequently updated to include the newly discovered chemical processes in the mechanism generator: recently, the mechanism generator has been extended to treat the chemistry of aromatic organic compounds based on the latest SAR developments (Jenkin et al., 2018a, b, 2019, 2020). Other recent updates include the chemistry of cyclic structures (Valorso et al., 2011), the gas-aqueous phase partitioning (Mouchel-Vallon et al., 2013) and the gas-to-chamber-wall partitioning (La et al., 2016). The chemical mechanism does not account for gas-phase chemistry of species with vapor pressure below 10⁻¹³ atm, which are assumed to be of enough low volatility to partition to the condensed phase. In our simulations, species vapor pressures are estimated via the approach by Nannoolal et al. (2008).

The mechanism generator follows a reduction protocol optimizing the number of generated reactions and chemical species (Valorso et al., 2011). This protocol has been established to reduce the computational cost of box model simulations. In the gas phase, isomer substitution is not allowed for position isomers when the overall yield for the product is > 5 %; otherwise only the highest-yield positional isomer is considered. For each reaction, pathways with branching ratios < 5 % are not accounted for, and the main reaction branches are scaled proportionally (Valorso et al., 2011). The chemical mechanism accounts for chemistry of peroxy and alkoxy radicals generated during precursors oxidation. SOA chemistry is a multi-generational process, and the generated chemical schemes follow the chemistry of reaction products up to the fifth generation in our simulations. The GECKO-A version used in this study has been recently enriched with SAR estimations of alkoxy radical decomposition and H-migration reaction rates (Vereecken and Peeters, 2009; La et al., 2016). Rate coefficients and branching ratios for gas-phase reactions of OH with aliphatic compounds and peroxy radicals are estimated based on recent SAR investigations (Jenkin et al., 2018a, b). Our simulations do not include autoxidation and condensed-phase chemistry, which represent limitations of GECKO-A modeling.

α-pinene and isoprene oxidation mechanisms have been investigated extensively both in model and experimental studies. Their oxidation mechanisms are well known, providing good estimations of resulting SOA composition. The SARs for both α-pinene and isoprene oxidation are implemented into GECKO-A with state-of-the-art protocols estimating radical reactions and respective rate constants. The base mechanism for α-pinene photooxidation and ozonolysis includes a detailed description of reactions with OH radicals (McVay et al., 2016) and branch reactions identified via quantum and theoretical chemical calculations (Peeters et al., 2001; Vereecken et al., 2007). Stabilized Criegee intermediate (SCI) chemistry is not treated explicitly, and stable products are directly assigned via a rule-base method. This scheme accounts for a direct reaction route to account for pinonaldehyde formation during the reaction of SCIs with water (McVay et al., 2016).

α-pinene GECKO-A mechanisms were evaluated in previous studies by comparisons with chamber experiments of SOA formation from photo-oxidation (Valorso et al., 2011; McVay et al., 2016) and ozonolysis (Denjean et al., 2015). The model captures the qualitative features of SOA formation with the variation of NO_x levels during photolysis experiments and with temperature for ozonolysis experiments. The measured SOA mass was overestimated by the model in Valorso et al. (2011) and Denjean et al. (2015), which could be in part explained by the loss of low-volatility organic compounds onto the chamber wall (McVay et al., 2016; La et al., 2016). In addition, McVay et al. (2016) highlighted a possible overestimation of the simulated contribution of later-generation oxidation products to SOA mass at high OH levels and an underestimation of the SOA growth at low OH levels that could be due to a lack of autoxidation processes in the mechanisms generated with GECKO-A. The influence of these limitations on the results simulated here is discussed along the paper. The GECKO-A isoprene oxidation scheme follows the main protocol developed by Aumont et al. (2005), and branching reactions for isoprene + OH are also considered (Paulot et al., 2009). Isoprene mechanisms generated with GECKO-A are for the first time indirectly evaluated in this study. The chemical schemes for α-pinene and isoprene oxidation are composed of 84 000 and 1900 species and of ∼ 700 000 and ∼ 17 000 reactions, respectively.

2.2 Box model and simulations setup

Each chemical mechanism is coupled to a box model representing SOA formation from α-pinene or isoprene oxidation in specific chamber experiments. The time evolution of species concentration is computed through a two-step method that explicitly solves stiff ordinary differential equations (Verwer, 1994, 1996). The experimental conditions are summarized in Table 1. All laboratory experiments were performed in absence of seed particles, but the box model used in our simulations does not treat nucleation. Thus, in our simulations the first steps of nucleation are approximated by the addition of seed particles with a particle radius of 5 nm and a concentration of 10⁴ particles cm⁻³ (McVay et al., 2016). The particle number concentration is assumed to remain constant during simulations (coagulation is not treated), while the particle radius evolves following the partitioning of organics. The density of SOA particles is assumed to be 1.2 g cm⁻³ (Kuwata et al., 2012). All simulations are conducted with a NO_x mixing ratio of 3 ppt to reflect low-NO_x conditions of chamber and flow tube experiments. For all simulations we consider reversible gas–wall partitioning of organic species with the vapor wall loss rate of of 10⁻³ s⁻¹ based on experimental observations (Lim and Ziemann, 2009; McVay et al., 2016). Photolysis rates are computed using the cross sections and the quantum yields from Aumont et al. (2005). Photolysis conditions are set to average daylight conditions of a mid-latitude spring day for isoprene photooxidation simulations and to zero for α-pinene dark ozonolysis experiments, coherent with the experimental conditions.

Table 1 Experimental conditions simulated in the box model to explore α-pinene and isoprene SOA viscosity, as well as simulated glass transition temperatures of dry SOA (T_g,org) and the O:C ratio.

^* Poke flow bead mobility (PF-BM), differential mobility analyzer (DMA), particle impactor (PI), and rebound impactor (R).

For continuous flow reactor experiments, the concentrations of VOC precursors and oxidants were stable by continuous inflow and outflow of gases from the tube reactor; in the box model simulations, precursor concentrations are fixed at steady-state values as measured in corresponding experiments. In chamber batch experiments, the reactants were injected once at the beginning of the experiment and the concentrations evolved; in box model simulations the VOC precursor and ozone concentrations are initially set to the respective experimental values, and species concentrations evolve due to gas oxidation, gas–particle partitioning and wall deposition. SOA particles were formed under dry conditions and then exposed to water vapor at different RH for viscosity measurements; the same procedure is applied in our simulations by forming SOA at RH = 0.5 %, and then viscosity estimations are conducted by considering hygroscopic growth at elevated RH. For each set of simulations, the measured RH and temperature were used to constrain the simulation conditions. Note that partitioning of water is not considered in the box model as particle water uptake should be very minor under the very low RH values at which SOA particles were formed.

The box model treats the mass transfer of gaseous organic species to particles and to chamber walls. Gas–particle partitioning is assumed to follow Raoult's law at equilibrium (La et al., 2016) (i.e., depending on the saturation vapor pressure of the organic compounds assuming ideal mixing). Gas-to-particle partitioning is described in the box model using the gas–particle mass transfer coefficient of a compound (k_gp in s⁻¹) with the Fuchs–Sutugin approach in the transition regime as (Seinfeld and Pandis, 2016)

$$\begin{array}{}\text{(1)}& {\displaystyle }{k}_{\mathrm{gp}}=\mathrm{4}\mathit{\pi }{D}_{\mathrm{g}}{r}_{\mathrm{P}}{N}_{\mathrm{P}}\mathit{\beta },\text{(2)}& {\displaystyle }\mathit{\beta }={\displaystyle \frac{\mathrm{0.75}\mathit{\alpha }(\mathrm{1}+Kn)}{K{n}^{\mathrm{2}}+Kn+\mathrm{0.283}Kn\mathit{\alpha }+\mathrm{0.75}\mathit{\alpha }}},\end{array}$$

where D_g (cm² s⁻¹) is the gas diffusivity, r_P(cm) is the particle radius, N_p (cm⁻³) is the particle number concentration, Kn is the Knudsen number and α is the mass accommodation coefficient. α, also termed as the bulk accommodation coefficient, represents the probability for a gas molecule colliding with surface to enter the bulk of the particle. Based on recent experiments and molecular dynamics simulations, α is assumed to be unity for the base case scenario of our simulations (Liu et al., 2019; Julin et al., 2014). This approach does not account for potential kinetic limitations caused by bulk diffusion in viscous particles.

A recent study has introduced an effective mass accommodation coefficient α_eff, which effectively considers the kinetics of accommodation at the surface, transfer across the gas–particle interface and further transport into the particle bulk (Shiraiwa and Pöschl, 2021):

$$\begin{array}{}\text{(3)}& {\mathit{\alpha }}_{\mathrm{eff}}={\mathit{\alpha }}_{\mathrm{s}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{{\mathit{\alpha }}_{\mathrm{s}}\mathit{\omega }{C}^{\mathrm{0}}}{\mathrm{4}{D}_b{\mathit{\rho }}_{\mathrm{p}}}\cdot \frac{r_{\mathrm{p}}}{\mathrm{5}}\times {\mathrm{10}}^{-\mathrm{12}}\frac{\mathrm{g}}{\mathrm{\mu }\mathrm{g}}\frac{{\mathrm{cm}}^{-\mathrm{3}}}{{\mathrm{m}}^{-\mathrm{3}}}}},\end{array}$$

where α_s is the surface accommodation coefficient (i.e., the probability for a gas molecule colliding with surface to adsorb on the particle surface), ω (cm s⁻¹) is the mean thermal velocity of the organic compound in the gas phase, D_b (cm² s⁻¹) is diffusivity in the condensed phase, ρ_p (g cm⁻³) is the particle density and C⁰ (µg m⁻³) is the pure compound saturation mass concentration. Application of α_eff in the Fuchs–Sutugin approach (e.g., α=α_eff) yields SOA partitioning with effective consideration of kinetic limitations induced by slow bulk diffusion in viscous particles. This approach can yield consistent results with a detailed kinetic multilayer model (KM-GAP, Shiraiwa et al., 2012) and two-film model solutions (Zaveri et al., 2014). We implement this method to explore the effects of mass accommodation and kinetic partitioning on predicted viscosity by comparing with the base case scenario (e.g., α fixed to 1).

2.3 Viscosity prediction implementation

The glass transition temperature of an organic compound i (T_g,i) can be estimated with the following recently developed parameterization based on its elemental composition (DeRieux et al., 2018):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{T}_{\mathrm{g},i}=& (n_{\mathrm{C}}^{\mathrm{0}}+\ln (n_{\mathrm{C}}\left)\right)b_{\mathrm{C}}+\ln \left(n_{\mathrm{H}}\right)b_{\mathrm{H}}+\ln \left(n_{\mathrm{C}}\right)\ln \left(n_{\mathrm{H}}\right)b_{\mathrm{CH}}\\ & +\ln \left(n_{\mathrm{O}}\right)b_{\mathrm{O}}+\ln \left(n_{\mathrm{C}}\right)\ln \left(n_{\mathrm{O}}\right)b_{\mathrm{CO}},\end{array}\end{array}$$

where the coefficient values for [$n_{\mathrm{C}}^{\mathrm{0}}$, b_C, b_H, b_CH, b_O, and b_CO] are [1.96, 61.99, −113.33, 28.74, 0, 0] for CH compounds and [12.13, 10.95, −41.82, 21.61, 118.96, −24.38] for CHO compounds. We have also included the recently developed parameterization for CHON compounds by Li et al. (2020):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{T}_{\mathrm{g},i}=& (n_{\mathrm{C}}^{\mathrm{0}}+\ln (n_{\mathrm{C}}\left)\right)b_{\mathrm{C}}+\ln \left(n_{\mathrm{O}}\right)b_{\mathrm{O}}+\ln \left(n_{\mathrm{N}}\right)b_{\mathrm{N}}\\ & +\ln \left(n_{\mathrm{C}}\right)\ln \left(n_{\mathrm{O}}\right)b_{\mathrm{CO}}+\ln \left(n_{\mathrm{C}}\right)\ln \left(n_{\mathrm{N}}\right)b_{\mathrm{CN}}\\ & +\ln \left(n_{\mathrm{O}}\right)\ln \left(n_{\mathrm{N}}\right)b_{\mathrm{ON}},\end{array}\end{array}$$

where the coefficient values for [$n_{\mathrm{C}}^{\mathrm{0}}$, b_C, b_O, b_N,b_CO, b_CN, and b_ON] are [5.34, 31.53, −7.06, 134.96, 6.54, −34.36, and −15.35]. These parameterizations can predict T_g,i with an uncertainty of about ±30 K; note that, however, in multicomponent SOA mixtures this uncertainty would be much smaller for ideal mixing conditions (Shiraiwa et al., 2017; DeRieux et al., 2018; Li et al., 2020). The T_g of a mixture of organic compounds (e.g. dry SOA) can be estimated using the Gordon–Taylor equation with a Gordon–Taylor constant (k_GT) of 1: $T_{\mathrm{g},\mathrm{org}}={\sum }_i{\mathit{\omega }}_i{T}_{\mathrm{g},i}$, where ω_i is the mass fraction of compound i (Dette et al., 2014).

The Gordon–Taylor equation can also be used to estimate the glass transition temperature of mixtures of organic compounds and water (T_g(w_org)):

$$\begin{array}{}\text{(6)}& T_{\mathrm{g}}\left(w_{\mathrm{org}}\right)={\displaystyle \frac{(\mathrm{1}-w_{\mathrm{org}})T_{\mathrm{g},\mathrm{w}}+\frac{\mathrm{1}}{k_{\mathrm{GT}}}{w}_{\mathrm{org}}{T}_{\mathrm{g},\mathrm{org}}}{(\mathrm{1}-w_{\mathrm{org}})+\frac{\mathrm{1}}{k_{\mathrm{GT}}}{w}_{\mathrm{org}}}},\end{array}$$

where w_org is the mass fraction of organics, T_g,w is the glass transition temperature of pure water (136 K), and k_GT is the Gordon–Taylor constant which is assumed to be 2.5 (Koop et al., 2011; Zobrist et al., 2008). The w_org is calculated using the mass concentration of water (mH₂O) and of the organics in SOA (mSOA) as w_org=mSOA$/(m$SOA +mH₂O). mH₂O can be calculated using the Kohler theory with effective hygroscopicity parameter (κ) (Petters and Kreidenweis, 2007; Prenni et al., 2007).

Once T_g(w_org) has been computed, the viscosity (η) can be derived with the following equation based on the Vogel–Tammann–Fulcher approach (DeRieux et al., 2018):

$$\begin{array}{}\text{(7)}& {\displaystyle }\log \mathit{\eta }=-\mathrm{5}+\mathrm{0.434}{\displaystyle \frac{T_{\mathrm{0}}D}{T-T_{\mathrm{0}}}},\text{(8)}& {\displaystyle }{T}_{\mathrm{0}}={\displaystyle \frac{\mathrm{39.17}{T}_{\mathrm{g}}\left(w_{\mathrm{org}}\right)}{D+\mathrm{39.17}}},\end{array}$$

where T₀ is the Vogel temperature and D is the fragility parameter. The D value is assumed to be 10 based on previous studies (Shiraiwa et al., 2017; DeRieux et al., 2018).

Viscosity can then be converted into bulk diffusivity using the Stokes–Einstein equation:

$$\begin{array}{}\text{(9)}& D_b={\displaystyle \frac{kT}{\mathrm{6}\mathit{\pi }a\mathit{\eta }}},\end{array}$$

where k is the Boltzmann constant (1.38 × 10⁻²³ J K⁻¹), T is the temperature (K) and a is the effective molecular radius (m). This relation, however, may not be accurate for very high viscosities and needs to be corrected with the fractional Stokes–Einstein (FSE) relation (Evoy et al., 2019):

$$\begin{array}{}\text{(10)}& D_b=D_{b,\mathrm{c}}{\left({\displaystyle \frac{{\mathit{\eta }}_{\mathrm{c}}}{\mathit{\eta }}}\right)}^{\mathit{\xi }},\end{array}$$

where ξ is an empirical fit parameter with the value of 0.93 (Evoy et al., 2019), η_c is the crossover viscosity with the value of 10⁻³ Pa s and D_b,c is the crossover diffusion coefficient at ${\mathit{\eta }}_{\mathrm{c}}={\mathrm{10}}^{-\mathrm{3}}$ Pa s as determined by the Stokes–Einstein equation (Eq. 9).

For validation and applicability of T_g parameterizations and the viscosity prediction method, they have been applied to high-resolution mass spectrometry data of toluene SOA and biomass burning aerosols (DeRieux et al., 2018), SOA generated by diesel fuels (Song et al., 2019), β-caryophyllene SOA (Maclean et al., 2021), and surrogate VOC mixtures by healthy and stressed plants (Smith et al., 2021), agreeing well with viscosity measurements.

3 Results and discussion

3.1 Viscosity predictions of isoprene and α-pinene SOA

Figure 1 shows the measured and simulated viscosities of isoprene SOA, showing a very good agreement with experimental observations. SOA generated by isoprene ozonolysis adopt a semi-solid state for RH < 40 % and a liquid state for RH > 60 %. The model predicted T_g,_org is 260 K, which is in agreement with a previous estimation of 255 K (Berkemeier et al., 2014). All simulated viscosities with κ = 0.1 fall within uncertainties of experimental measurements except at 70 % RH, which is slightly below the lower bound of the error bar of experimental measurements. The simulated O:C ratio is 1.0, which is within the range of their estimation of 0.64–1.1 based on previous chamber experiments (Song et al., 2015). To examine the effects of κ on the calculated viscosity of isoprene SOA, sensitivity studies are conducted by varying κ within the range of 0.05–0.15, showing that lower and higher κ would lead to higher and lower viscosity, respectively.

Figure 1 Viscosity of isoprene SOA as a function of RH. Circle markers with error bars represent measurements and associated uncertainties by Song et al. (2015), and square markers represent simulations with the GECKO-A box model with different κ values.

Isoprene photooxidation in chamber experiments was conducted with very high precursor concentrations and a short reaction time such that a traditional gas-phase oxidation scheme dominates, which can be captured well in GECKO-A (Aumont et al., 2005). Note that the formation of isoprene epoxydiols (IEPOX) (Paulot et al., 2009; Bates et al., 2014) and subsequent uptake of IEPOX into acidic particles followed by a series of multiphase reactions including the formation of oligomers and organosulfates (Riva et al., 2019) are found to be important in the ambient atmosphere (Wennberg et al., 2018). The uptake of IEPOX can be limited by bulk diffusion in a viscous matrix (Zhang et al., 2018), and IEPOX-SOA is reported to have high viscosity (Riva et al., 2019). For viscosity simulations of atmospherically relevant IEPOX-derived SOA, IEPOX multiphase processes would need to be treated in future studies.

Figure 2 shows viscosity of α-pinene SOA measured in a number of studies and the corresponding simulated results by the GECKO-A box model for each experimental condition. Based on previous laboratory measurements κ was assumed to be 0.1 (Lambe et al., 2011; Pajunoja et al., 2014; Petters et al., 2019). Overall, the model simulations can reproduce the RH-dependent viscosity of α-pinene SOA, showing that α-pinene SOA is semi-solid between 0 %–65 % RH and liquid for RH > 65 % at the room temperature. Simulated viscosity values of α-pinene SOA are by a few orders of magnitude higher than the ones simulated for isoprene SOA, which is consistent with experimental observations. Simulation results for Renbaum-Wolff's and Grayson's experiments fall within the relatively large uncertainties of experimental measurements for the 40 %–60 % RH range. In simulations for Renbaum-Wolff's experiments at RH > 60 %, the predicted viscosity values are about 1 order of magnitude lower than the measured values. Larger deviations are observed for RH < 40 %, where simulated viscosities are up to 4 orders of magnitude lower than experimental measurements. The simulated viscosities are about 1 order of magnitude higher than Grayson's measurements for 20 < RH < 40 %. The simulated viscosity for Bateman's experiment agrees well with the observed experimental value within uncertainties.

Figure 2 Comparison of measured (solid markers) and simulated (open markers) viscosity of α-pinene SOA as a function of RH. For each set of laboratory experiments the measured viscosities are reported with uncertainties. The solid and dashed lines represent the polynomial regression for measurements and model simulations, respectively, to guide eyes.

Kidd et al. (2014) provided only one data point at 85 % RH, which also represents inferred viscosity based on their impactor measurement; the predicted viscosity is 3 orders of magnitude below the inferred value but only 1 order of magnitude below its lower bound. For simulations of Zhang et al. (2015) our predictions are about 3 to 4 orders of magnitude lower than measurements.

The simulated T_g,org values are within 255–285 K (see Table 1), which is in good agreement with a previous estimation of 278 K (Berkemeier et al., 2014) and a measurement of 272 K (Petters et al., 2019). To assess overall performance of our viscosity predictions, measured and simulated viscosities were fitted respectively with a second-order polynomial regression. While the RH dependence of viscosity is very well captured by our simulations, the simulated values underestimate viscosity by about 2 or 3 orders of magnitude on average. The variance of model simulations and experimental measurements calculated for the logarithm of the viscosity are very similar (R²=0.85), indicating that the GECKO-A box model can reproduce the variability observed in experimental measurements with different experimental conditions.

Grayson et al. (2016) investigated the impact of SOA mass loadings on viscosity, observing lower viscosity with higher SOA mass concentrations in the range of 10²–10⁴ µg m⁻³. We conducted four different simulations where precursor concentrations were varied in the range of the experimental values to match the final SOA mass loadings. As shown in Fig. 3, our simulated viscosities are within the experimental uncertainties, capturing this dependence very well. At higher mass loadings, compounds with relatively high volatility and low T_g can condense on particles, leading to a reduction of SOA viscosity; for lower mass loadings, condensation of lower-volatility compounds with higher T_g would be dominant, resulting in higher viscosity.

Figure 3 Impacts of SOA mass loadings on viscosity of α-pinene SOA. The red and green lines represent the upper and lower bounds of measured viscosity in Grayson et al. (2016). Orange markers represent simulated viscosities by the GECKO-A box model.

3.2 Chemical composition and functional group distributions of SOA

The GECKO-A box model tracks the concentrations of species both in the gas and particle phases, while retaining information on molecular properties including molar mass, O:C ratio, vapor pressure, and volatility or pure compound saturation mass concentration (C₀). Figure 4 shows most abundant particle-phase 500 compounds in simulated isoprene and α-pinene SOA products in the two-dimensional volatility basis set framework of O:C ratio vs. log C₀ (Donahue et al., 2011). The markers are color-coded with T_g, and the marker size is scaled with particle-phase concentration in each simulation. Compounds with lower volatility tend to have higher O:C ratio and higher T_g, in good agreement with recent experimental and model studies (Zhang et al., 2019; Li et al., 2020). Isoprene oxidation products are found to have higher O:C ratio, while α-pinene oxidation products have lower C₀ and higher T_g. These results are in line with higher viscosity of α-pinene SOA compared to isoprene SOA, as measured and modeled in Figs. 1 and 2.

Figure 4 Chemical composition of SOA derived from oxidation of (a) isoprene and (b) α-pinene of all simulated experimental conditions reported in Table 1 plotted in the 2D-VBS framework of volatility and O:C ratio. Markers represent 500 most abundant particle-phase compounds from each simulation. The markers are color-coded with T_g. The marker size is scaled with particle-phase concentration in each simulation.

None of these experimental studies measured average O:C ratio of SOA particles, but some of them reported estimated O:C ratios based on previous experiments with similar experimental conditions. Song et al. (2015) estimated the O:C ratio of isoprene SOA to be 0.64–1.1, which is consistent with our simulated O:C ratio of 1.0. Renbaum-Wolff et al. (2013) estimated the O:C ratio of α-pinene SOA to be 0.3–0.4, which is slightly lower than our simulated value of 0.49. Valorso et al. (2011) reported that the GECKO-A box model tends to overestimate the O:C ratio of SOA generated by α-pinene photooxidation. Denjean et al. (2015) showed that the O:C ratio of SOA simulated with GECKO-A from α-pinene ozonolysis was coherent with the one simulated with the Master Chemical Mechanism and within the uncertainty range of the experimental value. A recent study by Gervasi et al. (2020), which selected 14–21 representative reaction products based on the Master Chemical Mechanism to simulate viscosity of the same measurement dataset, estimated the O:C ratio to be 1.1 for isoprene SOA and 0.51 for α-pinene SOA. Overall, the O:C ratios simulated by the GECKO-A box model are reasonable and in line with available measurements and modeling studies, while further studies are warranted for simultaneous measurements of O:C ratio and viscosity along with model applications.

Our methods to estimate T_g of organic compounds and SOA viscosity are based on elemental composition without accounting for molecular structure and specific intramolecular interactions. Some organic compounds with reactive functional groups may undergo particle-phase reactions, which are not treated in our simulations. To explore these aspects, we investigate the functional groups distributions computed by GECKO-A and box modeling. The simulated functionality group distributions of particle-phase SOA compounds are shown in Fig. 5. Isoprene SOA is characterized by high concentrations of alcohols (–ROH) and hydroperoxides (–ROOH), followed by lower concentrations of ketones (–RC(O)R) and aldehydes (–RCHO). Our results are consistent with experimental measurements showing that polyols and organic peroxides are the primary species formed upon isoprene photooxidation under low-NO_x conditions (Surratt et al., 2006). Most experiments for α-pinene led to similar chemical composition with high concentrations of –RC(O)R and –ROOH, followed by slightly lower concentrations of –RCHO, –ROH and carboxylic acids (–RC(O)OH). The simulations for Grayson's experiments are characterized by the highest concentrations of –ROOH and –ROH and by lower fractions of –RCHO. The simulations for Kidd's experiments show the second highest abundance of –ROH and noticeable amounts of –RC(O)OR. The simulation for Bateman's experiment is characterized by high –RC(O)R abundance. A previous study suggested that the hydrogen-bond formation among organic compounds in the SOA organic matrix may influence the viscosity of α-pinene SOA significantly (Kidd et al., 2014). It has been further shown that alcohol and carboxylic acid groups can increase viscosity due to the formation of hydrogen bonding and ionic interactions if carboxylic acids are dissociated (Grayson et al., 2017; Rothfuss and Petters, 2017). These effects are not explicitly considered in our viscosity estimation method, which may be one of the reasons for the discrepancies between measurements and simulations of SOA viscosity.

Figure 5 Simulated functional group distributions in particle-phase compounds derived from isoprene photooxidation (Song et al., 2015) and α-pinene ozonolysis (Bateman et al., 2015; Kidd et al., 2014; Grayson et al., 2016; Zhang et al., 2015; Renbaum-Wolff et al., 2013).

Chemical composition simulated for Zhang's experiments is characterized by the highest –RCHO fraction among all experiments, which may explain the largest difference between predicted and measured viscosities. Aldehydes are known to be highly reactive in the condensed phase by reacting with alcohols and hydroperoxides to form peroxyhemiacetals and oligomers (Ziemann and Atkinson, 2012). Such multiphase reactions could be accelerated under acidic conditions in the presence of carboxylic acids (Bakker-Arkema and Ziemann, 2020; Shiraiwa et al., 2013). In addition, the GECKO-A model treats neither autoxidation to form highly oxygenated organic molecules (HOMs) nor gas-phase dimerization reactions of peroxy radicals, which are recently found to play a significant role in SOA formation and growth (Bianchi et al., 2019). These higher-molar-mass compounds with very low volatility have higher T_g, hence leading to higher SOA viscosity (Koop et al., 2011; Zhang et al., 2019; Champion et al., 2019; Li et al., 2020). While the extent of effects of these processes may vary among reaction conditions applied in each experiment, the lack of treatment of these processes may be another plausible reason of lower simulated viscosity compared to experimental measurements.

3.3 Effects of mass accommodation

All the above simulations were conducted with mass accommodation coefficient equal to unity. Hereby, we investigate potential kinetic limitations of bulk diffusion on SOA partitioning and resulting viscosity by simulating two distinct experimental conditions of flow tube experiments by Renbaum-Wolff et al. (2013) and batch experiments by Kidd et al. (2014). The effective mass accommodation coefficient (α_eff) is computed with Eq. (3), which is implemented into Eqs. (1) and (2) to simulate partitioning that effectively accounts for kinetic limitations of bulk diffusion.

Figure 6 shows the simulated viscosity estimations for SOA formation under dry conditions followed by water exposure at different RH. Final viscosity values are higher for simulations with α_eff compared to those with α fixed to 1. For simulations of Renbaum-Wolff's experiments there is very small difference in viscosity between the two simulations. For simulations of Kidd's experiments, there is a large difference by several orders of magnitude in viscosity between two simulations of different α when RH < 40 %; for example, at 40 % RH the simulated SOA viscosity is ∼ 10⁴ Pa s with α=1, while it is ∼ 10⁸ Pa s with α=α_eff. The contrasting response between the simulations of Kidd's and Renbaum-Wolff's experiments can be explained by stark differences in SOA composition (see Fig. S1 in the Supplement). This effect is caused by different experimental setups and oxidant and precursor concentrations: Kidd et al. (2014) formed SOA in a batch chamber with a low O${}_{\mathrm{3}}/\mathit{\alpha }$-pinene ratio of 0.8, while Renbaum-Wolff et al. (2013) used a continuous flow reactor with high O${}_{\mathrm{3}}/\mathit{\alpha }$-pinene ratio of 3.

Figure 6 Simulated viscosities of α-pinene SOA with experimental conditions of Kidd et al. (2014) (blue) and Renbaum-Wolff et al. (2013) (red). The simulations were conducted in two different scenarios: α = 1 (dashed lines) and α=α_eff (solid lines).

For the simulation of Kidd's experiments, ROOH and ROH compounds with relatively large molar mass, high O:C ratio and low volatility are formed at the early stage of reaction followed by rapid condensation into the particle phase (Figs. S1, S2). The relatively low ozone load limits further oxidation of ROOH and ROH compounds for fragmentation. Hence, SOA has a relatively high O:C ratio of 0.73 and low volatility, leading to very high viscosity and strong reduction of α_eff (Fig. S5) to retard SOA growth. Consequently, the simulated SOA mass loading is 15 µg m⁻³, which is much lower than the simulation with α= 1 (i.e., 7000 µg m⁻³), where higher-volatility compounds would condense without kinetic limitations, resulting in lower viscosity. The observed mass loadings in Kidd et al. (2014) were approximately ∼ 1000 µg m⁻³. This indicates that particle viscosity and resulting kinetic limitations may be overestimated in α_eff simulations or particle-phase reactions contribute substantially on SOA formation, as discussed in Sect. 3.2.

For the simulation of Renbaum-Wolff's experiments, ROH compounds are the major species formed in the gas phase during the first minutes of reaction (Fig. S1), and the continuous input of O₃ leads to the fragmentation and decomposition of ROH species into RC(O)R compounds with lower molar mass (Figs. S2, S3). Thus, SOA has relatively low O:C ratio of 0.49 and moderately low volatility (Fig. S4), resulting in lower viscosity compared to the case of simulation for Kidd et al. (2014) In this case, α_eff stays unity for low-volatile compounds (C⁰ < 10⁻² µg m⁻³), while it is reduced for semi-volatile compounds. The reduction of α of semi-volatile compounds in α-pinene SOA has been recently observed by Liu et al. (2019), as shown in Fig. S5. Note that it is difficult to make a direct comparison between simulated and measured α_eff due to the lack of detailed information on SOA formation conditions in Liu et al. (2019). The mass loadings are on the same order, being 15 and 79 µg m⁻³ for simulations with α_eff and α=1, respectively. Both results are close to the experimental measurement value of ∼ 50 µg m⁻³. As viscosity measurements were not conducted at lower RH in Kidd et al. (2014) and time evolution of SOA mass concentrations and viscosity are unavailable for both experimental studies, it is hard to fully resolve impacts of mass accommodation in this study; we intend to further investigate this aspect in a follow-up study.

4 Conclusions

We applied the T_g parameterizations and viscosity prediction method in the GECKO-A box model to simulate the evolution of viscosity and composition of SOA generated via α-pinene and isoprene oxidation. A range of simulations were performed with experimental conditions applied for viscosity measurements in order to explore various effects such as chemical composition, relative humidity, mass loadings and mass accommodation coefficient on SOA viscosity. Simulated viscosities are consistent with the observed RH dependence of viscosity for α-pinene and isoprene SOA, demonstrating the robustness of our viscosity prediction method. Simulated viscosity values for isoprene SOA are in good agreement with measurements, while those for α-pinene SOA were lower than experimental measurements by a few orders of magnitude. The simulated chemical composition and functional group distributions indicate that α-pinene SOA contain substantial amounts of aldehydes, ketones and carboxylic acids. These compounds are reactive and may undergo dimerization and oligomerization reactions. In addition, GECKO-A does not treat detailed chemistry involving stabilized Criegee intermediates, autoxidation and gas-phase dimerization by peroxy radicals. These processes are known to lead to the formation of high-molecular-mass compounds that would increase viscosity of SOA. Moreover, the model assumes that gas–particle partitioning of organics follows Raoult's law with ideal mixing conditions. The implementation of these processes is warranted in future studies for better representation of SOA chemical composition and viscosity.

Experiments were conducted to form SOA under dry conditions followed by water exposure for viscosity measurements, and it was unnecessary to consider uptake of water into the condensed phase upon SOA formation. Previous studies have suggested that water vapor can significantly affect the composition of SOA. At higher RH water vapor interacts with Criegee intermediates, leading to production of carboxylic acids or aldehydes and ketones via decomposition routes (Kristensen et al., 2014). Kidd et al. (2014) have reported that as the RH at which the SOA is formed increases, there is a decrease in viscosity, accompanied by an increasing contribution from carboxylic acids and a decreasing contribution from higher-molecular-mass products. These aspects including chemistry of Criegee intermediates and water uptake should be a subject of future studies for simulating particle viscosity under humid conditions.

We have also explored the effects of kinetic partitioning by accounting for reduction of mass accommodation coefficient in a semisolid or viscous phase. The simulation results suggest that kinetic limitations in the particle phase would result in a decrease in SOA mass loading and an increase in viscosity due to suppression of condensation of semi-volatile compounds. It is still challenging to accurately assess the extent of this effect on viscosity due to the lack of experimental data for comparison and to the absence of some chemical processes in our modeling method. Future experiments with simultaneous measurements of chemical composition in the gas and particle phases as well as particle viscosity under different RH conditions will be enormously helpful to further improve and constrain the model. While this study simulated viscosity of pure organic particles, future studies should also investigate viscosity of inorganic–organic particles and effects of non-ideal interactions including phase separation (Zuend and Seinfeld, 2012; You et al., 2014) and gel formation (Richards et al., 2020). The developed model in this study is a useful tool for further exploration of the phase state of atmospherically relevant SOA, and it should also be useful in refining the representation of SOA chemical evolution and phase state in regional and global air quality models.