Non-equilibrium interplay between gas–particle partitioning and multiphase chemical reactions of semi-volatile compounds: mechanistic insights and practical implications for atmospheric modeling of polycyclic aromatic hydrocarbons

Polycyclic aromatic hydrocarbons (PAHs) are carcinogenic air pollutants. The dispersion of PAHs in the atmosphere is influenced by gas–particle partitioning and chemical loss. These processes are closely interlinked and may occur at vastly differing timescales, which complicates their mathematical description in chemical transport models. Here, we use a kinetic model that explicitly resolves mass transport and chemical reactions in the gas and particle phases to describe and explore the dynamic and non-equilibrium interplay of gas–particle partitioning and chemical losses of PAHs on soot particles. We define the equilibration timescale τ_eq of gas–particle partitioning as the e-folding time for relaxation of the system to the partitioning equilibrium. We find this metric to span from seconds to hours depending on temperature, particle surface area, and the type of PAH. The equilibration time can be approximated using a time-independent equation, ${\mathit{\tau }}_{\mathrm{eq}}\approx \frac{\mathrm{1}}{k_{\mathrm{des}}+k_{\mathrm{ads}}}$, which depends on the desorption rate coefficient k_des and adsorption rate coefficient k_ads, both of which can be calculated from experimentally accessible parameters. The model reveals two regimes in which different physical processes control the equilibration timescale: a desorption-controlled and an adsorption-controlled regime. In a case study with the PAH pyrene, we illustrate how chemical loss can perturb the equilibrium particulate fraction at typical atmospheric concentrations of O₃ and OH. For the surface reaction with O₃, the perturbation is significant and increases with the gas-phase concentration of O₃. Conversely, perturbations are smaller for reaction with the OH radical, which reacts with pyrene on both the surface of particles and in the gas phase. Global and regional chemical transport models typically approximate gas–particle partitioning with instantaneous-equilibration approaches. We highlight scenarios in which these approximations deviate from the explicitly coupled treatment of gas–particle partitioning and chemistry presented in this study. We find that the discrepancy between solutions depends on the operator-splitting time step and the choice of time step can help to minimize the discrepancy. The findings and techniques presented in this work not only are relevant for PAHs but can also be applied to other semi-volatile substances that undergo chemical reactions and mass transport between the gas and particle phase.

1 Introduction

Polycyclic aromatic hydrocarbons (PAHs) are air pollutants that are structurally characterized by their fused aromatic ring systems (Keyte et al., 2013). Given their carcinogenic properties (Boström et al., 2002), developmental toxicity (Billiard et al., 2008), and abundance in the environment (Ravindra et al., 2008), PAHs pose a risk to human health (Kim et al., 2013).

PAHs are semi-volatile compounds that may exist in the gas phase, adsorbed on the surface of aerosol particles, or absorbed into the bulk of aerosol particles. We describe the suspension of nano- to micrometer-sized particles in outside air as atmospheric aerosols. Typical atmospheric aerosol particles include sea salt, mineral dust, sulfate, and organic particles (Pöschl, 2005). Two key types of organic particles include soot, formed during fossil-fuel combustion, and secondary organic aerosols (SOA), formed from condensation of organic vapors in the atmosphere. The mass transfer and distribution of PAHs between the gas phase and particle phase are referred to as gas–particle partitioning. For PAHs, an accurate model description of gas–particle partitioning is needed to interpret monitoring data, determine atmospheric burden and lifetime, and ultimately assess the hazards their emissions pose to human health. Moreover, after inhalation, the distribution of a semi-volatile compound between the gas phase and the particle phase can determine its bioaccessibility (Liu et al., 2017; Wei et al., 2020). While gas-phase PAHs can directly partition into the epithelial lining fluid of the lung, particle-phase PAHs first have to dissolve from a matrix and may hence be less bioaccessible (Lammel et al., 2020).

In equilibrium, the flux of a semi-volatile species from the gas phase to the particle surface is equal to the flux that desorbs back into the gas phase. The state of the system at equilibrium can be described mathematically with thermodynamic entities (Junge, 1977; Yamasaki et al., 1982; Pankow, 1994). The equilibrium gas–particle partitioning of semi-volatile compounds is determined by both adsorption onto the particle surfaces and absorption into the particle bulk. The relative contributions of these processes depend on the concentrations, composition, and phase state of particles. Physicochemical properties of a compound – such as the octanol–air partition coefficient K_OA (Finizio et al., 1997; Harner and Bidleman, 1998), the soot–air partition coefficient K_SA (Dachs and Eisenreich, 2000; Lohmann and Lammel, 2004), and Abraham descriptors (Arp et al., 2008; Shahpoury et al., 2016) – are typically used to predict the position of equilibrium. In terms of surface adsorption, soot or black carbon particles may be especially relevant for the gas–particle partitioning of PAHs as they exhibit large energies of desorption (Kubicki, 2006; Guilloteau et al., 2010) and are often co-emitted (Miguel et al., 1998).

Field observations of PAH gas–particle partitioning in the form of particulate fractions often differ from the predictions of equilibrium models (Lohmann et al., 2000; Mandalakis et al., 2002; Terzi and Samara, 2004; Callén et al., 2008; Akyüz and Çabuk, 2010; Lammel et al., 2010; Wei et al., 2015). These discrepancies can be induced by perturbations to the equilibrium – such as emissions, changes in temperature, dry and wet deposition, or chemical reaction (Pankow and Bidleman, 1992) – but only persist if the rate of the adsorption–desorption kinetics is slow compared to rate of the perturbing process. In other words, partitioning equilibrium will still be maintained if the equilibration timescale is shorter than the perturbation timescale. This means that accurate knowledge of equilibration timescales is vital to assess whether non-equilibrium effects may occur.

Gas–particle partitioning is an important process in chemical transport models (CTMs), which describe the long-distance transport and chemical degradation of atmospheric constituents such as PAHs (Shrivastava et al., 2017; Mu et al., 2018). In CTMs, PAHs are partitioned, transformed, and transported in discrete time steps, often using the method of operator splitting. With operator splitting, the partitioning equilibrium is restored at each model time step through instantaneous equilibration (IE) (Galarneau et al., 2014) rather than treating gas–particle partitioning continuously. PAH concentrations predicted by CTMs have been shown to depend on the employed treatment of gas–particle partitioning. For instance, Lammel et al. (2009) found that using different equilibrium partitioning models influenced atmospheric cycling, the total environmental fate, and long-range transport potential of PAHs. Friedman et al. (2014) found that implementing a partitioning scheme in which PAHs slowly evaporate from aerosol particles yielded the better agreement between observed and simulated concentration and partitioning data compared to the instantaneous-equilibration approach. Overall, CTMs that assume equilibrium partitioning tend to be more common than those accounting for mass-transfer limitations explicitly, as can be seen from a recent review of partitioning methods in regional-scale transport models of SOA (McFiggans et al., 2015).

Equilibration timescales of gas–partitioning may be estimated theoretically, using analytical equations or numerical models. By solving analytical transport equations, the equilibration timescales of partitioning for volatile inorganic compounds were found to depend on the size of aerosol particles (Meng and Seinfeld, 1996). More recently, there have been numerical simulations for SOA as a function of temperature and relative humidity (Shiraiwa and Seinfeld, 2012; Li and Shiraiwa, 2019). Alternatively, equilibration timescales may also be obtained experimentally. For example, Saleh et al. (2013) found the equilibration timescale of SOA formed by α-pinene ozonolysis to be less than 30 min following a perturbation in temperature. Furthermore, the interplay between partitioning and multiphase reaction of OH with alkanes was shown to influence the distribution of product isomers (Zhang et al., 2015).

For PAHs, several studies have investigated the timescales of gas–particle partitioning from the perspective of absorptive partitioning. Rounds and Pankow used a radial diffusion model to investigate the kinetic limitations of partitioning resulting from diffusion of a semi-volatile compound absorbed within a particle (Rounds and Pankow, 1990). Odum et al. (1994) additionally included a parameter to account for mass-transfer limitations at the surface. In chamber experiments, Kamens et al. (1995) examined the equilibration timescales of PAHs. However, an in-depth analysis of the important case of PAH adsorption onto the surface of soot remains elusive. In recent years, the desorption rate coefficients of PAHs from soot have been experimentally parameterized over a range of atmospherically relevant temperatures (Guilloteau et al., 2010). However, a systematic comparison between the equilibration timescales of partitioning and the timescales of loss processes has not been carried out.

In this study, we use a kinetic model to (1) examine the timescales of gas–particle partitioning for six PAHs and (2) investigate the chemical loss of PAHs by explicitly coupling the partitioning and oxidation chemistry of the PAH pyrene. The model uses the conventions of the Pöschl–Rudich–Ammann framework (Pöschl et al., 2007) and is based on the kinetic double-layer model for aerosol surface chemistry (K2-SURF; Shiraiwa et al., 2009). We quantify the equilibration timescales of six PAHs on the model surface of solid soot particles for different temperatures and particle number concentrations (Sect. 3.2). We illustrate how the combination of slow partitioning and chemical loss of PAHs can perturb the particulate fraction from equilibrium (Sect. 3.3.2) and alter chemical lifetime (Sect. 3.3.1) in the example of the PAH pyrene. We detail how a dominant loss of pyrene from the particle phase may decrease the particulate fraction. Likewise, in the case of dominant loss of pyrene from the gas phase, the particulate fraction would increase. Compared to instantaneous partitioning, which would conserve equilibrium particulate fractions, chemical lifetime may be affected through depletion of pyrene in the more reactive phase. We apply the knowledge gained from the kinetic model calculations to the description of gas–particle partitioning in CTMs by comparing the explicitly coupled (EC) solution to a method mimicking operator splitting with instantaneous equilibration and evaluate the performance of both methods in different scenarios.

2 Methods

2.1 Kinetic model

A modified version of the kinetic double-layer model K2-SURF is used for all simulations (Shiraiwa et al., 2009). The original K2-SURF model consists of a near-surface gas phase and surface layer, with gas diffusion from the far-surface gas phase represented by a correction factor (Fig. A1). In this study, we added explicit treatment of gas diffusion to K2-SURF to track gas-phase PAH concentrations. PAHs reversibly desorb and adsorb between the aerosol particle surface and the near-surface gas phase. The rate coefficient for PAH desorption from the particle surface $k_{\mathrm{des}}=A{e}^{-E_{\mathrm{A}}/RT}$ (s⁻¹) depends on temperature T and two parameters determined from experiment: the Arrhenius factor A and the activation energy of desorption E_A (Guilloteau et al., 2010). R is the gas constant. Aerosol particles are assumed to be monodisperse and to consist of a spherical, impenetrable solid carbon core. The system is closed with respect to aerosol particles and PAH in all simulations. In simulations involving chemical reactions, the system is open with respect to oxidants; i.e., gas-phase OH and O₃ concentrations are fixed.

The differential equations in Eqs. (1)–(3) describe the time evolution of [PAH]_g (cm⁻³), [PAH]_gs (cm⁻³), and [PAH]_s (cm⁻²), which are the number concentrations (i.e., a unitless count of molecules per unit volume or unit area) of PAH in the gas phase, in the near-surface gas phase, and on the surface of aerosol particles, respectively. J_des (${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), J_ads (${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), and J_diff (s⁻¹) are the desorption flux, adsorption flux, and gas diffusion flux. Each flux term is described in detail in Appendix A.

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\mathrm{g}}}{\mathrm{d}t}}=-J_{\mathrm{diff}}{N}_{\mathrm{p}}-L_{\mathrm{g}}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\text{gs}}}{\mathrm{d}t}}={\displaystyle \frac{(J_{\mathrm{des}}-J_{\mathrm{ads}})d_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }+J_{\mathrm{diff}}}{V_{\mathrm{gs}}}}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\text{s}}}{\mathrm{d}t}}=-J_{\mathrm{des}}+J_{\mathrm{ads}}-L_{\mathrm{s}}\end{array}$$

The surface area of a single aerosol particle with diameter d_p (cm) is $d_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }$, V_gs (cm³) is the volume of gas in the near-surface gas phase for a single aerosol particle, and N_p (cm⁻³) is the particle number concentration. L_g (${\mathrm{cm}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) is the rate of chemical loss in the gas phase, and L_s (${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the rate of chemical loss in the particle phase. Reactions of PAHs within the near-surface gas phase are considered to be negligible due to the small fraction of PAHs in this volume. Sources of PAHs are not considered in this study.

2.2 Chemical reactions

The surface reaction between pyrene and O₃ is modeled using a Langmuir–Hinshelwood mechanism, including reversible adsorption of O₃ onto the surface of aerosol particles and reaction of surface-adsorbed O₃ with surface-adsorbed pyrene.

$$\begin{array}{}\text{(R1)}& {\displaystyle }& {\displaystyle }{\mathrm{O}}_{\mathrm{3}\left(\mathrm{gs}\right)}\rightleftharpoons {\mathrm{O}}_{\mathrm{3}\left(\mathrm{s}\right)}\text{(R2)}& {\displaystyle }& {\displaystyle }{\mathrm{PAH}}_{\left(\mathrm{s}\right)}+{\mathrm{O}}_{\mathrm{3}\left(\mathrm{s}\right)}\to {\mathrm{product}}_{\left(\mathrm{s}\right)}\end{array}$$

The rate coefficient for reaction of pyrene with O₃, $k_{\mathrm{PAH}+{\mathrm{O}}_{\mathrm{3}}}=\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{17}}$ cm² s⁻¹, and the corresponding mass-transport parameters are taken from Shiraiwa et al. (2009) (Table A1). Reaction products are treated as inert and non-volatile. Note that the reaction between O₃ and benzo(a)pyrene on the surface of soot has been suggested to involve the formation of long-lived reactive oxygen intermediates (Shiraiwa et al., 2011). Such a detailed chemical mechanism is beyond the scope of this study, which instead focuses on the interaction of partitioning and chemistry, and has thus been omitted for simplicity. The desorption rate coefficients of both pyrene and O₃, which are temperature dependent and explicitly included in the model, are expected to be the main driver of sensitivity in the model with regards to temperature. It should be noted that the accommodation coefficient and surface layer reaction rate coefficient may also exhibit temperature dependence but without further quantitative parameters cannot be included in the model. The gas-phase reaction between O₃ and pyrene is considered negligible and is therefore not included (Keyte et al., 2013). The reaction between pyrene and OH is accounted for both in the gas phase and on the surface of particles.

$$\begin{array}{}\text{(R3)}& {\displaystyle }{\mathrm{PAH}}_{\left(\mathrm{g}\right)}+{\mathrm{OH}}_{\left(\mathrm{g}\right)}\to {\mathrm{product}}_{\left(\mathrm{g}\right)}\text{(R4)}& {\displaystyle }{\mathrm{PAH}}_{\left(\mathrm{s}\right)}+{\mathrm{OH}}_{\left(\mathrm{gs}\right)}\to {\mathrm{product}}_{\left(\mathrm{s}\right)}\end{array}$$

The gas-phase reaction between pyrene and OH is modeled with the rate coefficient $k_{\mathrm{PAH}+\mathrm{OH}}=\mathrm{6.58}\times {\mathrm{10}}^{-\mathrm{11}}$ cm³ s⁻¹ (theoretical calculation at 298 K; Zhang et al., 2014). The reaction between pyrene and OH on the surface of particles is treated with an Eley–Rideal-like mechanism using a surface reaction probability of 0.32 (obtained for pyrene; Bertram et al., 2001) and assuming an OH gas diffusion coefficient D_g of 0.21 cm² s⁻¹ (Tang et al., 2014). The temperature dependence of the gas-phase OH reaction of PAHs has been found experimentally to be “slight to nonexistent” (Brubaker and Hites, 1998). Likewise, the reaction probability of OH on a pyrene surface has been found to exhibit only a slight temperature dependence (Liu et al., 2012). We therefore do not include temperature dependence of chemical rate coefficients in this model. The uptake of OH onto the surface of particles is considered to be irreversible.

2.3 Particulate fraction

The measured distribution of PAHs (and other semi-volatiles) between the particle phase and the gas phase is commonly described with the particulate fraction Φ, i.e., the fraction of total PAHs associated with aerosol particles (Eq. 4).

$$\begin{array}{}\text{(4)}& \mathrm{\Phi }={\displaystyle \frac{[\text{PAH}{]}_{\mathrm{p}}}{[\text{PAH}{]}_{\mathrm{p}}+[\text{PAH}{]}_{\mathrm{g}}}}\end{array}$$

The total concentration of PAH adsorbed on the surface of aerosol particles [PAH]_p (cm⁻³) is the product of the surface area of a single particle $d_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }$ with diameter d_p, the particle number concentration N_p, and surface concentration of PAH [PAH]_s (cm⁻²; Eq. 5).

$$\begin{array}{}\text{(5)}& [\text{PAH}{]}_{\mathrm{p}}=d_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }{N}_{\mathrm{p}}[\text{PAH}{]}_{\mathrm{s}}\end{array}$$

2.4 Equilibration timescale τ_eq

To quantify the time for PAHs to reach their equilibrium distribution between the gas phase and the particle phase, we use the equilibration timescale (τ_eq), defined as e-folding time for the relaxation of the system to gas–particle partitioning equilibrium. Figure 1 shows results from a kinetic box model simulation with a concentration of pyrene in air of 5×10⁵ cm⁻³, temperature T=298 K, particle number concentration $N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}\mathrm{particles}{\mathrm{cm}}^{-\mathrm{3}}$, and particle diameter d_p=200 nm. No chemical loss of pyrene is included here. τ_eq is obtained numerically from model outputs by interpolating the time required by the system to achieve $\mathrm{1}-\frac{\mathrm{1}}{e}$ (i.e., ≈63.2 %) of the difference ΔΦ between an initial particulate fraction Φ₀ and the equilibrium particulate fraction Φ_eq.

Figure 1 The evolution of the particulate fraction Φ of pyrene with respect to time, starting from initial particulate fractions Φ₀=0.1 (red) and Φ₀=0.9 (blue). The equilibration timescale τ_eq is defined as the time required for the system to achieve 63.2 % of ΔΦ, the difference between Φ₀ and Φ_eq.

In this example, pyrene reaches Φ_eq after ≈2 min, and the equilibration timescale is independent of the initial particulate fraction; i.e., τ_eq is the same regardless of whether Φ₀=0.1 or Φ₀=0.9. In fact, τ_eq is found to be independent of the choice of Φ₀ for most conditions due to the first-order and hence mono-exponential nature of the adsorption and desorption processes. This allows for consistent intercomparison across different temperatures and particle number concentrations without changing starting distributions. Exceptions may occur in cases where surface adsorption is not strictly a first-order process, due to either surface saturation effects or gas-phase diffusion limitations. These conditions occur at very low particle number concentrations (typically $<\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{particles}{\mathrm{cm}}^{-\mathrm{3}}$), and further details are given in Appendix B.

3 Results and discussion

3.1 Extreme cases of multiphase chemistry and partitioning interaction

Three extreme cases can be formulated when partitioning, and chemical-loss processes of a semi-volatile compound take place at different relative timescales (Fig. 2).

Figure 2 Schematic on how the gas–particle partitioning equilibrium of a semi-volatile compound may be perturbed from an initial state (center) due to chemical loss, depending on equilibration timescales. If the timescales of partitioning are shorter than the timescales of chemical loss, the system is able to maintain equilibrium (a). However, the combination of rapid gas-phase loss and slow replenishment from the particle phase increases the particulate fraction above the equilibrium value (b). In turn, the combination of rapid particle-phase loss and slow replenishment by condensation decreases particulate fraction below the equilibrium value (c).

When the timescale of partitioning is short compared to the timescales of chemical loss, molecules are redistributed quickly between both phases (case a, Fig. 2). In this case, the relative amounts of gas- and particle-phase species will remain very close to their equilibrium values (Φ≈Φ_eq). This is independent of whether molecules are lost primarily from the gas phase or from the particle phase.

In contrast, if the timescale of the partitioning process is slow and the chemical loss rates from the gas and particle phase differ substantially, the particulate fraction will be perturbed from its equilibrium value (cases b and c, Fig. 2). When the loss rate in the gas phase exceeds the loss rate in the particle phase, the particulate fraction increases beyond its equilibrium value (Φ>Φ_eq; case b, Fig. 2). However, when the loss rate in the particle phase is greater than that in the gas phase, the particulate fraction decreases (Φ<Φ_eq; case c, Fig. 2).

Unlike these scenarios, chemical-loss and partitioning timescales may not differ substantially. Likewise, chemical losses are likely to take place in both phases simultaneously. Every real system must therefore be seen as superposition of these cases. The extent to which perturbation occurs depends upon the difference between partitioning and chemical reaction timescales. An in-depth discussion on the magnitude of perturbation is provided in Sect. 3.3.

Hence, two preconditions are required for the particulate fraction Φ of the system to be perturbed from the equilibrium particulate fraction predicted by equilibrium partitioning theory Φ_eq: (1) slow partitioning relative to the timescale of chemical loss and (2) an imbalance of chemical loss between the gas and particle phases.

If timescales of chemical loss and partitioning were known for all natural systems, they could be classified and mathematically treated in the respective limiting case. In this manuscript, we (1) estimate the partitioning timescales of PAH on soot as a function of atmospheric conditions, (2) compare these timescales to typical chemical loss rates in order to investigate whether perturbations from equilibrium exist, and (3) explore the implications of treating partitioning and chemistry separately in chemical transport models.

3.2 Partitioning equilibration timescales for PAHs on soot

τ_eq depends on the molecular structure of the PAH, particle number concentration, and temperature. This is explored in the following section with a series of simulations using a fixed total concentration of PAHs in air of 5×10⁵ cm⁻³ and particles with a diameter of 50 nm.

Figure 3 (a) Equilibration timescale τ_eq for different PAHs as a function of particle number concentration at a constant temperature (T=298 K). (b) Comparison of the equilibration timescale of pyrene to the timescales of adsorption τ_ads and desorption τ_des, highlighting the transition between adsorption-controlled and desorption-controlled behavior.

3.2.1 Particle number concentration

The effect of varying the particle number concentration N_p on the equilibration timescale shows a distinct behavior (Fig. 3a): τ_eq is particle number independent at lower N_p, while τ_eq is particle number dependent at higher N_p. The equilibration timescales of the less strongly adsorbed PAHs – including anthracene, fluoranthene, and pyrene – are not significantly affected by particle number concentration until a fairly high threshold particle number concentration is achieved (≈10⁵, 10⁴, and 10⁴ particles cm⁻³, respectively).

Once the threshold particle number concentration is reached, a linear relationship in the double logarithmic dependence of equilibration timescale and particle number concentration emerges. The more strongly adsorbed PAHs – chrysene, benzo(e)pyrene, and benzo(a)pyrene – reach this limit at a much lower N_p (≈10² particles cm⁻³). This can be understood when looking at Fig. 3b, in which the equilibration timescale of pyrene is shown together with the individual timescales of desorption τ_des (dashed gray line, calculated with Eq. 6) and adsorption τ_ads (dotted gray line, calculated with Eq. 7).

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\mathrm{des}}={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{des}}}}=A{e}^{-E_{\mathrm{A}}/RT}\end{array}$$

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_{\mathrm{ads}}={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{ads}}}}={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-{\mathit{\theta }}_{\mathrm{s}}){\mathit{\alpha }}_{\mathrm{s},\mathrm{0}}{d}_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }{N}_{\mathrm{p}}\mathit{\omega }/\mathrm{4}}}\approx {\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{s},\mathrm{0}}{d}_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }{N}_{\mathrm{p}}\mathit{\omega }/\mathrm{4}}}\end{array}$$

In the limit of an adsorbate-free surface, adsorption and desorption are first-order processes with respect to the near-surface gas and surface concentrations of PAH, respectively, and can therefore be described with rate coefficients k_ads (s⁻¹) and k_des (s⁻¹). The desorption timescale τ_des depends on the Arrhenius factor A and the activation energy for desorption from the aerosol particle surface (E_A), the gas constant R, and temperature T. The adsorption timescale τ_ads depends on the surface accommodation coefficients on an adsorbate-free substrate α_s,0; the particle number concentration N_p; the surface area of a single aerosol particle $d_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }$ with diameter d_p; and the mean thermal velocity ω. The surface coverage θ_s is very small for typical particle number concentrations and will therefore be neglected in the following. In general, τ_eq can be approximated as a function of both timescales according to Eq. (8) (see Appendix B for details of the terms and derivation). This approximation holds as long as gas diffusion is sufficiently fast and does not limit equilibration.

$$\begin{array}{}\text{(8)}& {\mathit{\tau }}_{\mathrm{eq}}\approx {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{des}}+k_{\mathrm{ads}}}}={\displaystyle \frac{\mathrm{1}}{A{e}^{-E_{\mathrm{A}}/RT}+{\mathit{\alpha }}_{\mathrm{s},\mathrm{0}}{d}_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }{N}_{\mathrm{p}}\mathit{\omega }/\mathrm{4}}}\end{array}$$

If one process (desorption or adsorption) dominates the behavior of τ_eq, the system can be said to fall into an adsorption-controlled regime (highlighted for pyrene with blue shading) or a desorption-controlled regime (highlighted with red shading). A multi-step process in which mass is lost and transferred in one direction can be described analogously to a series of resistors in an electrical circuit, and the term limiting can be used to describe the slowest step. In contrast, the gas–particle partitioning system is a reversible system in which mass is transferred in both directions, and the relative rates of these mass-transfer processes determine the position of equilibrium. We therefore observe in Fig. 3b (and also Fig. 4b) that the equilibration time is determined primarily by the fastest process (i.e., that with the shortest timescale). We therefore adopt the term controlled to characterize this behavior.

In the low particle number concentration limit, the system is in a desorption-controlled regime and the equilibration timescale is thus strongly influenced by strength of the PAH–soot interaction, which explains the large differences in equilibration timescale between PAHs in Fig. 3a. In the high particle number concentration limit, the equilibration timescale is determined primarily by the adsorption of PAH onto particles from the near-surface gas phase and is therefore independent of PAH type as can be seen from the convergence of curves in Fig. 3a. The equilibration timescale here coincides with the adsorption timescale τ_ads, and the system is in an adsorption-controlled regime. The transition between both regimes occurs where τ_ads intersects τ_des and coincides with the point Φ_eq=0.5. At this specific point, equal amounts of PAH are in the gas and particle phases, and the timescales of desorption and adsorption contribute equally to the equilibration time.

As surface coverages θ_s are very small and PAHs generally have surface accommodation coefficients on an adsorbate-free substrate of ${\mathit{\alpha }}_{\mathrm{s},\mathrm{0}}=\mathrm{1}$ (Julin et al., 2014), we find in this study a special case of the adsorption-controlled regime where molecular collision of gas molecules is the sole controller of partitioning. For adsorbates with lower α_s,0, the adsorption timescale would be longer and the system may be in the desorption-controlled regime.

3.2.2 Temperature

The effect of varying temperature T on the equilibration timescale shows a behavior similar to the one seen for the particle number concentration (Fig. 4a): τ_eq is temperature independent at low T, while τ_eq is temperature dependent and begins to decrease at higher T. For the most weakly adsorbed PAH anthracene, τ_eq begins decreasing at 240 K towards higher T and at 298 K is already less than 5 s. The equilibration timescales for fluoranthene and pyrene begin decreasing at ≈260 K and at 298 K are both less than 100 s. Strongly adsorbed PAHs – including chrysene, benzo(e)pyrene, and benzo(a)pyrene – do not undergo a significant change in equilibration timescale in the investigated temperature range.

Figure 4 (a) Equilibration timescale τ_eq for different PAHs as a function of temperature at a constant particle number concentration ($N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{particles}{\mathrm{cm}}^{-\mathrm{3}}$). (b) Comparison of the equilibration timescale of pyrene to the timescales of adsorption τ_ads and desorption τ_des, highlighting the transition between adsorption-controlled and desorption-controlled behavior.

Again, the adsorption-controlled and desorption-controlled regimes explain this behavior (Fig. 4b). Between 210 and 240 K, PAH molecules possess little kinetic energy and are prevented from escaping into the gas phase, thus exhibiting long desorption lifetimes (Fig. A3) and high equilibrium particulate fractions (Fig. A2b). As most PAH is adsorbed on the surface of aerosol particles, molecular collision determines equilibration time. The system is in the adsorption-controlled regime, highlighted for pyrene with blue shading and signified by the coincidence with the adsorption timescale τ_ads (dotted gray line). The number of collisions between gas-phase PAHs and particles slightly increases as the thermal velocity increases, but this effect is much smaller compared to the effect of temperature increase on desorption rates. Note that the surface accommodation coefficient is assumed to be temperature independent in this study. Overall, upon increase in temperature, the desorption process becomes increasingly important. At high temperature, the system is in the desorption-controlled regime, highlighted for pyrene with red shading in Fig. 4b and signified by the coincidence with the timescale of desorption τ_des (dashed gray line).

3.3 Interplay of multiphase chemistry and partitioning

Chemical reactions with O₃ and OH are important loss processes for PAHs. If the rate of chemical loss is fast relative to gas–particle partitioning, the gas–particle distribution may be perturbed from its equilibrium state (cf. Fig. 2, cases b and c). This effect is exemplified for pyrene by including surface chemistry with O₃ (0, 1, 10, and 100 ppb) or gas-phase and surface chemistry with OH (0, 0.01, 0.1, and 1 ppt) in the model. A concentration of 10 ppb is representative of surface background O₃ concentrations (Vingarzan, 2004), while 100 ppb O₃ is characteristic of concentrations at more polluted sites (Wang et al., 2017). An OH concentration of 0.01 ppt is representative of concentrations measured at night (Geyer et al., 2003), while 0.1 ppt is representative of daytime concentrations (Stone et al., 2012), and 1 ppt is an upper limit only encountered in highly polluted conditions (Hofzumahaus et al., 2009) and smoke plumes (Hobbs et al., 2003). We employ the following conditions in the pyrene–soot system: T=280 K, $N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$ particles cm⁻³, and d_p=50 nm. At the start of the simulation, the initial total concentration of pyrene (5×10⁵ cm⁻³) is distributed between the gas and particle phases according to the particulate fraction expected at equilibrium (i.e., ${\mathrm{\Phi }}_{\mathrm{0}}={\mathrm{\Phi }}_{\mathrm{eq}}=\mathrm{0.24}$). The kinetic rate coefficients for all reactions are given in Sect. 2.2. Oxidant concentrations are fixed during simulations.

3.3.1 Non-equilibrium effects on chemical lifetime and particulate fraction

Figure 5a shows the decrease in the total (gas + particle) concentration of pyrene at different O₃ concentrations. The chemical lifetime of pyrene is determined by calculating the time needed for the concentration to reach $\mathrm{1}/e$ (i.e., ≈36.8 %) of its initial value. As the concentration of O₃ increases from 1 to 10 and 100 ppb, the chemical lifetime of pyrene subsequently decreases to 23.9, 3.1, and 1.2 h, respectively. It is informative to compare these chemical lifetimes to those calculated by assuming partitioning is infinitely fast (i.e., the particulate fraction is fixed to Φ_eq during simulations). In this fictional scenario, the lifetimes of pyrene are significantly shorter at 23.1, 2.4, and 0.3 h and correspond to decreases of 3, 23, and 75 %, respectively. This comparison demonstrates that partitioning and chemical loss are closely interlinked. In this specific case, slow partitioning prolongs the lifetime of pyrene.

Figure 5 The total pyrene concentration (gas phase+particle phase) with respect to time at different concentrations of (a) O₃ and (b) OH. O₃ reacts only with pyrene on the surface of particles, whereas OH reacts with pyrene in simultaneous gas and surface reactions. The total concentration of pyrene determined in a fictional scenario with infinitely fast partitioning is shown with gray lines. For the same simulations, the particulate fraction of pyrene with respect to time is shown: (c) O₃ and (d) OH. The equilibrium particulate fraction Φ_eq is shown with a black line.

This effect can be understood by observing the change in particulate fraction over time during each of the simulations (Fig. 5b). As each simulation proceeds, the particulate fraction Φ drops below the equilibrium particulate fraction Φ_eq and eventually reaches a quasi-steady state Φ_qs. At O₃ concentrations of 10 and 100 ppb, the particulate fractions reach values of Φ_qs=0.18 and 0.05, respectively. This effect can be explained by slow partitioning: chemical loss reduces the surface concentration of pyrene faster than its replenishment from the gas phase (non-equilibrium case c in Fig. 2). In the quasi-steady state, chemical loss and repartitioning are balanced. Importantly, both values differ significantly from Φ_eq. In contrast, when the O₃ concentration is low enough (1 ppb) the particulate fraction remains approximately equal to its value at equilibrium ($\mathrm{\Phi }\approx {\mathrm{\Phi }}_{\mathrm{eq}}=\mathrm{0.24}$). At this O₃ concentration, the rate of partitioning is sufficiently high so that pyrene lost from the particle surface can be fully replenished from the gas phase (equilibrium case a in Fig. 2). Hence, non-equilibrium behavior increases with oxidant concentration.

Figure 5c shows the decrease in total concentration of pyrene due to the simultaneous gas and surface reactions with OH. The lifetimes of pyrene with OH concentrations of 0.01, 0.1, and 1 ppt are 18.9, 1.9, and 0.2 h, respectively. Nearly identical lifetimes are obtained if partitioning is assumed to be infinitely fast, thus indicating that non-equilibrium effects on chemical lifetime are insignificant for this system. Figure 5d shows that, in contrast to the behavior of the O₃ system, the highest concentration of OH perturbs the particulate fraction to a quasi-steady state above its equilibrium value. The particulate fraction reaches a quasi-steady state with a value of Φ_qs=0.37 at 1 ppt OH. Although chemical loss takes in both phases simultaneously, the turnover of pyrene is higher in the gas phase. The particulate fraction thus increases, characteristic of the non-equilibrium case b in Fig. 2. At lower concentrations of OH, the extent of the perturbation becomes only slight (Φ_qs=0.25 at 0.1 ppt) and eventually disappears (0.01 ppt in Fig. 5d). Hence, non-equilibrium effects on particulate fraction can be significant, even if they were insignificant for chemical lifetime. This is due to the short reaction timescale of the OH–pyrene system compared to its partitioning timescale: pyrene reaches $\mathrm{1}/e$ of its initial concentration before the quasi-steady state is established.

3.3.2 Visualization of non-equilibrium effects with phase portraits

The dynamic behavior of the system may be visualized as a trajectory in the phase space of gas-phase and particle-phase pyrene concentrations, [PAH]_g and [PAH]_p (Fig. 6). At every point in the phase portrait, a vector illustrates how the system would change with time. Here, the direction of each vector arrow indicates the extent to which pyrene is being lost or transferred between phases, and its length indicates the rate of change. The exact characteristics of the phase portrait depend on temperature, available particle surface area, the strength of the PAH–soot interaction, and the rate of the chemical reactions involved. For a system where pyrene partitions without chemical loss, all trajectories converge onto a central line at which the system stops changing, known in mathematics as the nullcline (Fig. 6a). This line represents the point of gas–particle partitioning equilibrium. The slope of the line represents the dimensionless gas–particle partitioning coefficient K_p, from which the equilibrium particulate fraction Φ_eq can be obtained (Eq. 9).

$$\begin{array}{}\text{(9)}& [\text{PAH}{]}_{\mathrm{p}}=K_{\mathrm{p}}[\text{PAH}{]}_{\mathrm{g}}={\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{eq}}}{\mathrm{1}-{\mathrm{\Phi }}_{\mathrm{eq}}}}[\text{PAH}{]}_{\mathrm{g}}\end{array}$$

Figure 6 The interplay of partitioning and chemistry illustrated as simulated trajectories in the phase space of gas- and particle-phase concentrations, [PAH]_g and [PAH]_p, of pyrene. (a) Partitioning only, (b) partitioning and chemical loss from the particle phase due to surface reaction with O₃ (100 ppb, red solid line), (c) partitioning and chemical loss from the gas and particle phase due to reaction with OH (1 ppt, blue solid line), and (d) comparison to schematics in Fig. 2. The direction of each small gray vector arrow indicates the extent to which pyrene is being lost or transferred between phases, and its length is proportional to the rate of change. The nullcline representing equilibrium particulate fraction in panel (a) (Φ_eq=0.24, dashed black line) is shown in panels (b, c), alongside the respective slow manifold representing the quasi-steady-state particulate fraction (Φ_qs, solid black line).

As seen previously, chemical reactions may cause perturbation of the partitioning equilibrium. Such a perturbation would be indicated by deviation of trajectories from the nullcline in the phase portrait. The difference between a perturbed and equilibrium system is depicted for [O₃]=100 ppb in Fig. 6b. The vector field fundamentally changes, and the trajectory of an exemplary simulated system (red solid line) does not converge to the nullcline obtained in Fig. 6a, despite starting at equilibrium conditions Φ_eq (shown as a dashed black line for reference). Instead, the trajectory converges onto a central trajectory termed the slow manifold (Fraser, 1988). All trajectories in this system (represented with dotted gray lines) converge towards this manifold, irrespective of initial conditions. After approaching the slow manifold, the trajectory proceeds towards the origin (i.e., full depletion of pyrene) with a constant slope. This constant slope indicates that a constant quasi-steady-state particulate fraction Φ_qs=0.05 has been reached. The deviation of the nullcline (Fig. 6a) and the slow manifold (Fig. 6b) can be used to indicate the extent of non-equilibrium effects in a multiphase chemical reaction system. For example, Fig. 6c shows that for the reaction with 1 ppt OH the discrepancy between the simulation trajectory and the partitioning nullcline is much smaller due to simultaneous loss in the gas and particle phases. The slow manifold here runs above the partitioning nullcline and is reached only just before all pyrene is consumed (compare solid blue lines in Fig. 5). Figure 6d shows how the nullcline and the slow manifolds above or below it can be interpreted using the diagrams in Fig. 2.

3.4 Implications for chemical transport models (CTMs)

In the previous sections, an explicitly coupled model of partitioning and chemistry is used. This means that mass-transport and chemical-loss processes are simultaneously evaluated in a set of differential equations. Hereinafter, this is referred to as the explicitly coupled approach (EC). As the explicitly coupled approach is computationally expensive, CTMs often treat the partitioning and chemical loss of PAHs separately using operator splitting (Brasseur and Jacob, 2017). Instantaneous equilibration (IE) is one type of operator-splitting approach: at each model time step (Δt), the gas–particle distribution of PAH is reset to the partitioning equilibrium (estimated by temperature, particle number concentration, PAH, and particle type), and chemical loss is then further integrated separately, starting from the newly established equilibrium. Time steps of global and regional transport models used to study PAH are typically around 15 min (Galarneau et al., 2014) or 30 min (Sehili and Lammel, 2007).

3.4.1 Influence of model time step length

The solution obtained using the IE approach can differ from the EC solution. Using the surface reaction of O₃ (100 ppb) with pyrene in the particle phase, we demonstrate that the magnitude and sign of this difference vary between Δt=4 min, 8 min, 30 min, and 1 h (Fig. 7a). The following conditions are used in the simulations: T=280 K, d_p=50 nm, and $N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{particles}{\mathrm{cm}}^{-\mathrm{3}}$.

Figure 7 Time evolution of total pyrene concentration in model calculations, comparing the explicitly coupled (EC) solution to the operator splitting with instantaneous-equilibration (IE) solution (Δt=4, 8, 30 min and 1 h) for systems with (a) surface reaction of particle-phase pyrene with 100 ppb O₃ and (b) reaction of gas- and surface-bound pyrene with 1 ppt OH. The operator-splitting time step leading to the smallest deviation is depicted as Δt_opt.

The lifetime of pyrene is underestimated when using time steps Δt of 4 and 8 min (Fig. 7a) but overestimated with Δt of 30 min and 1 h. In this example, an optimal time step exists for which the deviation from EC is minimized, Δt_opt=19.9 min (Fig. 7a). It is close to the equilibration timescale of gas–particle partitioning τ_eq of pyrene, which is around 15 min under these conditions. τ_eq could thus serve as an initial guess for Δt_opt. Δt_opt is determined using a golden-section search optimization algorithm (Kiefer, 1953) to minimize the absolute difference between EC and IE model outputs. Of note, a deviation between IE and EC and hence a dependence of the model result on the operator-splitting time step only arises if a significant departure from partitioning equilibrium occurs. Under equilibrium partitioning conditions, a range of sufficiently short Δt can describe the system accurately. In an example with OH (1 ppt) reacting with gas-phase and surface-bound pyrene, all IE calculations produce negligible errors, irrespective of the time step used (Fig. 7b). This is due to the particulate fraction being very close to Φ_eq until the majority of pyrene has reacted (cf. Fig. 5d).

3.4.2 Deviation from explicitly coupled (EC) approach

The discrepancy between the EC and IE solutions depends not only on the length of Δt but also on the relative rates of partitioning and chemical loss. In this section, this discrepancy is explored as a function of desorption rate and is therefore characteristic of a range of PAHs. The reaction rate coefficients of pyrene are used as a best guess for generic PAHs. The discrepancy can be quantified with an error metric, E_loss, which can be interpreted as the relative difference in loss rates (Eq. 10). Δ[PAH]_EC(t) and Δ[PAH]_IE(t) are the accumulated losses of PAH at each time point t using EC and IE, respectively (Eqs. 11 and 12). This metric is chosen as it detects discrepancies in model solutions independent of the absolute turnover, which is important when comparing scenarios at high and low oxidant concentrations. E_loss ranges between −1 and 1 and is evaluated until either t_90 % or t=24 h is reached. E_loss is positive when the IE solution overpredicts the loss of PAHs compared to the reference EC solution and negative when loss is underestimated.

$$\begin{array}{}\text{(10)}& {\displaystyle }{E}_{\mathrm{loss}}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{t=\mathrm{1}}^n{\displaystyle \frac{\mathrm{\Delta }\left[\text{PAH}{]}_{\mathrm{IE}}\right(t)-\mathrm{\Delta }[\text{PAH}{]}_{\mathrm{EC}}\left(t\right)}{\mathrm{\Delta }\left[\text{PAH}{]}_{\mathrm{EC}}\right(t)+\mathrm{\Delta }[\text{PAH}{]}_{\mathrm{IE}}\left(t\right)}}\text{(11)}& {\displaystyle }\mathrm{\Delta }\left[\text{PAH}{]}_{\mathrm{EC}}\right(t)=[\text{PAH}{]}_{\mathrm{EC}}\left(\text{0}\right)-\left[\text{PAH}{]}_{\mathrm{EC}}\right(t)\text{(12)}& {\displaystyle }\mathrm{\Delta }\left[\text{PAH}{]}_{\mathrm{IE}}\right(t)=[\text{PAH}{]}_{\mathrm{IE}}\left(\text{0}\right)-\left[\text{PAH}{]}_{\mathrm{IE}}\right(t)\end{array}$$

Figure 8 Discrepancy between the explicitly coupled (EC) solution and the operator-splitting (IE) solution for surface reaction of particle-phase PAHs with O₃, determined with the error metric E_loss. Regions where IE overestimates the loss of PAH are indicated in red ($\mathrm{0}<E_{\mathrm{loss}}\le \mathrm{1}$), and regions where it underestimates PAH loss in blue ($-\mathrm{1}\ge E_{\mathrm{loss}}<\mathrm{0}$). The result is shown using two different operator-splitting time step lengths: (a) Δt=1 min and (b) Δt=30 min. The dark gray area indicates the conditions under which the time taken for total PAH to deplete by 90 % (t_90 %) is less than a single operator-splitting time step. The desorption rates for the PAHs chrysene, pyrene, and anthracene at 280 K are indicated on the top axis.

Figure 8 shows the extent and direction of deviation of IE from EC in a case study of PAH surface chemistry in which the desorption rate coefficient k_des is varied between $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{8}}$ and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{1}}$ s⁻¹ and the concentration of O₃ between 0 and 120 ppb for IE time steps of Δt=1 min and Δt=30 min. When Δt=1 min (Fig. 8a), IE overestimates PAH loss compared to EC, indicated by red coloring. Deviation is largest when k_des is between $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹. Here, the IE time step of 1 min causes PAH to transfer onto particles at an artificially high rate. This increases the particle-phase concentration of PAH and results in faster chemical loss. The IE solution hence shows fewer non-equilibrium effects of slow partitioning on multiphase chemistry compared to the reference EC model. At the highest k_des ($>\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹), non-equilibrium effects of slow partitioning still occur, but in the desorption-controlled regime (cf. Fig. 4 at 280 K) an increase in k_des leads to a reduction in equilibration timescale. This leads not only to weaker non-equilibrium effects of slow partitioning in the EC solution but also to a better match between EC equilibration timescale and IE time step. Hence, the discrepancy between the IE and EC approach is reduced, as evident by the more faint red coloring. At low k_des ($<\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹), most PAH is located on the surface of particles at all times, and re-partitioning of gas-phase PAH after depletion of particle-phase PAH is negligible. Thus, no deviation of the IE from the EC approach occurs.

In contrast, with a time step of Δt=30 min (Fig. 8b), the IE approach generally underestimates the loss of PAH compared to the EC approach, indicated by blue coloring. The largest underestimations are found at high k_des and high [O₃]_g. Underestimation of PAH loss occurs because the re-partitioning induced by the longer IE time step of 30 min is slower than the true equilibration rate in the EC model. In this scenario, the IE approach thus leads to stronger non-equilibrium effects of slow partitioning compared to the EC model. When the equilibration timescale becomes shorter, at high k_des ($>\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹), the discrepancy between the IE and the EC solution further increases, especially at high [O₃]_g. Notably, at $k_{\mathrm{des}}\approx \mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, the IE approach slightly overestimates loss of PAH at high [O₃]_g. This is due to the EC equilibration timescale dropping below the equivalent of Δt=30 min just before non-equilibrium effects of slow partitioning vanish at the lowest k_des. In between, a zero contour (labeled “0”), and hence no deviation between both methods, occurs when $k_{\mathrm{des}}\approx \mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ s⁻¹. Here, the IE approach matches the EC equilibration timescale by alternatingly underestimating and overestimating the concentration of PAH at different points of the simulation, and a cancellation of errors occurs (cf. Fig. 7a). This case is distinct from the left zero contour at $k_{\mathrm{des}}\approx \mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹, in which departure from equilibrium does not occur and both methods return truly identical results (compare panels a and c of Fig. A5). Another region exists: for $k_{\mathrm{des}}<\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{s}}^{-\mathrm{1}}$, the simulation proceeds for less than 30 min, and therefore fewer than one IE time step is evaluated (Fig. A5). In this region partitioning effectively does not take place, and we chose not to report the numerical value of E_loss.

E_loss is negligibly small when the concentration of [OH]_g is varied between 0 and 0.5 ppt with reaction both in the gas phase and on the surface of particles (Fig. A4a). Unlike the example with O₃, PAH loss rates due to reaction with OH in each phase are similar enough that Φ is not perturbed far from Φ_eq under these conditions. Often, in chemical transport models only the gas-phase reaction of PAH with OH is included. Varying the concentration of OH between 0 and 0.5 ppt with reaction in the gas phase only causes IE to overestimate the loss of PAH (red area, Fig. A4b).

To estimate errors for global models, it is also informative to present the discrepancy between the EC and IE approaches as a percentage difference. In Fig. 8 for instance, when $E_{\mathrm{loss}}=-\mathrm{0.5}$, a model using the IE assumption underestimates the loss of PAH by 34 % compared to the EC solution (Fig. 8a). Likewise, if $E_{\mathrm{loss}}=-\mathrm{0.1}$, a model using the IE assumption overestimates the loss of PAH by 23 %. In order to fully quantify the error introduced by the instantaneous-equilibration assumption, it would be necessary to implement a non-equilibrium partitioning scheme directly into a CTM.

It should be noted that, alongside the gas-phase and particle-phase concentrations of PAHs, CTMs are often evaluated by comparing the predicted particulate fraction and partitioning coefficient to observational data. Both of these metrics depend on the relative concentrations of gas-phase and particle-phase PAHs. Therefore, due to error propagation, the particulate fraction and partitioning coefficient may be more sensitive than absolute concentration to the effects of a non-equilibrium partitioning scheme.

4 Atmospheric implications

This study shows that the chemical loss of polycyclic aromatic hydrocarbons (PAHs) and their partitioning between the gas and particle phases are closely interlinked. The equilibration timescales of adsorptive partitioning are quantified for six PAHs on the surface of soot. Our model predicts that equilibration timescales range from seconds to hours depending on temperature, available particle surface area, and molecular structure of the PAH. We highlight the molecular processes governing this timescale with two regimes: adsorption-controlled and desorption-controlled partitioning.

Soot constitutes only a fraction of total ambient aerosol particles (Pöschl, 2005). Thus, a logical next step will be to investigate how equilibration timescales vary for other types of particle surfaces. For example, given the weaker desorption energies of PAHs such as anthracene on the surface of NaCl (75.3 kJ mol⁻¹) compared to soot (87.9 kJ mol⁻¹; Chu et al., 2010), one would expect equilibration timescales on NaCl to be shorter.

On other particle types, PAH molecules can undergo absorptive partitioning by diffusing through surface layers into the bulk of the particle. For secondary organic aerosol (SOA), the particle phase state influences the rates of condensation and evaporation (Shiraiwa and Seinfeld, 2012). Equilibration timescales for PAHs are therefore also expected to be dependent on particle composition and humidity. For absorptive partitioning, the equilibration timescales of PAHs are expected to be even longer than the equilibration timescales for adsorptive partitioning. Alongside the contributions to the equilibration timescale from the adsorption- and the desorption-controlled regimes, absorptive partitioning is also controlled by the diffusion of PAHs through the bulk of aerosol particles. The complex interplay of partitioning and reaction in the gas and particle phases plays a critical role in the growth of SOA particles (Shiraiwa et al., 2013; Berkemeier et al., 2020), and departure from partitioning equilibrium adds to this complexity (Cappa and Wilson, 2011). However, the role of bulk diffusion in determining equilibration timescales is beyond the scope of this study and will be investigated in a follow-up study that builds on the framework provided here.

Chemical reaction of pyrene with O₃ on the surface of particles perturbs the particulate fraction from partitioning equilibrium at atmospherically relevant oxidant concentrations. As the extent of this perturbation increases with concentration of O₃, the largest deviations from equilibrium particulate fraction are likely to take place in the most polluted air. The reaction of pyrene with OH in both phases results in much smaller perturbations. In general, the biggest deviations from equilibrium particulate fraction are expected for low-volatility PAHs when atmospheric conditions induce slow partitioning (i.e., cold temperatures and low particle number concentrations). Other chemical-loss processes may also be important for PAHs, such as the reaction with NO₃ in both the gas phase (Zhang et al., 2014) and the particle phase (Gross and Bertram, 2008), as well as aqueous-phase photodegradation processes (Fasnacht and Blough, 2002). These reactions must eventually be studied simultaneously in order to establish whether loss in both phases balances out or whether perturbation from equilibrium takes place.

Using existing observational datasets, it may be possible to establish how the size of the deviation in particulate fraction (between observed values and the predictions of equilibrium partitioning models) depends on the concentrations of OH, O₃, NO₃, and other perturbing variables. This could help identify and compare key perturbing variables in a real-world setting.

It should also be noted that, while in this study simulations involving chemical loss are initialized at the point of equilibrium (Φ₀=Φ_eq), in reality PAHs may be emitted in a state far from partitioning equilibrium. Depending on the prevailing loss processes, such an effect could either enhance or inhibit perturbations from equilibrium.

Non-chemical-loss processes, such as dry and wet deposition, remove PAHs from the gas and the particle phases at different rates and may also cause perturbations from partitioning equilibrium. The fastest loss processes – i.e., those operating at the shortest timescales – will cause the greatest perturbations. In the case of polybrominated diphenyl ethers (PBDEs), Li et al. (2015) attempted to include the effect of loss by deposition on partitioning and derived an equation for the partitioning coefficient assuming a steady state rather than equilibrium. However, given that for PAHs the estimated lifetimes due to dry deposition (1 to 14 d) and wet deposition (5 to 15 months; Škrdlíková et al., 2011) are much longer than lifetimes due to chemical loss (in this study less than 24 h), chemical loss is expected to be the loss process that is most likely to perturb the partitioning equilibrium.

The methodology described in this study is universally applicable to semi-volatile compounds on solid surfaces if mass-transfer parameters and chemical reaction rate coefficients are available. In some cases it may be necessary to estimate these parameters. In quantum mechanical calculations, graphene surfaces could be used as a model for soot and desorption energies. Such values are already available for PBDEs (Ding et al., 2014) and other small organic molecules (Lazar et al., 2013).

It has to be noted that an explicitly coupled solution of partitioning and chemical loss is computationally too expensive for inclusion in typical regional and global CTMs. Hence, alternative algorithms would be highly anticipated. Knowledge about the position of the partitioning steady state in the presence of chemical loss (as indicated by the slow manifold that can be visualized in a phase portrait of gas- and particle-phase concentrations) could be used to develop such a method for global and regional models.

Appendix A: Kinetic model

The flux of PAH molecules from the gas phase to the near-surface gas phase J_diff with concentrations [PAH]_g and [PAH]_gs, respectively, is calculated with Eq. (A1).

$$\begin{array}{}\text{(A1)}& J_{\mathrm{diff}}=\mathrm{2}\mathit{\pi }(d_{\mathrm{p}}+\mathrm{2}\mathit{\lambda })D_{\mathrm{g}}\left(\right[\text{PAH}{]}_{\mathrm{g}}-\left[\text{PAH}{]}_{\mathrm{gs}}\right)\end{array}$$

The gas-phase diffusion coefficient D_g is fixed at 0.06 cm² s⁻¹ for all PAH compounds, based on measurements for anthracene and pyrene in nitrogen (Siddiqi et al., 2009). d_p is the diameter of aerosol particles. The mean free path λ is defined in Eq. (A2).

$$\begin{array}{}\text{(A2)}& \mathit{\lambda }={\displaystyle \frac{\mathrm{3}{D}_{\mathrm{g}}}{\mathit{\omega }}}\end{array}$$

The mean thermal velocity of a molecule ω depends on temperature T and its molar mass M (Eq. A3).

$$\begin{array}{}\text{(A3)}& \mathit{\omega }=\sqrt{{\displaystyle \frac{\mathrm{8}RT}{\mathit{\pi }M}}}\end{array}$$

The adsorption flux J_ads of molecules from the near-surface gas phase to the particle phase is described using Eq. (A4).

$$\begin{array}{}\text{(A4)}& J_{\mathrm{ads}}={\mathit{\alpha }}_{\mathrm{s},\mathrm{0}}(\mathrm{1}-{\mathit{\theta }}_{\mathrm{s}})J_{\mathrm{coll}}\end{array}$$

The surface accommodation coefficient on an adsorbate-free substrate α_s,0 describes the probability that a molecule adsorbs upon collision with an adsorbate-free aerosol particle and for PAH molecules is assumed to be ${\mathit{\alpha }}_{\mathrm{s},\mathrm{0}}=\mathrm{1}$ (Julin et al., 2014). The total sorption layer coverage θ_s is calculated as the sum of the fractional coverages of PAH and O₃, θ_PAH and ${\mathit{\theta }}_{{\mathrm{O}}_{\mathrm{3}}}$, respectively (Eq. A5). θ_PAH and ${\mathit{\theta }}_{{\mathrm{O}}_{\mathrm{3}}}$ are calculated using the surface concentrations [PAH]_s and [O₃]_s and molecular cross sections σ_PAH and ${\mathit{\sigma }}_{{\mathrm{O}}_{\mathrm{3}}}$ of PAH and O₃, respectively (Eq. A6). In order to estimate σ, each benzene-like ring of a PAH molecule is assumed to occupy $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{15}}$ cm². The collision flux J_coll – i.e., the flux of molecules colliding with the surface – is defined in Eq. (A7).

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{s}}={\mathit{\theta }}_{\mathrm{PAH}}+{\mathit{\theta }}_{{\mathrm{O}}_{\mathrm{3}}}\text{(A6)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{PAH}}=\mathit{\sigma }[\text{PAH}{]}_{\mathrm{s}}\text{(A7)}& {\displaystyle }& {\displaystyle }{J}_{\mathrm{coll}}={\displaystyle \frac{[\text{PAH}{]}_{\mathrm{gs}}\mathit{\omega }}{\mathrm{4}}}\end{array}$$

The temperature-dependent desorption flux (J_des) due to PAH molecules evaporating from the surface of an aerosol particle depends on the rate coefficient for desorption (k_des) and the surface concentration of PAH [PAH]_s (Eq. A8).

$$\begin{array}{}\text{(A8)}& J_{\mathrm{des}}=k_{\mathrm{des}}[\text{PAH}{]}_{\mathrm{s}}\end{array}$$

k_des depends on the Arrhenius factor (A) and the activation energy for desorption from the aerosol particle surface (E_A; Eq. A9).

$$\begin{array}{}\text{(A9)}& k_{\mathrm{des}}=A{e}^{-E_{\mathrm{A}}/RT}\end{array}$$

The temperature dependence of the desorption rate coefficient was previously determined for seven PAHs on fresh kerosene soot (Guilloteau et al., 2008, 2010), and the obtained parameters are implemented in this model (see Table A1). These activation energies of desorption for PAHs on soot are consistent with those obtained theoretically on pure graphene (Lechner and Sax, 2014) and coronene (Kubicki, 2006). It should also be noted that different types of soot can have different effects on gas–particle partitioning (Mader and Pankow, 2002) and more aged soot may have a reduced affinity for PAH. Despite the simplifications of this model, we aim to provide a basis to which further complexity can be added.

Table A1 Mass-transport parameters to determine the desorption rate coefficients of PAHs and O₃ ($k_{\mathrm{des}}=A{e}^{-E_{\mathrm{A}}/RT}$): Arrhenius parameter A and activation energy of desorption E_A (Guilloteau et al., 2008, 2010).

Irreversible reactions between pyrene and either O₃ on the surface of aerosol particles or OH in the gas phase and on the surface of aerosol particles are investigated with the model. The equations of mass transport for O₃ are identical to those for PAH, and the corresponding parameters are reported in Tables A1 and A2. As the uptake of OH is considered to proceed via an Eley–Rideal mechanism, the diffusion of OH from the gas phase to the near-surface gas phase is treated using a gas-phase diffusion correction factor C_g,OH (Eq. A10). The full equation for C_g,OH can be found in Eq. (14) of Pöschl et al. (2007).

$$\begin{array}{}\text{(A10)}& [\mathrm{OH}{]}_{\mathrm{gs}}={\mathrm{C}}_{\mathrm{g},\mathrm{OH}}[\mathrm{OH}{]}_{\mathrm{g}}\end{array}$$

Table A2 Mass-transport parameters for adsorption of PAHs and O₃: molar mass M; surface accommodation coefficient on an adsorbate-free substrate α_s,0 (Julin et al., 2014); and molecular cross section σ. The gas-phase diffusion coefficient D_g used for PAHs is 0.06 cm² s⁻¹ (Siddiqi et al., 2009), and that for O₃ is 0.14 cm² s⁻¹ (Massman, 1998).

The rate of PAH loss from the particle surface due to chemical reaction with OH L_s,OH depends on probability γ_OH,PAH that reaction occurs following collision of OH with PAH on the particle surface (Eq. A11). The rate of gas-phase PAH loss by OH L_g,OH and the rate of PAH loss from the surface due to reaction with O₃ $L_{\mathrm{s},{\mathrm{O}}_{\mathrm{3}}}$ are defined by Eqs. (A12) and (A13), respectively. Further details of these reactions and their parameters can be found in Sect. R1.

$$\begin{array}{}\text{(A11)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{s},\mathrm{OH}}={\mathit{\gamma }}_{\mathrm{OH},\text{PAH}}{\mathit{\theta }}_{\mathrm{PAH}}{J}_{\text{coll},\mathrm{OH}}\text{(A12)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{g},\mathrm{OH}}=k_{\mathrm{g}}[\text{PAH}{]}_{\mathrm{g}}[\mathrm{OH}{]}_{\mathrm{g}}\text{(A13)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{s},{\mathrm{O}}_{\mathrm{3}}}=k_{\mathrm{s}}\left[\text{PAH}{]}_{\mathrm{s}}\right[{\mathrm{O}}_{\mathrm{3}}{]}_{\mathrm{s}}\end{array}$$

Figure A1 The model used in this study consists of three layers: gas phase, near-surface gas phase, and sorption layer. The surface-layer reaction (SLR) with O₃ and the gas-surface reactions (GSR) and gas-phase reactions (GR) with OH are modeled alongside mass transport and partitioning. J_diff, J_des, and J_ads represent the gas diffusion flux between the gas phase and the near-surface gas phase, the desorption flux from the particle surface to the near-surface gas phase, and the adsorption flux from the near-surface gas phase to the particle surface, respectively. Note that, for the gas diffusion flux of OH, gas diffusion was implemented through a correction factor C_g,OH.

Figure A2 Equilibrium particulate fraction Φ_eq of selected PAHs as a function of (a) particle number concentration (T=298 K) and (b) temperature ($N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$) for particles of size d_p=50 nm.

Figure A3 Desorption lifetime (${\mathit{\tau }}_{\mathrm{des}}=\mathrm{1}/k_{\mathrm{des}}$) of selected PAHs as a function of temperature.

Figure A4 Discrepancy, determined with the error metric E_loss, between the explicitly coupled (EC) solution and the operator-splitting (IE) solution with a time step of Δt=30 min for the reaction of OH (1 ppt) with (a) both gas- and particle-phase pyrene and (b) only gas-phase pyrene. Regions where EC overestimates the loss of pyrene are indicated in red ($\mathrm{0}<E_{\mathrm{loss}}\le \mathrm{1}$) and underestimated in blue ($-\mathrm{1}\ge E_{\mathrm{loss}}<\mathrm{0}$). The conditions are 280 K and 1×10³ particles cm⁻³.

Figure A5 The explicitly coupled (EC) solution and instantaneous-equilibration (IE) solutions with time steps Δt=1 and 30 min at specific points for the evaluation of the error metric E_loss (Fig. 8a and b). These points represent a horizontal cross section of Fig. 8 at $\left[{\mathrm{O}}_{\mathrm{3}}\right]=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{12}}{\mathrm{cm}}^{-\mathrm{3}}$ and the following k_des: (a) 10^−6.3, (b) 10^−4.7, (c) 10^−3.1, and (d) 10^−1.5 s⁻¹.

Figure A6 The explicitly coupled (EC) solution and instantaneous-equilibration (IE) solutions with time steps Δt=1 and 30 min at specific points for the evaluation of the error metric E_loss (Fig. A4b). These points represent a horizontal cross section of Fig. A4b at $\left[\mathrm{OH}\right]=\mathrm{1}\times {\mathrm{10}}^{\mathrm{7}}{\mathrm{cm}}^{-\mathrm{3}}$ and the following k_des: (a) 10^−6.3, (b) 10^−4.7, (c) 10^−3.1, and (d) 10^−1.5 s⁻¹.

Table A3 Symbols, definitions and units.

Appendix B: Derivation of equation for equilibration time

An approximate equation for equilibration time τ_eq (s) is obtained analytically using the relaxation time of a simple reversible reaction (Bernasconi, 1976). The equation approximates the numerically obtained results from the kinetic model and is derived by assuming that gas–particle partitioning can be described as two first-order processes, adsorption and desorption, with rate coefficients k_des (s⁻¹) and k_ads (s⁻¹). We find there is good agreement between the numerically obtained results and the approximate equation as long as the gas diffusion flux J_diff does not significantly affect gas–particle partitioning (i.e., [PAH]_g≈[PAH]_gs) and surface crowding effects do not significantly inhibit adsorption of PAH onto the surface (i.e., θ_s is small). The equilibration timescale of PAHs predicted by the equation are within 10 % of the numerically obtained timescale among our test conditions ($\ge {\mathrm{10}}^{\mathrm{3}}\mathrm{particles}{\mathrm{cm}}^{-\mathrm{3}}$, $\mathrm{210}\mathrm{K}\ge T\ge \mathrm{310}\mathrm{K}$). For a gas-phase [PAH]_g (cm⁻³) and particle-phase concentration [PAH]_p (cm⁻³), the rate equation can be expressed in the following way:

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\mathrm{g}}}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\mathrm{p}}}{\mathrm{d}t}}=k_{\mathrm{des}}[\text{PAH}{]}_{\mathrm{p}}-k_{\mathrm{ads}}[\text{PAH}{]}_{\mathrm{g}}.\end{array}$$

At equilibrium, the concentration of PAH in the gas and particle phases is [PAH]_g,eq and [PAH]_p,eq. As there is no net transfer in mass,

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\mathrm{g}}}{\mathrm{d}t}}=k_{\mathrm{des}}[\text{PAH}{]}_{\mathrm{p},\text{eq}}-k_{\mathrm{ads}}[\text{PAH}{]}_{\mathrm{g},\mathrm{eq}}=\mathrm{0}.\end{array}$$

The concentrations of PAH in the gas and particle phases can therefore be expressed as a displacement from their values at equilibrium, Δ[PAH]_g and Δ[PAH]_p (Eqs. B3 and B4).

$$\begin{array}{}\text{(B3)}& {\displaystyle }& {\displaystyle }[\text{PAH}{]}_{\mathrm{g}}=[\text{PAH}{]}_{\mathrm{g},\mathrm{eq}}+\mathrm{\Delta }[\text{PAH}{]}_{\mathrm{g}}\text{(B4)}& {\displaystyle }& {\displaystyle }[\text{PAH}{]}_{\mathrm{p}}=[\text{PAH}{]}_{\mathrm{p},\text{eq}}+\mathrm{\Delta }[\text{PAH}{]}_{\mathrm{p}}\end{array}$$

Substitution of $\mathrm{\Delta }[\text{PAH}{]}_{\mathrm{g}}=-\mathrm{\Delta }[\text{PAH}{]}_{\mathrm{p}}=x$ into Eqs. (B3) and (B4) allows the displacements from equilibrium to be expressed in terms of a single variable x (Eqs. B5 and B6).

$$\begin{array}{}\text{(B5)}& {\displaystyle }& {\displaystyle }[\text{PAH}{]}_{\mathrm{g}}=[\text{PAH}{]}_{\mathrm{g},\mathrm{eq}}+x\text{(B6)}& {\displaystyle }& {\displaystyle }[\text{PAH}{]}_{\mathrm{p}}=[\text{PAH}{]}_{\mathrm{p},\text{eq}}-x\end{array}$$

Insertion of Eqs. (B5) and (B6) into Eq. (B2) gives

$$\begin{array}{}\text{(B7)}& {\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\mathrm{g}}}{\mathrm{d}t}}=+k_{\mathrm{des}}\left(\right[\text{PAH}{]}_{\mathrm{p},\text{eq}}-x)-k_{\mathrm{ads}}([\text{PAH}{]}_{\mathrm{g},\mathrm{eq}}+x)=\mathrm{0}.\end{array}$$

Expansion of brackets, followed by insertion of Eq. (B2), gives

$$\begin{array}{}\text{(B8)}& {\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\mathrm{g}}}{\mathrm{d}t}}=-(k_{\mathrm{des}}+k_{\mathrm{ads}})x=\mathrm{0}.\end{array}$$

Insertion of Eq. (B5) into Eq. (B8) gives Eq. (B9); then noting that $\mathrm{d}[\text{PAH}{]}_{\mathrm{g},\mathrm{eq}}/\mathrm{d}t=\mathrm{0}$ gives Eq. (B10).

$$\begin{array}{ll}\text{(B9)}& {\displaystyle \frac{\mathrm{d}\left(\right[\text{PAH}{]}_{\mathrm{g},\mathrm{eq}}+x)}{\mathrm{d}t}}& {\displaystyle }={\displaystyle \frac{\mathrm{d}[\text{PAH}{]}_{\mathrm{g},\mathrm{eq}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}{\displaystyle }& {\displaystyle }=-(k_{\mathrm{des}}+k_{\mathrm{ads}})x\\ {\displaystyle }& {\displaystyle }=\mathrm{0}\\ \text{(B10)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}& {\displaystyle }=-(k_{\mathrm{des}}+k_{\mathrm{ads}})x\end{array}$$

Both sides are integrated (Eq. B11), and the equation rearranged (Eq. B12).

$$\begin{array}{}\text{(B11)}& {\displaystyle }& {\displaystyle }\underset{x_{\mathrm{0}}}{\overset{x}{\int }}{\displaystyle \frac{\mathrm{1}}{x}}\mathrm{d}t=\underset{t_{\mathrm{0}}}{\overset{t}{\int }}-(k_{\mathrm{des}}+k_{\mathrm{ads}})\mathrm{d}t\text{(B12)}& {\displaystyle }& {\displaystyle }x=x_{\mathrm{0}}{e}^{-(k_{\mathrm{des}}+k_{\mathrm{ads}})t}\end{array}$$

The equilibration time τ_eq can therefore be considered the e-folding time required for the displacement x to decrease to $\mathrm{1}/e$ of its initial value (Eq. B13).

$$\begin{array}{}\text{(B13)}& {\mathit{\tau }}_{\mathrm{eq}}\approx {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{des}}+k_{\mathrm{ads}}}}={\displaystyle \frac{\mathrm{1}}{A{e}^{-E_{\mathrm{A}}/RT}+{\mathit{\alpha }}_{\mathrm{s},\mathrm{0}}{d}_{\mathrm{p}}^{\mathrm{2}}\mathit{\pi }{N}_{\mathrm{p}}\mathit{\omega }/\mathrm{4}}}\end{array}$$