Exploration of the influence of environmental conditions on secondary organic aerosol formation and organic species properties using explicit simulations: development of the VBS-GECKO parameterization

Atmospheric chambers have been widely used to study secondary organic aerosol (SOA) properties and formation from various precursors under different controlled environmental conditions and to develop parameterization to represent SOA formation in chemical transport models (CTMs). Chamber experiments are however limited in number, performed under conditions that differ from the atmosphere and can be subject to potential artefacts from chamber walls. Here, the Generator for Explicit Chemistry and Kinetics of Organics in the Atmosphere (GECKO-A) modelling tool has been used in a box model under various environmental conditions to (i) explore the sensitivity of SOA formation and properties to changes on physical and chemical conditions and (ii) develop a volatility basis set (VBS)-type parameterization. The set of parent hydrocarbons includes n-alkanes and 1-alkenes with 10, 14, 18, 22 and 26 carbon atoms, α-pinene, β-pinene and limonene, benzene, toluene, o-xylene, m-xylene and p-xylene. Simulated SOA yields and their dependences on the precursor structure, organic aerosol load, temperature and NO_x levels are consistent with the literature. GECKO-A was used to explore the distribution of molar mass, vaporization enthalpy, OH reaction rate and Henry's law coefficient of the millions of secondary organic compounds formed during the oxidation of the different precursors and under various conditions. From these explicit simulations, a VBS-GECKO parameterization designed to be implemented in 3-D air quality models has been tuned to represent SOA formation from the 18 precursors using GECKO-A as a reference. In evaluating the ability of VBS-GECKO to capture the temporal evolution of SOA mass, the mean relative error is less than 20 % compared to GECKO-A. The optimization procedure has been automated to facilitate the update of the VBS-GECKO on the basis of the future GECKO-A versions, its extension to other precursors and/or its modification to carry additional information.

1 Introduction

Fine particulate matter impacts visibility (e.g. Han et al., 2012), human health (e.g. Lim et al., 2012; Malley et al., 2017) and climate (e.g. Boucher et al., 2013). A large fraction of fine particles is organic, representing between 20 % and 90 % of the total mass (e.g. Jimenez et al., 2009). This organic fraction can be directly emitted into the atmosphere (primary organic aerosol, POA) or formed by gas–particle partitioning of low-volatility species produced during the oxidation of gaseous organic compounds (secondary organic aerosol, SOA) (e.g. Mader et al., 1952; O'Brien et al., 1975; Grosjean, 1992). The secondary origin of organic aerosol (OA) dominates the primary fraction in most environments (e.g. Gelencsér et al., 2007). SOA is composed of a large number of compounds having various chemical structures and spanning therefore a wide range of physicochemical properties (reactivity, volatility, molar mass, solubility, optical absorption). The amount of SOA formed from a gaseous precursor depends mainly on its structure (carbon chain length, degree of unsaturation, number and type of functional groups, etc.) and on environmental conditions that influence (i) the concentration and the structure of organic compounds produced during gas-phase oxidation (oxidant concentrations and NO_x levels, temperature, photolysis, humidity, etc.) and (ii) their partitioning between the gas and the condensed phase (pre-existing OA mass (C_oa), temperature, etc.) (e.g. Kroll and Seinfeld, 2008; Ng et al., 2017; Shrivastava et al., 2017).

SOA formation is represented in chemical transport models (CTMs) using simplified parameterizations. These parameterizations describe the formation of secondary organic species from the oxidation of gaseous precursors and their consumptions through the use of surrogate species having different volatilities (e.g. Odum Jay et al., 1996; Pun et al., 2002, 2003, 2006; Tsigaridis and Kanakidou, 2003; Couvidat et al., 2012; Couvidat and Sartelet, 2015), sometimes discretized into decadal volatility bins (e.g. Robinson et al., 2007; Donahue et al., 2006, 2012; Jimenez et al., 2009; Zhang and Seinfeld, 2013). The gas–particle partitioning of a surrogate species or a given volatility bin is usually described considering an absorption process according to the theory described by Pankow (1994). Comparisons with field observations show that current CTMs have difficulties in reproducing the spatial and temporal variations of PM_2.5 mass concentrations (e.g. de Gouw et al., 2005; Johnson et al., 2006; Volkamer et al., 2006; Heald et al., 2005; Hodzic et al., 2010; Shrivastava et al., 2011; Couvidat et al., 2012; Solazzo et al., 2012; Im et al., 2015), mainly because of the difficulties to represent the secondary organic fraction of the aerosol (Solazzo et al., 2012).

SOA parameterizations have typically been developed and optimized based on atmospheric chamber data. Atmospheric chambers have been widely used to study SOA formation from various precursors under different controlled environmental conditions. Experiments, being limited in number, do not cover the diversity of atmospheric chemical and physical conditions that influence SOA formation. Smog chamber experiments are also usually performed under conditions that differ from the atmosphere (high level of oxidants and precursors, different light spectrum and intensity, low humidity, short times of ageing, etc.). In addition, SOA formation experiments can be subject to potential artefacts from chamber wall surfaces, such as heterogeneous reactions or aerosol and vapour wall losses (e.g. McMurry and Grosjean, 1985; Matsunaga and Ziemann, 2010; La et al., 2016).

The objectives of this study are (i) to explore the influence of environmental conditions on SOA formation and the distribution of physicochemical properties (volatility, enthalpy of vaporization, molar mass, OH rate constant and water solubility) of the evolving molecular mixture and (ii) to optimize a volatility basis set (VBS)-type parameterization (e.g. Donahue et al., 2006) on the basis of a deterministic and explicit modelling of SOA formation. Deterministic and explicit simulations allow (i) taking into account the influence of various environmental conditions encountered in the atmosphere, (ii) simulating SOA formation under representative concentration regimes and timescales and (iii) accessing detailed properties of organic compounds in both gas and condensed phases. Nearly explicit representation of SOA formation under various environmental conditions is possible with the Generator for Explicit Chemistry and Kinetics of Organics in the Atmosphere (GECKO-A) (Aumont et al., 2005; Camredon et al., 2007). The environmental scenarios and the explicit representation of SOA formation are described in Sect. 2. Explicit simulations are used to explore the influence of environmental conditions on SOA formation and the distribution of the physicochemical properties of species produced during organic compound oxidation in Sect. 3. A VBS-type parameterization is optimized and evaluated using explicit simulations as a benchmark in Sect. 4.

2 Explicit simulations of SOA formation

2.1 The environmental scenarios

Scenarios were developed in a box model to represent the variability of SOA formation within the range of environmental conditions encountered in the atmosphere (temperature, relative humidity, solar zenith angle, aerosol mass concentration and chemical concentrations). Scenarios were built with various NO_x concentrations to cover chemical regimes from low NO_x to high NO_x conditions (Table 1).

Table 1 Constant environmental conditions used in the various scenarios. For each scenario, the concentrations of OH, HO₂, NO and NO₂ are at the stationary state induced by the constant conditions. Values in bold represent the conditions used in the reference scenario.

To examine the sensitivity of SOA formation to each parameter, chemical and physical conditions were fixed in each scenario. Physical parameters set in the simulations were temperature (270 and 298 K), relative humidity (70 %), zenith angle (30, 50 and 70^∘) and C_oa (0.1, 1 and 10 µg m⁻³). The chemistry considered was the reactivity of inorganic species (i.e. the ${\mathrm{O}}_x/{\mathrm{NO}}_x/{\mathrm{HO}}_x$ chemistry) and the oxidation of CO and CH₄. The influence of non-methane volatile organic compound (VOC) oxidation on the HO_x and NO_x cycles was included using a surrogate species. The OH lifetime with respect to reaction with the surrogate VOC was set to 1 s, in the range of what is observed in situ (Lou et al., 2010; Hansen et al., 2014). The oxidation of this surrogate VOC with OH leads to the formation of a surrogate peroxy radical, with a chemistry assumed to be similar to CH₃O₂. Dry deposition was considered for O₃, H₂O₂ and HNO₃ (with deposition rates of $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{3}}$ s⁻¹, respectively). Chemical concentrations were fixed for CO (120 ppb), methane (1750 ppb), O₃ (10, 40 and 70 ppb) and NO_x (ranging from 10⁻² to 10² ppb). Photolysis frequencies were computed for each fixed zenith angle at the surface using the Tropospheric Ultraviolet and Visible (TUV) radiation model (Madronich and Flocke, 1999) with rules for organic chromophores described in Aumont et al. (2005). Each simulation was run in a box model until the OH, HO₂, NO and NO₂ stationary state concentrations were reached. These stationary states were used as inputs for our SOA formation scenarios. In each scenario, the initial concentration of a given SOA precursor is set to a low enough value (10 pptC) to not affect significantly both the stationary state conditions of the simulations and the prescribed amount of pre-existing aerosol (C_oa).

Figure 1 Simulated stationary state concentration of OH and HO₂ and corresponding reaction rate ratio (RRR) as a function of NO_x concentrations for the various scenarios. The blue points are the simulated values for the scenario performed at a temperature of 298 K, a zenith angle of 50^∘ and 40 ppb of ozone. The grey zone corresponds to the range of values simulated in all the scenarios.

Figure 1 shows the simulated stationary state concentrations of OH and HO₂ as a function of NO_x concentrations for the various scenarios. Simulated evolutions and concentrations of OH and HO₂ with NO_x are in agreement with the observations (e.g. Stone et al., 2012). The distribution of secondary species depends on the evolution of RO₂ radicals. The branching ratio of RO₂ reacting with NO, called hereafter RRR (for reaction rate ratio), is often used to characterize the chemical environment (e.g. Lane et al., 2008). In the scenarios, the loss of RO₂ with CH₃O₂ (the major RO₂ radical) represents less than 10 %. The RRR was therefore calculated here without considering the RO₂ + RO₂ reactions as

$$\begin{array}{}\text{(1)}& \mathrm{RRR}={\displaystyle \frac{k_{{\mathrm{RO}}_{\mathrm{2}}+\mathrm{NO}}\left[\mathrm{NO}\right]}{k_{{\mathrm{RO}}_{\mathrm{2}}+\mathrm{NO}}\left[\mathrm{NO}\right]+k_{{\mathrm{RO}}_{\mathrm{2}}+{\mathrm{HO}}_{\mathrm{2}}}\left[{\mathrm{HO}}_{\mathrm{2}}\right]}},\end{array}$$

where $k_{{\mathrm{RO}}_{\mathrm{2}}+\mathrm{NO}}$ and $k_{{\mathrm{RO}}_{\mathrm{2}}+{\mathrm{HO}}_{\mathrm{2}}}$ are the rate constants for the reactions of the peroxy radicals with NO and HO₂, respectively, and [NO] and [HO₂] the concentrations of the different radicals. Figure 1 also presents the variation of RRR with NO_x. Rate constants of RO₂ with NO and HO₂ used to compute RRR are calculated at 298 K: $k_{{\mathrm{RO}}_{\mathrm{2}}+\mathrm{NO}}$ = $\mathrm{9.0}\times {\mathrm{10}}^{-\mathrm{12}}$ cm³ molec⁻¹ s⁻¹ (Jenkin et al., 1997) and $k_{{\mathrm{RO}}_{\mathrm{2}}+{\mathrm{HO}}_{\mathrm{2}}}$ = $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{11}}$ cm³ molec⁻¹ s⁻¹ (e.g. Boyd et al., 2003, assuming a large carbon skeleton for RO₂). Figure 1 shows that most of the RO₂ reacts with NO for NO_x concentrations higher than 1 ppb. The NO_x and HO_x concentrations used in the various scenarios are consistent with typical chemical characteristics of low- to high-NO_x environments. These stationary state scenarios are representative of the evolution of the tropospheric ${\mathrm{O}}_x/{\mathrm{NO}}_x/{\mathrm{HO}}_x$ chemical system. They were used here to explore explicitly the formation of SOA under the ranges of NO_x levels, oxidant concentrations, temperature, solar radiation and pre-existing OA mass encountered in the troposphere.

2.2 The explicit representation of SOA formation

Aromatic and long-chain aliphatic hydrocarbons emitted by anthropogenic activities, and isoprene and terpenic compounds emitted by biogenic sources, are known to be major SOA precursors (e.g. Carlton et al., 2009; Kroll and Seinfeld, 2008; Nordin et al., 2013; Tsigaridis et al., 2014; Zhao et al., 2016). Explicit oxidation chemical schemes were generated with the GECKO-A modelling tool for five linear alkanes (decane, tetradecane, octadecane, docosane and hexacosane), five linear 1-alkenes (decene, tetradecene, octadecene, docosene and hexacosene), five aromatic compounds (benzene, toluene, m-xylene, o-xylene and p-xylene) and three terpenic compounds (α-pinene, β-pinene and limonene). The GECKO-A modelling tool does not handle the polycyclic aromatic structures, and the current oxidation protocol is not up to date to represent SOA formation from isoprene oxidation. SOA formation from polycyclic aromatic hydrocarbons and isoprene oxidation has thus not been studied here.

The GECKO-A tool is a computer programme that automatically writes chemical schemes on the basis of a prescribed protocol. Elementary data are taken from laboratory studies if available and structure/property relationships if not. The protocol implemented in GECKO-A is described by Aumont et al. (2005), with chemistry updates for gas-phase oxidation performed by Valorso et al. (2011) and Aumont et al. (2012, 2013). The evolution of alkoxy radicals through decomposition and isomerization is estimated using the structure–property relationship developed by Vereecken et al. (2009, 2010), as described by La et al. (2016). In the current version, GECKO-A generates explicit chemical schemes for non-aromatic species and uses the Master Chemical Mechanism (MCM, version 3.3.1) to describe the chemistry of species including an aromatic structure (Jenkin et al., 2003; Bloss et al., 2005).

The explicit description of the oxidation of organic precursors involves millions of secondary species and reactions. Explicit chemical schemes cannot be solved for precursors having more than eight atoms of carbon, even in a box model. Simplifications are therefore needed. In this study, the gaseous reactivity of species having a saturation vapour pressure lower than 10⁻¹³ atm is not considered as these species are dominantly in the condensed phase. Lumping of position isomers is also allowed if the production yield of a species is lower than 10⁻³ (Valorso et al., 2011). RO₂ mainly reacts with NO, HO₂ and other RO₂. In our scenario, the RO₂ family is dominated by CH₃O₂. Among the self RO₂ reactions, only the RO₂+CH₃O₂ are considered when chemical schemes are generated. Furthermore, writing of the oxidation scheme is stopped after formation of the 15th-generation stable products.

The gas–particle partitioning of each stable organic compound was implemented as described in Camredon et al. (2007). Phase partitioning is described by an absorption process following Raoult's law (e.g. Pankow, 1994):

$$\begin{array}{}\text{(2)}& C_{\mathrm{p}}^i={{\displaystyle \frac{RT{C}_{\mathrm{OA}}}{{\mathit{\gamma }}_{i}P}}}^{\mathrm{sat}}{\stackrel{\mathrm{‾}}{M}}_{\mathrm{OA}}\times C_{\mathrm{g}}^i,\end{array}$$

where $C_{\mathrm{p}}^i$ and $C_{\mathrm{g}}^i$ are the concentration of the species i in the particulate and the gas phase, respectively, R is the gas constant, T the temperature, ${\stackrel{\mathrm{‾}}{M}}_{\mathrm{OA}}$ the mean molar weight of the OA (250 g mol⁻¹), $P_i^{\mathrm{sat}}$ the saturation vapour pressure of the species i and γ_i the activity coefficient of the species i in the condensed phase. The condensed phase is considered homogeneous, ideal (γ_i=1), inert and at equilibrium with the gas phase. Saturation vapour pressure of each organic compound is estimated using the Nannoolal method (Nannoolal et al., 2004, 2008). Equation (2) is solved for the set of organic species using a simple iterative method (Pankow, 1994). For the gas-phase oxidation, time integration is solved using the two-step solver (Verwer and VanLoon, 1994; Verwer et al., 1996).

Figure 2 Temporal evolution of the simulated carbon distribution between the precursor, CO and CO₂, and gaseous and condensed organic compounds. Simulated results are shown for o-xylene, α-pinene and n-octadecane for the scenario performed at a temperature of 298 K, a zenith angle of 50^∘, a C_oa of 1 µg m⁻³, 40 ppb of ozone and a RRR of 90 %.

3 Exploration of SOA formation and organic species properties

3.1 Carbon distribution and SOA formation

Figure 2 presents the temporal evolution of the carbon during the oxidation of a precursor in the scenario performed at a temperature of 298 K, a zenith angle of 50^∘, a C_oa of 1 µg m⁻³, 40 ppb of ozone and a RRR of 90 %. This scenario is called hereafter the “reference” scenario. Examples are shown for o-xylene, α-pinene and n-octadecane oxidation. Under these average conditions, the precursors are oxidized in 10, 6 and 2.5 h for o-xylene, n-octadecane and α-pinene, respectively. The carbon is first dominantly in the form of secondary gaseous organic compounds (around 55 %, 85 % and 40 % at the gaseous maximum for o-xylene, α-pinene and n-octadecane, respectively). For o-xylene, these gas-phase intermediates are mainly transformed into CO and CO₂, with less than 1 % of the carbon present as secondary OA. After 5 days of oxidation, the carbon mainly ends in the form of CO and CO₂ (∼85 %), the remaining fraction being mostly gas-phase organic carbon. The complete oxidation of α-pinene into CO and CO₂ takes longer than for o-xylene. During the oxidation, up to 7 % of the carbon is in SOA. After 5 days, the carbon initially present as α-pinene mainly ends in the form of CO and CO₂ (∼85 %); the remaining fraction of the carbon is gaseous at 10 % and condensed at 5 %. For n-octadecane, the full oxidation of the carbon into CO and CO₂ is slow. Gas-phase intermediates produce dominantly low-volatility species which transfer to the condensed phase. A maximum carbon fraction of around 70 % is in the particle phase after 10 h of oxidation. The carbon fraction in SOA is slowly transferred back to the gas phase when gaseous oxidation proceeds, leading ultimately to the formation of CO and CO₂. The oxidation of the carbon into CO and CO₂ occurs only by gas-phase oxidation here as the particle phase is represented as inert in this model configuration. After 5 days of oxidation, the carbon is mainly in the aerosol (∼62 %) and in the form of CO and CO₂ (∼38 %).

Figure 3 (a) Dependence of Y_max on the precursor's carbon number for the scenario performed at a temperature of 298 K, a zenith angle of 50^∘, a C_oa of 1 µg m⁻³, 40 ppb of ozone and a RRR of 90 %. (b) Dependence of $Y_{\text{max}}^{\mathrm{270}\mathrm{K}}/Y_{\text{max}}^{\text{ref}}$ on the precursor's carbon number. (c) Dependence of $Y_{\text{max}}^{\mathrm{0.1}\mathrm{\mu }\mathrm{g}\mathrm{m}-\mathrm{3}}/Y_{\text{max}}^{\text{ref}}$ (open symbols) and $Y_{\text{max}}^{\mathrm{10}\mathrm{\mu }\mathrm{g}\mathrm{m}-\mathrm{3}}/Y_{\text{max}}^{\text{ref}}$ (filled symbols) on the precursor's carbon number.

3.2 Influence of environmental parameters on SOA formation

Figure 3a shows the evolution of the maximum SOA yield (Y_max, calculated as the mass ratio of the produced SOA to the reacted precursor quantity) as a function of the number of carbons in the parent compound. Results are shown for the reference scenario; the maximum SOA yield is noted $Y_{\text{max}}^{\text{ref}}$. For a given 1-alkane or 1-alkene series, simulated SOA yield increases with the length of the precursor's carbon skeleton. This behaviour, as well as the simulated yields, is consistent with observations from chamber experiments during n-alkane and 1-alkene oxidation at ambient temperature and high NO_x (e.g. Lim and Ziemann, 2009; Matsunaga et al., 2009). Y_max reaches a plateau of about 1.2 for precursor chain lengths higher than 18 carbons, independently of their carbon number. These precursors lead to the formation of first-generation compounds being dominantly in the condensed phase (e.g. Aumont et al., 2012). For aromatic compounds, simulated Y_max decreases with increasing number of methyl groups. This Y_max evolution with the aromatic precursor's structure is consistent with experimental trends observed in chambers (e.g. Ng et al., 2007; Li et al., 2016). For terpenes, simulated Y_max increases in the following precursor's order: α-pinene, β-pinene and limonene, in agreement with observed behaviours (e.g. Zhao et al., 2015).

The influence of temperature and pre-existing OA on Y_max is presented in Fig. 3b and c. Results are shown for a scenario in which one environmental parameter is changed from the reference scenario: the temperature from 298 to 270 K or the C_oa from 1 to 0.1 or 10 µg m⁻³ (the maximum SOA yields are noted hereafter $Y_{\text{max}}^{\mathrm{270}\mathrm{K}}$, $Y_{\text{max}}^{\mathrm{0.1}\mathrm{\mu }\mathrm{g}\mathrm{m}-\mathrm{3}}$ and $Y_{\text{max}}^{\mathrm{10}\mathrm{\mu }\mathrm{g}\mathrm{m}-\mathrm{3}}$, respectively). As expected, simulation results show that aerosol yields increase when C_OA increases (e.g. Odum et al., 1996) and when the temperature decreases (e.g. Takekawa et al., 2003). In both cases, the highest sensitivity is simulated for species bearing between 6 and 14 carbon atoms. For these precursors, the secondary organic species contributing to SOA are semi-volatile and their partitioning depends largely on the temperature and C_oa. A low sensitivity of Y_max is simulated for the species with the longest carbon skeleton (i.e. with n_C>18). These low-volatility precursors rapidly form low-volatility compounds that partition dominantly to the condensed phase, regardless of environmental conditions.

Figure 4 Dependence of Y_max on RRR (here in percent) for n-alkanes (a), 1-alkene (b), terpenes (c) and aromatic species (d) for the scenarios performed at a temperature of 298 K, a zenith angle of 50^∘, a C_oa of 1 µg m⁻³ and 40 ppb of ozone. Points are connected with dashed lines to help in seeing trends of Y_max with RRR.

The influence of NO_x on Y_max is presented in Fig. 4. Results are shown for the reference scenario (RRR of 90 %) and scenarios into which the RRR is changed from the reference scenario to 0 %, 10 %, 50 % and 100 %. Simulated results show that SOA yields generally decrease with the increase of NO_x. These trends are in agreement with observations (e.g. Ng et al., 2007; Donahue et al., 2005). Again, a weak sensitivity to NO_x is simulated for largest chain length species (i.e. with n_C>18). These larger chain length compounds produce first oxidation products able to partition dominantly in the condensed phase at all NO_x levels. The explicit simulations reproduce the observed tendencies of SOA yields with the precursor's structure, NO_x levels, temperature and pre-existing OA mass. These agreements support their use to explore the physicochemical properties of organic compounds and as a reference for the development of reduced chemical schemes for SOA formation.

Figure 5 Distribution of organic species (gaseous and particulate summed) in terms of volatility ($P_{\mathrm{298}\mathrm{K}}^{\text{sat}})$ and (a) molar mass (Mw), (b) vaporization enthalpy (ΔH_vap), (c) OH reactivity (k_OH) and (d) solubility ($H_{\mathrm{298}\mathrm{K}}^{\text{eff}})$. Simulated results are shown for o-xylene, α-pinene and n-octadecane after 6 h of oxidation for the scenario performed at a temperature of 298 K, a zenith angle of 50^∘, a C_oa of 1 µg m⁻³, 40 ppb of ozone and a RRR of 90 %. Each bubble is a species and the volume of the bubble is proportional to the carbon ratio represented by the species. The bubbles are coloured by volatility classes (see Table 2 for the bin bounds). The distribution of the properties weighted by the concentration in carbon atoms of the species in each volatility class at 6 h is represented as box plots showing the 25th, 50th and 75th percentiles, and the minimal and maximal values. Red points are the mean properties weighted by the concentration in carbon atoms of the species in each volatility class.

3.3 Physicochemical properties of organic compounds

The explicit simulations were used to explore the distribution of organic species physicochemical properties that influence OA formation, in particular

1. the molar mass (Mw), needed to convert the aerosol concentration of the species from molecular to mass units;

2. the volatility (i.e. saturation vapour pressure) that influences directly the gas–aerosol partitioning of a species. The saturation vapour pressure at a given temperature is calculated in GECKO-A using the Clausius–Clapeyron expression with P^sat at 298 K and vaporization enthalpy (ΔH_vap) estimated from the Nannoolal et al. method (2004, 2008);

3. the solubility (i.e. Henry's law coefficient) that influences the dry deposition of organic gases but also potentially its gas–aerosol partitioning for water-soluble species. The effective Henry's law coefficient (H^eff) at 298 K is calculated in GECKO-A with the GROHME method (Raventos-Duran et al., 2010);

4. the reaction rate with OH (k_OH) that drives the chemical lifetime of organic species. The OH reaction rate at 298 K of each organic compound is estimated in GECKO-A using the Kwok and Atkinson (1995) structure activity relationship with updates described by Aumont et al. (2012) if no experimental data are available.

Figure 5 presents the simulated distribution of (a) molar mass, (b) vaporization enthalpy, (c) OH reactivity and (d) solubility ($H_{\mathrm{298}\mathrm{K}}^{\text{eff}})$ as a function of the volatility ($P_{\mathrm{298}\mathrm{K}}^{\text{sat}})$ of organic species. Results are presented for a time step close to the maximum of secondary organic compound, after 6 h of oxidation, for the reference scenario. Examples are shown for o-xylene, α-pinene and n-octadecane oxidation. Thousands of secondary organic species spanning several orders of magnitude in properties are formed during the oxidation of a given precursor. The volatility distributions range from 10⁻¹⁵ to 1 atm, independently of the precursor structure. Molar mass increases with the decrease of the volatility of organic compounds, varying from around 50 to 400 g mol⁻¹ for the highest and lowest volatile species. Simulated results show that species formed at high NO_x have a molar mass higher than at low NO_x (see Fig. S1 in the Supplement). Explicit simulation shows a ΔH_vap increase of 12.5 kJ mol⁻¹ per logarithm unit decrease of $P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$, independently of NO_x levels (see Fig. S1). Rate constants of organic species with OH at 298 K vary between 10⁻¹³ and 10⁻¹⁰ molec⁻¹ cm³ s⁻¹. Simulations do not show a clear trend of k_OH with the volatility of the species, neither with the level of NO_x, k_OH being scattered over 3 orders of magnitude for a given P^sat and RRR. The effective Henry's law constant increases with the decrease of organic compound volatility, varying between around 1 and 10¹⁶ mol L⁻¹ atm⁻¹ for the highest and lowest volatile species. NO_x levels do not have a significant influence on the H^eff distribution of species with the volatility (see Fig. S1).

4 Development and evaluation of the VBS-GECKO parameterization

4.1 Structure of the parameterization

A VBS-type parameterization (e.g. Donahue et al., 2006) was developed on the basis of the explicit GECKO-A simulations. The parameterization, called VBS-GECKO, lumps the secondary organic compounds (SOCs) into bins of saturation vapour pressures at 298 K ($P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$) and considers the following:

• The formation of n volatility bins (VB_k,i with i varying from 1 to n) from the gas-phase oxidation for each precursor k (precu_k). The VB_k,i values are considered to be formed from OH oxidation for all the precursors and also from O₃ and NO₃ oxidation for alkenes and terpenes. The reaction rate constants for the oxidation of a given precursor k (k_precuk+OH, $k_{\text{precuk}+{\mathrm{O}}_{\mathrm{3}}}$ and $k_{\text{precuk}+{\mathrm{NO}}_{\mathrm{3}}})$ are taken from the GECKO-A database (see Table S1). Each VB_k,i is formed with a molecular stoichiometric coefficient that depends on NO_x via the RRR ($a_{k,\text{RRR},i}$, $b_{k,\text{RRR},i}$ and $c_{k,\text{RRR},i}$ for OH, O₃ and NO₃ reaction, respectively):

  $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{precu}}_k^{\left(g\right)}+\mathrm{OH}\stackrel{k_{\text{precuk + OH}}}{⟶}\to a_{k,\text{RRR},\mathrm{1}}{\text{VB}}_{k,\mathrm{1}}+a_{k,\text{RRR},\mathrm{2}}\\ \text{(R1)}& {\displaystyle }& {\displaystyle }{\text{VB}}_{k,\mathrm{2}}+\mathrm{\dots }+a_{k,\text{RRR},n}{\text{VB}}_{k,n}{\displaystyle }& {\displaystyle }{\text{precu}}_k^{\left(g\right)}+{\mathrm{O}}_{\mathrm{3}}\stackrel{k_{\text{precuk + O3}}}{⟶}\to b_{k,\text{RRR},\mathrm{1}}{\text{VB}}_{k,\mathrm{1}}+b_{k,\text{RRR},\mathrm{2}}\\ \text{(R2)}& {\displaystyle }& {\displaystyle }{\text{VB}}_{k,\mathrm{2}}+\mathrm{\dots }+b_{k,\text{RRR},n}{\text{VB}}_{k,n}{\displaystyle }& {\displaystyle }{\text{precu}}_k^{\left(g\right)}+{\mathrm{NO}}_{\mathrm{3}}\stackrel{k_{\text{precuk + NO3}}}{⟶}{c}_{k,\text{RRR},\mathrm{1}}{\text{VB}}_{k,\mathrm{1}}+c_{k,\text{RRR},\mathrm{2}}\\ \text{(R3)}& {\displaystyle }& {\displaystyle }{\text{VB}}_{k,\mathrm{2}}+\mathrm{\dots }+c_{k,\text{RRR},n}{\text{VB}}_{k,n}.\end{array}$$

• The ageing of the VB_k,i from gaseous OH oxidation and photolysis redistributes the matter between the various VB_k,i values.. The OH and photolysis gaseous ageing is considered for all the VB_k,i values except for the lowest volatility bin, the gas-phase fraction of that bin being expected to be negligible under atmospheric conditions. The OH reaction rate constant, k_OH, was set to the same value for each VB_k,i. These VB${}_{k,i}+\mathrm{OH}$ reactions lead to the formation of the VB_k,j with a stoichiometric coefficient $d_{k,\text{RRR},i,j}$ depending on RRR ($\forall i\ne n$):

  $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{VB}}_{k,i}^{\left(g\right)}+\mathrm{OH}\stackrel{k_{\text{OH}}}{⟶}{d}_{k,\text{RRR},i,\mathrm{1}}{\text{VB}}_{k,\mathrm{1}}+d_{k,\text{RRR},i,\mathrm{2}}\\ \text{(R4)}& {\displaystyle }& {\displaystyle }+{\text{VB}}_{k,\mathrm{2}}+\mathrm{\dots }{d}_{k,\text{RRR},i,n}{\text{VB}}_{k,n}.\end{array}$$

  Each VB_k,i (i≠n) is photolysed with a photolysis constant being proportional to the acetone one (J_acetone). The proportionality factor, ϕ_k, is considered to be the same for the VB_k,i of a given precursor k. The VB_k,i photolysis leads to a loss of matter ($\forall i\ne n$):

  $$\begin{array}{}\text{(R5)}& {\displaystyle }{\text{VB}}_{k,i}^{\left(g\right)}+\mathrm{OH}\stackrel{{\mathit{\phi }}_k{J}_{\mathrm{acetone}}}{⟶}\text{carbon lost}.\end{array}$$

• The partitioning of the precursor k and the VB_k,i between the gas and the particle phase. The partitioning is described, as in the reference GECKO-A simulations, by an absorption process following Raoult's law and considering a homogeneous, ideal, inert condensed phase at equilibrium with the gas phase:

  $$\begin{array}{}\text{(R6)}& {\displaystyle }{\text{precu}}_k^{\left(g\right)}\leftrightarrow {\text{precu}}_k^{\left(p\right)}\text{(R7)}& {\displaystyle }{\text{VB}}_{k,i}^{\left(g\right)}\leftrightarrow {\text{VB}}_{k,i}^{\left(p\right)}.\end{array}$$

Table 2 Properties of the seven bins of the VBS-GECKO parameterization.

4.2 VBS-GECKO properties

Explicit simulations were used to select the number and the range of the volatility bins and their properties (molar weight, k_OH, vaporization enthalpy and Henry's law constant). Tests were performed to establish the best number and range of volatility bins to get a compromise between the reliability and the size of the VBS-GECKO parameterization. The partitioning of organic species having a $P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$ lower than 10⁻¹² atm is almost exclusively in the particulate phase under typical atmospheric conditions. Therefore, the lower volatility bin in the VBS-GECKO parameterization lumps all species having a $P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$ lower than 10^−12.5 atm. Similarly, organic species having a $P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$ above 10⁻⁶ atm are dominantly in the gas phase under typical atmospheric conditions and volatility bin is included for species having a $P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$ above 10^−5.5 atm. Between these thresholds, the gas–particle partitioning of a species depends sensitively on temperature and OA load, so that a finer volatility discretization is desirable. A total of seven volatility bins were selected for the VBS-GECKO parameterization. The boundaries of the k volatility bins i, VB_k,i, are the same for all the precursors k. Bin intervals at 298 K are shown in Table 2 and in Fig. 5 by the different colours of the SOC bubbles.

The properties of the k volatility bins i, VB_k,i, are set to the same value for the various precursors k. Each of the VB_k,i values has for assigned saturation vapour pressure the central value of its logarithmic interval (see Table 2). Molar weight, vaporization enthalpy and Henry's law constant for each bin i are set to the mean properties of the SOCs in this volatility bin (i.e. the SOC property weighted by its concentration) during the alkanes' explicit simulations (C10, C14, C18, C22 and C26). The same k_OH value of 4.10⁻¹¹ cm³ molec⁻¹ s⁻¹ is used for all the VB_k,i values (except for the lowest volatility bin which is considered inert). This value is in the upper range of the SOC rate constants seen in Fig. 5 and is similar to that used earlier in VBS parameterization (e.g. Robinson et al., 2007; Grieshop et al., 2009; Hodzic et al., 2010). Note that the performance of the VBS-GECKO parameterization should only be weakly sensitive to the k_OH value assigned to each VB_k,i. Indeed, compensations occur during the optimization of the $d_{k,\text{RRR},i,j}$ stoichiometric coefficients (i.e. an upper k_OH value for the reactivity of the VB_k,i will lead to high values of the $d_{k,\text{RRR},i,i}$ associated with the production of the reacting VB_k,i bin and low values of the $d_{k,\text{RRR},i,j}$ associated with the production of the other bins). The properties assigned to each bin are summarized in Table 2.

In the published VBS-type parameterizations, molar weights are usually set to a value of 250 g mol⁻¹ (e.g. Tsimpidi et al., 2010; Robinson et al., 2007). This value is smaller than the mean values provided by the explicit simulations for the less volatile bins. The vaporization enthalpy used for each bin is also generally smaller in the VBS parameterization, especially for the less volatile bins (e.g. Robinson et al., 2007; Grieshop et al., 2009; Hodzic et al., 2010; Tsimpidi et al., 2010). The difference reaches 40 kJ mol⁻¹ for the less volatile bins. The vaporization enthalpies recently used by Hodzic et al. (2016) and derived from experimental data are however consistent with the VBS-GECKO values obtained here. Dry deposition of low-volatility organic compounds has been generally ignored in VBS models and no Henry's law coefficient was therefore assigned to volatility bins until Hodzic et al. (2014). Hodzic et al. (2014) provided Henry's law coefficient for the various volatility bins based on GECKO-A explicit simulations for anthropogenic and biogenic precursors. As expected, Henry's law coefficients provided in Table 2 agree with most of the values given in Hodzic et al. (2014), both being derived from similar types of simulations.

Table 3 Simulation conditions of the learning scenarios.

4.3 VBS-GECKO optimization

4.3.1 Optimization scenario and method

Explicit simulations were used to optimize stoichiometric coefficients and photolysis factors of the VBS-GECKO parameterization. A set of stoichiometric coefficients was optimized for each of the 18 precursors and for each RRR. The entire range of RRR is covered by linear interpolation of the stoichiometric coefficients. Each set (k, RRR) of stoichiometric coefficients is optimized using the explicit simulations presented in Sect. 2, with some reactions turned successively off in the GECKO-A chemical schemes to specifically evaluate the coefficients related to the reaction with each oxidant (OH, O₃, NO₃) and photolysis, as described hereafter and summarized in Table 3.

For each n-alkane and aromatic k and for each RRR, the seven $a_{k,\text{RRR},i}$ ($i=\mathrm{1},\mathrm{7}$) and 42 $d_{k,\text{RRR},i,j}$ ($i=\mathrm{1},\mathrm{6}$; $j=\mathrm{1},\mathrm{7}$) coefficients were optimized simultaneously using six explicit simulations performed with a zenith angle of 50^∘, 40 ppb of ozone, at two temperatures (270 and 298 K) and for three values of C_oa (0.1, 1 and 10 µg m⁻³). For these optimizations, the photolysis reactions of secondary organic compounds were removed in both the GECKO-A and the VBS-GECKO chemical schemes.

For each 1-alkene and terpene k, and for each RRR, the seven $a_{k,\text{RRR},i}$ ($i=\mathrm{1},\mathrm{7}$), the seven $b_{k,\text{RRR},I}$ ($i=\mathrm{1},\mathrm{7}$), the seven $c_{k,\text{RRR},i}$ ($i=\mathrm{1},\mathrm{7}$) and the 42 $d_{k,\text{RRR},i,j}$ ($i=\mathrm{1},\mathrm{6}$; $j=\mathrm{1},\mathrm{7}$) coefficients were optimized simultaneously on 11 explicit simulations and considering in each simulation the oxidation of the precursor k with a single oxidant only (either OH, O₃ or NO₃). For the simulations taking into account the precursor reaction with OH only, the same set of six simulations described above for alkanes and aromatics was used, i.e. a zenith angle of 50^∘, 40 ppb of ozone, at two temperatures (270 and 298 K) and for three values of C_oa (0.1, 1 and 10 µg m⁻³). For the simulations accounting only for the precursor reaction with O₃, a set of three simulations was used with three O₃ concentrations (10, 40 and 70 ppb), at 298 K and a C_oa of 1 µg m⁻³. For the simulations accounting only for the precursor reaction with NO₃, a set of two simulations was used with 40 ppb of ozone, at 298 K, a C_oa of 1 µg m⁻³ and for two NO₃ concentrations (0.04 and 0.4 ppt for the terpenes, 0.4 and 40 ppt for the 1-alkenes). For each oxidant case, the reactions of the precursor with the other oxidants and the photolysis reaction of secondary organic compounds were turned off in both the GECKO-A and the VBS-GECKO chemical schemes.

The set of 49 stoichiometric coefficients for the alkane and aromatic precursors and the 63 stoichiometric coefficients for the 1-alkenes and terpenes were optimized for five RRRs (0, 0.1, 0.5, 0.9 and 1). After the optimization of these coefficients, the photolysis factor ϕ_k was finally optimized for each of the 18 precursors k. Each photolysis factor ϕ_k does not depend on RRR and was optimized simultaneously using the 15 simulations performed for a precursor k with three values of zenith angle (30, 50 and 70^∘) for the five RRRs, 40 ppb of ozone, at 298 K and a C_oa of 1 µg m⁻³. The optimization was performed for different durations depending on the lifetime of the studied precursor with OH. For the less reactive species (aromatic compounds and n-alkanes), the optimization was realized for a duration of 5 days. For 1-alkenes and terpenes, the optimization duration was 2 and 1 days, respectively. The time step used for the optimization was set to 20 min.

The optimizations were performed minimizing the root mean square error (RMSE) between the GECKO-A and VBS-GECKO simulated evolutions of the total mass concentration of each of the seven volatility bins (particulate and gaseous summed). The RMSE was calculated on the entire duration of a set of reference simulations as

$$\begin{array}{}\text{(3)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\sum _{i=\mathrm{1}}^n\left\{\frac{\sum _{j=\mathrm{1}}^{n_{\text{time}}}{\left(m_{i,j}^{\text{RED}}-m_{i,j}^{\text{REF}}\right)}^{\mathrm{2}}}{n_{\text{time}}}\right\}}{n}}},\end{array}$$

where $m_{i,j}^{\text{REF}}$ and $m_{i,j}^{\text{RED}}$ are the mass of the volatility bin i at the time step number j simulated with the GECKO-A and the VBS-GECKO, respectively, n_time the number of time steps over the reference simulations and n the total number of bins. The optimization was run following a bound optimization by quadratic approximation (BObyQA) method (Powell, M., 2009) under the R software (R Core Team, 2017). The fitted stoichiometric coefficients were bounded between 0 and 1, and the photolysis factors between 0 and 100. The iterative process to minimize the RMSE stops when none of the parameters vary more than 0.2 % between two successive iterations.

For each n-alkane or aromatic precursor, the 246 parameters (49 stoichiometric coefficients for each of the five RRRs plus one photolysis coefficient) were adjusted on a total of 45 learning simulations (six simulations for each of the five RRRs plus 15 simulations for the photolysis coefficient), leading to a total of 16 200 simulated time steps (45 simulation of 120 h each with three time steps per hour) providing distinct mass concentrations for the seven volatility bins. Similarly, for each 1-alkene or terpene, the 316 parameters (63 stoichiometric coefficients for each of the five RRRs plus one photolysis coefficient) were optimized on 70 learning simulations (11 simulations for each of the five RRRs plus 15 simulations for the photolysis coefficient) leading to a total 10 080 simulated time steps for 1-alkene (70 simulations of 48 h each with three time steps per hour) and 5040 for terpene (70 simulations of 24 h each with three time steps per hour). The fitted parameters of the VBS-GECKO parameterization are provided in Table S1.

Figure 6 Temporal evolution of volatility bin mass concentration (a) and organic aerosol mass (b) simulated with GECKO-A (continuous lines) and VBS-GECKO (dashed lines). Simulated results are shown for o-xylene, α-pinene and n-octadecane for the scenario performed at a temperature of 298 K, a zenith angle of 50^∘, a C_oa of 1 µg m⁻³, 40 ppb of ozone and a RRR of 90 %. Parameters have been optimized on duration represented by white background on the figure.

4.3.2 Optimization results

Figure 6 presents the temporal evolution of the mass concentration of each volatility bin (gas and particle summed) and of OA in the GECKO-A and VBS-GECKO simulations during the oxidation of a precursor. Results are shown for the simulations performed at a temperature of 298 K, a zenith angle of 50^∘, a C_oa of 1 µg m⁻³, 40 ppb of ozone and a RRR of 90 % (i.e. in the scenario previously called the reference scenario in Sect. 3.1). Examples are given for o-xylene, α-pinene and n-octadecane oxidations. Note that, in Fig. 6, organic species with $P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$ values greater than the VBS-GECKO volatility bin upper boundary (i.e. 10^−5.5 atm; see Fig. 5) are not reported. At the beginning of the oxidation, GECKO-A explicit simulations show that first oxidation species are mainly in the highest volatility bins of the VBS-GECKO volatility range: (1) in VB1 and VB2 for o-xylene, (2) mainly in VB1 but also in VB2, VB3 and VB4 for α-pinene and (3) in VB2 and VB3 for n-octadecane. The oxidation of these first oxidation species leads to the formation of further generation species with a lower volatility, increasing the VB5, VB6 and VB7 mass concentration. As oxidation proceeds, the concentration of the volatility bins generally increases to reach a maximum before decreasing. Nevertheless, this behaviour is not observed for

• the VB7 temporal evolution. In this version of the GECKO-A tool, VB7 has no chemical sink (see above) and the concentration of the VB7 therefore increases to reach a plateau.

• the VB2 temporal evolution of o-xylene. Two species contribute to most of this bin in the explicit simulations: a dinitrocresol and a dinitro-dimethyl-phenol. Under the conditions of this simulation, these second-generation species have long lifetimes with OH and NO₃ (hundreds of hours), and therefore the chemical loss of VB2 cannot be observed on the timescale of the simulation.

Figure 6 shows that the structure of the VBS-GECKO parameterization reproduces the temporal evolution of the seven volatility bins' mass concentration during this learning simulation. The VBS-GECKO parameterization shows however discrepancies in representing the evolution of the highest volatility bins (VB1 and VB2). The temporal evolution of organic species in VB1 and VB2 is driven by the reactivity of species having $P_{\mathrm{298}\mathrm{K}}^{\text{sat}}$ values greater than 10^−5.5 atm. As no volatility bin is included for these species in VBS-GECKO, the VB1 and VB2 mass evolution is less well reproduced than the bins of lower volatility. Species in VB1 and VB2 do not, however, partition substantially into the condensed phase except under heavy aerosol loading and low temperature for VB2. The reliability of the optimization was evaluated comparing the temporal evolution of OA mass (which is not directly optimized during the fitting process) between GECKO-A and VBS-GECKO simulations (see Fig. 6). The relative root mean square error (RRMSE) was calculated a posteriori for the aerosol mass concentration as

$$\begin{array}{}\text{(4)}& \text{RRMSE}=\sqrt{{\displaystyle \frac{\sum _{j=\mathrm{1}}^{n_{\text{time}}}{\left\{\frac{\left(M_j^{\text{RED}}-M_j^{\text{REF}}\right)}{M_j^{\text{REF}}}\right\}}^{\mathrm{2}}}{n_{\text{time}}}}},\end{array}$$

with $M_j^{\text{REF}}$ and $M_j^{\text{RED}}$ are the mass concentration of OA formed at each time step j with the GECKO-A and the VBS-GECKO, respectively, and n_time the total number of time steps. The temporal evolution of OA mass is well reproduced by the parameterization for o-xylene, α-pinene and n-octadecane, with a RRMSE lower than 20 % (10.2 %, 17.1 % and 1.8 % for o-xylene, α-pinene and n-octadecane, respectively).

Figure 7 Distributions of the RRMSE on organic aerosol mass concentration calculated between VBS-GECKO and GECKO-A on the learning simulations. The bottom and top of the boxes are the first and third quartiles; the bands inside the boxes represent the medians. The ends of the whiskers represent the lowest RRMSE still within 1.5 times the interquartile range (IQR) of the lower quartiles and the highest RRMSE still within 1.5 IQR of the upper quartiles. Red points represent the mean values and empty points the outliers. There are 45 simulations for n-alkane or aromatics and 70 for 1-alkene or terpene. Time steps for which the OA mass simulated with GECKO-A is under a threshold of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ µg m⁻³ (i.e. a SOA yield lower than around 1 %) are not considered in RRMSE calculation.

Figure 7 shows the distributions of the RRMSE for the OA mass concentration for the various learning simulations and for each precursor. For the n-alkanes and 1-alkenes, the RRMSE increases when the carbon chain length decreases. The mean RRMSE on OA mass evolutions is less than 5 % for n-alkanes and 1-alkenes bearing more than 14 carbon atoms. It reaches 23 % and 31 % for n-decane and 1-decene. For terpenes, the mean of RRMSEs are of 24 %, 23 % and 17 % for α-pinene, β-pinene and limonene, respectively. For the aromatic precursors, the mean of RRMSEs are of 23 %, 45 %, 32 %, 36 % and 20 % for benzene, toluene, o-, m- and p-xylene, respectively. The VBS-GECKO parameterization reproduces SOA masses with a RRMSE lower than 50 % for more than 90 % of the learning simulations.

4.4 VBS-GECKO evaluation

4.4.1 Evaluation scenario

Results of the GECKO-A simulations performed in the previous sections were used as a training dataset to fit the coefficients of the parameterization. Here, new simulations are conducted with GECKO-A to produce a validation dataset to assess VBS-GECKO. New box model scenarios were therefore designed for that purpose. These scenarios use a mixture of precursors and environmental conditions not encountered during the training phase, including variations of environmental parameters with time (C_oa, temperature, photolysis and NO_x levels). Scenarios were run for two constant C_OA values (2 and 20 µg m⁻³), two types of meteorological conditions for temperature, relative humidity and zenith angle (a “summer” and a “winter” case, called hereafter SUM and WIN, respectively) and two NO_x conditions (a “high NO_x” and a “low NO_x” case, called hereafter HNOx and LNOx, respectively). A total of eight scenarios were run.

Diurnal temperature profiles were represented by a sinusoid function with a 5 K amplitude, a maximum reached at 14:00 LT and average values of 290 and 275 K for the SUM and WIN scenarios, respectively. Diurnal photolysis frequencies were computed for midlatitudes using the TUV model (Madronich and Flocke, 1999) with zenith angle variations corresponding to 1 July for the SUM case and 1 January for the WIN scenario. The relative humidity was set to 60 % and 85 % for the SUM and WIN scenarios, respectively. For the HNOx and LNOx cases, the initial concentration of NO_x was set at 50 and 0.5 ppb, respectively. The OH lifetime with respect to the reaction with the surrogate VOC was again set to 1 s. Chemical concentrations were also fixed for CO, methane and O₃ at 120, 1750 and 40 ppb, respectively. The evaluation scenarios started at 12:00 LT and lasted 5 days. Simulated evolutions of NO, NO₂, NO₃, OH, HO₂ and O₃ concentrations for the HNOx_WIN, HNOx_SUM, LNOx_WIN and LNOx_SUM scenarios are presented in Figs. S2–S5.

For each scenario, three mixtures of precursors summing up to an initial concentration of 10 pptC were considered: a “biogenic”, an “intermediate” and an “anthropogenic” case (called hereafter BIO, INT and ANT, respectively). In the BIO case, biogenic compounds (i.e. α-pinene, β-pinene and limonene) represent 90 % in carbon of the initial mixture and anthropogenic compounds (i.e. the five alkanes, the five alkenes and the five aromatic compounds) the remaining 10 %. In the INT and ANT cases, the biogenic∕anthropogenic ratio in carbon is 50 %∕50 % and 10 %∕90 %, respectively. In the biogenic or anthropogenic family, each species was arbitrarily set to the same level (in pptC units).

A total of 24 simulations (eight environmental scenarios using each three sets of initial mixture of precursors) were run with the GECKO-A chemical scheme and VBS-GECKO parameterization.

Figure 8 Simulated organic aerosol (OA) mass with GECKO-A (continuous lines) and VBS-GECKO (dashed lines), for the different evaluation scenarios using C_oa=2 µg m⁻³. Panel (a) presents results for the ANT cases, (b) for the INT cases and (c) for the BIO cases. The different colours represent the different scenarios (see text).

Figure 9 Distribution of the relative errors (RE_i) on organic aerosol mass concentration calculated between VBS-GECKO and GECKO-A for the 24 evaluation scenarios. The bottom and top of the boxes are the first and third quartiles; the bands inside the boxes represent the medians. The ends of the whiskers represent the lowest and the highest RE_i values. There are 361 time steps for each evaluation scenario.

4.4.2 Evaluation results

Figure 8 shows the temporal evolution of the OA mass concentrations formed according to the GECKO-A and the VBS-GECKO, for the different evaluation scenarios performed with a C_oa of 2 µg m⁻³. Temporal evolutions simulated with a C_oa of 20 µg m⁻³ are presented in Fig. S6. Figure 9 presents the distribution of the relative errors on the OA mass concentrations for the 24 evaluation scenarios. The relative error (RE) is calculated at each time step i as

$$\begin{array}{}\text{(5)}& {\text{RE}}_i={\displaystyle \frac{\left(M_i^{\text{RED}}-M_i^{\text{REF}}\right)}{M_i^{\text{REF}}}},\end{array}$$

where $M_i^{\text{REF}}$ and $M_i^{\text{RED}}$ are the mass concentrations of OA formed at each time step i with the GECKO-A and the VBS-GECKO, respectively.

Reference GECKO-A simulations show generally (1) a fast increase of OA mass concentration during the first 10 h of oxidation, (2) a clear diurnal cycle of OA concentration linked to the temperature variations, with a larger amplitude for the WIN than for the SUM scenarios, (3) higher simulated OA concentrations for the LNOx scenarios with maximum concentration obtained during the summer, (4) a slight increase of OA concentrations with the increase of anthropogenic precursors in the HNOx scenarios and with the increase of biogenic precursors in the LNOx scenarios. These trends are well represented by the VBS-GECKO parameterization as well as the orders of magnitude of the OA concentrations (mean errors between ±20 %). Low errors are simulated for the ANT mixture (with $\left|{\text{RE}}_i\right|<\mathrm{15}$ %). The errors generally increase with the increase of biogenic precursors in the mixture (with $\left|{\text{RE}}_i\right|<\mathrm{30}$ % for the BIO mixture). In particular, for the summer LNOx scenarios, VBS-GECKO overestimates OA mass from the second day. This overestimation increases with the length of the simulation.

The origin of these errors is not entirely clear, as the role of many important driving variables is masked by the training and fitting procedures. It is, however, possible to list some of the major simplifications and approximations that were made in developing the VBS-GECKO parameterization: (i) discretization into seven decadal bins (rather than continuous vapour pressure distribution, for example); (ii) use of the same matrix of VB plus OH stoichiometric coefficients to represent the reactivity of organic species formed from the oxidation of the precursor with OH, O₃ and NO₃; (iii) use of a single photolytic scale factor (relative to acetone photolysis) to represent the overall photolysis of a complex mixture containing many photo-labile species; (iv) discretization into five RRRs with separate optimization and linear interpolation to cover the entire range; and (v) selected simulation duration used for the optimization that is expected to influence the weight given at first or further generation products. More detailed consideration of these important parameters presents a clear opportunity to further improve the parameterization.

5 Conclusions

The GECKO-A modelling tool was used to explore SOA formation from the oxidation of various hydrocarbons in a box model under stationary state scenarios representative of environmental conditions. The set of parent hydrocarbons includes n-alkanes and 1-alkenes with 10, 14, 18, 22 and 26 carbon atoms, α-pinene, β-pinene and limonene, benzene, toluene, o-xylene, m-xylene and p-xylene. The developed environmental scenarios allow investigation of the sensitivity of SOA formation to changes of physical (photolysis, C_oa, T, relative humidity) and chemical (${\mathrm{O}}_x/{\mathrm{NO}}_x/{\mathrm{HO}}_x$) conditions. Simulated trends of maximum SOA yields are consistent with the literature data. In particular, Y_max increases with the carbon skeleton size of the parent the n-alkane or 1-alkenes, from 0.2 for the C₁₀ to a plateau value of 1.2 for the C_>18. For the aromatic compounds, Y_max slightly decrease with number of methyl groups with values below 0.15. For terpenes, simulated Y_max increases in the following precursor's order: α-pinene, β-pinene and limonene. As expected, Y_max increases when temperature decreases and C_OA increases. NO_x conditions drive the distributions of secondary organic species and changes substantially the volatility distribution of the species. Simulated SOA yields generally decrease with the increase of NO_x (or RRR). SOA yields from terpene are particularly sensitive to NO_x levels, with values in the 1.1–1.75 range under low NO_x conditions (RRR close to 0) to values in the 0.05–0.2 range under high NO_x conditions (RRR close to 1). SOA yields for species having a large carbon skeleton (i.e. with n_C>18) show a low sensitivity to temperature, C_OA and NO_x regimes as they rapidly form compounds that partition dominantly to the condensed phase. Explicit simulations were used to explore the distribution of the millions of secondary organic species formed during the oxidation of the precursors in term of molar mass, vaporization enthalpy, OH reaction rate and Henry's law coefficient.

Simulation results obtained with the GECKO-A modelling tool were used as a reference to fit a VBS-type parameterization. The VBS-GECKO parameterization is designed for 18 selected hydrocarbons and includes OH oxidation for all the hydrocarbons, as well as NO₃ and O₃ oxidation for the terpenes and 1-alkenes. For each hydrocarbon, seven volatility bins are considered, ageing by reaction with OH, photolysis and partitioning according to Raoult's law. The dependence of SOA formation on NO_x levels is represented by a linear interpolation of different sets of parameters optimized for various RRRs. Properties of the bins, i.e. molar weight, vaporization enthalpy, OH rate constants and Henry's law coefficient, are selected from the reference explicit simulations. The evaluation of the parameterization shows that VBS-GECKO captures the dynamic of SOA formation for a large range of conditions with a mean relative error on OA mass temporal evolution lesser than 20 % compared to explicit simulations.

Explicit simulations of organic compound oxidation provided by GECKO-A have been previously used to develop a reduced chemical scheme for the ${\mathrm{HO}}_x-{\mathrm{NO}}_x-\mathrm{VOC}$ chemistry (Szopa et al., 2005) as well as to develop reduced parameterizations for SOA formation from semi-volatile and intermediate volatility n-alkanes (Hodzic et al., 2016). For the first time in this study, a reduction protocol has been developed and evaluated to optimize SOA parameterizations from explicit simulations for a large number of precursor species and under various environmental conditions. The optimization procedure has been automated to facilitate the update of the VBS-GECKO on the basis of the future GECKO-A versions, its extension to other precursors and/or its modification to carry additional information (e.g. optical properties, hygroscopicity, N∕C ratio). The implementation and the evaluation of VBS-GECKO in a 3-D model are the purpose of a companion paper.