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

multiphase chemical reactions of semi-volatile compounds: mechanistic insights and practical implications for atmospheric modeling of PAHs Jake Wilson1, Ulrich Pöschl1, Manabu Shiraiwa2,*, and Thomas Berkemeier1,* 1Multiphase Chemistry Department, Max Planck Institute for Chemistry, Mainz, Germany 2Department of Chemistry, University of California, Irvine, CA, USA Correspondence: Thomas Berkemeier (t.berkemeier@mpic.de) and Manabu Shiraiwa (m.shiraiwa@uci.edu)


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). 25 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. As atmospheric aerosols, we describe the suspension of nano-to micrometer sized particles in outside air. 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 30 phase and particle phase is 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 35 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 40 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 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 (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 60 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 70 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 75 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 80 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) to 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 PRA (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 85 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 90 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 solution to a method mimicking operator splitting with instantaneous equilibration and evaluate the performance of both methods in different scenarios.

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 100 gas phase. The rate coefficient for PAH desorption from the particle surface k des = Ae −EA/RT (s −1 ) 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 3 concentrations are fixed.

105
The differential equations in Eqs. 1-3 describe the time evolution of which are the number concentrations (i.e. a unitless count of molecules per unit volume or unit area) of PAH in the gas phase, near-surface gas phase and on the surface of aerosol particles, respectively. J des (cm −2 s −1 ), J ads (cm −2 s −1 ) and J diff (s −1 ) are the desorption flux, adsorption flux and gas diffusion flux. Each flux term is described in detail in section A of the appendix.
The surface area of a single aerosol particle with diameter d p (cm) is d 2 p π, V gs (cm 3 ) is the volume of gas in the nearsurface gas phase for a single aerosol particle and N p (cm −3 ) is the particle number concentration. L g (cm −3 s −1 ) is the rate of chemical loss in the gas phase and L s (cm −2 s −1 ) is the rate of chemical loss in the particle phase. Reactions of PAHs within 115 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.

Chemical reactions
The surface reaction between pyrene and O 3 is modeled using a Langmuir-Hinshelwood mechanism, including reversible adsorption of O 3 onto the surface of aerosol particles and reaction of surface-adsorbed O 3 with surface-adsorbed pyrene.
The rate coefficient for reaction of pyrene with O 3 , k PAH+O3 = 2.7 × 10 −17 cm 2 s −1 , 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 125 the reaction between O 3 and benzo(a)pyrene on the surface of soot has been suggested to involve the formation of long-lived reactive oxygen intermediates (ROI; 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 3 , 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 130 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 3 and pyrene is considered negligible and is therefore not included (Keyte et al., 2013). The reaction between pyrene and OH is accounted for in both the gas phase and on the surface of particles.
The gas-phase reaction between pyrene and OH is modeled with the rate coefficient k PAH+OH = 6.58 × 10 −11 cm 3 s −1 (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.21cm 2 s −1 (Tang et al., 2014). The temperature dependence of the gas-phase 140 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.

Particulate fraction 145
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).
The total concentration of PAH adsorbed on the surface of aerosol particles [PAH] p (cm −3 ) is the product of the surface area of a single particle d 2 p π with diameter d p , the particle number concentration N p , and surface concentration of PAH, [PAH] s 150 (cm −2 ; Eq. 5).

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.
155 Figure 1 shows results from a kinetic box model simulation with a concentration of pyrene in air of 5×10 5 cm −3 , temperature T = 298 K, particle number concentration N p = 1×10 4 particles cm −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 1 − 1 e (i.e. ≈ 63.2 %) of the difference ∆Φ between an initial particulate fraction Φ 0 and the equilibrium particulate fraction Φ eq .

160
In this example, pyrene reaches Φ eq after ≈ 2 minutes and the equilibration timescale is independent of the initial particulate fraction, i.e. τ eq is the same regardless of whether Φ 0 = 0.1 or Φ 0 = 0.9. In fact, τ eq is found to be independent of the choice of Φ 0 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, either due to 165 surface saturation effects or gas phase diffusion limitations. These conditions occur at very low particle number concentrations (typically < 1 × 10 3 particles cm −3 ) and further details are given in Appendix B.
3 Results and discussion

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 170 at different relative timescales (Fig. 2).  and Φ0 = 0.9 (blue). The equilibration timescale τeq is defined as the time required for the system to achieve 63.2 % of ∆Φ, the difference between Φ0 and Φeq.
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.

175
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).

180
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 reactions 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 185 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 Gas-phase loss is faster than particlephase loss Particle-phase loss is faster than gas-phase loss Particulate fraction larger than equilibrium value

Initial State
Particulate fraction at equilibrium value Particle-phase loss is faster than gas-phase loss Gas-phase loss is faster than particlephase loss

A B C
Particulate fraction smaller than equilibrium value 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).
as a function of atmospheric conditions, 2) compare these timescales to typical chemical loss rates in order to investigate 190 whether perturbations from equilibrium exist, and 3) explore the implications of treating partitioning and chemistry separately in chemical transport models.

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 5 cm −3 and particles 195 with a diameter of 50 nm.

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 num-  ber concentration until a fairly high threshold particle number concentration is achieved (≈ 10 5 , 10 4 and 10 4 particles cm −3 , 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 2 particles cm −3 ). This can be understood when looking at 205 Fig. 3b, in which the equilibration timescale of pyrene is shown together with the individual timescales of desorption τ des (gray dashed line, calculated with Eq. 6) and adsorption τ ads (gray dotted line, calculated with Eq. 7).
In the limit of an adsorbate-free surface, adsorption and desorption are first-order processes with respect to the near-surface 210 gas and surface concentrations of PAH respectively and can therefore be described with rate coefficients k ads (s −1 ) and k des (s −1 ). 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 2 p π with diameter d p and the mean thermal velocity ω. The surface coverage θ s is very small for typical 215 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 for details of the terms and derivation). This approximation holds as long as gas diffusion is sufficiently fast and does not limit equilibration.
If one process (desorption or adsorption) dominates the behavior of τ eq , the system can be said to fall into an adsorption-220 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 analagously 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 225 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 230 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.

235
As surface coverages θ s are very small and PAHs generally have surface accommodation coefficients on an adsorbate-free substrate of α s,0 = 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.  highlighting the transition between adsorption-controlled and desorption-controlled behavior.

240
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 245 undergo a significant change in equilibration timescale in the investigated temperature range.
Again, the adsorption-controlled and desorption-controlled regimes explain this behavior (Fig. 4b). Between 210 K 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 (gray dotted 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, 255 highlighted for pyrene with red shading in Fig. 4b and signified by the coincidence with the timescale of desorption τ des (gray dashed line).

Interplay of multiphase chemistry and partitioning
Chemical reactions with O 3 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).

260
This effect is exemplified for pyrene by including surface chemistry with O 3 (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. 10 ppb is representative of surface background O 3 concentrations (Vingarzan, 2004), while 100 ppb O 3 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 265 (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 p = 1×10 3 particles cm −3 , d p = 50 nm. At the start of the simulation, the initial total concentration of pyrene (5×10 5 cm −3 ) is distributed between the gas and particle phases according to the particulate fraction expected at equilibrium (i.e. Φ 0 = Φ eq = 0.24).  This effect can be understood by observing the change in particulate fraction over time during each of the simulations 280 (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 3 concentrations of 10 ppb 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 285 the O 3 concentration is low enough (1 ppb) the particulate fraction remains approximately equal to its value at equilibrium (Φ ≈ Φ eq = 0.24). At this O 3 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.

290
The lifetimes of pyrene with OH concentrations of 0.01 to 0.1 and 1 ppt, are 18.9 h, 1.9 h 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 3 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 295 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 Figure 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 OHpyrene system compared to its partitioning timescale: pyrene reaches 1/e of its initial concentration before the quasi-steady 300 state is established.

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 trans-305 ferred 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 310 equilibrium particulate fraction Φ eq can be obtained (Eq. 9).
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 perturbed and equilibrium system is depicted for [O 3 ] = 100 ppb in Fig. 6b. The vector field fundamentally changes and the trajectory of an 315 exemplary simulated system (red solid line) does not converge to the nullcline obtained in Fig. 6a, despite starting at equilibrium conditions Φ eq (shown as black dashed line for reference). Instead, the trajectory converges onto a central trajectory termed the slow manifold (Fraser, 1988). All trajectories in this system (represented with gray dotted 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 320 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 cm cm cm cm cm cm reached only just before all pyrene is consumed (compare solid blue lines in Fig. 5). Fig. 6d shows how the nullcline and the 325 slow manifolds above or below it can be interpreted using the diagrams in Fig. 2.

Implications for chemical transport models (CTMs)
In the previous sections, an explicit, 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. Hereon, this is referred to as the explicitcoupled approach (EC). As the explicit-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).

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 3 (100 ppb) with pyrene in the particle phase, we demonstrate that the magnitude and sign of this difference varies 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, N p = 1 × 10 3 particles cm −3 .

340
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 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 345 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).

Deviation from explicit-coupled (EC) approach
The discrepancy between the EC and IE solutions not only depends 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 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 355 (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. 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 5×10 −8 and 5 × 10 −1 s −1 , and the concentration of O 3 between 0 and 365 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 1×10 −4 and 1×10 −2 s −1 . 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 less non-equilibrium effects of slow partitioning on multiphase chemistry compared to the reference EC model. At the highest k des (> 1 × 10 −2 s −1 ), non-equilibrium effects of 370 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 not only leads 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 IE and EC approach is reduced, as evident by the more faint red coloring. At low k des (< 1×10 −5 s −1 ), 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 375 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 3 ] 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 380 partitioning compared to the EC model. When the equilibration timescale becomes shorter, at high k des (> 1 × 10 −2 s −1 ), the discrepancy between the IE and the EC solution further increases, especially at high [O 3 ] g . Notably, at k des ≈ 1 × 10 −4 s −1 , the IE approach slightly overestimates loss of PAH at high [O 3 ] 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 des ≈ 1 × 10 −3 s −1 . 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 des ≈ 5 × 10 −5 s −1 , in which departure from equilibrium does not occur and both methods return truly identical results (compare Figs. A5a and A5c). Another region exists: for k des < 1 × 10 −4 s −1 , the simulation proceeds for less than 30 min and therefore less than one IE time step is evaluated (Fig. A5). In this region partitioning effectively does not take place and 390 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 in both the gas phase and on the surface of particles (Fig. A4a). Unlike the example with O 3 , 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 395 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 loss = -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 loss = -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 400 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 405 scheme.

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 tempera-410 ture, 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 −1 ) compared to soot (87.9 kJ mol −1 ; Chu et al., 2010), one 415 would expect equilibration timescales on NaCl to be shorter.
On other particle types, PAHs 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 420 longer than the equilibration timescales for adsorptive partitioning. Alongside the contributions to the equilibration timescale from the adsorption and the desorption controlling 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., in review, 2020) and departure from partitioning equilibrium adds to this complexity (Cappa and Wilson, 2011). However, the role of bulk diffusion in determining 425 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 3 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 3 , the largest deviations from equilibrium particulate fraction are likely to take place in the most polluted air. The reaction of pyrene 430 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 3 in both the gas phase  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 435 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) depend on the concentrations of OH, O 3 , NO 3 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 440 (Φ 0 = Φ eq ), in reality PAHs may be emitted in a state far from partitioning equilibrium. Depending on the prevailing loss processes, such an effect could both 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 445 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 days) 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.

450
The methodology described in this study is universally applicable to semi-volatile compounds on solid surfaces if masstransfer 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 455 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.
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.
The gas-phase diffusion coefficient D g is fixed at 0.06 cm 2 s −1 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 465 in Eq. A2.
The mean thermal velocity of a molecule ω depends on temperature T and its molar mass M (Eq. A3).
The adsorption flux J ads of molecules from the near-surface gas phase to the particle phase is described using Eq. A4.

470
J ads = α s,0 (1 − θ s )J coll (A4) 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 α s,0 = 1 (Julin et al., 2014 and σ O3 of PAH and O 3 , respectively (Eq. A6). In order to estimate σ, each benzene-like ring of a PAH molecule is assumed to occupy 2 × 10 −15 cm 2 . The collision flux J coll , i.e. the flux of molecules colliding with the surface, is defined in Eq. A7.
The temperature dependent desorption flux (J des ) due 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).
k des depends on the Arrhenius factor (A) and the activation energy for desorption from the aerosol particle surface (E A ; Eq. A9).
The temperature dependence of the desorption rate coefficient was previously determined for seven PAHs on fresh kerosene soot (Guilloteau et al., 2008(Guilloteau et al., , 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 soots can have different effects 490 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.
Irreversible reactions between pyrene and either O 3 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 3 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 495 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).
[OH] gs = C g,OH [OH] g (A10) 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 500 by OH L g,OH and the rate of PAH loss from the surface due to reaction with O 3 L s,O3 are defined by Eq. A12 and Eq. A13, respectively. Further details of these reactions and their parameters can be found in section 2.2.
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 −1 ) and k ads (s −1 ). We find there is good agreement between the numerically-obtained results and the 510 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 Both sides are integrated (Eq. B11) and the equation rearranged (Eq. B12).