Kinetic modeling of formation and evaporation of secondary organic aerosol from NO₃ oxidation of pure and mixed monoterpenes

Organic aerosol constitutes a major fraction of the global aerosol burden and is predominantly formed as secondary organic aerosol (SOA). Environmental chambers have been used extensively to study aerosol formation and evolution under controlled conditions similar to the atmosphere, but quantitative prediction of the outcome of these experiments is generally not achieved, which signifies our lack in understanding of these results and limits their portability to large-scale models. In general, kinetic models employing state-of-the-art explicit chemical mechanisms fail to describe the mass concentration and composition of SOA obtained from chamber experiments. Specifically, chemical reactions including the nitrate radical (NO₃) are a source of major uncertainty for assessing the chemical and physical properties of oxidation products. Here, we introduce a kinetic model that treats gas-phase chemistry, gas–particle partitioning, particle-phase oligomerization, and chamber vapor wall loss and use it to describe the oxidation of the monoterpenes α-pinene and limonene with NO₃. The model can reproduce aerosol mass and nitration degrees in experiments using either pure precursors or their mixtures and infers volatility distributions of products, branching ratios of reactive intermediates and particle-phase reaction rates. The gas-phase chemistry in the model is based on the Master Chemical Mechanism (MCM) but trades speciation of single compounds for the overall ability of quantitatively describing SOA formation by using a lumped chemical mechanism. The complex branching into a multitude of individual products in MCM is replaced in this model with product volatility distributions and detailed peroxy (RO₂) and alkoxy (RO) radical chemistry as well as amended by a particle-phase oligomerization scheme. The kinetic parameters obtained in this study are constrained by a set of SOA formation and evaporation experiments conducted in the Georgia Tech Environmental Chamber (GTEC) facility. For both precursors, we present volatility distributions of nitrated and non-nitrated reaction products that are obtained by fitting the kinetic model systematically to the experimental data using a global optimization method, the Monte Carlo genetic algorithm (MCGA). The results presented here provide new mechanistic insight into the processes leading to formation and evaporation of SOA. Most notably, the model suggests that the observed slow evaporation of SOA could be due to reversible oligomerization reactions in the particle phase. However, the observed non-linear behavior of precursor mixtures points towards a complex interplay of reversible oligomerization and kinetic limitations of mass transport in the particle phase, which is explored in a model sensitivity study. The methodologies described in this work provide a basis for quantitative analysis of multi-source data from environmental chamber experiments but also show that a large data pool is needed to fully resolve uncertainties in model parameters.

1 Introduction

Atmospheric aerosol particles play an important role in the Earth system by influencing weather and climate, enabling long-range transport of chemical compounds and negatively affecting public health (Pöschl, 2005; Fuzzi et al., 2006). A major contributor to the global aerosol burden is the oxidation of volatile organic compounds (VOCs) to condensable organic species, which leads to formation of secondary organic aerosol (SOA; Kanakidou et al., 2005). Important classes of SOA precursors include alkanes and aromatic compounds, which are often emitted from anthropogenic sources, as well as alkenes such as isoprene, monoterpenes, and sesquiterpenes, which are predominantly emitted by trees (Hallquist et al., 2009). The monoterpenes α-pinene and limonene are among the most abundant and well-studied SOA precursors (Seinfeld and Pandis, 2016). Atmospheric oxidation of alkenes occurs mainly through three oxidants: the hydroxyl radical (OH), which is produced in daylight and is short lived; the abundant but comparatively slow reacting ozone (O₃); and the nitrate radical (NO₃), which is the major source of SOA at nighttime but also contributes to SOA formation during daytime, despite its quick photolysis (Liebmann et al., 2019). The oxidation of VOCs by NO₃ results in the formation of high yields of various nitrated organic compounds, alkyl nitrates and peroxy acyl nitrates, which are produced in lower quantities through other atmospheric oxidation channels such as reaction of organic peroxy radicals (RO₂) or hydroperoxy radicals (HO₂) with nitric oxide (NO) (Perring et al., 2013; Ng et al., 2017). These organic nitrates (ON) play an important role in the atmospheric nitrogen budget by serving as temporary or permanent sinks for highly reactive nitrogen oxides (NO_x = NO + NO₂).

Due to their sufficiently low volatility, ON can be taken up into atmospheric aerosol particles, where they are shielded from gas-phase chemical decomposition, causing NO_x to be temporarily removed from atmospheric oxidation cycling. While NO_x can be recycled back into the atmosphere via photolysis (Müller et al., 2014), photooxidation (Nah et al., 2016b), and thermal decomposition of ON, permanent removal can occur through ON hydrolysis (Takeuchi and Ng, 2019) and deposition processes (Nguyen et al., 2015).

Furthermore, the presence of ON affects the formation and persistence of organic aerosol (OA) (Ng et al., 2017). The contribution of particulate ON mass (pON) to total organic aerosol (pON ∕ OA) has been investigated previously in laboratory studies by mass-spectrometric methods (Fry et al., 2009, 2011, 2014; Boyd et al., 2015; Nah et al., 2016b; Boyd et al., 2017; Faxon et al., 2018; Takeuchi and Ng, 2019) and a radioactive tracer method (Berkemeier et al., 2016), revealing that pON ∕ OA can reach up to 0.9 in the particle phase under certain conditions. Although ambient pON ∕ OA varies strongly temporally and regionally, measured values of the ratio of organic mass in ON to the total organic mass have been shown to reach up to 0.77 (Ng et al., 2017, and references therein).

Despite the importance of ON to the dynamics of SOA formation, the chemical mechanism for their formation in the gas and particle phases is still under discussion (Kurtén et al., 2017; Claflin and Ziemann, 2018; Draper et al., 2019). The Master Chemical Mechanism (MCM) provides a resource of the gas-phase degradation chemistry of typical SOA precursors with atmospheric oxidants (Saunders et al., 2003; Jenkin et al., 2003). However, application of MCM to the oxidation of monoterpenes with NO₃ leads to a significant underestimation of particle mass and pON ∕ OA as this mechanism is missing several important chemical reactions, for example, oxidation of the second double bond of limonene (Boyd et al., 2017; Faxon et al., 2018).

It has been hypothesized and shown recently that a majority of SOA might exist in oligomerized form (Kalberer et al., 2004; Gao et al., 2010), which might alter their evaporation behavior (Baltensperger et al., 2005; D’Ambro et al., 2018). In that case, the evaporation timescale is determined by chemical decomposition instead of equilibrium partitioning due to volatility (Pankow, 1994). Additionally, organic aerosol particles can exhibit a highly viscous phase state (Virtanen et al., 2010; Koop et al., 2011; Reid et al., 2018), which leads to kinetic limitations in evaporation (Vaden et al., 2011), slowing of particle-phase chemistry (Gatzsche et al., 2017), and non-equilibrium partitioning (Cappa and Wilson, 2011).

To describe kinetic limitations in mass transport, a number of kinetic multi-layer models have been developed recently to describe aerosol particles and cloud droplets, including KM-SUB (Shiraiwa et al., 2010), KM-GAP (Shiraiwa et al., 2012), ADCHAM (Roldin et al., 2014), and MOSAIC (Zaveri et al., 2008, 2014). These models are capable of explicitly resolving mass transport and chemical reactions within aerosol particles. Using these models, Shiraiwa et al. (2013) and Zaveri et al. (2018) were able to find evidence for diffusion limitation affecting SOA formation dynamics by inspection of the evolution of particle size distributions. Yli-Juuti et al. (2017) and Tikkanen et al. (2019) used an evaporation model based on KM-GAP to describe the interaction of volatility and viscosity during isothermal dilution as a function of different environmental conditions. However, to the best of our knowledge, no model has been presented that describes all aspects of gas-phase chemistry, particle-phase chemistry, gas–particle partitioning, and bulk diffusion of SOA.

A model capable of describing all these aspects of SOA formation must rely on a large set of kinetic parameters, which are often not readily accessible. However, model parameters can be systematically altered so the model matches experimental data, an approach often referred to as inverse modeling. Simultaneously optimizing multiple model parameters can often be unfeasible via manual optimization and prompts the use of global optimization methods (Berkemeier et al., 2013, 2017). As opposed to local optimization methods, global optimization algorithms are not as easily stuck in local minima and are able to reliably find solutions of difficult optimization problems. In conjunction with a kinetic model, global optimization algorithms represent a powerful tool that allows inference of molecular-level information from macroscopic data. Thus, global optimization algorithms based on differential evolution, such as the Monte Carlo genetic algorithm (MCGA), have become increasingly popular in the modeling of complex multiphase chemical systems (Berkemeier et al., 2017; Marshall et al., 2018; Tikkanen et al., 2019).

In a previous study, Boyd et al. (2017) showed that the retained aerosol mass from oxidation of limonene with NO₃ after heating from 25 to 40 ^∘C is significantly different than the mass obtained from oxidizing limonene at 40 ^∘C. They further showed that the evaporation behavior of mixtures of limonene SOA and β-pinene SOA crucially depends on the order in which oxidation occurred. Oxidation of limonene followed by subsequent oxidation of β-pinene led to an aerosol that exhibits much slower evaporation of limonene compared to an aerosol produced by simultaneous oxidation of the two precursors. At the time, it was only postulated that diffusion limitations and/or oligomerization reactions could have led to these observations. In this work, we conduct new environmental chamber experiments and apply a novel kinetic modeling framework to investigate whether gas-phase chemistry, equilibrium partitioning, and particle-phase chemistry can accurately describe the formation and evaporation of monoterpene SOA from oxidation of α-pinene, limonene, and mixtures of both precursors with NO₃. α-Pinene is chosen over β-pinene since it shows a more distinct evaporation behavior to limonene SOA and is the overall better-understood SOA precursor. We perform experiments at a lower initial temperature compared to Boyd et al. (2017) to include a second heating stage in the experiments. We focus the modeling efforts on the experimental observables aerosol mass and organic nitrogen content (pON ∕ OA) as a function of time in the reaction chamber. The model uses a simplified, lumped kinetic mechanism based on MCM (Berkemeier et al., 2016) but modifies some of the branching ratios in RO₂ chemistry and adds chemical reactivity in the particle phase. Building on the observations of Boyd et al. (2017) in their mixed precursor experiments, we investigate the linearity of these two observables by quantitative comparison of formation and evaporation of SOA from pure and mixed monoterpene precursors. We first test the hypothesis of whether particle-phase oligomerization in a well-mixed liquid phase can explain the observed behavior. Then, we use the kinetic model to perform a sensitivity analysis on the potential effect of retarded bulk diffusion due to a viscous phase state. The kinetic modeling framework consisting of a kinetic multi-layer model based on KM-GAP and the MCGA algorithm is used as an analysis tool to explore the mechanistic interactions between reactive intermediates and oxidation products that can lead to non-additivity of the investigated reaction systems.

2 Experimental and theoretical methods

2.1 Georgia Tech Environmental Chamber (GTEC)

The aerosol formation and evaporation experiments are performed as batch reactions in the GTEC facility, which consists of two separate 12 m³ Teflon chambers in a temperature- and humidity-controlled enclosure (Boyd et al., 2015). A consistent experimental routine is maintained for all experiments presented in this study and resembles the method used by Boyd et al. (2017) with small updates. Concentrations of O₃ and NO_x are determined with a UV absorption O₃ analyzer (Teledyne T400) and a chemiluminescence NO_x monitor (Teledyne 200 EU), respectively. Aerosol particle number and volume concentrations are measured using a scanning mobility particle sizer (SMPS, TSI), which consists of a differential mobility analyzer (DMA, TSI 3040) and a condensation particle counter (CPC, TSI 3775). Bulk aerosol composition is measured using a high-resolution time-of-flight aerosol mass spectrometer (HR-ToF-AMS, DeCarlo et al., 2006).

The Teflon chamber is flushed with zero air for at least 24 h, and the chamber enclosure is cooled to 5 ^∘C several hours prior to each experiment, to ensure full equilibration with regard to temperature, pressure, and humidity. Monoterpene oxidation is initiated at 5 ^∘C and under dry conditions (RH < 5 %). All experiments are conducted using ammonium sulfate seed particles. Seed particles are generated by atomizing a 15 mM ammonium sulfate solution into the chamber for 20 min, which typically results in particle number concentrations around 20 000 cm⁻³ and mass concentrations of 28–41 µg∕m³. Simultaneously, monoterpene precursors are injected into the chamber. Injection volumes of the precursors are chosen to achieve consistent total aerosol mass concentrations around 100 µg∕m³ in all experiments, based on knowledge about aerosol yields in trial experiments for this study. For α-pinene, we use a microsyringe to inject a known volume of liquid into a mildly heated glass bulb from which a 5 L∕min zero air flow carries the evaporating fumes into the chamber. For limonene, the required liquid volume is so low that the use of microsyringes is a source of non-negligible uncertainty, and hence a gas cylinder filled with 0.85 ppm limonene, calibrated and confirmed using gas chromatography with flame ionization detection (GC-FID), is used to inject a known volume of gas into the chamber over the course of several minutes. NO₃ is produced by oxidation of NO₂ with O₃ (generated by passing zero air through a photochemical ozone generator) in a 1.5 L flow tube (0.9 L∕min flow, 100 s residence time). The reaction mixture is optimized so NO₃ and N₂O₅ are produced in high yields, with no significant amount of O₃ entering the chamber. This is achieved by using a 2:1 ratio of NO₂ and O₃. N₂O₅ decomposes in the chamber to release NO₃ over time. Injection of NO₃ ∕ N₂O₅ marks the beginning of the reaction.

When peak SOA mass is reached, which is typically achieved in under 4 h, the chamber enclosure temperature is raised to 25 ^∘C and, after another waiting period, to 42 ^∘C. The temperature changes take approximately 90 min in both cases. Temperature profiles are reported alongside the experimental results in Fig. 2.

In total, four experiments are conducted, either with a single monoterpene precursor, pure α-pinene (APN) and pure limonene (LIM), or with a mixture of both precursors. In the case where both precursors are used, the oxidation occurred in one of two variants: simultaneous (MIX) or sequential oxidation (SEQ). In case of the MIX experiment, both precursors are injected simultaneously into the chamber prior to NO₃ ∕ N₂O₅ injection. In case of the SEQ experiment, peak SOA mass after the first precursor oxidation is first awaited. Then, a second NO₃ ∕ N₂O₅ injection and injection of the second VOC follow in sequence. An 8-fold excess of N₂O₅ is used for pure limonene experiments, and a 4-fold excess is used for pure α-pinene experiments. In the mixed precursor experiments, the amount of injected NO₃ ∕ N₂O₅ is determined using the same ratios proportionately. A summary of all experimental conditions, including injected precursor amounts, seed mass, total aerosol mass, SOA mass excluding seed, and SOA yields, can be found in Table 1. It is noted that we refer to the total aerosol mass concentration (sum of inorganic seed mass concentration and SOA mass concentration) in the chamber simply as “aerosol mass” in our discussions. We use the term “SOA yield” to refer to the ratio of produced organic aerosol mass concentration to the initial VOC mass concentration.

2.2 Kinetic model

The kinetic model calculations in this study are performed with a multi-compartmental model akin to the KM-SUB/KM-GAP model family (Shiraiwa et al., 2010, 2012). The model code is set up as a generator script that uses an input chemical mechanism to generate a system of differential equations that is able to describe the key physical and chemical processes in the GTEC chamber. The model compartments include the chamber wall, the wall near-surface gas phase, the chamber gas phase, the particle near-surface gas phase, the particle surface, and the particle bulk. The processes explicitly described in the model include injection of chemical compounds, irreversible loss of wall-adsorbed species, temperature change, gas diffusion to the chamber wall, gas diffusion to particles, condensation and evaporation at the wall and particle surfaces, and chemical reactions in the gas and particle phases. Wall loss of particles is implicitly accounted for in this study by using wall-loss-corrected SMPS data (Keywood et al., 2004; Nah et al., 2017).

All product molecules with vapor pressures lower than 1 Pa at 298 K are allowed to partition into the topmost layer of the particles, according to their volatility. Gas–particle partitioning is explicitly treated in the model, and equilibration between the particle near-surface gas phase and the particle surface is achieved by balancing surface adsorption and desorption rates. This way, evaporation and condensation kinetics are treated more realistically than in a model assuming instantaneous equilibrium partitioning. The adsorption flux J_ads, X of a molecule X is calculated from the collision flux from the particle near-surface gas phase to the particle surface, which in turn is calculated from the mean thermal velocity ω_X, and the accommodation coefficient α_s, X. α_s, X is assumed to be 0.5 for all organic species in this study, in line with previous investigations (Julin et al., 2013; von Domaros et al., 2020). A sensitivity study on the effect of α_s, X on model output can be found in Fig. S1 in the Supplement.

$$\begin{array}{}\text{(1)}& J_{\mathrm{ads},X}={\mathit{\alpha }}_{\mathrm{s},X}\cdot {\displaystyle \frac{{\mathit{\omega }}_X}{\mathrm{4}}}\cdot [X{]}_{\mathrm{gs}}\end{array}$$

The desorption flux from the particle surface to the gas phase J_des, X is dependent on the vapor pressure p_vap, X and the ratio of the concentration of X in the particle near-surface bulk layer [X]_b1 (in unit cm⁻³) and the sum of all other species Y_j in that layer.

$$\begin{array}{}\text{(2)}& J_{\mathrm{des},X}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{s},X}\cdot {\mathit{\omega }}_X\cdot p_{\mathrm{vap},X}\cdot N_{\mathrm{A}}\cdot [X{]}_{\mathrm{b}\mathrm{1}}}{\mathrm{4}\cdot R\cdot T\cdot \sum [Y_j{]}_{\mathrm{b}\mathrm{1}}}}\end{array}$$

Here, R is the universal gas constant, T is the temperature in kelvin (K), and N_A is Avogadro's number. The vapor pressure of product compounds is assumed to be temperature dependent with a precursor-dependent effective enthalpy of volatilization, ΔH_vap, Z in kilojoules per mole (kJ∕mol), where Z is the precursor of X. We assume this single effective enthalpy to be representative for the entire product spectrum and hence independent of C^*.

$$\begin{array}{}\text{(3)}& p_{\mathrm{vap},X}\left(T\right)=p_{\mathrm{vap},X}\left(\mathrm{298}\mathrm{K}\right)\cdot \exp {\displaystyle \frac{-\mathrm{\Delta }{H}_{\mathrm{vap},Z}}{R\cdot (T-\mathrm{298})}}\end{array}$$

Note that, while the employed model is inherently a multi-layer model, only a single well-mixed layer is used to describe the aerosol phase in the default calculations in this study. Multiple layers were used for the calculations in Sect. 3.5.3, leading to Fig. 6. New particle formation from low-volatility vapors is not treated in this model, so seed particles have to be pre-defined. Seed particles are initialized as covered with a very small amount of non-volatile organics ($\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ ppb gas-phase mixing ratio) to aid in computation of gas–particle partitioning. The model can be run in two modes: lumped mode, in which only vapor pressure bins are defined; and explicit mode, in which vapor pressures must be pre-supplied for all participating species. In the following, we will describe the specific lumped mode used in this study.

Reversible and irreversible vapor wall loss is described following Huang et al. (2018), with slight modifications to fit into the KM-SUB/KM-GAP model structure. The Teflon wall is described using two layers: a surface layer, to which vapor molecules partition reversibly; and an inner layer, into which vapor molecules diffuse irreversibly on the timescale of the experiment. The wall adsorption flux $J_{\mathrm{ads},X,\mathrm{wall}}$ is parameterized according to Eq. 4.

$$\begin{array}{}\text{(4)}& J_{\mathrm{ads},X,\mathrm{wall}}={\mathit{\alpha }}_{\mathrm{wall}}\cdot {\displaystyle \frac{{\mathit{\omega }}_X}{\mathrm{4}}}\cdot [X{]}_{\mathrm{ws}}\end{array}$$

[X]_ws is the wall near-surface gas-phase concentration of X. The wall accommodation coefficient α_wall, X is parameterized according to Eq. 5.

$$\begin{array}{}\text{(5)}& {\mathit{\alpha }}_{\mathrm{wall},X}={\mathrm{10}}^{-\mathrm{2.744}}\cdot {C_X^{*}}^{-\mathrm{0.6566}}\end{array}$$

$C_X^{*}$ is the saturation mass concentration of X, which indicates the organic aerosol mass at which a semi-volatile organic substance would be in the gas and particle phase in equal parts. The wall desorption flux $J_{\mathrm{des},X,\mathrm{wall}}$ is parameterized according to Eq. (6).

$$\begin{array}{}\text{(6)}& J_{\mathrm{des},X,\mathrm{wall}}={\mathit{\alpha }}_{\mathrm{wall}}\cdot {\displaystyle \frac{{\mathit{\omega }}_X}{\mathrm{4}}}\cdot {\displaystyle \frac{{\mathit{\gamma }}^{\mathrm{\infty }}{C}_X^{*}{M}_{\mathrm{W},\mathrm{wall}}}{{\mathrm{10}}^{\mathrm{3}}{C}_{\mathrm{wall}}{M}_{\mathrm{W},X}}}\end{array}$$

Here, γ^∞ is the activity coefficient in Teflon and C_wall the effective organic mass concentration of the wall itself and is set to be 32.2 mg m⁻³ (Huang et al., 2018). M_W, X and M_W, wall denote the molecular weight of X and the effective molecular weight of the Teflon wall, respectively.

$$\begin{array}{}\text{(7)}& {\mathit{\gamma }}^{\mathrm{\infty }}={\mathrm{10}}^{\mathrm{3.299}}\cdot {C_X^{*}}^{-\mathrm{0.6407}}\end{array}$$

The gas diffusion flux from the chamber interior to the wall near-surface gas phase $J_{\mathrm{dif},X,\mathrm{ws}}$ is described using the Fickian gas diffusivity coefficient D_g, X and an additional Eddy diffusivity coefficient k_e, which was estimated to be 0.03 s⁻¹ for the GTEC chamber in a previous study (Nah et al., 2016a).

$$\begin{array}{}\text{(8)}& J_{\mathrm{dif},X,\mathrm{ws}}={\displaystyle \frac{\mathrm{2}}{\mathrm{pi}}}\sqrt{D_{\mathrm{g},X}{k}_e}[X{]}_{\mathrm{g}}\end{array}$$

Note that the explicit treatment of a near-surface gas phase at the wall constitutes a slight variation from the framework of Huang et al. (2018), who treated gas diffusion and adsorption simultaneously in a resistor-style approach. The two resistor terms were split into the separate fluxes $J_{\mathrm{ads},X,\mathrm{wall}}$ and $J_{\mathrm{dif},X,\mathrm{ws}}$ in this study. The thickness of the near-wall gas phase had only little impact on calculation results in the range of 0.1 mm–1 cm and was set to the higher limit of 1 cm for numerical stability. Irreversible transport from the Teflon surface layer to the inner Teflon layer is assumed to occur at a first-order rate l_w, i and treated as independent of volatility of the organic molecule. l_w,i is obtained by fitting the model to experimental data and typically falls around 10⁻⁴ s⁻¹ or 0.3 h⁻¹ (cf. Table 1).

2.3 Lumped chemical mechanism

The gas-phase chemical mechanism, summarized in Fig. 1a, is modeled after the initial reaction steps in the MCM but does not assume specific sum or structural formulas of product molecules. The validity of this approach has been shown in previous work (Berkemeier et al., 2016). For limonene SOA, we apply the same general chemistry but consider the oxidation of both double bonds individually, which leads to the more complex reaction scheme shown in Fig. S2. Note that oxidation of the second double bond of limonene with NO₃ is not considered in MCM. However, we have shown previously that including oxidation of the second double bond leads to a significantly improved correlation between a kinetic model and chamber experiments (Boyd et al., 2017).

Figure 1 (a) Schematic representation of the lumped chemical mechanism for oxidation of monoterpenes with one double bond (e.g., α-pinene). The asterisk stands for chemical reaction with NO, NO₃, and RO₂. (b) The stable products are divided into six product bins each with a different volatility (grey arrows; bin 1–bin 6), according to a probability distribution (example graphs on the right). (c) Oligomerization occurs in equilibrium reactions in the particle phase under conservation of precursor origin and volatility bin.

To account for chemical identity, the major product classes, nitrated and non-nitrated organic molecules, are subdivided into logarithmically spaced volatility bins (Fig. 1b) following the concept of a volatility basis set (VBS; Donahue et al., 2011). To minimize the number of model parameters, six volatility bins are chosen with higher resolution in and around the experimental range (1−1000 µg∕m³) to achieve high sensitivity. To also cover a wide range of volatilities, a very low volatility and a very high volatility bin are included at the ends of the spectrum: (1) $\mathrm{9.91}\times {\mathrm{10}}^{-\mathrm{8}}$ Pa (C^* = 0.01 µg∕m³), (2) 9.91$\times {\mathrm{10}}^{-\mathrm{6}}$ Pa (C^* = 1 µg∕m³), (3) 9.91$\times {\mathrm{10}}^{-\mathrm{5}}$ Pa (C^* = 10 µg∕m³), (4) 9.91$\times {\mathrm{10}}^{-\mathrm{4}}$ Pa (C^* = 100 µg∕m³), (5) 9.91$\times {\mathrm{10}}^{-\mathrm{3}}$ Pa (C^* = 1000 µg∕m³), and (6) 9.91$\times {\mathrm{10}}^{-\mathrm{1}}$ Pa (C^* = 100 000 µg∕m³) at 298 K. Oligomeric species are chosen to be fully non-volatile and hence technically form a seventh volatility bin. The average molar mass of molecules in the organic aerosol phase is assumed to be 250 g∕mol, which is similar to assumptions in previous publications (Berkemeier et al., 2016) and consistent with our measurements using chemical ionization high-resolution time-of-flight mass spectrometry with a special filter inlet that samples both the aerosol and gas phase (FIGAERO-HR-ToF-CIMS; Lopez-Hilfiker et al., 2014) that were conducted alongside this study (Takeuchi and Ng, 2019).

A specific aim of this study is the mechanistic analysis of ON formation. Therefore, the gas-phase formation of ON is treated in detail and has been expanded from the MCM template, which is detailed in Fig. S3. We assume that chemical reaction of NO₃ with the terpenic precursor yields a nitrated peroxy radical (R^NO₂). The fate of the nitrate group (–ONO₂) in this radical is dependent on its radical branching ratios. Following MCM, we assume that the reaction of R^NO₂ with HO₂ yields a stable organic nitrate product, whereas reaction with NO, NO₃, RO₂, or unimolecular decay leads to formation of a nitrated alkoxy radical (R^NO), which can further stabilize under elimination of the nitrate group. Reaction of two RO₂ may also yield dimers. Another channel of ON formation is the reaction of a non-nitrated peroxy radical (RO${}_{\mathrm{2}}^{\mathrm{II}}$) with NO. Following MCM, we assume that only RO${}_{\mathrm{2}}^{\mathrm{II}}$, which is the main intermediate in monoterpene OH oxidation and a secondary intermediate of monoterpene ozonolysis, can undergo this reaction and is in that regard distinct from RO${}_{\mathrm{2}}^{\mathrm{I}}$, which is the main intermediate in monoterpene ozonolysis. However, this RO${}_{\mathrm{2}}^{\mathrm{II}}$ + NO reaction channel has only minor implications in this study due to the low prevalence of NO under the employed reaction conditions, i.e., injection of NO₃ ∕ N₂O₅ as well as no irradiance with UV lights.

Particle-phase chemistry is included as formation and decomposition of oligomers from monoterpene oxidation products. Possible reaction pathways for oligomerization include the formation of esters, aldols, hemiacetals, acetals, peroxyhemiacetals, and peroxyacetals from alcohol, aldehyde, hydroperoxide, and carboxylic acid moieties in the monoterpene oxidation products (Ziemann and Atkinson, 2012) but are lumped into a single reaction for simplicity. These oligomers are assumed to be non-volatile, which is in line with recent investigations (DePalma et al., 2013; Barsanti et al., 2017), but can re-partition back to the gas phase after decomposition into the monomeric building blocks. Oligomer decomposition is treated as temperature dependent with a precursor-specific activation energy $E_{\mathrm{A},\mathrm{decom},Z}$ of precursor Z to be used in an Arrhenius equation. The information about volatility and nitration degree of monomers is retained during oligomerization and reinstated after their decomposition. This process is outlined in Fig. 1c. A discussion of the oligomerization scheme is provided in the Supplement, Sect. S1. An overview of all reactions of the lumped model in the gas and particle phases is given in Table S1.

2.4 Global optimization

The Monte Carlo genetic algorithm (MCGA; Berkemeier et al., 2017) is applied for inverse fitting of the kinetic model to the experimental data and determining the non-prescribed kinetic parameters listed in Table 1. The MCGA method consists of two steps: a Monte Carlo step and a genetic algorithm step. During the Monte Carlo step, kinetic parameter sets are randomly sampled from a defined parameter range, and the residual between the model result and the experimental data is determined for each parameter set through evaluation of the kinetic model. During the genetic algorithm step, the parameter sets are optimized mimicking processes known from natural evolution: a survival mechanism retains best-fitting parameter sets, the recombination mechanic generates new parameter sets by combing parameters of high scoring sets, and the mutation step prevents early homogenization of the sample of parameter sets. To determine the model–experiment correlation, we use a least-squares approach that minimizes the sum of the squares of the residuals, Eq. (9). The estimator is normalized to the magnitude of the largest data point in a given sample, max (Y_i), and the number of data points n_i of data set i. Additionally, optional weighting factors w_i can be used to guide the optimization process. In this study, pure precursor experiments are each weighted twice as high as the mixed precursor experiments to ensure that any non-linearity in the mixed precursor experiments is detected as a deviation between the model and the experiment for those experiments. pON ∕ OA data are weighted by a factor of 4 less than SOA mass data as the focus of this paper is the formation and evaporation behavior of SOA and more assumptions go into the determination of pON ∕ OA.

$$\begin{array}{}\text{(9)}& f_i=w_i\sqrt{{\displaystyle \frac{\mathrm{1}}{n_i}}\sum \left({\displaystyle \frac{Y_{\mathrm{model}}-Y_{\mathrm{data},i}}{\max \left(Y_{\mathrm{data},i}\right)}}\right)}\end{array}$$

After an optimization result is returned, a one-dimensional golden-section search (Press et al., 2007, Sect. 10.2) is used to ensure conversion into a minimum of the optimization hypersurface. The simplex method (Press et al., 2007, Sect. 10.5) is used to find other combinations of parameters that lead to equivalent model results (test of uniqueness). Weighting factors w_i can be used to assign a lower importance to data sets that, e.g., exhibit large scatter due to experimental noise, represent experimental artifacts, or are deemed only supplementary for the purpose of the optimization.

Table 1 Fit parameters of the kinetic model. Error estimates for the volatility distribution (parameters f_apin and f_lim) can be found in Fig. S4 in the Supplement; error estimates for all other parameters are ranges in which a parameter can be varied until the model–experiment correlation decreases by 10 %. For a full list of kinetic parameters, see Table S1.

Note that for the experiments discussed in this paper, multiple model solutions can be obtained, dependent not only on the choice of data sets that is optimized to, but also on the choice of weighting factors. In the following sections, we focus our discussion on one fit of the model to experimental data as it scored best in our choice of model–experiment correlation estimator (fit 1, f=0.88 according to Eq. 9). Multiple evaluations of MCGA typically give similar results to fit 1 but sometimes get stuck in local minima that are significantly worse. This is a direct consequence of undersampling with MCGA, given the large amount of model input parameters. Typically, about 150 000 parameter sets were sampled during a MCGA run, which is not sufficient given the number of input parameters but marks an upper achievable range for this study as it takes about 3 d to complete on an 80 CPU computer cluster. Among the inferior fits that were obtained, we also found a distinct fit that scores worse overall (fit 2, f=0.097) but scores better in some aspects of the data set and will be discussed alongside fit 1. We will discuss the dependence of the best fit on weighting factors and the uniqueness of the obtained model solution in Sects. 3.5 and 4.

3 Results and discussion

3.1 Pure limonene oxidation (LIM)

3.1.1 Experimental observations (LIM)

Figure 2a shows the total aerosol mass concentration (denoted as “aerosol mass”) during an experiment of limonene oxidation with NO₃ in the presence of ammonium sulfate seed particles, as well as subsequent evaporation in the GTEC chamber, here referred to as the “LIM” experiment. Oxidation at 5 ^∘C initially causes a fast increase in aerosol mass (black open markers, left axis) from 29 µg∕m³ of seed mass to about 70 µg∕m³ of aerosol mass within the first 20 min of the experiment. Afterwards, aerosol growth slows down considerably, so that the peak aerosol mass of 110 µg∕m³ is reached only after 5 h. The slow increase in SOA mass is likely an important feature of the experimental data for determination of mass transfer and chemical reaction rates.

Figure 2 (a) Comparison of experimental and modeling results for oxidation of limonene with NO₃. Open black markers are experimental aerosol mass obtained using an SMPS. The red solid line represents the best-fit model result, the red dotted line represents an alternative model fit, and the grey dashed line corresponds to the experimental temperature profile. (b) Analysis of the oligomerization state of particle-phase products in the model according to fit 1. (c) Analysis of the occupation of volatility bins of all products according to fit 1 and at peak SOA mass. Shadings in the bar plot denote where molecules of a certain volatility bin reside: gas phase (grey) or particle phase (colored). Products in the particle phase are further distinguished as organic nitrates (green) and non-nitrated organics (orange).

The produced SOA mass corresponds to a SOA yield of 130 % (Table 2) and is observed to be constant in the chamber for several hours at 5 ^∘C. Note that this observation is different from previous experiments conducted at 25 ^∘C and 40 ^∘C (Boyd et al., 2017), where peak aerosol mass was achieved swiftly and SOA yields at aerosol mass loadings similar to this study were determined to be 174 % and 94 %, respectively. While the lower SOA yield at 40 ^∘C compared to 25 ^∘C can be explained with the equilibrium partitioning theory, the slightly lower mass yield observed at 5 ^∘C in this study cannot be.

Table 2 Experimental conditions for environmental chamber experiments presented in this study alongside aerosol masses and SOA yields at 5 ^∘C.

^* Aerosol masses are calculated from aerosol volume concentrations using a density of (NH₄)₂SO₄ seed particles of 1.75 g∕cm³, the organic phase of 1.64 g∕cm³ for limonene SOA (Boyd et al., 2017), 1.46 g∕cm³ for α-pinene SOA (Nah et al., 2016b), and 1.55 g∕cm³ for the mixtures. SOA mass is calculated as the difference between peak total aerosol mass and pre-growth seed mass. All the reported masses are wall loss corrected.

After 7 h of total experiment time, the temperature set point of the chamber enclosure is increased to 25 ^∘C. The new temperature plateau is reached inside the Teflon chamber 90 min later (grey dashed line, right axis). The temperature change causes a slight reduction in aerosol mass from 110 to about 104 µg∕m³. At the new temperature set point, aerosol mass is not constant but rather decays at a constant rate. After about 19 h, the temperature set point is increased to 42 ^∘C, which again causes an immediate slight reduction in aerosol mass from 90 to about 83 µg∕m³. At the new temperature plateau of 42 ^∘C, aerosol mass once again decays at a constant rate that is comparable to the one previously observed.

3.1.2 Kinetic modeling results (LIM)

In the following, kinetic modeling results are discussed in terms of a best fit (fit 1) that is obtained using the Monte Carlo genetic algorithm (MCGA). An alternative fit (fit 2) was obtained but is indistinguishable from fit 1 for the LIM experiment. The uniqueness of these fits and potential pitfalls of the optimization process are discussed in Sects. 3.5 and 4.

Under the conditions employed in this study, limonene precursor oxidation is dominated by NO₃ oxidation. RO₂ fate is dominated by reaction with NO₃ and RO₂ as very little NO and HO₂ are present in the chamber. The kinetic model (red solid and dotted line in Fig. 2a) is able to reproduce the observed aerosol formation and evaporation behavior. In the model run at hand, the initial quick increase in aerosol mass is due to condensation of dimers formed in the gas phase through the RO₂ + RO₂ channel (from now on referred to as “gas-phase dimers”), making up about 50 % of condensing material in the initial seconds. Subsequent growth is due to condensation of monomeric oxidation products (from now on referred to as “monomers”) of sufficiently low volatility (Fig. 2b). When half of the peak SOA mass is reached, the particle phase is to a large extent comprised of monomeric compounds, about 40 % of which still contain a C–C double bond (Fig. S5). These mono-unsaturated oxidation products either partition back into the gas phase where they can be oxidized further or co-oligomerize in the particle phase with other oxidation products.

The vapor pressures of the non-nitrated and nitrated mono-unsaturated oxidation products were fitted during the MCGA optimization and determined to have saturation mass concentrations C^* of 5.5×10⁵ and 7.43×10³ µg∕m³ at 298 K, respectively. This means that the non-nitrated intermediate is fully volatile and the non-nitrated intermediate partitions to some extent into the particle phase. At peak SOA mass, 33 % of oxidation products still contain a double bond in this model run, all of which are nitrated and present as oligomers. Note that this is possible because we do not consider the oxidation of unsaturated compounds in the particle phase.

The volatility distribution key determined by global optimization can be found in Fig. S4. A large fraction of limonene oxidation products in this model run occupies the sixth and highest volatility bin ($C^{*}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ µg∕m³ at 298 K), which is mostly present in the gas phase under these reaction conditions. Figure 2c shows the resulting volatility distribution of organics in the particle phase according to the model at peak SOA mass, which lacks organic material from the highest volatility bin. In the model, the slow increase in aerosol mass from about 80 to 110 µg∕m³ is due to oligomerization of monomers forming higher molecular weight structures through accretion reactions in the particle phase (from now on referred to as “oligomers”). According to the model fit, oligomerization occurs at a rate of $k_{\mathrm{form},\lim }=\mathrm{17.2}$ h⁻¹. Barsanti et al. (2017) compiled accretion rate coefficients with relevance to SOA formation and report rate coefficients for hemiacetal formation under neutral conditions in methanol of 0.1 ${\mathrm{M}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ and peroxyhemiacetal formation of 0.5–70 ${\mathrm{M}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$. Assuming that every limonene oxidation product has two reactive sites to undergo oligomer formation, k_form,lim can be translated into a second-order reaction rate coefficient of 1.4 ${\mathrm{M}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ and thus lies in close proximity to literature values.

Oligomerization slowly removes semi-volatile species in the particle phase from the partitioning equilibrium, which in turn causes a flux of semi-volatile molecules from the gas phase into the particle phase. The highest volatility components partition into the particle oligomer phase slowest, causing the slow increase in limonene SOA mass over 5 h. Quantum chemical and mechanistic studies have previously predicted such pronounced differences between the volatility of typical oxidation products of monoterpenes and their oligomers of several orders of magnitude (DePalma et al., 2013; Barsanti et al., 2017).

At peak SOA mass, the model predicts most of the organic material in the particle phase to exist in an oligomeric state (Fig. 2b), which explains the lack of initial evaporation caused by an increase in chamber temperature. A model fit to the LIM experimental data was attempted without inclusion of particle-phase oligomerization reactions. The model output of this simulation run shows an overall low correlation to the experimental data as it cannot explain the long time to reach peak SOA mass and the slow mass decrease at 42 ^∘C (Fig. S6).

The slow decay of aerosol mass between 6 and 24 h of the experiment is attributed in the model to a slow unimolecular decay of oligomeric material with a rate of 0.09 h⁻¹ and subsequent evaporation of monomers at elevated temperatures, followed by deposition and irreversible loss of vapors on the chamber walls. The decomposition rate suggested by the model agrees well with the rate of 0.06–0.2 h⁻¹ reported by D’Ambro et al. (2018) for SOA formed from ozonolysis of α-pinene. Following Le Chatelier's principle, removal of monomers from the equilibrium causes a constant flux of organic matter from oligomeric to monomeric state. Since the volatility of the monomeric subunit is retained in the model, this process is faster for monomers that have higher volatilities because they partition into the gas phase more quickly and readily, causing an enrichment of low-volatility monomeric subunits in the particle phase. The (meta-)stability of organic material in the particle phase can hence be attributed not only to the stability of the oligomer bond, but also the volatility of the monomeric building blocks at that temperature.

Monomers are removed from the system by deposition onto and diffusion into the chamber walls, which is the main driver of loss of organic mass. The irreversible loss rate of wall-adsorbed molecules into the chamber wall is determined to be $l_{\mathrm{w},\mathrm{i}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, which is within the range of values reported as re-evaluation from literature data in Fig. 5 of Huang et al. (2018). Figure S7 shows the distribution of organic molecules between wall, particle, and gas phase in the model for all experiments conducted in this study. The dependence of model output on l_w, i is explored in Fig. S8, indicating that the model output simulating the LIM experiment is more sensitive to changes in l_w, i than the simulation of the APN experiment described below, which can be attributed to the slow uptake and oligomerization process of semi-volatile molecules that stands in competition with irreversible wall loss.

The global optimization returned a value of 164 kJ∕mol for the effective enthalpy of vaporization ΔH_vap of limonene oxidation products. This number stands in contrast to values used for monoterpenes in SOA models such as ECHAM-HAM (59 kJ∕mol; Saathoff et al., 2009), GEOS-Chem (42 kJ∕mol; Chung and Seinfeld, 2002), or GISS modelE (72.9 kJ∕mol; Tsigaridis et al., 2006) but agrees with the value of about 160 kJ∕mol obtained in Boyd et al. (2017) at a similar mass loading. Boyd et al. compared SOA yields at two different temperatures for a range of initial precursor concentrations and determined ΔH_vap based on the Clausius–Clapeyron equation. A sensitivity study on the effect of ΔH_vap on model output is shown in Fig. S9.

The results obtained in this study can be compared to and used to interpret results in a previous study by Boyd et al. (2017). This study observes a lower SOA yield at 5 ^∘C (130 %) compared to the previous experiments performed at 25 ^∘C (174 %). This finding cannot be explained by gas–particle partitioning alone, as lower temperatures should give rise to higher SOA yields. A probable cause could be the temperature dependence of the gas-phase oxidation chemistry; however, test calculations using the temperature-dependent rate coefficients reported in the MCM mechanism showed hardly any effect of temperature on SOA yield. Thus, another promising explanation is the temperature dependence of the oligomerization rate coefficient. As the model calculations highlight, condensation of vapors onto the suspended particles stands in competition with loss to the chamber walls, which should not be strongly temperature dependent. When oligomerization occurs more slowly, oxidation products from higher volatility bins are increasingly lost to the walls before they can be incorporated into the particle oligomer phase. This is confirmed by a sensitivity study that shows a strong influence of oligomer formation rate k_form, lim on model output (Fig. S10). In addition to temperature dependence of the rate coefficient itself, oligomerization turnover might be effectively depressed by a semi-solid phase state at 5 ^∘C as discussed in Sect. 3.5.3.

Another observation in Boyd et al. (2017) was a lower SOA mass when directly forming limonene SOA at 40 ^∘C compared to first forming limonene SOA at 25 ^∘C and then heating to 40 ^∘C. This observation could also be explained by the successive condensation and oligomerization of semi-volatile vapors suggested by the model in this study. The fraction of chemical species from the higher volatility bins that partitions into the particle phase is much smaller at 40 ^∘C compared to 25 ^∘C. This may prevent the additional slow mass accumulation through oligomerization of semi-volatile oxidation products at 40 ^∘C and result in a lower SOA yield.

3.2 Pure α-pinene oxidation (APN)

3.2.1 Experimental observations (APN)

Figure 3a shows the aerosol mass during the corresponding experiment of α-pinene oxidation with NO₃, here referred to as the “APN” experiment. Similar to the LIM experiment described above, oxidation at 5 ^∘C initially causes a fast increase in aerosol mass (black open markers); however, peak aerosol mass is already reached after 3 h of oxidation at 109 µg∕m³ and a corresponding SOA yield of 25.2 % (Table 2). At a comparable organic mass, this yield is significantly lower than observed in the limonene oxidation experiment. Note that, in order to achieve similar SOA mass loadings among all experiments in this study, a larger amount of precursor is added in the α-pinene oxidation experiment.

The SOA yield in this study appears to be larger than previously reported for the oxidation of α-pinene with NO₃: Hallquist et al. (1999) measured a 7 % yield (corresponding to 52.9 µg∕m³ organic aerosol) at 15 ^∘C. Nah et al. (2016b) measured a yield of 3.6 % (corresponding to 2.4 µg∕m³ organic aerosol) at room temperature. Fry et al. (2014) reported no significant aerosol growth at room temperature. This is indicative of the low temperature employed in the experiments having a significant impact on SOA yield.

Figure 3 (a) Comparison of experimental and modeling results of aerosol mass for oxidation of α-pinene with NO₃. Open black markers are experimental aerosol masses obtained using an SMPS. The blue solid line represents the best-fit model result, the blue dotted line represents an alternative model fit, and the grey dashed line corresponds to the experimental temperature profile. (b) Analysis of the oligomerization state of particle-phase products in the model according to fit 1. (c) Analysis of the occupation of volatility bins of all products according to fit 1 and during peak SOA mass. Shadings in the bar plot denote where molecules of a certain volatility bin reside: gas phase (grey) or particle phase (colored). Products in the particle phase are further distinguished as organic nitrates (green) and non-nitrated organics (orange).

After about 4 h of total experiment time, the temperature set point of the chamber enclosure is increased to 25 ^∘C, leading to a sharp and significant evaporation of organic material from aerosol particles. When the new temperature plateau is reached after 7 h, aerosol mass has decreased to 80 µg∕m³. Since evaporation has hardly slowed down by that time, heating to the new temperature set point of 42 ^∘C is initiated after 8 h of experiment time (i.e., without long waiting time at the 25 ^∘C temperature plateau) to avoid losing too much SOA mass from evaporation. After a chamber temperature of 42 ^∘C is reached after 10 h, evaporation slows down considerably and continues at a slow rate until the end of the experiment, where a minimum aerosol mass of 57 µg∕m³ is observed. With a seed mass of 37.3 µg∕m³, this corresponds to a retained SOA mass of about 20 µg∕m³ (cf. Table 2).

3.2.2 Kinetic modeling results (APN)

The kinetic model (blue solid line in Fig. 3a) shows a reasonable correlation to the experimental data. The detailed model analysis in Fig. 3b reveals that at peak SOA mass, the aerosol is composed of about 73 % of monomers, 5 % oligomers, and 22 % gas-phase dimers (Fig. 3b). These monomers mostly occupy the $C^{*}=\text{1–100}$ µg∕m³ volatility bins (Fig. 3c). Note that in Fig. 3c, a large fraction of α-pinene oxidation products occupies the $C^{*}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ µg∕m³ volatility bin, which explains the overall low SOA yield.

Upon increase in chamber temperature, evaporation of monomers in volatility bins $C^{*}=\mathrm{10}$–100 µg∕m³ and decomposition of oligomers lead to a decrease in the monomer and oligomer mass, respectively. As a result, the gas-phase dimers represent a greater fraction of the total condensed mass, and their mass fraction increases from 22 % to 74 %. Hence, the slowing of evaporation of organic material toward the end of the experiment can be attributed to the fact that the remaining organic aerosol is only comprised of gas-phase dimers ($C^{*}=\mathrm{0.01}$ µg∕m³), low-volatility monomers ($C^{*}=\mathrm{0.01}$–1 µg∕m³ volatility bins), and oligomers composed of low-volatility monomer building blocks.

Compared to the LIM experiment, peak aerosol mass is reached more quickly in the APN experiment, which is even exaggerated in the model solution. In the model, the oligomer formation rate is low at 0.124 h⁻¹, which is 2 orders of magnitude slower than determined for the LIM experiment. On the other hand, the oligomer decomposition rate is determined to be 19.0 h⁻¹, which is 2 orders of magnitude quicker than that determined for the LIM experiment and the rates reported by D’Ambro et al. (2018) for SOA from α-pinene ozonolysis. This leads to an overall lower, more labile oligomer content for the APN experiment according to the model. The higher gas-phase dimer concentration can be explained by the higher initial precursor concentration used in the APN experiment that leads to a higher momentary RO₂ concentration (cf. Fig. S11) and hence a more pronounced RO₂ + RO₂ gas-phase chemistry compared to the LIM experiment. The branching coefficient c₁ for dimer formation (cf. Fig. S3) is not included in the original MCM mechanism but was determined here from the inverse modeling to be $\mathrm{1.96}\times {\mathrm{10}}^{-\mathrm{2}}$.

The effective enthalpy of vaporization ΔH_vap of α-pinene oxidation products is determined to 81.3 kJ∕mol, which is only slightly larger than values used in the SOA models ECHAM-HAM (59 kJ∕mol; Saathoff et al., 2009), GEOS-Chem (42 kJ∕mol; Chung and Seinfeld, 2002), or GISS modelE (72.9 kJ∕mol; Tsigaridis et al., 2006).

Figure 3a also shows an alternative fit to the experimental data (fit 2). While this fit scored overall lower in our metric for model–experiment correlation, mostly due to misrepresentation of particulate organic nitrate content (pON ∕ OA, cf. Sect. 3.4), it leads to a better representation of SOA mass for the APN experiment. In fit 2, oligomer fraction is overall higher, leading to a slower increase and slower decline of SOA mass, which is more in line with experimental data. This is achieved by a much faster oligomerization rate of 9.0 h⁻¹ and a slower oligomer decomposition rate of 5.8 h⁻¹ (Table S2). The oligomerization state of SOA according to fit 2 in analogy to Fig. 3b is shown in Fig. S12, and the fit parameters for both, fits 1 and 2, are compared in Table S2. A discussion of the inability of the model to fit both SOA mass and pON ∕ OA for α-pinene at the same time can be found in Sect. 3.5.3.

3.3 Simultaneous and sequential oxidation experiments (MIX and SEQ)

In addition to oxidation experiments with single precursors, experiments are performed where α-pinene and limonene are oxidized simultaneously (MIX) or in sequence (SEQ) to investigate whether their co-existence affects growth or evaporation of SOA. In Fig. 4a (MIX) and 4b (SEQ), aerosol mass is displayed for these two scenarios alongside kinetic modeling results. The experiments are set up in a way that the produced aerosol mass is comparable in magnitude to the pure precursor experiments and both precursors contribute to the produced mass in equal parts. Table 2 lists the experimental SOA yields along with injected precursor amounts.

3.3.1 Experimental observations (MIX and SEQ)

In the MIX experiment (Fig. 4a), most of the initial increase in aerosol mass (black open markers) is rapid, and peak SOA mass is reached after about 3 h, which is comparable to the pure α-pinene oxidation experiment. The evaporation pattern upon chamber heating shows a less pronounced decrease in particle mass compared to the APN experiment but is more pronounced than observed in the LIM experiment. Overall, the mass loss during the 5 to 25 ^∘C evaporation step is more pronounced than mass loss during the 25 to 42 ^∘C step.

Figure 4 Overview of experimental and modeling results of aerosol mass for experiments with mixed monoterpene precursors. The experiments in the two panels differ in the way the precursors were added: (a) simultaneous oxidation of a mixture of α-pinene and limonene; (b) sequential oxidation of firstly α-pinene and secondly limonene with NO₃. Open black markers are experimental aerosol mass obtained using an SMPS. The colored solid and dotted lines represent model results from two different fits to the experimental data. The grey dashed line indicates the experimental temperature profile.

In the SEQ experiment (Fig. 4b), initial growth of α-pinene SOA onto the inorganic seed particles is rapid. After subsequent injection of limonene precursor, the second increase in aerosol mass is more gradual, as would be expected from the pure LIM experiment. This might be due to slow formation of oligomers but also simply because the lower amount of limonene precursor and proportionately lower injected NO₃ leads to a longer reaction time. However, the modeled reaction times for α-pinene and limonene to reach 5 % of their initial concentration after precursor injection were both about 15 min (cf. Fig. S11), which is a short time frame in comparison to the slow increase in limonene mass. The evaporation pattern in the SEQ experiment is less pronounced than the one in the MIX experiment during the 5 to 25 ^∘C temperature increase and equally marginal from 25 to 42 ^∘C.

3.3.2 Kinetic modeling results (MIX and SEQ)

The model result of the best-fit modeling scenario (fit 1, solid green and purple lines) shows fair correlation to the experimental data in the MIX experiment (Fig. 4a) but lacks in correlation in the SEQ experiment (Fig. 4b). The alternative modeling scenario (fit 2, dotted green and purple lines) shows very similar behavior. Strikingly, the mass at peak aerosol growth is overestimated by the model in both scenarios. Furthermore, initial evaporation is overestimated such that aerosol mass in the middle and late stages of the experiments agrees between the model and the experiment for the MIX experiment. Towards the end of the experiment, evaporation is further overestimated in the SEQ experiment, such that predicted aerosol mass becomes lower than the experimentally observed mass.

We note that, while peak mass does not coincide between the model and the experiment for the MIX and SEQ experiment, it is possible to obtain model fits in which this is the case. It is however not possible to match both the peak mass and the experimentally observed evaporation pattern. Slight overestimation of peak mass in the fits at hand can hence be seen as a consequence of the optimization algorithm trying to minimize the least-squares error when in reality the evaporation pattern could not be reproduced.

Figure S13 shows the time evolution of α-pinene- and limonene-derived oxidation products over time in the MIX and SEQ experiments. More α-pinene than limonene oxidation products evaporate from the particles in these model simulations, as would be expected from the pure precursor experiments. However, the fact that model–experiment correlation in the MIX and SEQ experiments is worse than in the APN and LIM experiments indicates non-linear behavior of the mixed precursor experiments. Because evaporation is overestimated by the model, especially in the SEQ experiment, effects not treated in the current model must lead to a slowing of evaporation speed in the mixed precursor experiments.

These results are similar to the findings of Boyd et al. (2017), who showed less evaporation of limonene SOA and more evaporation of β-pinene SOA in a SEQ-type experiment (β-pinene SOA condensing on preformed limonene SOA) compared to their MIX-type experiment. The study postulated a core–shell morphology of a limonene SOA core and a β-pinene SOA shell that is sustained due to incomplete mixing, though oligomerization between limonene and β-pinene oxidation products could also play a role. Here, we show in a proof of concept that oligomerization mechanics alone cannot fully explain the evaporation of monoterpene SOA mixtures. In Sect. 3.5, we will take a closer look at possible explanations.

3.4 Organic nitrate fractions

In this study, the organic nitrate fraction (pON ∕ OA) is presented as a ratio of the total mass concentration of particulate ON (which includes the organic part and nitrate part of the ON compounds) to the total mass concentration of organic aerosol (which includes both ON and non-nitrated organics) (Takeuchi and Ng, 2019). It can be inferred from AMS data using Eq. (10). In this formula, it is assumed that all organic aerosol mass is found in the organic and nitrate signals of the AMS (AMS_ORG and AMS_NO3) and all AMS nitrate is ON. When MW_pON is the average molar mass of the ON (i.e., 250 g∕mol in this study) and ${\mathrm{MW}}_{{\mathrm{NO}}_{\mathrm{3}}}$ the molar mass of the nitrate group (i.e., 62 g∕mol), the pON mass can be determined by scaling the AMS signal with the ratio of these molar masses.

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{pON}}{\mathrm{OA}}}={\displaystyle \frac{{\mathrm{AMS}}_{{\mathrm{NO}}_{\mathrm{3}}}\cdot \frac{{\mathrm{MW}}_{\mathrm{pON}}}{{\mathrm{MW}}_{{\mathrm{NO}}_{\mathrm{3}}}}}{{\mathrm{AMS}}_{{\mathrm{NO}}_{\mathrm{3}}}+{\mathrm{AMS}}_{\mathrm{ORG}}}}\approx {\displaystyle \frac{\mathrm{4.03}}{\mathrm{1}+\frac{{\mathrm{AMS}}_{\mathrm{ORG}}}{{\mathrm{AMS}}_{{\mathrm{NO}}_{\mathrm{3}}}}}}\end{array}$$

Figure 5 depicts measured and modeled values of pON ∕ OA for all four experiments. Panel a shows that in the LIM experiment, pON ∕ OA is high, with a mass ratio of about 0.8 in the particle phase, and only slightly increases over time, which is reproduced in the model. An alternative representation, showing the contribution of dinitrated, mononitrated, and non-nitrated organics to SOA mass, is shown in Fig. S14 in the Supplement and reveals that, despite the high pON ∕ OA, a significant fraction of products remains non-nitrated and the high pON ∕ OA is caused by the presence of dinitrated oxidation products. Note that the average molar mass of ON might change during the experiment, e.g., by evaporation of lower molecular weight components, which is not considered in our calculation. In the model, the slow evaporation of limonene SOA is caused by oligomer decomposition followed by evaporation of volatile monomers. The fact that pON ∕ OA is rather constant over time thus gives no evidence that decomposition rates of oligomers consisting of nitrated or non-nitrated monomeric building blocks might differ, and we use the same oligomer decomposition rate irrespective of nitration state of the product bin. Note that in the absence of oligomerization, a constant pON ∕ OA could only be obtained if nitrated and non-nitrated organics were evenly distributed across the evaporating volatility bins.

Figure 5 Experimental and modeling results of particulate organic nitrate content (pON ∕ OA) for four different types of chamber-generated SOA. (a) Only limonene, (b) only α-pinene, (c) a mixture of α-pinene and limonene, and (d) sequential oxidation of firstly α-pinene and secondly limonene. Cross markers are experimental nitration degrees inferred using a high-resolution time-of-flight aerosol mass spectrometer (HR-ToF-AMS). The colored solid and dotted lines represent results of the kinetic model. The grey dashed line indicates the experimental temperature profile.

Panel b shows pON ∕ OA in the APN experiment. The initial nitrate content is lower than in the LIM experiment with a value of about 0.45. During the first temperature increase in the APN experiment, ON content increases with the reduction in organic mass, indicating predominant evaporation of non-nitrated oxidation products. During the second evaporation step, ON content decreases, indicating predominant evaporation of nitrated oxidation product. The best-fit model run (solid blue line) captures the ON content very well. As Fig. 3b highlights, the model suggests the higher-volatility monomers to be non-nitrated and the lower-volatility monomers to be nitrated, which causes the distinct trend of pON ∕ OA. The alternative model run (dotted blue line), however, fails to capture the ON time dependence. This is due to the high oligomer content of fit 2 and due to the model not distinguishing between nitrated oligomer and non-nitrated oligomer decomposition rates. Hence, while the oligomer-heavy fit 2 shows a better correlation to α-pinene SOA mass, it fails at describing pON ∕ OA.

The measured and simulated ON contents for the experiments with multiple precursors are shown in panels c and d of Fig. 5 for the MIX and SEQ experiment, respectively. While both experiments use approximately the same concentrations of α-pinene and limonene, the measured pON ∕ OA are slightly different. Simultaneous oxidation (MIX) leads to an initial pON ∕ OA of 0.53, which is surprisingly low and closer to the value measured for pure α-pinene SOA. Sequential oxidation (SEQ) leads to an initial pON ∕ OA of 0.52 after α-pinene oxidation and increases to 0.6 after oxidation of limonene has concluded. This value in the SEQ experiment is closer to the expected value when assuming linear additivity of ON content. The unexpectedly low ON content in the MIX experiment points either towards non-linear effects in chemistry that are not captured by the model or towards uncertainties in the pON ∕ OA measurements. For the latter, there are two major sources of uncertainty. First, a default value of relative ionization efficiency (RIE) of 1.1 is used for AMS nitrate in this study (Canagaratna et al., 2007). This value is typically associated with inorganic nitrate as the RIE of nitrate derived from pON has not yet been experimentally measured to the knowledge of the authors. It is thus not clear how this value depends on chemical composition or if exposure to higher temperature may lead to variation of RIE over the course of an experiment. Second, a constant molecular weight of pON (250 g∕mol) is assumed for calculation of pON ∕ OA. However, it is possible that changes in chemical composition result in changes in the average molecular weight during an experiment. However, qualitatively, the time and temperature dependence of the ON fraction is overall captured well by the model for the mixed precursor experiments. In both cases, predominant evaporation of α-pinene oxidation products, which are the more-volatile and less-nitrated components of the mixture, leads to an increase in pON ∕ OA until the highest temperature.

The model parameters that mainly the determine pON ∕ OA are the volatility distributions of the nitrated and non-nitrated oxidation products but also the branching coefficients of the gas-phase chemical mechanism (cf. Fig. S3). The chemical mechanism presented in this study deviates from the MCM template in that it allows nitrated alkoxy radicals (R^NO) to stabilize without elimination of the nitrate function. This is realized in the model using a branching coefficient c₄ that determines the fraction of R^NO that loses its nitrate group during the conversion to a stable oxidation product. c₄ is determined to be 0 for the α-pinene system and 0.52 for the limonene system, both indicating a significant retrieval of stable organic nitrates from nitrated alkoxy radicals. A small value of c₄ stands in contrast to the findings of Kurtén et al. (2017), who ascribed the low organic nitrate yield in the oxidation of α-pinene with NO₃ to a predominant stabilization of R^NO to the volatile and non-nitrated pinonaldehyde. Note that these calculations were performed at 25 ^∘C, while α-pinene oxidation occurred at 5 °C in our experiments and model. c₄ itself is unlikely to have a positive temperature dependence, as the reaction pathway with the lower activation barrier should be even more favored at lower temperature. However, it may be possible that the fraction of alkyl radicals that undergo rearrangement (Vereecken et al., 2007) is enhanced at low temperature. The peroxy and alkoxy radicals resulting from such a rearrangement do not lose NO₂ upon stabilization. In addition, oxidation products with aldehyde moieties might be nitrated in a secondary reaction with NO₃ (Atkinson and Arey, 2003). This represents another channel of increasing pON ∕ OA and is not considered in our model. Thus, the simple gas-phase chemistry branching coefficients c₂–c₄ obtained through inverse modeling may be seen as effective parameters that represent gas-phase radical chemistry in the context of a certain experiment and volatility distribution, but their numerical values should not be evaluated in isolation.

A notable observation from modeling is that dimers from the gas-phase reaction of RO₂ + RO₂ are mainly nitrates because most RO₂ radicals originate from the reaction of alkene with NO₃ and are hence nitrated. This is especially significant for the α-pinene + NO₃ reaction system since the high momentary RO₂ radical concentrations in these experiments lead to a high estimated contribution of gas-phase dimers to SOA mass of 22 % at peak SOA mass and 74 % after heating to 42 ^∘C (cf. Fig. 3).

In summary, the experimental and modeling results in this study confirm previous studies and report a high efficiency of nitration in the reaction of monoterpenes with NO₃, with a nitrated SOA fraction larger than 50 % under most experimental conditions studies (Ng et al., 2017, and references therein). Limonene SOA shows overall higher nitration degrees than α-pinene SOA, which can be understood by the higher number of double bonds of the VOC precursor compound itself and hence more possibilities to introduce a nitrate group during oxidation. The temporal evolution of limonene SOA pON ∕ OA was constant, which can be explained with a particle phase mostly consisting of oligomers whose decomposition rates do not differ for nitrated and non-nitrated building blocks. The temporal evolution of α-pinene SOA pON ∕ OA can only be retrieved if the particle phase is predominantly comprised of monomers: sequential evaporation of nitrated and non-nitrated monomers with different vapor pressure leads to modulation of pON ∕ OA.

3.5 Deviation between the model and the experiment

From Sect. 3.3, we can conclude that while peak aerosol mass can be reconciled between the four simulated experiments with the kinetic model, the evaporation pattern in experiments MIX and SEQ cannot be brought fully into agreement with the pure precursor experiments LIM and APN. Hence, the kinetic model must lack a process that leads to resistance in evaporation in the mixed precursor scenarios compared to the pure precursor experiments. Possible mechanisms introducing such non-linearity include

1. non-linear gas-phase chemistry,

2. augmented particle-phase oligomerization chemistry, and

3. mass transfer limitations.

In general, none of these points can be fully excluded based on the results presented in this paper. However, in the following, we will go through the obtained evidence and evaluate these points to make an informed guess on how likely they are to affect aerosol formation and evaporation.

3.5.1 Gas-phase chemistry

Non-linear effects in gas-phase chemistry branching ratios could lead to a mixture of oxidation products that is more readily oxidized or dimerized and hence would show a reduced evaporation rate upon increase in chamber temperature. One possible mechanism for this is an increased yield of gas-phase dimers due to bimolecular reaction of two RO₂ radicals from different precursors, forming hetero-dimers of oxidation products. Formation of hetero-dimers is considered in the model; however, the branching ratio is assumed to be similar for limonene- and α-pinene-derived molecules, and hence self-reactions are of the same speed as cross-reactions. Berndt et al. (2018) showed that cross-reactions of two different α-pinene-derived RO₂ radicals can be faster than the respective self-reaction rates. If such an effect existed for hetero-dimers of α-pinene and limonene oxidation products, this would cause a higher dimer fraction in the product spectrum, which in turn would lead to reduced evaporation of SOA from precursor mixtures due to overall lower volatility. Since in precursor mixtures the number of RO₂ radicals is diversified, more cross-reactions will occur naturally, which would lead to more gas-phase dimers and in turn explain the slower evaporation in the MIX experiment. The SEQ experiment, however, also shows slow evaporation compared to the pure precursor experiment. Since oxidation occurred separately and cross-reactions are not enhanced by diversification of RO₂ radicals, formation of hetero-dimers in the gas phase cannot be the cause for reduced product volatility in the SEQ experiment.

Of note, any explanation for a decreased volatility of oxidation products due to gas-phase chemistry would not only change the evaporation behavior of the SOA mixture, but also likely alter the SOA yield. However, while we observe slower evaporation in mixed precursor experiments compared to single-precursor experiments, elevated SOA yields are not observed. It is hence unlikely that altered gas-phase chemistry leads to the observed reduced evaporation rates of the SOA mixtures.

3.5.2 Oligomerization

Augmented oligomerization in the particle phase is a possible explanation of reduced evaporation rates in case mixtures of oxidation products from different precursors oligomerize more readily together than the pure components in isolation. Unlike the gas-phase chemistry scenarios described above, these effects could be observed in both MIX and SEQ experiments, since particle-phase oligomerization may occur retroactively after the second oxidation step in the sequential oxidation experiment. Moreover, oligomerization of already low-volatility products would not alter SOA yields as strongly as gas-phase chemical effects would but could have a pronounced influence on evaporation rates.

In general, an augmentation effect leading to a higher oligomerization degree in mixtures could be achieved when the hetero-oligomers were formed more efficiently than a linear combination of formation rates of both homo-oligomers. A similar effect would be achieved when oxidation products of one of the two precursors were such efficient oligomer-formers that they would cause the oxidation products of the other precursors to oligomerize more readily and pull them into the oligomer phase. Therefore, during development of the model, we tested an implementation of the oligomerization scheme where formation of hetero-oligomers occurs at a combined rate using their logarithmic mean value, but first-order decomposition rates remain unaffected by the precursor type. The model solution exhibited a large discrepancy in oligomerization rates of a few orders of magnitudes, with limonene oxidation products oligomerizing quickly and readily and α-pinene oxidation products hardly oligomerizing in isolation. As a result, mixtures of oxidation products still oligomerized significantly, driven by the high individual oligomer formation rate of limonene oxidation products. Equilibrium oligomerization degree is governed by both oligomer formation and decomposition rates but is also naturally capped to a value of 100 %. Hence, in conclusion, mixing a strong oligomer former that reaches this cap in isolation with a weak oligomer former can lead to a higher combined oligomerization degree of the mixture. However, this pure theoretical result seems unphysical as it requires a very high oligomerization degree of pure limonene SOA and a very small degree of oligomerization in pure α-pinene SOA, which has not been observed in experimental studies (Faxon et al., 2018; Takeuchi and Ng, 2019).

3.5.3 Mass transfer limitations

Increased mass transfer limitations caused by high viscosity can cause a reduction of volatilization. This is due to surface concentrations of the evaporating components being depleted when the mixing timescale in the particle is longer than the evaporation timescale. Mass transfer limitation is not treated in the model runs previously shown in this study. Instead, a well-mixed bulk phase is assumed, and any resistance in evaporation is explained with oligomerization reactions. The slow evaporation of limonene SOA is hence solely caused by significant oligomerization in the model runs previously presented but could also be caused by mass transfer limitations induced by a high bulk-phase viscosity, especially if a high fraction of particle-phase oligomers would have formed that depresses mobility of molecules in the condensed phase (Baltensperger et al., 2005; D’Ambro et al., 2018). Hence, limonene SOA might exhibit a more viscous phase state than α-pinene SOA. The high viscosity caused by limonene oxidation products might in turn affect evaporation in the mixed precursor experiments and cause the observed non-linear effects.

In a first approximation, viscosities of mixtures can be assumed to be a linear combination of the individual viscosities and follow a logarithmic mixing rule (Gervasi et al., 2019). This entails that the change in the rate of mass transport between pure compounds and their mixtures can reach orders of magnitudes. This would be in line with volatilization rates observed in the mixed precursor experiments being more similar to the pure LIM experiment, which was observed in this and a previous study (Boyd et al., 2017). Notably, while evaporation steps immediately following a change in chamber temperature are overall similar between the MIX and SEQ experiments, the slope of the aerosol mass versus time curve is steeper in the MIX experiments. This might suggest that in the SEQ experiment, limonene SOA might be covering the preformed α-pinene oxidation products in a core–shell morphology and thus hampering their volatilization.

To test the effect of impeded bulk diffusivity on the evaporation of SOA, we perform a sensitivity study in which we increase viscosity in the model to evaluate whether the evaporation rates in the MIX and SEQ experiments can be brought into agreement with observations. We use the best case fitting scenario shown in Fig. 2 and raise the viscosity in the simulation to 1×10⁷ Pa s (Fig. 6). This viscosity is in the typical range for SOA under dry conditions and falls into the semi-solid phase state region (Koop et al., 2011; Shiraiwa et al., 2011; Abramson et al., 2013; Zhang et al., 2015; Grayson et al., 2016; Gervasi et al., 2019). Using the Stokes–Einstein relation (Einstein, 1905) and an effective molecular radius of 2 nm, this viscosity corresponds to a bulk diffusion coefficient of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{16}}$ cm²∕s at 298 K. The effective radius is approximated from geometric considerations assuming spherical molecular shape, a molar mass of 250 g∕mol, and density of 1.55 g∕cm³. The temperature dependence of this diffusion coefficient is approximated with a constant activation enthalpy of diffusion ΔH_dif=50 kJ∕mol according to Eq. (11).

$$\begin{array}{}\text{(11)}& D_{\mathrm{b}}\left(T\right)=D_{\mathrm{b}}\left(\mathrm{298}\mathrm{K}\right)\cdot \exp {\displaystyle \frac{-\mathrm{\Delta }{H}_{\mathrm{dif}}}{R\left(\frac{\mathrm{1}}{T}-\frac{\mathrm{1}}{\mathrm{298}}\right)}}\end{array}$$

Figure 6 shows that a reduced bulk diffusivity leads to a reduction in peak SOA mass and a shallower evaporation profile. The reduction in SOA yield is caused by a reduction in particle-phase oligomerization: monomer building blocks cannot freely diffuse into the particle but rather partition to the near-surface layers predominantly. This effectively lowers their uptake coefficient, which stands in competition to uptake by the chamber walls. Despite the lower oligomer fraction in these calculations, evaporation is significantly slowed down compared to the well-mixed case. This model result insinuates that the co-presence of limonene SOA and α-pinene SOA might strongly reduce the mobility of α-pinene oxidation products so that the fast evaporation of α-pinene oxidation products observed in the pure α-pinene oxidation experiment does not take place.

Figure 6 Sensitivity study on the influence of viscosity on model simulation results based on the best case fitting scenario in the (a) MIX and (b) SEQ experiments. Model simulations were performed for the default well-mixed case (solid lines) and at a diffusivity coefficient of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{16}}$ cm²/s, which corresponds to a bulk viscosity of 1×10⁷ Pa s according to the Stokes–Einstein relation and falls within the semi-solid phase state range.

The outcome of this 1-D sensitivity study has to be treated with caution since introducing slow diffusion of oxidation products also causes a shift in all other optimization parameters. For example, with the default parameter set, the slow evaporation of limonene SOA in the model is purely attributed to oligomer formation. The slow-down in evaporation in this sensitivity study hence suggests that the high oligomerization degree suggested by the model for limonene SOA in the previous best-fit solutions might have been overestimated. In fact, a particularly high oligomer content was not observed for limonene SOA from oxidation with NO₃ in measurements using FIGAERO-CIMS (Faxon et al., 2018). Distinction of these two effects (oligomerization vs. mass transfer limitation of slow evaporation) could be possible with the model and the MCGA but is not attempted in this study due to the prohibitive computational cost of model calculations at low diffusivities and will be subject of future studies. Furthermore, the slow growth of particles in the pure limonene oxidation experiment is attributed in the well-mixed model (Sect. 3.1.2) to dissolution and subsequent oligomerization of high-volatility compounds in the fourth to sixth volatility bins. In a viscous particle model, the volatility of these bins might shift down, while maintaining the same particle growth velocity.

We have seen in Sect. 3.2.2 and 3.4 that experimental α-pinene SOA mass can only be matched with a model run that ascribes a high oligomer content to α-pinene SOA (fit 2), which is typically not reported in the literature (Romonosky et al., 2017). In return, a high oligomer content cannot describe the time evolution of α-pinene pON ∕ OA properly. Hence if, hypothetically, a semi-solid phase state of α-pinene SOA were to slow down evaporation so that SOA evaporation is reconciled between the model and the experiment with a particle phase mostly comprised of oligomers, the distinct temporal evolution of pON ∕ OA could still be matched.

Taken together, it is possible that increased mass transfer limitation led to the observed reduced evaporation rates of the SOA mixtures as postulated in Boyd et al. (2017). However, there are still large uncertainties and a high computational expense associated with a model treatment of highly viscous SOA systems. While frameworks for the determination of viscosity of mixtures have recently been developed (Gervasi et al., 2019), these rely on structural information about individual compounds. Furthermore, while the Stokes–Einstein relation seems to hold for similar systems at viscosities of up to 10⁴ Pa s (Ullmann et al., 2019), it is not clear whether it also holds for viscosities of  10⁷ Pa s derived in this study (Evoy et al., 2019).

Additionally, treatment of slow particle-phase diffusion requires many model layers to describe the steep concentrations gradients arising at the particle surface upon evaporation. In combination with the multitude of tracked species in the particle phase, computational costs quickly reach unfeasible ranges. Ideally, the spatial resolution model layers would have to be generated upon model runtime by an algorithm that detects steep concentration gradients. This detailed description will be presented in a forthcoming publication.

4 Conclusions and outlook

In this study, an inverse modeling approach is utilized alongside laboratory chamber experiments to gain insights into the molecular-level processes which occur during the formation and evaporation of SOA from the oxidation of α-pinene, limonene, and mixtures of both precursors with NO₃. We find α-pinene SOA to form and evaporate rather quickly and limonene SOA to form and evaporate more slowly. Both SOA types, however, show retardation in evaporation compared to instantaneous equilibration of a specified volatility basis set, which can in part be explained by the presence of particle-phase oligomers. A mixed oxidation and a sequential oxidation of both precursors show the expected linear additivity of SOA yields but a non-linear reduction in evaporation behavior, which could not be fully explained without including diffusion limitations in the particle phase into the model calculations. Since it is computationally difficult to treat the effects of slow mass transport fully in these models, this paper focuses first on oligomerization and tries to make cases for and against oligomerization as the sole cause for our observations.

The oxidation products of both SOA types are found to be heavily nitrated. The results highlight the significance of NO₃ as oxidant in SOA formation and the importance of ON as products of monoterpene oxidation. The study finds evidence for non-equilibrium partitioning caused by slow particle-phase chemistry and slow diffusion, which is currently not considered in global models and may lead to underestimation of SOA persistence and hence underestimated global SOA burdens.

The modeling approach applied in this study comprises a combination of the kinetic model based on KM-GAP (Shiraiwa et al., 2012) with the automated global optimization suite MCGA (Berkemeier et al., 2017) and details the full chemistry and physics of SOA particle growth and shrinkage. The underlying SOA formation and evaporation mechanism uses a simplified and lumped version of the Master Chemical Mechanism (MCM; Jenkin et al., 2003; Saunders et al., 2003; Berkemeier et al., 2016), extends it with a reversible particle-phase oligomerization and gas-phase dimerization scheme, and treats gas–particle partitioning with a volatility basis set approach (Donahue et al., 2006, 2011) for each product bin. The study focuses on NO₃ oxidation of monoterpenes and their mixtures, but the model framework can be ported to other chemical systems. The depth resolution capabilities of the model allow for a sensitivity study of the influence of particle phase state on the evaporation of these particles. A full treatment of composition-dependent, depth-resolved viscosity as a global optimization parameter is ultimately needed to disentangle the interactions of particle-phase diffusion and particle-phase chemistry. Due to the computational expense of finely resolved computational layers and the general uncertainty in the physical and chemical parameters, this will be subject of follow-up studies. In such studies, offline analysis of the oligomerization degree of SOA material can help to constrain oligomerization and oligomer decomposition rates, and thermodynamic models can be used to provide estimates for composition dependence of viscosities and diffusivities (DeRieux et al., 2018; Gervasi et al., 2019).

While there is significance to the general conclusions drawn from the model analysis, the individual model parameters that are returned by the inverse modeling approach must be treated with caution and evaluated in the context of the model and experimental data that are employed. Given the large number of fitting parameters and the limited number of experimental data sets, it cannot be ensured that a true and correct global minimum is obtained in this isolated case study. With a simplified multi-parameter model and experimental data sets that are aggregate observables and subject to uncertainty, the concept of a single global minimum and multiple local minima on the optimization hypersurface becomes blurred, and several extended areas on the optimization hypersurface can exhibit a minimal function value. For example, Fig. S4 shows an estimate of the uncertainty in the volatility distributions obtained in this study. The error bars in Fig. S4 are standard deviations of individual re-fits of volatility distributions that all lead to a similar calculation outcome and hence quantify their uniqueness (or lack thereof). Figures S1, S8, S9, and S10 show sensitivity case studies of very influential model parameters, and Table 1 shows a local sensitivity analysis of the remaining input parameters, which gives an impression of their range within a single model fit. However, the true parameter ranges can be much larger than apparent from these local sensitivity analyses. For example, changes in branching ratios in the gas-phase chemical mechanism can in principle be offset with changes in the oxidation products' volatility distributions, thus forming a co-dependent parameter subset. The uniqueness of the obtained parameter set can be enhanced by inclusion of more experimental data at different conditions or by a priori determination of model parameters such as measurements of volatility distributions, oligomerization degrees, or particle viscosities, which will be an imperative task in follow-up studies.

However, despite the remaining uncertainties in derived model parameters, the modeling suite presented here constitutes a step forward in the computational, data-driven evaluation of SOA formation with kinetic models. In this work, only a small set of laboratory chamber data is utilized for optimization as proof of concept. We postulate that, by reconciling and cross-comparing large sets of experimental data, we will be able to significantly enhance our understanding of SOA and close the gap between our expanding theoretical knowledge about the detailed gas-phase chemistry, gas–particle partitioning, particle phase state of SOA, and the application of this knowledge in chemical transport models.