Advanced modeling of gas chemistry and aerosol dynamics with SSH-aerosol v2.0

SSH-aerosol is developed to represent the evolution of primary and secondary pollutants in the atmosphere by processes linked to gas-phase chemistry, aerosol dynamics (coagulation, condensation/evaporation and nucleation) and intra-particle reactions. The representation of process complexity can be adjusted based on the user's choices.

The model uses a sectional size distribution, and offers the capability to discretize chemical composition to account for the mixing state of particles. The algorithms are designed to represent the evolution of ultrafine particles: conservation of mass and number during numerical resolution, taking into account the Kelvin effect, the condensation dynamics of nonvolatile compounds, and nucleation. Different parameterizations are provided for nucleation: binary, ternary, heteromolecular and organic nucleation depending on the compounds involved.

For gas-phase chemistry, schemes of different complexities can be handled: from simple schemes to model ozone, oxidants and inorganic chemistry (e.g. CB05, RACM2, Melchior2), to more complex schemes, e.g. from the Master Chemical Mechanism (MCM). The complexity of the schemes used for secondary organic aerosol (SOA) formation may also be adjusted: from schemes built from chamber data to near-explicit schemes from MCM. SOA schemes reduced using the GENOA algorithm are also provided for several precursors (toluene, a sesquiterpene and three monoterpenes), together with their evaluation against chamber or flow-tube experiments. A wall-loss module has also been added for easier comparisons to chamber experiments.

Specific developments were made in version 2.0 to automatically link the chosen gas-phase mechanism to SOA formation by using the SMILES structure of organic compounds, allowing for the determination of their hydrophilic and hydrophobic properties and for the partitioning in both organic and aqueous phases. The gas/particle partitioning may also be represented with different complexities. For the organic phase, viscosity may be modelled, adapting the aerosol viscosity to its composition, and coupling organic and inorganic thermodynamics. The dynamic evolution of the partitioning may be computed explicitly or thermodynamic equilibrium may be assumed. Different options are also provided to simulate the chemistry of organic compounds inside the particles with different types of reactions: irreversible 1st order reactions, bulk oligomerization, hydratation of aldehydes and reactions of organic compounds with inorganic ions.

The SSH-aerosol model may be installed with a docker for standalone use. It has also been coupled to several 3D models to represent gas and aerosol concentrations: from the local scale with computational fluid dynamic and street network models to the regional scale with chemistry-transport models.

1 Introduction

The simulation of three-dimensional atmospheric pollutant concentrations is influenced by numerous factors, each associated with potentially significant uncertainties, including emissions, meteorological conditions, transport processes, and deposition. The formation and evolution of aerosols involve many chemical compounds and several complex processes, which are often very simplified in 3D air-quality and climate models, or even not taken into account, despite the fact that some of these processes may be critical for accurately representing aerosol properties and concentrations. The creation of aerosol box models allows for the independent evaluation of the various processes involved in aerosol dynamics and atmospheric chemistry, free from the uncertainties associated with 3D modeling and the computational limitations it imposes. They can therefore be used to simulate both aerosol formation and evolution in enclosed environment (like atmospheric smog chambers or flow tubes) with a detailed representation of processes, or inside 3D air quality models with generally more simplified approaches. Furthermore, the creation of box models enables the coupling of the same aerosol model with different 3D models across various scales (e.g. in regional and urban air quality models, Maison et al., 2024a), ensuring a consistent representation of aerosol chemistry and dynamics.

Challenges in the development of aerosol box models are numerous: representation of the mixing state (Riemer et al., 2019), gas-particle interactions (Pöschl et al., 2007), formation of secondary inorganic and organic aerosols and their properties, formation and growth of ultrafine particles (Lee et al., 2019). Consequently, the choice of the model type and the aerosol representation is guided by the need to adequately represent and resolve these distinct challenges. Several aerosol box models have been developed, and most represent aerosol dynamics using a sectional approach (SIREAM, Debry et al., 2007; SCRAM, Zhu et al., 2015; ATRAS2, Matsui, 2017; MAFOR, Karl et al., 2022; MOZAIC, SALSA, Kokkola et al., 2018; ARCA, Clusius et al., 2022), a modal approach (MAM, Sartelet et al., 2006; M7, Vignati et al., 2004) or a particle-based approach (Curtis et al., 2024).

In SSH-aerosol, the sectional approach was selected over the modal approach because it is numerically simpler and offers greater precision in representing the growth of ultrafine particles (Devilliers et al., 2013). The SCRAM module for aerosol dynamics is used because of its ability to represent the particle mixing state by discretizing both size and composition. Few models account for the mixing state; one example is ATRAS2, which specifically represents the mixing state of black carbon, and M7, which uses strong simplifying assumptions when mixing different modes. In SSH-aerosol, the mixing state of each compound can be represented explicitly, as in the SCRAM module, but modeling the mixing state may lead to longer computation times than in the classical case where aerosols are assumed to be internally mixed (Zhu et al., 2016).

For atmospheric chemistry, various mechanisms are commonly used to represent the formation of ozone and oxidants, using different strategies to lump volatile organic compounds (VOCs). Examples include CB05 (Yarwood et al., 2005), RACM2 (Goliff et al., 2013), SAPRC (Carter, 2010) and MELCHIOR2 (Derognat et al., 2003). For the formation of organic oxidation products that could be semi-volatile and partition to form secondary organic aerosols (SOAs), different gaseous oxidation mechanisms of varying complexity exist. Near-explicit mechanisms, such as the Master Chemical Mechanism (MCM) (Jenkin et al., 2012) and those generated by the Generator for Explicit Chemistry and Kinetics of Organics in the Atmosphere (GECKO-A) (Aumont et al., 2005; Camredon et al., 2007), involve a large number of compounds and reactions, whereas empirical mechanisms built from chamber experiments use only a limited number of compounds and reactions, e.g. the Volatile Basis Set VBS (Donahue et al., 2006), the 2-products (Odum et al., 1996) and surrogate (Couvidat et al., 2012) approaches. Recently, the development of the GENerator of reduced Organic Aerosol mechanisms (GENOA) has allowed the reduction of near-explicit scheme while retaining major formation pathways and key properties of the resulting secondary organic aerosol, such as functional groups (Wang et al., 2022, 2023b). Reduced mechanisms are built by reducing and lumping reactions and species using different strategies (e.g. lumping, replacing, jumping, removing) under a variety of environmental conditions. These reductions preserve the main chemical pathways and represent radical reactions, which may strongly impact the formation of organic particles (Wang et al., 2024). Conserving molecular properties of organic compounds allows the estimation of their partitioning in both aqueous and organic phases, the influence of non-ideality and organic-inorganic interactions (Kim et al., 2011). The partitioning between gas and particles may be strongly influenced by viscosity (Kim et al., 2019; Schervish et al., 2024) and intra-particle reactions (Pun and Seigneur, 2007), which are still very rarely taken into account in 3D models. In this paper, the term SOA is used in a broad sense to refer to the compounds involved in the formation of organic particles, regardless of whether they are in the gas or particle phase.

The SSH-aerosol model was built to address these different challenges, with the ability to adjust the model's complexity depending on the user's choice. Numerous test cases are provided to check the model accuracy and to illustrate its capacities. Section 2 summarizes the features from the previous version and describes the code structure as well as new features. Section 3 focuses on the use and evaluation of chemical schemes of different complexity for SOA formation. Section 4 describes the advance gas/particle partitioning taking into account viscosity and particle-phase reactions, and Sect. 5 addresses the modeling of ultrafine particles. Tools for coupling to 3D models and applications using SSH-aerosol coupled to 3D models are listed in Sect. 6.

2 Model structure and new features

As illustrated in Fig. 1, the model was initially based on the merge of three state-of-the-art models:

• SCRAM: The Size-Composition Resolved Aerosol Model (Zhu et al., 2015), written in Fortran 90, which simulates the dynamics (coagulation, condensation/evaporation and nucleation), classifying particles by both composition and size with the possibility to resolve the mixing state of atmospheric particles.

• SOAP: The Secondary Organic Aerosol Processor (Couvidat and Sartelet, 2015), a thermodynamic model that computes the partitioning of organic compounds, taking into account hygroscopicity, absorption into the aqueous and organic phases of particles, non-ideality, viscosity and phase separation. it is written in C++.

• H²O: The Hydrophilic/Hydrophobic Organics (Couvidat et al., 2012) mechanism, which uses a molecular surrogate approach to represent the myriad of formation of semivolatile organic compounds formed from the oxidation in the atmosphere of volatile organic compounds. Gas-phase chemistry is solved using algorithms written in Fortran 90.

The size distribution is represented using a sectional approach (Debry et al., 2007), in which the particle population is approximated by discrete size sections characterized by their mass, number, and mean diameter. Depending on user settings, the section boundaries may either remain fixed or evolve in time, and the mean diameter may be constrained to the geometric midpoint of the section or allowed to vary within its bounds. Particles are assumed to be spherical, so that the mass M, number N, and mean diameter d_m in each section satisfy $M=N\mathit{\rho }\mathit{\pi }/\mathrm{6}{d}_{\mathrm{m}}^{\mathrm{3}}$, with ρ the particle density, which can be fixed or vary with the particle composition.

Figure 1 Workflow of the SSH-aerosol model.

2.1 Model structure

The SSH-aerosol package contains different directories for source code, configuration files, input files, output files, and output visualization. The model structure is the same as previous versions: the main program file is ssh-aerosol.f90 where the configuration file is read, as well as the list of gaseous and particle compounds and their properties. The different variables, such as meteorological parameters and initial concentrations are read from files. Partition coefficients for coagulation (coefficients indicating how particles formed by the coagulation of two smaller particles should be redistributed into size sections) are computed once at the beginning of simulations. They can be written into files, which can then be used to avoid their recalculation. The gaseous chemistry (e.g. H²O module) and aerosol dynamics (SCRAM and SOAP modules) are solved sequentially with a time step determined by the user. Note that this time step differs from the time steps used to solve gaseous chemistry and aerosol dynamics, which are automatically adapted, as detailed in Sect. 2.3.

In SSH-aerosol version 1.1, gaseous chemical reactions were read by a pre-processor and interpreted before their use in SSH-aerosol. In version 2.0, this pre-processor is removed, and the chemical reactions are directly read and interpreted in SSH-aerosol, eliminating the need to recompile the program each time the chemical scheme changes.

Simulations require a namelist configuration file that details the chemical species and reactions, input and output conditions, as well as various simulation setups. The gas phase species are listed in the file gas_phase_species together with their molar mass, while aerosol species are detailed in a separate file aerosol_species_list_file. Each row represents an aerosol species and its specific properties: the species type (1 for dust, 2 for elemental carbon, 3 for inorganics and 4 for organics), the group to which the species belongs in case of mixing-state resolved particles, the name of corresponding gas-phase species, its molar weight (g mol⁻¹), the collision factor (probability that a particle–particle encounter leads to effective coalescence), molecular diameter (Angstrom), surface tension (N m⁻¹), accommodation coefficient (between 0 and 1), density in kg m⁻³ and whether or not the species is nonvolatile. The column “partitioning” indicates whether the species condense on only an organic phase (keyword: HPHO), only an aqueous phase (keyword: HPHI) or both (keyword: BOTH). In the next column, the SMILES structure of the species or its decomposition into the 60 UNIFAC functional groups (see Table 4) can be given, then the saturation vapor pressure (in Torr), the enthalpy of vaporization, the Henry constant, and the reference temperature at which it is specified. Note that the Henry constant is not needed if both the saturation vapor pressure and the SMILES structure (or decomposition into functional groups) are provided. The saturation vapor pressures are provided as model inputs and may be derived through different methods or parameterizations, including tools such as UManSysProp (Topping et al., 2016).

As a new feature in SSH-aerosol v2.0, outputs can now be written in a NetCDF format rather than a text format. Writing NetCDF output files is faster than text or binary outputs. Furthermore, SSH-aerosol v2.0 can be installed using Docker, allowing for a consistent setup across different platforms (macOS, Windows, and Linux/Unix).

2.2 Physical processes

The processes included in SSH-aerosol, as well as the main options, are illustrated in Fig. 2. New processes and options are highlighted in red, while those that have been improved since version 1.1 are highlighted in orange. A wide variety of configuration files, i.e. “test cases”, are presented in the supplementary material (user manual) to illustrate the model capabilities. These test cases are also used for model checks from one model version to the next. Only modifications that do not affect the reproduction of the test cases or that can be explained by model improvements or corrections are kept.

Figure 2 Main processes and options available in SSH-aerosol. New processes and options are highlighted in red, while those that have been improved since version 1.1 are highlighted in orange.

Several test cases were presented in Sartelet et al. (2020). They concern:

• Coagulation and condensation for internally-mixed particles or mixing-state resolved particles. The test cases consider an initial size distribution of internally-mixed particles following the urban/hazy conditions of Seigneur et al. (1986) and Zhang et al. (1999), which represent stringent conditions for coagulation/condensation respectively. For the condensation test case, an initial concentration of sulfuric acid that condenses onto particles is specified (9.9 µg m⁻³). The number and volume distribution after 12 h of simulation are the same as in Zhang et al. (1999). These test cases are also extended to the case of particles made of two populations of different composition to illustrate how coagulation/condensation affects the mixing state of particles, with the number and mass distribution being a function of the size and fraction of the chemical compounds.

• Condensation/evaporation: influence of redistribution. Because of condensation/evaporation, the bounds of diameter sections may evolve (due to growth or shrinking of particles). Because they need to be fixed if both condensation/evaporation and coagulation are solved, as well as in 3D Eulerian models, redistributing the mass and number concentrations among fixed sections is necessary. This redistribution leads to numerical diffusion. Different algorithms for redistribution are proposed and illustrated using the condensation test case (the moving-diameter approach or the Euler-coupled approach).

• Condensation/evaporation: thermodynamic equilibrium assumption versus dynamic. To save computational time, gas-phase concentrations at the particle surface are often assumed to be at equilibrium with the bulk gas concentrations, simplifying the calculation of condensation/evaporation. Bulk particle equilibrium concentrations are distributed between the different sections, using weights that depend on the condensation/evaporation kernel of the condensation/evaporation rate. A hybrid method is used where particles above a certain cutoff diameter are treated dynamically, while smaller particles are assumed to be in thermodynamic equilibrium. A simulation of Tokyo's pollution episode of the 25 June 2001 (Sartelet et al., 2006) shows that while dynamic modeling is more accurate, it is slower compared to using thermodynamic equilibrium, with the hybrid method offering a compromise.

• Condensation/evaporation: role of Kelvin effect. To highlight the significance of the Kelvin effect in the growth of ultra-fine particles, a simulation is conducted based on the emission of particle of nonadecane (C₁₉H₄₀) from a diesel engine exhaust, described in Devilliers et al. (2013). Modeling condensation/evaporation dynamically, two simulations are performed: one including the Kelvin effect and one without. The results show that the Kelvin effect is crucial for accurately modeling the evolution of ultra-fine particles made of a semi-volatile compound.

• Condensation/evaporation: influence of viscosity using a constant coefficient. The mass transfer of hydrophobic organic aerosols can be significantly influenced by particle viscosity. To explore this, the condensation of organic surrogates with varying volatilities is studied for different viscosities, represented by varying diffusion coefficients (D). The particles are modeled with the approach of Couvidat and Sartelet (2015), which implicitly represents the diffusion of organic aerosol inside aerosol, with five aerosol layers. The results show that the condensation of low-volatility surrogates is unaffected by viscosity, while higher-viscosity particles (low D) may inhibit the condensation of more volatile surrogates.

• Condensation/evaporation: absorption into organic or aqueous phases, influence of ideality. The oxidation of VOCs produces less volatile compounds that may condense onto particles, increasing particle mass. To study this, experiments on exhaust emissions from a Euro 5 gasoline car are compared to chamber measurements of organic aerosol concentrations from Platt et al. (2013). The simulation shows that more than half of the particle mass after 5 h comes from the condensation of aged intermediate/semivolatile VOCs (IVOCs/SVOCs). Additionally, the presence of ammonia (NH₃) emissions significantly increased inorganic aerosol concentrations. The study also explored how the hygroscopicity of compounds affects their phase, showing that isoprene oxidation products prefer to condense onto aqueous phase. The impact of ideality on SOA formation shows that activity coefficients greatly influence the SVOC gas-to-particle partitioning, particularly in the aqueous phase.

• Coupling between nucleation, coagulation and condensation. SSH-aerosol model v1.1 incorporates homogeneous binary and ternary nucleation schemes. A test case is run under the hazy conditions of Zhang et al. (1999) solving simultaneously nucleation, condensation, and coagulation. The results show that ultra-fine particles grow over time because of the combined influence of the different processes. Furthermore, the growth is faster when considering extremely-low-volatility compounds, such as the Monomer surrogate (an ELVOC from terpene oxidation).

2.3 Numerical implementation

The gaseous chemistry and aerosol modules of SSH-aerosol are invoked in the routine ssh-aerosol.F90. Gaseous chemistry (H²O module) and aerosol dynamics (SCRAM and SOAP modules) are solved sequentially. The outer loop time step is user-defined and corresponds to the output frequency of concentrations; however, the internal time steps for gaseous chemistry and aerosol dynamics differ from this output interval. Each module selects its own adaptive time step based on its numerical requirements.

Because the governing equations are stiff, dedicated stiff solvers are employed. Gaseous and organic chemistry are integrated using the two-step solver of Verwer (1994), which replaces the Rosenbrock scheme used in version 1. Aerosol dynamics are solved with the second-order explicit trapezoidal rule (ETR) (Ascher and Petzold, 1998; Sartelet et al., 2006). These algorithms ensure numerical stability and computational efficiency while conserving particle mass and number across the discrete size spectrum. For each process, the initial time step is set to the minimum value specified by the user, and it is subsequently adjusted adaptively according to the user-defined numerical tolerance. In SCRAM, coagulation may be solved either simultaneously with nucleation and condensation/evaporation or treated separately, depending on the user's choice. When the processes are split, coagulation is computed first because it is the slowest of the three. Nucleation and condensation/evaporation, which compete for the same gas-phase precursors, are always solved simultaneously. When bulk thermodynamic equilibrium is assumed, the equilibrium calculation is performed at the end of each main time step, after resolving coagulation, nucleation, and condensation/evaporation for the dynamic sections and the low-volatility species. Evaporation can lead to the complete depletion of certain size bins. To ensure numerical robustness in such cases, the model employs mass-threshold parameters (TINYM and TINYN) that control accuracy and prevent the accumulation of non-physical residual masses. For layers inside particles, the volume ratio of layers is assumed to remain constant. Therefore, a layer cannot vanish. Depending on the user's choice, mass and number concentrations can be remapped onto the fixed size sections after each time step. This remapping/redistribution becomes mandatory when coagulation is activated, because the coagulation kernels and partition coefficients are precomputed for fixed bin boundaries.

2.4 New features

Some of the new features of v2.0 concern gas-phase chemistry for SOA formation, with the possibility to use near-explicit chemical schemes for SOA formation. The properties of the organic species can now be detailed in the configuration file using their SMILES (simplified molecular input line entry system) code or functional group distribution according to the UNIFAC (Unified Functional Activity Coefficient) nomenclature, as well as their saturation vapor pressure and enthalpy of vaporization. When provided, SMILES code are automatically decomposited in UNIFAC functional groups. Chemical schemes can be obtained from MCM, and various near-explicit schemes are added in the software, namely for the oxidation of toluene, naphtalene, monoterpenes and sesquiterpenes. More specifically, the new functionalities related to gaseous chemistry concern:

• The coupling to the GENerator of Reduced Organic Aerosol Mechanisms (GENOA) and Master Chemical Mechanism

• The treatment of specific MCM and GECKO-A kinetic rates

• The replacement of the ROS2 numerical solver for gaseous chemistry by the two-step solver (Verwer et al., 1996) to reduce computational time.

For easier comparisons with chamber experimental results, a wall-loss module for vapors and particles and numerous comparisons to chamber experiments have been added.

New features also concern the partitioning between gas and particle, more specifically

• Addition of an intra-particle reaction module (oligomerization, hydrolysis, and reactions between organic and inorganic ions).

• For organic thermodynamic, the particle viscosity can now be estimated depending on the particle composition

• Coupled inorganic-organic thermodynamic can be handled in SOAP.

• Addition of nucleation and numerical schemes to simulate ultrafine particles:

Other new features are related to the dynamics of particles and the modelling of ultrafine particles:

• Consideration of the Kelvin effect to determine size distribution of condensing and evaporating species when assuming thermodynamic equilibrium

• Faster computation of coagulation partition coefficients

• New parameterizations of nucleation

These new features are detailed and illustrated with test cases in the following sections.

2.5 Runtime

As detailed in the appendix A, many options are available, allowing the representation of the different physical processes with a large range of complexity. It is illustrated with different test cases in this paper. Increasing complexity may often lead to a large increase in computational time (CPU), as illustrated in Table 1. Near-explicit chemical schemes may be 3 to 65 times more expensive than parameterized schemes (comparisons of simulation numbers 1, 2, 3; numbers 4 and 7). Taking into account ideality has a limited impact on CPU time, with an increase by a factor about 1.5–2 (comparisons of simulation numbers 4 and 5; numbers 4 and 6; numbers 9 and 10). The impact of oligomerisation on CPU time may be limited (comparison of simulation numbers 12 and 13), while the impact of viscosity may be high, with an increase by a factor almost 8 (comparison of simulation numbers 10 and 11). Dynamic gas/particle partitioning may lead to a large increase of CPU time, from a factor 2 to a factor 20 (comparisons of simulation numbers 13 and 14; numbers 15 and 16; numbers 17 and 18; numbers 19 and 20). The relations between CPU time and model complexity presented here are only illustrative, and they may vary depending on the environmental conditions specified in the model's input.

Table 1 Runtime (in s) for different simulations presented in the Figures of this paper. Rdc. stands for reduced and Expl. for quasi-explicit.

3 Chemical schemes for SOA formation

The oxidation of VOCs in the troposphere results in lower volatility organic compounds that can condense on particles under favorable conditions and form SOA. To address the influence of detailed gas-phase chemical processes of VOC oxidation, SOA formation may be modeled with schemes of different complexities. In SSH-aerosol v1.1, the H²O scheme (Couvidat et al., 2012; Sartelet et al., 2020) is available. This scheme was built from chamber experimental yields in order to represent SOA formation with only a few chemical surrogate species. This scheme is suitable for 3D simulations because it allows a limited number of compounds to be used, which means that calculation time is affordable. This scheme also has other advantages, such as differentiation of low-NO_x and high-NO_x pathways, and the allocation of a molecular structure to most surrogate species with molecules, allowing to estimate the partitioning between gaseous and organic or gaseous and aqueous phases. However, the chemical pathways are too simplified, which makes it difficult to accurately predict future air quality, for example in the event of significant changes in NO_x levels (Wang et al., 2024). In SSH-aerosol v2.0, semi-volatil organic compound formation can be represented using a range of chemical schemes, including highly detailed chemical mechanisms, such as those of the Master Chemical Mechanism (MCM) (Jenkin et al., 1997), as well as those presented in Lannuque et al. (2023), Wang et al. (2023b), and Lannuque and Sartelet (2024), or schemes reduced from near-explicit chemical schemes using GENOA (Wang et al., 2023b; Sartelet et al., 2024). The SOA module may be customized to use more or less complex SOA chemical scheme depending on the SOA precursor.

The inorganic reactions governing oxidant formation do not always need to be included in the reaction list. Depending on the user's settings, oxidant concentrations can either be prescribed using constant input profiles, or they can be computed using existing gas-phase chemistry mechanism, such as CB05, RACM2 or MELCHIOR2. These schemes can be complemented with SOA chemical schemes describing the degradation of VOCs leading to SOA formation. To offer flexibility, the model allows users to construct customized chemical mechanisms, covering both gas-phase and SOA chemistry, by combining a chosen gas-phase mechanism with one or more SOA schemes (typically one per SOA precursor class). For VOCs already represented explicitly in the selected gas-phase mechanisms (e.g. toluene, monoterpenes), the associated SOA scheme can be added without modifying oxidant concentrations. However, some VOCs are not included in these gas-phase mechanisms, in which case the inclusion of specific reactions for radicals, oxidants as well as SOA formation may be necessary. For example, the impact of naphthalene on ozone and radical production is not represented in the CB05, RACM2 or MELCHIOR2 mechanisms.

Instructions for downloading the MCM chemical schemes from the MCM website and a converter for preparing input files in the SSH-aerosol format are detailed in the supplementary material (user manual). Since aerosol properties are not provided in MCM chemical schemes, users also need to define variables such as the method used to compute saturation vapor pressures for SVOCs. The UManSysProp (Topping et al., 2016) is used by the converter to compute aerosol-related properties accordingly, often using the vapor pressure estimation methods (“v0”, Myrdal and Yalkowsky, 1997), with the boiling point estimation methods (“b0”, Nannoolal et al., 2004).

In the following section, a simulation with the MCM beta-caryophyllene scheme with fixed oxidant concentrations is first detailed. Then more general test cases illustrate the influence of the SOA scheme complexity on SOA formation depending on NO_x levels.

3.1 Example of use of the MCM chemical scheme: beta-caryophyllene

The simulation of the beta-caryophyllene oxidation scheme (BCARY) in MCM v3.3.1 (Jenkin et al., 2012) is illustrated by a test case. Users may choose to fix oxidant concentrations, in order to reduce SOA chemical schemes of a specific VOC using the GENOA tool (Wang et al., 2022, 2023b). For that application, the presence of RO₂ from other sources is taken into account by introducing the notion of background RO₂. Indeed, as MCM reactions may involve radical reactions, SOA may be formed as a result of RO₂ species reacting with hydroperoxy radicals (HO₂) and other RO₂ radicals. Hence, RO₂-RO₂ reactions can be treated using the concept of a “RO₂ pool” that include RO₂ from all precursors. For example, a RO₂-RO₂ reaction can be written as follows:

NBCO2 -> 3.000E-01 BCBNO3
+ 7.000E-01 BCAL + 7.000E-01 NO2
KINETIC RO2 1 ARR 9.200E-14

where ARR indicates it is an Arrhenius type of reaction, “RO2 1” indicates that the NBCO2 species reacts with the total concentration of all RO₂ species in the RO₂ pool with index “1”. The number following “RO2” is required and specifies which RO₂ pool the species reacts with. It can be adjusted to refer to a different RO₂ pool if multiple pools are considered (the numbering starts from 1). All RO₂ pool species are defined in an input RO₂ species file, where each line lists an RO₂ species name and its corresponding RO₂ pool number, separated by a space.

The RO₂-RO₂ reaction can be treated in four different ways by adopting different options (tag_RO2) in the namelist:

• tag_RO2=0, compute no RO₂-RO₂ reaction.

• tag_RO2=1, compute RO₂-RO₂ reactions with only the produced RO₂ concentrations. The total concentrations of all RO₂ species recorded in the RO₂ species list are added up to compute the kinetic rates.

• tag_RO2=2, compute RO₂-RO₂ reactions with only the background RO₂ concentrations. For this option to work, the background concentrations should be provided for the gas-phase species RO2pool in the input constant concentration file.

• tag_RO2=3, compute RO₂-RO₂ reactions with both the produced and background RO₂ concentrations. In this case, the kinetic rate is calculated based on a sum of the background concentrations from the constant concentration file and the produced RO₂ concentrations of all RO₂ species.

In order to record the output concentrations of the RO₂ pool, a specific gas-phase species, RO2pool, has been added to the gas-phase species list.

The oxidation of BCARY for different tag_RO2 options is illustrated in Fig. 3. The influence of the RO₂ pool option is low for BCARY, because the RO₂ + RO₂ pathways, i.e., the part influenced by the background RO₂ pool, is not the dominant pathway for BCARY SOA formation.

Figure 3 Time evolution of BCARY-SOA yields (the fraction of th BCARY's reacted mass that is converted into SOA) simulated with the BCARY MCM chemical mechanism and different RO₂-RO₂ reaction options.

3.2 Influence of SOA scheme complexity on SOA formation for different NO_x levels

In these test cases, chemical schemes of different complexities are used and compared with respect to their simulated SOA concentrations for different NO₂ levels. Initial gas concentrations are taken to be representative of atmospheric conditions simulated by Sartelet et al. (2022) over Greater Paris in July 2009. The initial NO₂ concentration is taken equal to 26 µg m⁻³, and the simulated OH concentration varies between 1.2 and $\mathrm{7.6}\times {\mathrm{10}}^{-\mathrm{4}}$ µg m⁻³. Simulations are performed with initial NO₂ concentrations multiplied by 2 (52 µg m⁻³) or divided by 2 (13 µg m⁻³). Low values are used for initial organic and inorganic aerosol concentrations, to serve as a condensation core. To keep the particle dynamic as simple as possible, only one size section (representing aerosol with diameters between 0.1585 and 0.4 µm) is used. The simulation is run for about 5 h with outputs every minute.

The different chemical schemes are combined with the CB05 mechanism for oxidant and O₃ formation. Four sets of configurations with SOA schemes of different complexity are studied here. In each setup, the scheme CB05 is supplemented with different SOA schemes for monoterpenes, sesquiterpene, toluene, and naphthalene. The SOA schemes are sumarized in Table 2. In the first setup, the SOA schemes represent near-explicit mechanisms. In the second setup, the SOA schemes of toluene and monoterpenes are schemes reduced using GENOA. The third and fourth setups involve SOA schemes from the H²O mechanism built from chamber experiments, with the third setup including the simple autoxidation scheme of Chrit et al. (2017) for monoterpenes. For monoterpenes, the near-explicit scheme is from MCM and PRAM (Peroxy Radical Autoxidation Mechanism) (Roldin et al., 2019), and the intermediate scheme is the GENOA reduced mechanism from Wang et al. (2023b). For toluene, the near-explicit scheme is from Lannuque et al. (2023). It is supplemented with the molecular rearrangement with ring opening of a bicyclic peroxy radical (BPR) with an O-O bridge (Iyer et al., 2023), as implemented in Sartelet et al. (2024), and the GENOA reduced scheme is from Sartelet et al. (2024). For naphthalene, the near-explicit scheme is from Lannuque and Sartelet (2024). For sesquiterpene (β-caryophyllene), the near-explicit scheme is from MCM. Differences in the simulated concentrations and their temporal evolution among the implicit, quasi-explicit and GENOA-reduced schemes arise from the richer formation pathways included in the quasi-explicit and GENOA-reduced mechanims.

Couvidat et al. (2013a)Lannuque and Sartelet (2024)Couvidat et al. (2012)Lannuque et al. (2023)Sartelet et al. (2024)Sartelet et al. (2024)Couvidat et al. (2012)Roldin et al. (2019)Wang et al. (2023b)Wang et al. (2023b)Couvidat et al. (2012)

Table 2 SOA chemical mechanisms used in the test cases of Sect. 3.2.

The toluene SOA concentrations simulated with the different configurations and with the different initial NO₂ conditions are shown in Fig. 4. The toluene SOA concentrations increase with time. Using the near-explicit schemes, the toluene SOA concentrations are lower when the initial NO₂ concentrations are decreased by a factor 2. If the initial NO₂ concentrations are increased by a factor 2, the toluene SOA concentrations are lower in the first 3 h of the simulations, but they are higher after. As shown in the left panel of Fig. 4, the toluene SOA concentrations of the reduced and near-explicit schemes are very similar. However, as shown in the right panel of Fig. 4, the toluene SOA concentrations of the H²O scheme are lower than those of the near-explicit schemes. This was already observed in Sartelet et al. (2024). This is partly attributed to the impact of molecular rearrangement with ring opening of a BPR with an O-O bridge present in the near-explicit scheme, but which is missing in the H²O scheme. The evolution of the toluene SOA concentrations is similar between H²O and the near-explicit schemes when initial NO₂ concentrations are reduced: the toluene SOA concentrations decrease. SOA concentrations also decrease with the H²O scheme when initial NO₂ concentrations are doubled. This decrease persists through the whole simulation, contrary to the near-explicit scheme.

Figure 4 Time evolution of toluene SOA mass with different chemical schemes and for different NO₂ initial concentrations (Ref. denotes the reference NO₂ levels, NO₂ × 2 the simulation with doubled NO₂, and NO₂ $/$ 2 the simulation with halved NO₂). In (a), the near-explicit chemical scheme of Lannuque et al. (2023) supplemented by molecular rearrangement with ring opening of a bicyclic peroxy radical with an O-O bridge (Iyer et al., 2023; Sartelet et al., 2024) is compared to the H²O mechanism. In (b), the near-explicit scheme (Expl.) is compared to the reduced scheme (Rdc.) of Sartelet et al. (2024).

The monoterpene SOA concentrations simulated with the different schemes and in the different initial NO₂ conditions are compared in Fig. 5. With all mechanisms, monoterpene SOA concentrations increase as initial NO₂ levels decrease, because of the NO_x regime, as already noted in Wang et al. (2024). The higher concentrations may be attributed to RO₂-HO₂ reaction rates, which tend to yield more highly oxidized monoterpene degradation products. As shown in the left panel of Fig. 5, the SOA concentrations simulated with the reduced schemes and their evolution with NO₂ levels are similar between the reduced and the near-explicit schemes. As shown in the lower panel of Fig. 5, the monoterpene SOA concentrations are significantly lower with the H²O mechanism that ignores autoxidation. Furthermore, monoterpene SOA concentrations decrease as NO₂ levels decrease in constrast to the simulation with the near-explicit mechanism. Taking into account autoxidation in the H²O mechanism leads to higher monoterpene SOA concentrations (upper right panel of Fig. 5), but they are still underestimated compared to the near-explicit scheme. Similarly to the near-explicit scheme, the H²O mechanism including autoxidation gives higher SOA concentrations with lower NO₂ levels.

Figure 5 Time evolution of monoterpene SOA mass with different chemical schemes and for different NO₂ initial concentrations (Ref. denotes the reference NO₂ levels, NO₂ × 2 the simulation with doubled NO₂, and NO₂ $/$ 2 the simulation with halved NO₂). In (a), the near-explicit chemical scheme (Expl.) of MCM and PRAM is compared to the reduced scheme (Rdc.) of Wang et al. (2023b). The near-explicit scheme is compared to the H²O mechanism in (b) and to the H²O mechanism without autoxidation in (c).

The naphthalene SOA concentrations simulated with the different configurations and in the different NO₂ initial conditions are shown in the left panel of Fig. 6. The naphtalene SOA concentrations increase as the initial NO₂ levels increase. The naphtalene SOA concentrations simulated with the H²O scheme are much lower than those simulated with the near-explicit scheme. The naphtalene SOA concentrations simulated with the H²O scheme can even be considered negligible, except when the NO₂ levels are doubled. The difference between the concentrations simulated with the near-explicit and the H²O scheme reflect large uncertainties in the H²O scheme, which was built by Couvidat et al. (2013a) from a specific set of experiments from Chan et al. (2009) that are not representative of the full range of atmospheric conditions. Even though these experiments were corrected from wall-loss effects of particles, they were not corrected for gas wall losses. Furthermore, the behaviour of the H²O scheme may not fully extrapolate to all experimental setups, especially flow tubes where residence times are much shorter than in smog chambers.

Figure 6 Time evolution of naphthalene (left panel, a) and β-caryophyllene (right panel, b) SOA mass with different chemical schemes and for different NO₂ initial concentrations (Ref. denotes the reference NO₂ levels, NO₂ × 2 the simulation with doubled NO₂, and NO₂ $/$ 2 the simulation with halved NO₂). The near-explicit chemical scheme (Expl.) of Lannuque and Sartelet (2024) is compared to the H²O mechanism.

The sesquiterpene SOA concentrations simulated with the different configurations and in the different NO₂ initial conditions are shown in the right panel of Fig. 6. During the first two hours of the simulations, the sesquiterpene SOA concentrations increase with increasing NO₂ levels. This is not the case in the simulation with the H²O mechanism, where the SOA concentrations are higher with the reference NO₂ levels.

These test cases highlight the importance of near-explicit chemical schemes, which incorporate various chemical pathways for SOA formation. In contrast, simplified schemes based on chamber experiments may simulate different evolution patterns under changing environmental conditions compared to near-explicit schemes. However, the reduced schemes generated by GENOA, as demonstrated in the examples presented in this study, successfully replicate the evolution predicted by near-explicit schemes.

3.3 Comparisons to flow-tube or chamber experiments

3.3.1 Wall losses

To allow comparison of concentrations to flow-tube or chamber experiments, a wall loss module accounting for losses of particles and organic vapors was added in v2.0.

Wall losses of organic vapors

A reversible loss of organic vapors is accounted for SVOC species defined in the species list of aerosol with the following reaction from Huang et al. (2018):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{g}}}{\mathrm{d}t}}=-k_{\mathrm{w}}^{\mathrm{g}}\left(C_{\mathrm{g}}-{\displaystyle \frac{C_{\mathrm{w}}}{K_{\mathrm{w}}{C}_{\mathrm{wall}}}}\right)\end{array}$$

with C_g the gas phase concentrations, $k_{\mathrm{w}}^{\mathrm{g}}$ the wall losses kinetic of organic vapors, C_w the equivalent wall concentration and K_w the wall partitioning constant. C_wall (in µg m⁻³) has to be defined in the namelist. If C_wall is not defined or equal to 0, wall deposition of organic vapors is not taken into account in the model. One possible method to estimate C_wall and K_w is to use the method of Huang et al. (2018):

$$\begin{array}{}\text{(2)}& C_{\mathrm{wall}}=\mathrm{10}\mathrm{800}{\displaystyle \frac{S}{V}}\end{array}$$

with $S/V$ the surface on volume ratio of the chamber, and

$$\begin{array}{}\text{(3)}& K_{\mathrm{w}}={\displaystyle \frac{K_{\mathrm{p}}}{{\mathit{\gamma }}^{\mathrm{wall}}}}\times {\displaystyle \frac{M_{\mathrm{m}}}{\mathrm{200}}}\end{array}$$

with M_m, the molar mass, γ^wall the wall activity coefficient, K_p the ideal partitioning constant:

$$\begin{array}{}\text{(4)}& K_{\mathrm{p}}={\displaystyle \frac{\mathrm{760}\times \mathrm{8.202}\times {\mathrm{10}}^{-\mathrm{5}}}{{\mathrm{10}}^{\mathrm{6}}{M}_{\mathrm{m}}{P}_i^{\mathrm{0}}\left(T\right)}}\end{array}$$

with $P_i^{\mathrm{0}}$(T) the saturation vapor pressure at temperature T (calculated from the properties defined in aerosol_species_list_file. The wall activity coefficient is calculated by this formula calculated for Teflon film (fluorinated ethylene propylene):

$$\begin{array}{}\text{(5)}& {\mathit{\gamma }}^{\mathrm{wall}}={\mathrm{10}}^{\mathrm{3.299}}{K}_{\mathrm{p}}^{\mathrm{0.6407}}\end{array}$$

$k_{\mathrm{w}}^{\mathrm{g}}$ can be provided directly by the user. Alternatively, $k_{\mathrm{w}}^{\mathrm{g}}$ can be calculated based on the eddy diffusion coefficient k_e and the $S/V$ ratio with the following formula:

$$\begin{array}{}\text{(6)}& k_{\mathrm{w}}^{\mathrm{g}}={\displaystyle \frac{S}{V}}{\left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}\sqrt{k_{\mathrm{e}}{D}_{\mathrm{g}}}}}+{\displaystyle \frac{\mathrm{4}}{{\mathit{\alpha }}_{\mathrm{w}}\stackrel{\mathrm{‾}}{c}}}\right)}^{-\mathrm{1}}\end{array}$$

with D_g the gas-phase diffusion coefficient, $\stackrel{\mathrm{‾}}{c}$ the mean speed velocity of the compound and α_w the accommodation coefficient for condensation onto the wall (Lannuque et al., 2023):

$$\begin{array}{}\text{(7)}& {\mathit{\alpha }}_{\mathrm{w}}={\mathrm{10}}^{-\mathrm{2.744}}{K}_{\mathrm{p}}^{-\mathrm{1.407}}.\end{array}$$

Wall losses of particles

The wall losses of particles are irreversible and are represented as:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{p}}}{\mathrm{d}t}}=-k_{\mathrm{w}}^{\mathrm{p}}{C}_{\mathrm{p}}\end{array}$$

with C_p the particle phase concentrations, $k_{\mathrm{w}}^{\mathrm{p}}$ the wall losses kinetic of particles.

$k_{\mathrm{w}}^{\mathrm{p}}$ can either be provided directly by the user by defining the value of the parameter or calculated for each size bin. In that case, the user has to provide the eddy diffusion coefficient k_e, the $S/V$ ratio, the radius of the chamber R_chamber and the minimal wall loss rate k_wp0 due to electrostatic forces. $k_{\mathrm{w}}^{\mathrm{p}}$ is calculated based on Pierce et al. (2008):

$$\begin{array}{}\text{(9)}& k_{\mathrm{w}}^{\mathrm{p}}=k_{\mathrm{wp}\mathrm{0}}+{\displaystyle \frac{\mathrm{6}\sqrt{k_{\mathrm{e}}D}}{\mathit{\pi }{R}_{\mathrm{chamber}}}}{D}_{\mathrm{1}}\left({\displaystyle \frac{\mathit{\pi }{v}_{\mathrm{s}}}{\mathrm{2}\sqrt{k_{\mathrm{e}}D}}}\right)+{\displaystyle \frac{v_{\mathrm{s}}}{\mathrm{4}{R}_{\mathrm{chamber}}/\mathrm{3}}}\end{array}$$

with D₁ is the Debye function, D the Brownian diffusivity of the particle (dependent on the aerosol diameter) and the settling velocity of the particle (dependent on the aerosol diameter).

This wall-loss module is used in the following test cases that describe comparison of SOA formation in flow-tube or chamber experiments.

3.3.2 SOA formation from α-pinene oxidation

This test case (illustrated in Fig. 7) is based on the work presented in Wang et al. (2023b). The two simulations represent the formation of SOA in an idealized oxidative flow reactor (OFR) (Xavier et al., 2019) from the ozonolysis of α-pinene. The oxidation of α-pinene is simulated using the near-explicit mechanism from MCM and PRAM, and the mechanism reduced with GENOA. The saturation vapor pressure is calculated using UManSysProp (Topping et al., 2016), with the method of Myrdal and Yalkowsky (1997) to compute saturation vapor pressures and the method of Nannoolal et al. (2004) for the boiling point. The SOA concentrations and their temporal evolution are similar with the near-explicit and the reduced schemes. The quantity of SOA formed at the end of the simulation corresponds well to that measured by Xavier et al. (2019) (error below 4 %).

Figure 7 Time evolution of monoterpene (α-pinene) SOA mass using the near-explicit mechanism (Ref) and the reduced mechanism (RdC). The results are compared to the measurements reported in Xavier et al. (2019).

3.3.3 SOA formation from toluene and naphthalene oxidation under experimental conditions

For toluene oxidation, the test case is illustrated in Fig. 8). It is based on the work of Lannuque et al. (2023). It represents the formation of SOA in an UV-lamp irradiated aerosol flow tube into which ammonium sulfate seeds and a high amount of gaseous toluene and isopropylnitrite (IPN) are injected. Toluene oxidation is here represented with the detailed mechanism of Lannuque et al. (2023), with the addition of the molecular rearrangement with ring opening of a bicyclic peroxy radical with an O-O bridge (ipso-BPR, Iyer et al., 2023; Sartelet et al., 2024). Wall losses for stable gaseous species and the irreversible condensation of glyoxal and methylglyoxal are taken into account. IPN chemistry and experimental conditions representation are described in Lannuque et al. (2021). The dynamic of the condensation/evaporation on both organic and aqueous phases is solved without assuming thermodynamic equilibrium. The differences between the simulations presented in Fig. 8a lie in the simulation of the interactions in the particulate phase: ideal case (activity coefficients equal to 1), or interactions between uncharged molecules only (short range activity coefficients computed with UNIFAC), or interactions between uncharged molecules and inorganic ions in aqueous phase (short, medium and long range activity coefficients computed with AIOMFAC). In the test-case conditions, considering interactions between uncharged molecules enhances the condensation while considering those with inorganic ions in the aqueous phase limits it, almost offsetting the first interactions. The SOA concentrations simulated after about 13 min compare reasonably well to the observed one, i.e. 16 µg m⁻³.

Figure 8 Time evolution of toluene SOA mass under experimental conditions taking into account or not interactions in condensed phases (a, left) or using different SOA mechanisms under ideal conditions (b, right).

The near-explicit toluene SOA mechanism of Lannuque et al. (2023) is replaced by reduced toluene SOA mechanisms generated with GENOA Sartelet et al. (2024). Three reduced toluene SOA mechanisms are presented in Fig. 8b. The first mechanism is reduced from MCM (“SART24-1”). The second one (“SART24-2”) is the reduced mechanism of Lannuque et al. (2023), to which is added the molecular rearrangement with ring opening of a bicyclic peroxy radical with an O-O bridge (ipso-BPR, Iyer et al., 2023). The third mechanism is the basic H²O mechanism built from chamber experiments. These three mechanisms give similar results, although the mechanism reduced using the MCM (“SART24-1”) tends to underestimate the concentrations. The H²O scheme significantly overestimates the observed concentrations.

A similar test case was also created for naphthalene SOA based on the work of Lannuque and Sartelet (2024) with conditions similar to those of toluene. In this test case, illustrated in Fig. 9, considering the interactions and/or the viscosity in the particle limits the condensation of semi-volatile compounds, leading to less SOA formation than observed in the experiment (6.6 µg m⁻³). When accounting for viscosity, diffusion inside the particles limits the condensation of semivolatile species because the particles are semisolid. However, both simulations with and without viscosity remain very close to the experimental data, one slightly overestimating and the other slightly underestimating, so this difference is not considered meaningful in regard of uncertainties on simulation parameters (such as the initial diameter of particles). For toluene SOA, viscosity has only a minor influence, likely because their oxidation products are less multifunctional. The calculation of viscosity is detailed in the next section.

Figure 9 Time evolution of naphthalene SOA mass under different experimental conditions taking into account or not interactions and viscosity in condensed phases.

4 Advance gas/particle partitioning taking into account viscosity and particle-phase reactions

As illustrated in the test cases of the previous sections, organic compounds can condense on both an organic and an aqueous phase. The condensation may be altered by the viscosity of the particles, as well by particle-phase reactions. The calculation of viscosity as a function of the particle composition, and the types of intra-particle reactions are now described.

4.1 Revision of the numerical solver

In order to represent diffusion inside viscous particles, those are discretized into layers. To explicitly represent this diffusion, a high number of layers would be necessary for the model to be accurate. To reduce the number of layers, alternative approaches involve the use of an effective mass accommodation coefficient (Shiraiwa and Pöschl, 2021; Lakey et al., 2023), or a methodology that represents implicitly the diffusion inside particles (Couvidat and Sartelet, 2015). The latter method is used in SSH-aerosol v2.0. As in Couvidat and Sartelet (2015) and in SSH-aerosol v1.1, it does not represent the exchange of compounds between layers but assumes a characteristic time to reach the interface. As illustrated by Fig. 10, diffusion is represented by exchange fluxes between the different layers inside the particle and the interface layer. The main consequence of this simplification is that the diffusion of a compound is computed assuming identical affinity across all layers, where this affinity represents the compound’s activity relative to the other particle-phase constituents. It may therefore miss some entrapment effect (impossibility for compounds to cross a layer because of a lack of affinity).

Figure 10 Illustration of implicit treatment of diffusion inside particle.

The flux of diffusion for a specific layer, bin and species ($J_{\mathrm{diff},i}^{\mathrm{bin},\mathrm{layer}}$ in µg m⁻³ of air s⁻¹) is calculated as:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& J_{\mathrm{diff},i}^{\mathrm{bin},\mathrm{layer}}=\\ & k_{\mathrm{diffusion}}^{\mathrm{bin},\mathrm{layer}}(A_{\mathrm{p},i}^{\mathrm{interface},\mathrm{layer}}{\displaystyle \frac{K_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}{M}_{\mathrm{o}}^{\mathrm{bin},\mathrm{layer}}}{K_{\mathrm{p},i}^{\mathrm{bin},\mathrm{interface}}{M}_{\mathrm{o}}^{\mathrm{bin},\mathrm{interface}}}}\\ & -A_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}})\end{array}\end{array}$$

with $k_{\mathrm{diffusion}}^{\mathrm{bin},\mathrm{layer}}$ the parameter of diffusion rate from the interface to the layer (in s⁻¹), “interface” the index of the layer at the gas/particle interface, $K_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}$ (in m³ µg⁻¹) the partitioning constant, $M_{\mathrm{o}}^{\mathrm{bin},\mathrm{layer}}$ (in µg m⁻³ of air) the mass of the organic phase in the layer and $A_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}$ the particle concentrations of the species i in the layer (in µg m⁻³ of air).

An hybrid explicit numerical method is used to solve condensation/evaporation/diffusion in Couvidat and Sartelet (2015) and SSH-aerosol v1.1. One issue with this method is the need to distinguish cases where a compound condenses or evaporates very rapidly (typically handled with an equilibrium approach) from cases that require a dynamic resolution. The numerical method is replaced by an iterative implicit numerical (backward Euler), which is capable of simultaneously representing both types of cases:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& B_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}(t+\mathrm{\Delta }t,\mathrm{iter}+\mathrm{1})=\\ & {\displaystyle \frac{A_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}\left(t\right)+{\mathrm{Prod}}_i^{\mathrm{bin},\mathrm{layer}}(t+\mathrm{\Delta }t,\mathrm{iter})\mathrm{\Delta }t}{\mathrm{1}+k_{\mathrm{loss},i}^{\mathrm{bin},\mathrm{layer}}\mathrm{\Delta }t}}\end{array}\end{array}$$

With iter the iteration index, $A_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}$ the initial concentration of compound i in a specific bin and layer, $B_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}(t+\mathrm{\Delta }t,\mathrm{iter}+\mathrm{1})$ the restimated concentrations, Prod${}_i^{\mathrm{bin},\mathrm{layer}}$ the production rate and $k_{\mathrm{loss},i}^{\mathrm{bin},\mathrm{layer}}$ the loss kinetic rate. Prod${}_i^{\mathrm{bin},\mathrm{layer}}$ and $k_{\mathrm{loss},i}^{\mathrm{bin},\mathrm{layer}}$ include the production and loss due to diffusion, chemical reactions as well condensation/evaporation for the layer at the interface with air.

This method is iterative, as Prod${}_i^{\mathrm{bin},\mathrm{layer}}$ is determined for concentrations at time t+Δt. Iteration are performed until convergence is reached with the iterative solver described in Couvidat and Sartelet (2015):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& A_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}(t+\mathrm{\Delta }t,\mathrm{iter}+\mathrm{1})=\\ & f\times B_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}(t+\mathrm{\Delta }t,\mathrm{iter}+\mathrm{1})+(\mathrm{1}-f)\\ & \times A_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}(t+\mathrm{\Delta }t,\mathrm{iter})\end{array}\end{array}$$

where $B_{\mathrm{p},i}^{\mathrm{bin},\mathrm{layer}}(t+\mathrm{\Delta }t)$ is the concentration estimated with Eq. (11) and f is a weighting factor ranging from 0 to 1 determined to ensure convergence. Following Couvidat and Sartelet (2015), f is divided by a factor 2 (therefore decreasing the variations of concentrations between two iterations) each time a non-convergence loop is detected to ensure convergence.

The implicit solver is found to be faster and more stable than the previous method based on the explicit solver. Note that the same solver is applied for aqueous-phase concentrations (but with only one layer and without representing diffusion inside the aqueous phase).

4.2 Viscosity

The gas/particle mass transfer is strongly affected by particle viscosity, modeled by diffusion coefficients that are used to compute $k_{\mathrm{diffusion}}^{\mathrm{bin},\mathrm{layer}}$ such as:

$$\begin{array}{}\text{(13)}& {\displaystyle }{k}_{\mathrm{diffusion}}^{\mathrm{bin},\mathrm{layer}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{dif}}{\mathit{\alpha }}_{\mathrm{layer}}}}\text{(14)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{dif}}={\displaystyle \frac{R_{\mathrm{p}}}{{\mathit{\pi }}^{\mathrm{2}}{D}_{\mathrm{org}}}}\end{array}$$

with R_p the radius of the particle, D_org the diffusion coefficient inside the particle and τ_dif the characteristic time of diffusion. The higher the viscosity is, the lower the diffusion coefficient is. To illustrate this effect, in SSH-aerosol v1.1, the condensation of hydrophobic organic surrogates of different volatility was studied for different values of viscosity.

In SSH-aerosol v2.0, the viscosity may be estimated as a function of the particle composition, with the AIOMFAC-VISC algorithm (Gervasi et al., 2020). In this algorithm, the UNIFAC parameters are used to compute the viscosity. Therefore, the viscosity of a mixture of organic compounds can be computed as long as the decomposition in functional groups of the compounds is given. It should be noted that the aqueous-phase is considered to be inviscid. In the organic phase, by allowing the diffusion coefficient to vary with viscosity, D_org is no longer constant and therefore Eq. (13) is rewritten, as a function of the apparent diffusion coefficient from the interface to the considered layer $D_{\mathrm{org},\mathrm{app},i}^{\mathrm{bin},\mathrm{layer}}$:

$$\begin{array}{}\text{(15)}& k_{\mathrm{diffusion}}^{\mathrm{bin},\mathrm{layer}}={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{2}}{D}_{\mathrm{org},\mathrm{app},i}^{\mathrm{bin},\mathrm{layer}}}{R_{\mathrm{p}}{\mathit{\alpha }}_{\mathrm{layer}}}}\end{array}$$

$D_{\mathrm{org},\mathrm{app},i}^{\mathrm{bin},\mathrm{layer}}$ is computed with AIOMFAC-VISC by using cumulated concentrations between the interface and the considered layer in order to provide an average diffusion coefficient.

Following DeRieux et al. (2018), in the AIOMFAC-VISC algorithm, a fragility parameter D_fragility is used to compute pure compound viscosity η₀:

$$\begin{array}{}\text{(16)}& {\log }_{\mathrm{10}}\left({\mathit{\eta }}_{\mathrm{0}}\right)=-\mathrm{5.0}+\mathrm{0.434}{\displaystyle \frac{D_{\mathrm{fragility}}\cdot T_{\mathrm{0}}}{T-T_{\mathrm{0}}}}\end{array}$$

with T₀ the Vogel temperature, or ideal glass transition temperature:

$$\begin{array}{}\text{(17)}& T_{\mathrm{0}}={\displaystyle \frac{T_{\mathrm{g}}\cdot \mathrm{39.17}}{D_{\mathrm{fragility}}+\mathrm{39.17}}}\end{array}$$

with T_g the glass transition temperature, a characteristic temperature (computed as a function of the number of N, C, and O atoms) at which the mixture changes from a glassy solid state to a semi-solid state. Although D_fragility is set to 10 in Gervasi et al. (2020), a parameterisation of D_fragility as a function of the $\mathrm{O}/\mathrm{C}$ ratio is used, following Schmedding et al. (2020). A

The diffusion coefficient $D_{{\mathrm{diff}}_c}$ of the compound c is deduced from the viscosity (η), the viscosity of pure water (η₀) and the diffusion coefficient in pure water D_diff,water, following the parameterization of Maclean et al. (2021):

$$\begin{array}{}\text{(18)}& D_{{\mathrm{diff}}_c}=D_{\mathrm{diff},\mathrm{water}}\times {\left({\displaystyle \frac{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{water}}}{\mathit{\eta }}}\right)}^{\mathit{\xi }}\end{array}$$

where ξ is a parameter used to compare the radius of the diffusing molecule (R_diff) to the average radius of the compound mixture:

$$\begin{array}{}\text{(19)}& \mathit{\xi }=\mathrm{1}-A\exp \left(B{\displaystyle \frac{R_{\mathrm{diff}}}{R_{\mathrm{mean}}}}\right)\end{array}$$

with A and B equal to 0.73 and 1.79 respectively (Maclean et al., 2021), and R_mean is calculated as the average value of R_diff weighted by molar fractions. R_diff is calculated from the density and molar mass of the compound:

$$\begin{array}{}\text{(20)}& R_{\mathrm{diff}}={\left({\displaystyle \frac{\mathrm{3}{M}_c}{\mathrm{4}{N}_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{c}}\cdot \mathrm{1000}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\end{array}$$

with M_c the molar mass (in g mol⁻¹), N_a the Avogadro number and ρ_c the compound density (in kg m⁻³).

This methodology allows to calculate the diffusion coefficients for each species and each particle layer and therefore to vary with particle composition. The effect of changing composition on the condensation of a semivolatile compound is illustrated in Fig. 11. The condensation of an organic compound (POAmP) is investigated for different particle compositions, and thus for different viscosities. POAmP is represented using the default molecular structure proposed by Couvidat et al. (2012) in the H²O mechanism for primary compounds and model species with unknown molecular structures. In this representation, the activity coefficient is set to unity, ensuring that the condensation process is not influenced by non-ideality effects. At the start of the simulation, a low concentration of POAmP in the gas phase (0.01 µg m⁻³) was set in order to not affect the viscosity of the organic phase. The time evolution of POAmP concentrations was computed for different relative humidity (RH) and for different compositions of the aerosol onto which POAmP condenses. The compositions differ by the number of alcohol groups (OH), of carbons, and of various types of functional groups (ketone, hydroxyl, acid, and nitrate). The simulations show that an aerosol composed of compounds with 16 carbons and 4 alcohol groups is very viscous at RH = 10 % and that condensation of POAmP hardly occurs. However, at RH = 70 %, the aerosol can be considered as inviscid and the condensation of POAmP is non-limited by viscosity. Adding just one OH group has a strong effect on viscosity at RH = 10 %. With 16 carbons and 2 OH groups, equilibrium between the gas and particle phases of POAmP is reached in about one hour. However, it takes more than 100 h to reach equilibrium with an aerosol composed of 16 carbons and 3 OH groups. Adding, acid, hydroperoxide or nitrate groups affects the viscosity even more significantly. Especially, just by adding one nitrate group, the aerosol switches from inviscid to extremely viscous. By comparison, the effect on viscosity of increasing the number of carbons is low. If the molecules have 2 OH groups, equilibrium is reached in less than 1 h with less than 16 carbons. Time to equilibrium increases progressively with the number of carbons. It takes about 100 h with 20 carbons.

Figure 11 Concentration of POAmP condensing onto 5 µg m⁻³ of a compound with 16 carbons C₁₆ as a function of time. The different graphs show the effect of changing the composition of the aerosol or humidity.

4.3 Particle-phase reactions

Reactions can occur inside the particles. These reactions can be activated or deactivated by modifying the file indicated by reaction_soap_file in the section &physic_condensation of the namelist. Different types of reactions can be accounted for. Examples on how to use this reaction file are provided in the supplementary material (user manual). Currently 5 types of particle-phase reactions are taken into account: 1st order irreversible reaction, 2nd order reversible oligomerization, hydratation of aldehydes, reactions of organic compounds with inorganic ions and a specific paramterization for pinonaldehyde oligomerization.

4.3.1 1st irreversible order reaction

An irreversible 1st order reaction is a reaction that transforms X into Y is represented with a kinetic parameter k such as:

$$\begin{array}{}\text{(21)}& X\stackrel{k}{\to }Y\end{array}$$

This type of reaction can account for catalysis by water and pH. It allows to set a lifetime of semi-volatile compounds in particles, e.g. a lifetime of 1 h corresponds to a kinetic rate of about $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹. This type of reaction could also be used to represent hydrolysis (for example the hydrolysis of nitrate groups into hydroxyl groups).

4.3.2 Type 1: Bulk oligomerization

In most modeling studies, a parameterization for oligomerization based on a simple first-order irreversible reaction is used (Carlton et al., 2008; Lemaire et al., 2016). However, most oligomerization pathways should involved second-order reversible reactions (like esterification, hemiacetalization, aldolization, peroxyhemiacetalization) catalyzed under acidic conditions and unfavored by humid conditions (Couvidat et al., 2018b).

In SSH-aerosol v2.0, oligomerization is based on the parameterization developped by Couvidat et al. (2018b). It is represented by a reversible process, unfavored by humid conditions, with a single reaction estimating if a compound is present as a monomer or as an oligomer. Oligomers from one species are therefore represented by a single oligomer species. The parameterization is termed “Bulk oligomerization” (as it does not explicitly represent all possible oligomers) and can be used to represent oligomerization in a complex mixture of compounds. It can be schematized with the following reaction (A representing the monomer and A_oligo representing any monomer block present inside the oligomers):

$$\begin{array}{}\text{(22)}& A\underset{k_{\mathrm{reverse}}}{\stackrel{k_{\mathrm{oligo}}}{\leftrightarrow }}{\displaystyle \frac{\mathrm{1}}{m_{\mathrm{oligo}}}}{A}_{\mathrm{oligo}}\end{array}$$

with A a monomer compound, A_oligo the monomer blocks of compound A inside oligomers, m_oligo the number of monomer blocks inside oligomers, k_oligo the kinetic rate parameter of oligomerization and k_reverse the kinetic rate parameter of the reverse reaction.

The net flux of oligomerization J_oligo is computed with the following equations:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{J}_{\mathrm{oligo}}& =-{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{A},\mathrm{monomer}}}{\mathrm{d}t}}\\ & =k_{\mathrm{oligo}}{a}_{\mathrm{A},\mathrm{monomer}}-k_{\mathrm{reverse}}{a}_{\mathrm{A},\mathrm{oligomer}}\end{array}\end{array}$$

with X_A,monomer the molar fraction of compound A, a_A,monomer the activity on a molar fraction basis of compound A and a_A,oligomer the activity on a molar fraction basis of the oligomer formed from compound A.

The kinetic rate of oligomerization k_oligo is computed as follows:

$$\begin{array}{}\text{(24)}& k_{\mathrm{oligo}}=k_{\mathrm{oligo}}^{\max }\sum _i{a}_{i,\mathrm{monomer}}\end{array}$$

with $k_{\mathrm{oligo}}^{\max }$ the maximum kinetic rate parameter for oligomerization, and ${\sum }_i{a}_{i,\mathrm{monomer}}$ the sum of monomer activities.

The reverse kinetic rate parameter is estimated from the equilibrium constant for oligomerization (due to esterification or a similar oligomerization mechanism), which may be written as

$$\begin{array}{}\text{(25)}& (K_{\mathrm{oligo}}^{\mathrm{eq}}{)}^{m_{\mathrm{oligo}}-\mathrm{1}}={\displaystyle \frac{a_{\mathrm{A},\mathrm{oligomer}}(a_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}{)}^{m_{\mathrm{oligo}}-\mathrm{1}}}{a_{\mathrm{A},\mathrm{monomer}}({\sum }_i{a}_{i,\mathrm{monomer}}{)}^{m_{\mathrm{oligo}}-\mathrm{1}}}}\end{array}$$

with $a_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ the activity of water on a molar basis, m_oligo is the number of oligomer blocks in oligomers, $K_{\mathrm{oligo}}^{\mathrm{eq}}$ is the the equilibrium oligomerization constant. $k_{\mathrm{oligo}}^{\max }$, m_oligo and $K_{\mathrm{oligo}}^{\mathrm{eq}}$ have to be provided by the user for each reactions.

To illustrate the formation of oligomers, the Platt test case, describing the SOA formation from exhaust emissions from a Euro 5 gasoline car in a chamber (Platt et al., 2013), modified by considering isoprene and NH₃ emissions (Sartelet et al., 2020), is run, allowing isoprene SOA to form oligomers based on the parameters proposed by Couvidat et al. (2018a) ($k_{\mathrm{oligo}}^{\max }=\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, m_oligo=3.35 and $K_{\mathrm{oligo}}^{\mathrm{eq}}=\mathrm{2.94}$). The evolution of the concentrations of isoprene SOA, monomers, and oligomers in the particles as a function of time for these simulations is shown in Fig. 12. Although, the SOA concentration does not change significantly when oligomerization is activated, the SOA composition is strongly impacted. After 1 h, almost all the SOA mass is constituted by monomer species. However, after 5 h, monomer species constitute only 50 % of isoprene SOA.

Figure 12 Time evolution of isoprene SOA, monomers, and oligomers in the particle phase.

4.3.3 Type 2: Hydratation of aldehydes

Aldehydes (and more generally carbonyl compounds) can undergo hydratation. In that case, the carbonyl group is transformed into a diol, and hydratation is represented as an equilibrium between two species, the aldehyde and the diol, using an hydratation constant K_hydratation:

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\left[\mathrm{diol}\right]}{\left[\mathrm{aldehyde}\right]}}=K_{\mathrm{hydratation}}{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{diol}}^{\mathrm{\infty }}{\mathit{\gamma }}_{\mathrm{aldehyde}}}{{\mathit{\gamma }}_{\mathrm{diol}}{\mathit{\gamma }}_{\mathrm{aldehyde}}^{\mathrm{\infty }}}}{a}_{\mathrm{w}}\end{array}$$

with a_w the activity of water (equal to RH at equilibrium). γ_diol and γ_aldehyde corresponds to the activity coefficient of the diol and the aldehyde. They are normalized by the activity coefficient at infinite dilution γ^∞, as the hydratation constant is generally measured in the aqueous phase.

The equilibrium constant between the diol and the aldehyde is generally quite low (0.57 for butanal according to GECKO-A, Camredon et al., 2007). However, glyoxal can be hydradated twice with a very high constant as illustrated in Fig. 13.

Figure 13 Illustration of glyoxal conversion into polyols due to hydration. The hydratation constants are taken from Ervens and Volkamer (2010).

The influence of hydratation on glyoxal absorption was tested by running the Platt isoprene test case including ammonia emission with and without hydratation. While the concentrations of glyoxal in the particle phase without hydratation are negligible (below 10⁻⁸ µg m⁻³), the total concentrations of glyoxal and its hydrate are higher and reach $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{3}}$ µg m⁻³.

4.3.4 Type 3: Reaction of organic compounds with inorganic ions

Organic compounds in the aqueous can react with inorganic ions to form other compounds. For example, the compound IEPOX (epoxide formed by the oxidation of isoprene) can react with SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and HSO${}_{\mathrm{4}}^{-}$ ions to form organosulfate (Couvidat et al., 2013b) or with NO${}_{\mathrm{3}}^{-}$ to form organonitrate. This type of reaction can account for catalysis by water (reaction rate multiplied by the activity of water) and pH (reaction rate multiplied by the activity of H⁺ ion). The Fig. 14 illustrates the formation of SOA from IEPOX based on Couvidat et al. (2013b) due to intra-particle reactions for different amounts of ammonia, and hence different pH of the particle. In the absence of ammonia (very acidic particles), the formation of organosulfate is very important. For partially neutralized sulfate, the formation of organosulfate decreases but is compensated by the formation of methyltetrols due to the hydrolysis of IEPOX. When sulfate is almost totally neutralized, the concentrations of SOA formed by the hydrolysis of IEPOX decreases significantly and is strongly dominated by the formation of methyltetrols.

Figure 14 Formation of SOA (top, a), organosulfate (top right, b) and methyltetrols (bottom left, c) due to the reaction of 10 µg m⁻³ IEPOX in the particle phase (in the presence of 10 µg m⁻³ of HNO₃ at the beginning of the reaction and 3 µg m⁻³ of sulfate) for a different ratio of ammonia to sulfate.

4.3.5 Type 4: Oligomerization of pinonaldehyde

Concerning organic reactions in the particles, oligomerization of pinonaldehyde (BiA0D, species in H²O) may be considered by the user. It corresponds to the parameterisation detailed in Pun and Seigneur (2007) and Couvidat et al. (2012), where BiA0D undergoes oligomerization via an effective Henry's law constant

$$\begin{array}{}\text{(27)}& H_{\mathrm{eff}}=H\left(\mathrm{1}+\mathrm{0.1}{\left({\displaystyle \frac{a\left(H^{+}\right)}{{\mathrm{10}}^{-\mathrm{6}}}}\right)}^{\mathrm{1.91}}\right)\end{array}$$

where H_eff is the effective Henry's law constant of BiA0D, H is the monomer Henry's law constant of BiA0D, and a(H⁺) is the activity of protons in the aqueous phase.

It should be noted that this parameterization (treating oligomerization with an equilibrium constant) has been initially developed to increase SOA formation from pinonaldehyde (Pun and Seigneur, 2007), and was present in SSH-aerosol v1.1. However, Couvidat et al. (2013b) showed that the uptake of pinonaldehyde on acidic particle is very slow and should not be treated with an equilibrium constant.

4.4 Irreversible condensation of methylglyoxal and glyoxal

Once uptaken into the aqueous phase of the particles, dicarbonyls can be oxidized by OH radicals present in this aqueous phase to irreversibly form highly soluble and low-volatility compounds (Hu et al., 2022). Empirical reactions representing this irreversible condensation pathway for glyoxal and methylglyoxal has been added in SSH-aerosol v2.0. The term “irreversible condensation pathway” is used in contrast to the classical “reversible condensation pathway” influenced by hydration and isomerization. The irreversible condensation pathway is represented by the following reaction:

$$\begin{array}{}\text{(28)}& \mathrm{DICARB}\stackrel{k_{\mathrm{irdi}}}{\to }\mathrm{IRDICARB}\end{array}$$

where DICARB is the semi-volatile species studied (glyoxal or methylglyoxal) and IRDICARB its non-volatile counterpart via irreversible partitioning. The kinetics k_irdi of this reaction follows the parameterization of Curry et al. (2018):

$$\begin{array}{}\text{(29)}& k_{\mathrm{irdi}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}v\times \mathit{\gamma }\times A_{\mathrm{surf}}\end{array}$$

with:

$$\begin{array}{}\text{(30)}& {\displaystyle }v=\sqrt{{\displaystyle \frac{\mathrm{8}RT}{\mathit{\pi }{M}_{\mathrm{dicarb}}}}}\text{(31)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\gamma }}}={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}+{\displaystyle \frac{v}{\mathrm{4}RT{H}_{\mathrm{dicarb}}^{\mathrm{eff}}\sqrt{k^{\mathrm{l}}{D}_{\mathrm{aq}}}}}\times {\displaystyle \frac{\mathrm{1}}{\coth q-\mathrm{1}/q}}\text{(32)}& {\displaystyle }q={\displaystyle \frac{R_{\mathrm{p}}}{\sqrt{\frac{D_{\mathrm{aq}}}{k^{\mathrm{l}}}}}}\end{array}$$

where γ is the reactive uptake coefficient, v and M_dicarb are the gas-phase thermal velocity and the molar mass of the DICARB species, A_surf is the aerosol surface, R is the universal gas constant, T is the temperature in Kelvins, $H_{\mathrm{dicarb}}^{\mathrm{eff}}$ is the effective Henry's law constant, q is a parameter that account for in-particle diffusion limitation, R_p is the particle radius, k^l is the first order loss rate due to reaction with OH radicals in the aqueous phase. The aqueous-phase diffusion coefficient D_aq is fixed to 10⁻⁹ m² s⁻¹, which is typical for small organics (Bird et al., 2006), and the mass accommodation coefficient α to 0.02, similar to that of formaldehyde uptake to water (Jayne et al., 1992). Henry constants of DICARB are adapted from Sander (2015) to match with Hu et al. (2022) observations and the Henry constant of OH used to estimate OH aqueous concentration from gaseous ones is adapted from Sander (2015) leading to aqueous concentration in agreement with typical mean value in the atmosphere from Herrmann et al. (2010). Further details and description of the parameters are available in Curry et al. (2018).

4.5 Coupling between oganics and inorganics

The standard configuration of SSH-aerosol is based on two separate thermodynamic modules: SOAP and ISORROPIA. ISORROPIA is called first to compute water, pH, and inorganics. Then, SOAP is called to compute the concentrations of organics. In this system, organics do not affect the computation of pH or the concentrations of inorganics. This requires solve organics and inorganics together.

SOAP was modified in order to be able to solve the equilibrium and the dynamics of inorganic and organic aerosol simultaneously. Organics can affect the formation of inorganics by:

1. Influencing the condensation of water, which in turn affects the absorption of HNO₃ and NH₃.

2. Influencing the computation of pH as the dissociation of acid organics is accounted for.

3. Influencing the computation of molalities (concentrations of ions in mol kg⁻¹). In ISORROPIA, molalities are computed as mol kg⁻¹ of water. In SOAP, all solvent molecules are accounted for (water but also organics).

SOAP is based on the AIOMFAC thermodynamic module that accounts for numerous ions, H⁺, OH⁻, NO${}_{\mathrm{3}}^{-}$, HSO${}_{\mathrm{4}}^{-}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, HCO${}_{\mathrm{3}}^{-}$, CO${}_{\mathrm{3}}^{\mathrm{2}-}$, Cl⁻, Na⁺, K⁺, Mg²⁺, Ca²⁺. Inorganic aerosols are assumed to be thermodynamically metastable (inorganic compounds are always present in the aqueous phase of particles and no solid salts are formed) with the exception of CaCO₃ that is assumed to have a very low solubility. Henry's law constant are taken from Fountoukis and Nenes (2007). The iterative solver in Eq. (12) is used to solve iteratively the system of equations. f in Eq. (12) is constrained in order to ensure that activity coefficients of inorganics and pH do not change significantly between two iterations in order to ensure stability and convergence.

Several test cases, such as the Platt test case with isoprene and NH₃, have been run using ISORROPIA or SOAP for inorganics. Using SOAP, the simulated concentrations are very similar to the results obtained using ISORROPIA, as shown in the Supplement (user manual). To illustrate some of the capabilities of the SOAP model for both inorganic and organic thermodynamic, a test case simulates the evolution of CaCO₃ in the presence or in the absence of HNO₃ with an equilibrium approach or a dynamic approach. As shown in Fig. 15, when CaCO₃ is not exposed to HNO₃ the particle does not evolve. Especially, it does not absorb water. This is due to the fact that CaCO₃ is currently the only particle for which delisquescence is accounted for due to its very low saturation value that would correspond to a delisquescence relative humidity close to 100 %. When [Ca²⁺ (mol kg⁻¹ of solvents)][CO${}_{\mathrm{3}}^{\mathrm{2}-}$ (mol kg⁻¹ of solvents)] > K (with K the equilibrium constant), a solid CaCO₃ is formed. When exposed to HNO₃, around 2 molec. of NO${}_{\mathrm{3}}^{-}$ will replace 1 molec. of CO${}_{\mathrm{3}}^{\mathrm{2}-}$ to maintain electroneutrality. It leads to the presence of the nitrate and water inside the particle, and part of Ca²⁺ that was previously associated with CO${}_{\mathrm{3}}^{\mathrm{2}+}$ is transferred to the aqueous phase. Therefore, due to the condensation of HNO₃, the particle grows bigger in size and becomes hygroscopic.

Figure 15 Time evolution of CaCO₃ particles in presence or in absence of HNO₃ with an equilibrium approach (a) or a dynamic approach (b).

5 Representation of ultrafine particle concentrations

Ultrafine particle (UFP) concentrations are strongly influenced by nucleation, coagulation and condensation of extremely-low volatile compounds. For coagulation, the discretized equations governing aerosol dynamics involve partition coefficients that account for the possibility that the coagulation of particles from two given size sections may produce particle sizes spanning multiple sections. The size mesh is assumed to remain fixed over time, so that these partition coefficients can be pre-computed. In previous model version, the coagulation partition coefficients were calculated using a Monte-Carlo approach. In v2.0, to gain CPU time when calculating coagulation partition coefficients, an explicit formulation was derived for particles that are internally mixed, as detailed in Jacquot and Sartelet (2025). If the mixing state is resolved, the coagulation partition coefficients are solved, as in v1.1, i.e. using a Monte-Carlo approach. Because the explicit formulation has negligible computational cost compared with the Monte Carlo approach, only the explicit formulation is retained for internally-mixed particles in v2.0.

Given the strong impact of coagulation on the growth of ultrafine particles, it should always be included in their modelling. Consequently, the particle-diameter limits are held constant over the course of the simulations. Because nucleation corresponds to the formation of nanometric particles from gas precursors, the lowest particle diameter bound of the size distribution should be about 1 nm. Different parameterizations are implemented to represent nucleation, as detailed in Sect. 5.1. Numerically, nucleation and condensation are always solved simultaneously, because nucleation and condensation are competing processes. However, the user can choose to split coagulation from other processes (variable splitting set to 1 in the namelist). To accurately represent the growth of UFP, the condensation of non-volatile compounds can be solved dynamically for all size sections. The influence of non-volatile compounds on the growth of UFP was illustrated in a test case in v1.1 and in the supplementary material (user manual). Semi-volatile inorganic and organic compounds are less likely to condense on UFP, because of the Kelvin effect. Therefore, to gain computational time, the condensation/evaporation of semi-volatile compounds may be solved by assuming bulk thermodynamic equilibrium between the gas and particle phases. An algorithm is then applied to redistribute the bulk aerosol concentrations across the size sections in a thermodynamically consistent way, explicitly accounting for the Kelvin effect, which increases the effective vapor pressure over smaller particles.

5.1 Nucleation

Several nucleation parameterizations are implemented. They can be used separately or together. An example of their use in 3D is illustrated in Sartelet et al. (2022). The nucleation parametrizations differ in the gas precursors (sulfuric acid, ammonia, extremely-low volatile organic compounds) and/or in the formulation of the parameterization.

• binary: water and sulfuric acid with the parameterisations of Vehkamaki et al. (2002) or Kuang et al. (2008);

• ternary: water, sulfuric acid and ammonia with the parameterisation of Napari et al. (2002) or Merikanto et al. (2007, 2009). To avoid artificially large nucleation rates in the parameterisation of Napari et al. (2002), a maximum nucleation rate of 10⁶ # particles cm⁻³ is set. A scaling factor of the nucleation rate can be applied (parameter scal_ternary in the namelist).

• heteromolecular: sulfuric acid, and extremely-low volatile organic compounds from monoterpene autoxidation with the parameterisation of Riccobono et al. (2014).

• organics: the user can choose the list of extremely-low volatile organic compounds for which nucleation is taken into account. The nucleation rate can be ajusted by specifying a exponent of the power law (parameter nexp_org) and scaling factor (parameter scal_org):

  $$\begin{array}{}\text{(33)}& J={\mathrm{scal}}_{\mathrm{org}}[\mathrm{Organics}{]}^{{\mathrm{nexp}}_{\mathrm{org}}}.\end{array}$$

Nucleated particles are assigned to the first size bin, whose lower bound should be set to about 1 nm to represent freshly formed clusters. The composition of the nucleation-size bin is determined by the precursor species involved in nucleation, and it is estimated based on the properties of these species.

5.2 Growth of ultrafine particles

To illustrate the growth of nucleated nanoparticles, the nucleation test case presented in Sartelet et al. (2006) is simulated. The initial distribution corresponds to the hazy conditions of Seigneur et al. (1986). The sulfuric acid production rate is 0.825 µg m⁻³ h⁻¹, the temperature is 288.15 K and the relative humidity is 60 %. The particles are initially assumed to be made of 70 % sulfate and 30 % ammonium. The initial gas phase ammonia concentration is taken to be 8 µg m⁻³. The concentrations of gas-phase ammonia and particulate-phase ammonium evolve with time due to both condensation/evaporation and ternary nucleation. The simulation is run for 1 h with output every 60 s. Figure 16 shows the size distribution of particles at the initial time and after 1 h. The nucleated particles grow under the effect of condensation of sulfuric acid and ammonia and coagulation. This growth can be accelerated by the presence of extremely low-volatility organic compounds (ELVOCs), such as the SSH-aerosol surrogate “Monomer”, underlying the importance of an accurate modelling of the formation of ELVOCs for the growth of ultra-fine particles.

Figure 16 Size distribution of particles in the nucleation test case. Under hazy conditions, nanoparticles nucleated from water, sulfuric acid and ammonia grow to larger size. A simulation adding ELVOCs (Monomer) to the initial conditions is performed to illustrate their effect on the growth of ultrafine particles. Number (left panel, a) and volume (right panel, b) concentrations.

5.3 Bulk equilibrium and Kelvin effect

To decrease CPU time associated with the treatment of condensation and evaporation, thermodynamic equilibrium between the gas and particle phases can be assumed. In that case, the bulk equilibrium is first computed, and then the bulk aerosol concentrations need to be distributed among the different size sections.

The redistribution was changed in version 2.0 to take into account the Kelvin effect, and for inorganics, the neutralization of acids by positive inorganic ions (sodium and ammonium). The amount of ammonium that condenses to neutralize sulfate is determined separately from the amount of ammonium that condenses to neutralize other compounds, such as nitric acid. Then the amount of ammonium that neutralizes sulfate is redistributed amongst the size sections, depending on the sulfate concentration of each section. The amount of ammonium that condenses with other semi-volatile compounds such as nitric acid is redistributed amonst the size section using a formulation that depends on the Kelvin effect and the condensation rate (depending on the gas-phase diffusivity of species and the diameters of particles).

To illustrate the role of the Kelvin effect versus the sulfate neutralization, the urban condition test case of Seigneur et al. (1986) is modified by adding NH₃ in the initial conditions, and by taking into account only condensation/evaporation processes (nucleation and coagulation are not considered). The initial particle concentration is assumed to be made of sulfate only, and the initial concentration of NH₃ is assumed to be 100 µg m⁻³. The simulation is run for 12 h.

Four simulations are run: in 2 simulations condensation/evaporation is computed dynamically with or without taking into account the Kelvin effect, in the other 2 simulations thermodynamic equilibrium is assumed. The number and volume size distribution simulated with the four numerical algorithms are compared in Fig. 17. The four numerical algorithms give similar number concentrations, because NH₃ condenses to neutralize the sulfate independently of the Kelvin effect.

Figure 17 Condensation and thermodynamic equilibrium test case for the condensation of ammonia neutralizing sulfate. Number (left panel, a) and volume (right panel, b) concentrations.

In the previous test cases, about 40 µg m⁻³ of ammonia is used to neutralize sulfate and the remaining concentration stays in the gas phase. To compare the effect of the Kelvin effect between dynamical calculation and thermodynamic assumption, nitric acid is added in the initial gas-phase concentration. Nitric acid condenses with the remaining ammonia to form ammonium nitrate, as shown in Fig. 18. As seen in the comparison of the number size distribution, the particles of diameters below about 0.01 µm do not grow by the condensation of ammonium nitrate if the Kelvin effect is taken into account, while these particles grow to larger diameters by condensation of ammonium nitrate if the Kelvin effect is not taken into account. The algorithm used to redistribute ammonium nitrate amongst size section is effective at limiting this condensation on particles of small diameters. This test case is also presented for organic aerosols in the supplementary material (user manual).

Figure 18 Condensation and thermodynamic equilibrium test case for the condensation of ammonia neutralizing sulfate, and for the condensation of ammonium nitrate. Number (left panel, a) and volume (right panel, b) concentrations.

6 Coupling to 3D models

SSH-aerosol can be coupled with external (3D) tools using the shared library (libssh-aerosol.so) available after compilation. A prototype of the typical workflow is described in the supplementary material (user manual). A “tool” repertory is added in v2.0. It contains routines that may be useful for the numerical representation of aerosols and the setup of simulations. For example, a routine is defined to partition, i.e. split, a size section into several size sections, while conserving both mass and number concentrations. This may be used to refine the number of size sections for initial, boundary conditions or emissions in 3D-simulations. Also a routine is provided to discretize into sections an aerosol distribution known as the sum of lognormal modes. This routine was used to build some of the test cases, and can be used in 3D simulations to link a lognormal representation of particles to the sectional one of SSH-aerosol.

SSH-aerosol has been coupled with various types of 3D model to simulate air quality. At the regional scale, it was coupled with the transport chemistry models Polyphemus/Polair3D (Sartelet et al., 2022) and CHIMERE (Menut et al., 2024) to study, in particular, the formation of ultrafine particles (Sartelet et al., 2022; Park et al., 2025; Jacquot and Sartelet, 2025) and the influence of organic aerosol patterns of different complexity (Wang et al., 2024; Sartelet et al., 2024). On a local scale, it was coupled to the MUNICH street network model (Kim et al., 2022), as well as to the OpenFoam and Code_Saturne computational fluid dynamic models (Lin et al., 2023), enabling the study of the secondary formation of pollutants in streets, in particular the role of ammonia emissions from vehicles (Lugon et al., 2021; Lin et al., 2023), the formation of condensables (Wang et al., 2023a) and secondary aerosols (Lin et al., 2023, 2024a, b; Sartelet et al., 2024), ultrafine particles (Park et al., 2025). The use of SSH-aerosol in models of different scales has made it possible to model primary and secondary pollutants over a whole city, from the regional down to the street scale (Lugon et al., 2021, 2022; Sarica et al., 2023a; Park et al., 2025; Sartelet et al., 2024; Squarcioni et al., 2025) and to study at both the regional and local scales the impact of certain processes that are often missing from modelling studies, such as asphalt emissions (Sarica et al., 2023b) or trees (Maison et al., 2024a, b), as well as to assess the impact of scenarios for changes in mobility and traffic (Lugon et al., 2022; Sarica et al., 2024).

One complexity of the use of SSH-aerosol inside 3D air quality models is due to the type of numerical solver used to calculate concentrations, especially if the model does not rely on a splitting approach where aerosol formation is not separated from other processes. However, Couvidat et al. (2025) proposed different methods to couple SSH-aerosol to 3D models depending on the type of numerical solvers.

7 Conclusions

The advancements in the SSH-aerosol v2.0 model represent a significant leap forward in the modeling of gas-phase chemistry and aerosol dynamics. By integrating state-of-the-art modules such as SCRAM, SOAP, and H²O, the model offers unparalleled flexibility and precision in simulating the formation and aging of aerosols, describing the properties of particles from the size and mixing-state distribution to the detailed particle composition. New features enhances the model's ability to capture complex atmospheric processes with high fidelity: incorporation of near-explicit chemical schemes, reduced mechanisms generated using GENOA, possibility of defining organic compounds according to their SMILES code or decomposition into functional groups in order to model the partitioning between aqueous and organic phases and the role of non-ideality, particle-phase reactions, composition-dependent viscosity, and fast algorithms to simulate the evolution of ultrafine particles.

The test cases included in the study underscore the robustness and versatility of SSH-aerosol v2.0. Whether modeling the effects of viscosity on gas-particle partitioning, the impact of radical chemistry on SOA formation, the role of NO_x levels in influencing secondary organic aerosol yields, or the formation and evolution of ultrafine particles, the model demonstrates its capability to handle a wide range of atmospheric scenarios. The ability to couple with 3D models further extends its applicability, making it suitable for both localized studies and regional-scale simulations. By bridging the gap between simplified and near-explicit chemical representations, it provides the means to explore and predict the behavior of aerosols under varying environmental conditions with greater accuracy.

Moreover, the improvements in computational efficiency, such as the direct interpretation of chemical reactions, the revision of numerical solvers and faster coagulation algorithms, ensure that the model remains practical for extensive applications in 3D. The addition of a wall-loss module, the modularity of the input data, such as the possibility to setup parameters for nucleation rates, enhance its usability for comparison to chamber experiments.

Overall, SSH-aerosol v2.0 stands out as a comprehensive tool for advancing our understanding of atmospheric chemistry and aerosol dynamics. Future developments could further enhance the model's capabilities, particularly coupling to molecular dynamic model for nucleation.

Appendix A

The table below lists all the options available in the namelist. Some parameters are not compulsory. If the parameter is not provided, a default value is used. Some parameters are indicated as “Optional”, they are not used unless a specific option is activated. In that case, the parameter must be provided. For example, parameter N_groups is only used when tag_external=1.

Vehkamaki et al. (2002)Kuang et al. (2008)Napari et al. (2002)Merikanto et al. (2007, 2009)

Table A1 Options available in the namelist of SSH-aerosol. Some parameters are not compulsory. If the parameter is not provided, a default value is used. Some parameters are indicated as “Optional”, they are not used unless a specific option is activated. In that case, the parameter must be provided. For example, parameter N_groups is only used when tag_external=1.