The Secondary Organic Aerosol Processor ( SOAP v 1 . 0 ) model : a unified model with different ranges of complexity based on the molecular surrogate approach

The Secondary Organic Aerosol Processor (SOAP v1.0) model is presented. This model is designed to be modular with different user options depending on the computing time and the complexity required by the user. This model is based on the molecular surrogate approach, in which each surrogate compound is associated with a molecular structure to estimate some prop5 erties and parameters (hygroscopicity, absorption into the aqueous phase of particles, activity coefficients and phase separation). Each surrogate can be hydrophilic (condenses only into the aqueous phase of particles), hydrophobic (condenses only into the organic phase of particles) or both (condenses into both the aqueous and the organic phases of particles). Activity coefficients are computed with the 10 UNIFAC thermodynamic model for short-range interactions and with the AIOMFAC parameterization for medium and long-range interactions between electrolytes and organic compounds. Phase separation is determined by Gibbs energy minimization. The user can choose between an equilibrium representation and a dynamic representation of the organic aerosol. In the equilibrium representation, compounds in the particle phase are 15 assumed to be at equilibrium with the gas phase. However, recent studies show that the organic aerosol (OA) is not at equilibrium with the gas phase because the organic phase could be semisolid (very viscous liquid phase). The condensation or evaporation of organic compounds could then be limited by the diffusion in the organic phase due to the high viscosity. A simplified dynamic representation of secondary organic aerosols (SOA) is used with OA divided into layers, 20 the first layer being at the center of the particle (slowly reaches equilibrium) and the final layer being near the interface with the gas phase (quickly reaches equilibrium). According to preliminary results, the non-ideality of the aerosol (including interactions with inorganic ions and separation of the organic phase into several phases) should be taken into account in models as it may affect SOA formation. Moreover, some compound should not be 25 assumed only hydrophilic or only hydrophobic as they can condense on both the organic and the aqueous phases.

Abstract.In this paper the Secondary Organic Aerosol Processor (SOAP v1.0) model is presented.This model determines the partitioning of organic compounds between the gas and particle phases.It is designed to be modular with different user options depending on the computation time and the complexity required by the user.This model is based on the molecular surrogate approach, in which each surrogate compound is associated with a molecular structure to estimate some properties and parameters (hygroscopicity, absorption into the aqueous phase of particles, activity coefficients and phase separation).
Each surrogate can be hydrophilic (condenses only into the aqueous phase of particles), hydrophobic (condenses only into the organic phases of particles) or both (condenses into both the aqueous and the organic phases of particles).Activity coefficients are computed with the UNIFAC (UNIversal Functional group Activity Coefficient; Fredenslund et al., 1975) thermodynamic model for short-range interactions and with the Aerosol Inorganic-Organic Mixtures Functional groups Activity Coefficients (AIOMFAC) parameterization for medium-and long-range interactions between electrolytes and organic compounds.Phase separation is determined by Gibbs energy minimization.
The user can choose between an equilibrium representation and a dynamic representation of organic aerosols (OAs).In the equilibrium representation, compounds in the particle phase are assumed to be at equilibrium with the gas phase.However, recent studies show that the organic aerosol is not at equilibrium with the gas phase because the organic phases could be semi-solid (very viscous liquid phase).The condensation-evaporation of organic compounds could then be limited by the diffusion in the organic phases due to the high viscosity.An implicit dynamic representation of secondary organic aerosols (SOAs) is available in SOAP with OAs divided into layers, the first layer being at the center of the particle (slowly reaches equilibrium) and the final layer being near the interface with the gas phase (quickly reaches equilibrium).Although this dynamic implicit representation is a simplified approach to model condensation-evaporation with a low number of layers and short CPU (central processing unit) time, it shows good agreements with an explicit representation of condensation-evaporation (no significant differences after a few hours of condensation).

Introduction
Fine particles are regulated because of their impact on human health (WHO, 2003).Furthermore, they degrade atmospheric visibility (Larson et al., 1989) and influence climate change (Kanakidou et al., 2005).It is therefore necessary to develop models able to predict particle formation, which can be used to predict their impact on health and environment and evaluate emission mitigation policies.Particulate organic matter (OM) represents a large fraction of the total fine particulate mass, typically between 20 and 60 % (Kanakidou et al., 2005;Yu et al., 2007;Zhang et al., 2007).Therefore, efforts have to be made to represent OM as accurately as possible in models.Three-dimensional (3-D) air quality models, which estimate particle concentrations, need to have an implicit rep-F.Couvidat and K. Sartelet: The SOAP model resentation of organic species.Because of the large number of organic species involved originating from diverse anthropogenic and biogenic sources, species need to be lumped according to their properties (e.g., by lumping species with similar saturation vapor pressures).In the surrogate-based methodologies, molecular structures are attached to surrogate compounds representing a large number of secondary organic aerosol (SOA) species to estimate several properties (e.g., condensation into an aqueous phase, oligomerization, hygroscopicity, non-ideality).
In 3-D air quality models, several assumptions are made on the thermodynamics of organic aerosols (OAs) such as equilibrium between the gas phase and the particle phase, ideality or no phase separation.However, these assumptions could strongly impact simulated OA formation.For example, some recent experimental studies emphasize the need to take into account dynamical aspects of the formation of OAs rather than assuming thermodynamic equilibrium with the gas phase because OAs can be highly viscous (Virtanen et al., 2010;Cappa and Wilson, 2011;Pfrang et al., 2011;Shiraiwa et al., 2011;Vaden et al., 2011;Shiraiwa and Seinfeld, 2012;Abramson et al., 2013).
Some OA models already represent the formation and the condensation of organic compounds using the surrogate approach: the Audio Electric Research GmbH (AER)/Electric Power Research Institute (EPRI)/Caltech (AEC) model (Pun et al., 2002(Pun et al., , 2006)), the Hydrophilic/Hydrophobic Organic model (H 2 O) (Couvidat et al., 2012(Couvidat et al., , 2013) ) and the Model to Predict the Multiphase Partitioning of Organics (MPMPO) (Griffin et al., 2003).However, none of these models take into account the dynamics of the condensation of organic compounds, the influence of interactions between organic and inorganic compounds on activity coefficients or the phase separation of particulate OM into several organic phases, whereas the interplay of these phenomena should be taken into account in models (Shiraiwa et al., 2013).Moreover, contrary to MPMPO AEC and H 2 O assume that each compound may condense into only one phase (the organic or the aqueous phase).Computation of activity coefficients and phase separation at equilibrium has been extensively developed in the thermodynamic model Aerosol Inorganic-Organic Mixtures Functional groups Activity Coefficients (AIOMFAC) (Zuend et al., 2008(Zuend et al., , 2010(Zuend et al., , 2011;;Zuend and Seinfeld, 2012).Shiraiwa et al. (2012) developed a multilayer model KM-GAP (kinetic multi-layer model of gasparticle interactions in aerosols and clouds) which explicitly treats the condensation and particle diffusion of organic compounds as well as heat transfer and particle-phase reactions.
To represent organic aerosol formation and take into account non-ideality, phase separation and the viscous state of OAs, the Secondary Organic Aerosol Processor (SOAP), destined to be implemented in 3-D air quality models, is developed and presented here.This model is designed to be modular with different user options depending on the computation time and the complexity required by the user.The user can choose between an equilibrium representation and a dynamic representation of organic aerosols, between ideality and non-ideality (with or without phase separation and interactions with inorganic ions) and on which phases each surrogate compound can condense (the organic phases, the aqueous phase or all phases).This paper describes the SOAP model and the results of several test cases.As the dynamic representation of SOAP is implicit in order to work with low numbers of particle layers and short CPU times, comparisons to an explicit representation of condensation-evaporation are also presented.

Overview
SOAP is based on the surrogate approach in which SOA species are divided into three types: hydrophilic species (which condense only into an aqueous phase when an aqueous phase is present, i.e., when the concentration of water in aerosol is higher than 0.01 µg m −3 ), hydrophobic species (which condense only into organic phases due to their low affinity with water) or both (which condense into both phases) as chosen by the user.The model can represent the non-ideality of aerosols (interactions between organic compounds and interactions between organic and inorganic compounds represented via activity coefficients), hygroscopicity, phase separation and formation of SOA following an equilibrium approach or a dynamic approach.
Semi-volatile organic compounds (SVOCs) are represented by surrogate compounds.To represent the nonideality of aerosols, molecular structures are assigned by the user to each surrogate compound.A default structure is provided for each surrogate.This structure is used to compute the non-ideality of aerosols via activity coefficients.If the user specifies that a compound is both hydrophilic and hydrophobic, the repartition between phases is done according to the value of activity coefficients.However, the user can decide that a compound is only hydrophilic or only hydrophobic (for example alkane or lowly oxidized compounds are probably not absorbed by the aqueous phase of particles).Moreover, if there is no compound that is both hydrophilic and hydrophobic, the condensation into organic phases can be solved separately from the condensation into the aqueous phase.The system is then uncoupled.On the other hand, if there is at least one compound which is both hydrophilic and hydrophobic, the condensation into organic phases and the condensation into the aqueous phase must be solved simultaneously.The system is then coupled and consumes more CPU time.The user may prefer an uncoupled system for 3-D application due to higher time efficiency.
To compute SOA formation, two different approaches can be used to model gas/particle partitioning: the equilibrium approach and the dynamic approach.In the equilibrium ap-proach, aerosols are assumed to be at thermodynamic equilibrium with the gas phase.The model then uses the Newton-Raphson method to efficiently compute the partitioning of each compound between the gas and particle phases.In the dynamic approach, condensation and particle diffusion of organic compounds are treated with a multi-layer representation of OAs (the organic mass is divided into several layers having different characteristic times to reach equilibrium).In this method, the particle size distribution is divided into sections (inside a section/bin, all particles are assumed to have the same diameter).Inside a bin, compounds condense into the aqueous phase or/and the organic phases.Furthermore, each bin can be separated into several layers (the number of layers being the same for each bin) and several phases (the number of phases can change from one bin/layer to another bin/layer).
SOAP is based on the methodology of the H 2 O and AEC model.However, as described above, several processes were added.The model can have species that can condense into the organic phases and the aqueous phase (like in MPMPO).It can compute the effects of the interactions between inorganic and organic compounds on the condensation of organics, phase separation of the organic aerosol, and dynamic formation of SOA.The dynamic approach developed in this study is an implicit approach to take into account particlephase diffusion with a low number of layers and to keep a computation time as low as possible in order to be used by 3-D air quality models.It does not describe particle-phase diffusion as thoroughly as KM-GAP (Shiraiwa et al., 2012).

Organic aerosol formation at equilibrium
The fundamental equations used in SOAP to represent the partitioning between gas and particles under the equilibrium assumption are described below.

Equilibrium between the gas phase and one organic phase
The equilibrium between the gas phase and one organic phase is described by Raoult's law: where P i is the partial pressure of the compound i, γ i,org the activity coefficient of i in the organic phase, X i,org the molar fraction of i in the organic phase and P 0 i the saturation vapor pressure of i. Pankow (1994) rewrote Eq. (1) (see Eq. 2) to compute the absorption of organic compounds by an organic phase: where A p,i is the concentration of i in the organic phase (in µg m −3 ), A g,i the concentration of i in the gas phase (in µg m −3 ), M o the concentration of the organic phase (in µg m −3 ) and K p,i the organic-phase partitioning coefficient (in m 3 µg −1 ) which is computed using Eq.(3) (Pankow, 1994).
where T is the temperature (in K), M ow the mean molar mass of the organic phase (in g mol −1 ) and P 0 i the saturation vapor pressure (in torr).In SOAP, activity coefficients are computed with the thermodynamic model UNIFAC (UNIversal Functional group Activity Coefficient; Fredenslund et al., 1975).Moreover, P 0 i or partitioning constants K p,i are the same as those used in H 2 O (they are generally determined by fitting to experimental results obtained in environmental chambers at a temperature T ref ).The partitioning coefficient is extrapolated from T ref to T by using the enthalpy of vaporization H i (in J mol −1 ) according to the Clausius-Clapeyron equation.
The absorption of water by the organic phase is computed with Eq. ( 5) derived from Eq. ( 1) applied for water.
where RH is the relative humidity, M water the molar mass of water (in g mol −1 ) and γ water,org the activity coefficient of water in the organic phase.

Equilibrium between the gas phase and an aqueous phase
SOAP does not currently take into account the formation of inorganic aerosols.An inorganic aerosol model like ISOR-ROPIA (Nenes et al., 1998) must be called separately and prior to the call of SOAP to provide inputs to SOAP: pH, concentrations of inorganic ions, ionic strength and the liquid water content of aerosols.
The equilibrium between the gas and aqueous phases can be described not only by Raoult's law (Eq.6) but also by Henry's law (Eq.7) if infinite dilution is assumed: where γ i,aq is the activity coefficient of i in the aqueous phase and X i,aq the molar fraction of i in the aqueous phase.
where C i is the concentration (in M) of species i in the aqueous phase, P i in atmosphere and H i the Henry's law constant (in M atm −1 ).

F. Couvidat and K. Sartelet: The SOAP model
As the Henry's law is often used to express the partitioning between the gas and aqueous phases, a modified Henry's law is used to extrapolate infinite dilution conditions to all conditions using an aqueous-phase partitioning coefficient K aq,i : where A aq,i is the concentration of i in the organic phase (in µg m −3 ), A g,i the concentration of i in the gas phase (in µg m −3 ), AQ the total concentration (organics + inorganics including water) of the aqueous phase (in µg m −3 ) and K aq,i the aqueous-phase partitioning coefficient (in m 3 µg −1 ) which is computed with Eq. ( 9): where ρ water is the density of the aqueous phase (in kg m −3 ), M aq the molar mass of the aqueous phase (in g mol −1 ), which can be slightly different from the molar mass of water due to the presence of other compounds, and ζ i the activity coefficient by reference to infinite dilution.ζ i is computed with Eq. ( 10): where γ ∞ i,aq is the activity coefficient at infinite dilution in water, which is computed with UNIFAC.However, the original UNIFAC (Fredenslund et al., 1975) only computes the activity coefficients due to short-range interactions and does not take into account medium-and long-range interactions due to the presence of electrolytes in the aqueous phase.In the aqueous phase, activity coefficients are computed from Eq. ( 11) (Zuend et al., 2008) γ LR , γ MR and γ SR are respectively the activity coefficients at long-, medium-and short-range interactions.γ SR is computed with UNIFAC, whereas γ LR and γ MR are computed with the AIOMFAC method.The last two parameters model the influence of inorganic ions on the partitioning of organic compounds (Zuend et al., 2008(Zuend et al., , 2011;;Zuend and Seinfeld, 2012).
Similarly to the case of condensation into an organic phase (Eq.4), the partitioning coefficient is extrapolated from T ref to T by using the enthalpy of vaporization H i (in J mol −1 ) as described in Eq. (12): However, some compounds are acids that can dissociate in the aqueous phase.Therefore, partitioning coefficients are modified to take into account acidic dissociation as done by Pun et al. (2006).
The absorption of water by the aqueous phase is computed with Eq. ( 13) derived from Eq. (1) A aq,water = M water AQ × RH γ water,aq M aq (13) where γ water,aq is the activity coefficient of water in the aqueous phase.As the amount of water absorbed by inorganics is already given by the inorganic model (i.e., ISORROPIA) and is used as input for SOAP, we assumed that the total amount of water (from inorganics and organics) computed by SOAP should at least be equal to the amount of water given by the inorganic model.Therefore, if the amount of water computed by SOAP is lower than the amount computed by the inorganic model, it is replaced by the amount computed by the inorganic model.This problem arises because inorganics are not computed within SOAP (in that case, there would be no numerical issue) and water concentrations can be lower in SOAP than in the inorganic model ISORROPIA for several reasons.First, there can be numerical differences between SOAP and ISORROPIA because different parameterizations are used in SOAP and ISORROPIA to compute water concentrations: SOAP computes the amount of water using the water activity computed by the thermodynamic models UNIFAC and AIOMFAC, whereas ISORROPIA computes the amount of water using the Zdanovskii-Stokes-Robinson (ZSR) correlation (Robinson and Stokes, 1965).Second, the amount of water absorbed by the aerosol can be different from the sum of the water amount absorbed by inorganics and the water amount absorbed by organics.Choi and Chan (2002) found that organic species can either reduce or enhance the water absorption of inorganic compounds.Depending on the conditions, the amount of water computed by SOAP could be higher without organics than with.Here we chose to assume that the amount of water should at least be equal to the amount of water given by ISORROPIA, although it would be possible to keep the amount of water computed by SOAP even when lower than the amount of water computed by ISORROPIA.However, if the amount of water given by SOAP is significantly lower than ISORROPIA, this would induce changes in the amount of inorganics in the aerosol.As the amount of inorganics may not be recomputed after SOAP, we considered that the amount of water is at least equal to the amount of water given by ISORROPIA; that is, the amount of water being recomputed only to provide a better estimate of the amount of water when water absorption are mainly due to organics.The best way to deal with this dilemma would be to fully couple inorganic and organic aerosol formation.

Equilibrium between the gas phase and several particulate phases
SOAP can compute the partitioning of compounds between the gas phase and several particulate phases.The user can specify for each compound if it is hydrophilic (condense into the aqueous phase of particles) or hydrophobic (condense into the organic phase of particles) or both hydrophilic and hydrophobic (condense into both phases).For an uncoupled system (compounds cannot be both hydrophobic and hydrophilic), for hydrophobic compounds, the total concentration of i in all phases A tot,i is which when combined with Eq. (2) results in A method of Newton-Raphson is then used to solve Eq. ( 15) and to minimize the error (an accuracy threshold is provided by the user): For an uncoupled system, for hydrophilic compounds, the total concentration of i in all phases A tot,i is which when combined with Eq. ( 8) results in A method of Newton-Raphson is then used to solve Eq. ( 18) and to minimize the error (an accuracy threshold is provided by the user): where AQ inorg is the concentration of inorganic compounds in the aqueous phase.For a coupled system, the total concentration of i in all phases A tot,i is which when combined with Eqs. ( 2) and (8) results in Similarly, a method of Newton-Raphson is used to simultaneously solve Eqs. ( 21) and ( 22) and to minimize the errors (error1 p and error aq ).

Saturation and separation of organic phases
If compounds having a low affinity with each other coexist inside a single organic phase, the organic phase may become saturated by some compounds and may become unstable.In this case, the separation of the organic phase into several organic phases (having different composition) may occur.To determine whether separating the organic phase into several organic phases makes the system more stable, the Gibbs energy G is computed as in Erdakos and Pankow (2004), Zuend et al. (2010), and Zuend and Seinfeld (2012) using different system configurations (different number of phases): where φ is an index of the various phases (gaseous and liquid), n φ i the number of moles of species i in phase φ and µ φ i the chemical potential of i in phase φ.
The most stable configuration of the system has the lowest Gibbs energy.Therefore, by adding one organic phase the Gibbs energy decreases, the system is more stable and phase separation takes place.If the Gibbs energy does not decrease, the previous solution (before adding one organic phase) is more stable and is kept.The number of organic phases is determined iteratively: one phase is added until the Gibbs energy increases.
The partitioning of organic compounds into an aerosol with several organic phases are determined with Eqs.21 and 22 generalized to several organic phases:

Dynamic gas uptake by organic particles
The dynamic approach, which is presented hereafter, is an implicit representation to take into account particle-phase diffusion with a low number of layers, and have a computation time as low as possible to be used in 3-D air quality models.The main assumptions are described here.
The first assumption is that the organic-phase diffusion coefficient is constant over the entire particle.It does not depend on the distance to the center of the particle.Although this assumption may not be valid, it is reasonable because, currently, to our knowledge, there is no parameterization to evaluate the order of magnitude of organic-phase diffusion coefficients.
The second assumption is that in the model there is no direct exchange of compounds between layers, and that the compounds condense directly from the gas phase to the layer or that they evaporate directly from the layer to the gas phase by taking into account an equilibration time specific of the www.geosci-model-dev.net/8/1111/2015/F. Couvidat and K. Sartelet: The SOAP model layer.Compounds condense into a layer or evaporate from a layer as if the other layers had the same affinity with compounds.Effects of entrapment of compounds inside central layers by the layers closer to the gas-phase interface (compounds inside the central layers having a low affinity with the compounds of the layers at the interface will not be able to evaporate, whereas compounds having a high affinity with the central layer but having a low affinity with the layers at the interface will not be able to condense into the central layer) are not taken into account.The model should give, however, a good estimation of the capacity of the organic aerosol to absorb compounds (by taking into account the time for each layer to reach equilibrium due to diffusion in the particle).
The third assumption is that the organic phases and the aqueous phase evolve separately, i.e., there is no kinetic transfer of compounds between the organic phases and the aqueous phase, due to the complexity of representing properly those transfers, which should strongly depend on the morphology of particles.If a compound tends to go from the aqueous phase to the organic phases, it has first to go from the aqueous phase to the gas phase and then from the gas phase to the organic phases.For example, in case of evaporation of the aqueous phase (that can be due to a strong change of the RH), this assumption can create some evaporation/recondensation issue in the model (compounds evaporate and recondense after some time into the organic phases according to condensation-evaporation fluxes).Actually, a part of organic compounds should go directly from the aqueous phase to the organic phases.It could also be possible that if an aqueous phase coexists with an organic phase into the same particle, an organic compounds will not condense directly from the gas phase into the organic phases (because the kinetic is too slow) but condense first into the aqueous phase and then go from the aqueous phase to the organic phases (if it is quicker for a compound to condense into the organic phases by this pathway).However, it can also be argued that if there is an aqueous phase (the RH is high), the organic phases may not be significantly viscous and therefore a high organic-phase diffusion coefficient should be used.
Finally, the model assumes that there is no gradient of the gas-phase concentrations near the interface with the particle.If the particle is divided into two separated regions (one aqueous and one organic), an angular gradient of gas-phase concentrations could influence the condensation of compounds into the two regions.To address properly this phenomenon, the particle and the gas phase at the vicinity of the particle should be discretized as a function of the angular gradient.However, due to the high diffusivity of organic compounds in the gas phase, the diffusion of organic compounds around the particle should be very quickly compared to the kinetic of condensation and therefore this assumption should have a low effect on condensation-evaporation.

Diffusion of organic compounds in spherical organic particles
The diffusion of organic compounds at a radius r at time t inside a spherical particle is governed by the following equation (Seinfeld and Pandis, 1998) where C(r, t) is the molar concentration at radius r at time t, D org the organic-phase diffusivity and R org the organicphase reaction rate.By assuming D org constant, the solution of this equation (with C(r, 0) = 0 and C(R p , t) = C s , R p being the radius of the particle) without organic-phase reaction is according to Seinfeld and Pandis (1998): By integrating Eq. ( 29) over the volume of a spherical particle, the following equation is found for the concentration in the particle phase A p : where τ dif is the characteristic time for diffusion in the center of the particle: where A eq is the organic-phase concentration at equilibrium with the gas phase (A eq = K p M o A g ).In Eq. ( 30), A p can be interpreted as the sum of an infinity number of layers of concentration A layer p : where is the concentration (in µg m −3 ) in a layer of volume fraction V layer determined by the fraction of the volume of the particle constituted by the layer and α layer the ratio between the characteristic diffusion time at the center of the particle τ dif to the characteristic diffusion time τ layer dif of the layer such as Numerically, Eq. ( 30) can be approached by a finite number of diffusion layers N dif layer (number of layers into which diffusion occurs) and by fitting the parameters V layer and α layer such as First, for a given concentration A eq , a given D org or R p , we estimate A p (t) using Eq. ( 32).Then V 1 , V 2 , V 3 , α 1 , α 2 , α 3 are fitted such as satisfying Eq. ( 35).The values of these parameters do not depend of the choice of A eq , D org and R p .
The evolution of concentrations A bin,layer p,i in a bin and a layer due to the condensation of species i limited by the diffusion of organic compounds in the organic phase is described by the Eq. ( 36), which describes the flux of diffusion J bin,layer diff,i , i.e., the evolution of the concentration A ) (with the assumption that the concentrations in one layer can be described independently from the concentrations in the other layers).
where M bin,layer o is the total mass of the organic phase in the layer computed with Eq. ( 37) by assuming that the density of the OM is constant over layers and k diffusion the kinetic coefficient for diffusion (in s −1 ) computed with Eq. ( 38) as the inverse for the characteristic diffusion time in a layer.

Diffusion of organic compounds in more complex particles
The previous equations correspond to the diffusion of organic compounds into an entirely organic spherical particle.However, in the atmosphere, particles generally have more complex geometries and can also be constituted by solid and/or inorganic phases.The morphology affects the time for an organic phase to be diffused in the particle.For example, for a same particle diameter, a particle entirely organic needs more time to reach equilibrium with the gas phase than a particle constituted by a solid core in the center surrounded by an organic phase (because the organic compounds do not penetrate the particle all the way to the center, they are only diffused inside the coating).As diffusion of organic compounds is affected by the presence of solid phases, morphology factors are designed to take into account this solid phase where diffusion of organics may not occur.We assume that this solid phase is located at the core of particles based on Katrinak et al. (1993), Sachdeva and Attri (2007), and Wang et al. (2014).Models, such as KM-GAP, treating explicitly the diffusion of compounds inside particles do not need those factors.As in SOAP, the volume of layer has to be constant, whereas only the characteristic time for diffusion can be affected by the morphology.We propose here to parameterize the effects of a solid core by defining morphology factors f layer morph , which affect the time to reach equilibrium with the gas phase.The characteristic time for diffusion of a layer τ layer dif expressed in Eq. ( 34) becomes Equation ( 36) then becomes with k diffusion defined as in Eq. ( 38).
The morphology factors can be determined numerically (at least for some simple case).We determined here morphology factors in the case of a spherical particle with a solid core at www.geosci-model-dev.net/8/1111/2015/ the center of the particle.The differential equation for diffusion (Eq.28) is solved and morphology factors are fitted to minimize the differences between Eqs. ( 28) and ( 40) for various volume fractions of the solid phase f s (volume of the solid phase in the particle over the volume of the particle).This fitting procedure is described in more detail in the Supplement.The variations of the morphology factors with the volume fraction f s are shown in Fig. 2 for three and four diffusion layers.The morphology factors for a solid core particle can be represented by a polynomial expression such as The values of the polynomial parameters are shown in Table 2. Typically, the presence of a solid phase would result in a morphology factor greater than 1, reducing the characteristic diffusion time of organics in the particle.This parameterization only takes into account the simple case of a spherical solid core at the center and not the wide range of morphologies present in the atmosphere.A similar methodology or a methodology based on observations could be applied to other morphologies which could greatly affect the characteristic times for diffusion.For example, morphology factors for the effect of the presence of an aqueous phase on characteristic time for diffusion could be determined (at least as a function of the mass of the aqueous phase and the affinity of compounds with water).Currently, the model assumes there is no effect of the aqueous phase on the characteristic times.

Condensation-evaporation of organic compounds into a viscous particle
In this part, the methodology to compute condensationevaporation-diffusion inside a viscous particle is described.
The algorithm to compute condensation-evaporationdiffusion fluxes is divided into two steps shown in Figs. 3   and 4. In the first step, the total flux of condensationevaporation-diffusion over the whole particle is computed whereas in the second step this flux is redistributed among layers.The complete methodology to compute fluxes is described hereafter.
The evolution equation of concentration of species i by diffusion in each bin and layer is described by Eq. ( 40).To properly take into account the condensation-evaporation of low volatile compounds, a thin surface layer, which is not limited by diffusion, is added to the model.In case of very low diffusivity, this layer allows non-volatile and low volatility compounds to condense at the surface of particles without any limitation due to diffusion.Let us denote the number of layers N layer into which the dynamic of condensationevaporation-diffusion is solved.We have A volume fraction V layer for the interface layer of 0.01 is used (using a lower volume fraction does not significantly change the results).Adding an interface layer between the gas and the particle slightly reduces the volume of the other layers, as described in Sect.2.3.1.In this paper, simulations are performed with four layers (three diffusion layers: V 1 = 0.6, V 2 = 0.26, V 3 = 0.13 and V 4 = 0.01) and with five layers (four diffusion layers: The flux of diffusion in the upper layer (the interface) is zero, as diffusion should not limit absorption at the interface: The evolution of concentrations by condensationevaporation at the interface, i.e., the flux of condensationevaporation at the interface, is where A eq,bin,interface g,i is the gas-phase concentration at equilibrium with the interface layer and k absorption the kinetic rate Table 1.Algorithm to compute the partitioning of compounds at equilibrium in the dynamic approach.
While the system has not converged (or has not reached a maximum number of iterations) 1. n = 1. 2. Factor = 1/n.3. Compute the new estimation of the concentrations A bin,layer,phase p,i,new according to Eq. ( 81).  of absorption, which is (Seinfeld and Pandis, 1998) where D p is the diameter of the particle, D air the diffusivity of compound i in air and f (Kn, α) the transition regime formula of Fuchs and Sutugin (1971).
The gas-phase concentration at equilibrium with the interface layer A eq,bin,interface g,i is computed as follows: where K bin,interface p,i is the partitioning coefficient into the interface layer for compound i and M bin,interface o the total organic mass concentrations of the interface layer and K effect the Kelvin effect correction coefficient: where σ is the surface tension, ρ organic the density of the organic phase, M ow the mean molar mass of the organic phase at the surface of the particle (molar mass of the layer at the interface) and R p the radius of the particle.A surface tension of 24 mN m −1 is chosen, which is roughly the surface tension of organic compounds according to Seinfeld and Pandis (1998) and a density of the organic phase of 1300 kg m −3 is used.
In a comprehensive dynamic model, diffusion would be represented explicitly.However, to have a model with a low number of layers and gain CPU time, the evolution of the concentrations in each layer is solved by considering the combined effect of condensation-evaporation and diffusion in the particle.The flux of condensation-evaporationdiffusion of a compound i inside a bin and a layer is noted J bin,layer tot,i and is used to compute the evolution of concentrations of a compound i inside a layer: The evolution of the total mass of the particle J bin tot,i can be computed by assuming that the characteristic time of the combined effect of condensation-evaporation and diffusion is equal to the sum of the characteristic times of condensation-evaporation and diffusion: 1 where J bin cond/evap,i is the flux limited by condensationevaporation and J bin diff,i the flux limited by diffusion.J bin diff,i is the sum of the diffusion fluxes over all bins (each of them being computed by Eq. 40): J bin cond/evap,i is obtained similarly to the interface flux of condensation-evaporation where A eq,bin g,i is the gas-phase concentration at equilibrium with the whole particle.
The gas-phase concentration at equilibrium with the particle is assumed to correspond to the sum of the gas-phase concentrations at equilibrium with each layer weighted by their kinetic of diffusion: where A eq,bin,layer g,i is the gas-phase concentration at equilibrium with the layer and with f bin,layer i a weighting factor that depends on the kinetic of diffusion: where K effect is the Kelvin effect correction coefficient computed with Eq. ( 47).
The coefficient f bin,layer i corresponds to the fraction of the compound i that arrives by diffusion in the layer if compound i condenses or the fraction of the compound i that departs from the layer if compound i evaporates.The weighting factor is computed with the following algorithm, which separates positive fluxes from negative fluxes.In the case of condensation of compound i (J bin tot,i > 0), the compound condenses preferentially on layers with positive diffusion fluxes J bin,layer diff,i .Therefore, if J bin,layer diff,i is negative inside a layer, this layer is not taken into account in the condensation flux.On the other hand, in the case of evaporation of compound i (J bin tot,i < 0), the compound evaporates from layers with negative J Once the weighting factors are computed, the flux of condensation-evaporation-diffusion over all layers J bin tot,i can be computed with Eqs. ( 49)-( 51).These weighting factors are used not only to redistribute the gas-phase concentration at equilibrium between layers from a bulk gas-phase concentration (Eq.52), but also to redistribute the total flux of condensation-evaporation-diffusion between layers from the bulk total flux of condensation-evaporation-diffusion.The computation of the total flux per layer J bin,layer tot,i is now explained.The total flux of condensation-evaporationdiffusion J bin tot,i is separated into three fluxes J + i (equal to the sum of diffusion fluxes which are positive), J − i (equal to the sum of diffusion fluxes which are negative) and the flux at the interface J bin,interface i : By definition, J + i and J − i are equal to where J bin,layer tot,i is the total flux of condensation-evaporationdiffusion of a compound i inside a bin and a layer.
If the compound i condenses into the particle (J bin tot,i > 0), we also have J bin cond/evap,i > 0 for Eq. ( 51).In this case, if the diffusion flux J bin,layer diff,i is negative, then the compound moves by diffusion oppositely to condensation and therefore it should not be affected by the kinetic of condensation-evaporation.More generally, compounds and layers for which diffusion reacts in a different direction as condensation-evaporation (fluxes of different signs) should not be affected by the kinetic of condensation-evaporation (because the compounds are transferred from one phase to another without exchange with the gas phase).In other words, The layers with J bin cond/evap,i × J bin,layer diff,i < 0 are the layers with a diffusion flux sign opposite to J bin tot,i .Therefore, the sum over these layers allows us to compute J − i with Eq. (58) if J bin tot,i > 0 and J + i with Eq. (57) if J bin tot,i < 0. If J bin tot,i > 0 (resp.J bin tot,i < 0), then J + i (resp.J − i ) corresponds to the part of the total flux that is affected by both condensationevaporation and diffusion and J − i (resp.J + i ) corresponds to the part of the total flux that is only affected by diffusion.J − i is computed by Eqs. ( 59) and (58) (resp.J + i is computed by Eqs.59 and 57), and J + i is deduced by Eq. ( 56).
The total positive (resp.negative) flux J + i (resp.J − i ) is redistributed onto layers by using the weighting factors f bin,layer i :

Characteristic time to reach equilibrium with the gas phase
The system of differential to solve Eq. ( 48) is stiff, as in the same layer/bin, some species reach equilibrium with the gas phase much quicker than others (Capaldo et al., 2000).To solve it efficiently, it is necessary to separately solve cases at equilibrium from other cases, which are solved dynamically.To determine the characteristic time to reach equilibrium with the gas phase, Eq. ( 48) is rewritten using the total concentrations and defining F where A eq is the concentration at equilibrium with the gas phase and τ eq the characteristic time to reach equilibrium.
τ eq is used to estimate the time necessary to reach equilibrium with the gas phase and therefore to identify cases that should be assumed at equilibrium when solving the system of Eq. ( 48).

Generalization to several organic phases
The OM can be separated into several organic phases.
Whereas the evolution of condensation-evaporation is dynamically modeled, phase separation and the repartition of compounds between organic phases are assumed to be at equilibrium.They are assumed to occur instantaneously: if an organic phase is saturated, it is divided instantaneously into several organic phases.The dynamic evolution due to condensation-evaporation in viscous particles is described by the following equation derived from Eq. ( 48) by taking into account the phase of the components: dA is the total organic concentrations in the specified organic phase in a specified bin and layer (in µg m −3 ), J bin,layer,phase diff , A bin,layer,phase p,i and K bin,layer,phase p,i the flux of condensation-evaporation-diffusion, concentration and partitioning coefficient of the compound i for a bin, a layer and a phase, respectively.The flux at the interface is computed via Eq.( 44) using the gas-phase concentration at equilibrium with the interface: The flux of condensation-evaporation over the whole particle is computed via Eq.( 51) by using the gas-phase concentration at equilibrium with a layer: At each time step, thermodynamic evolution is first computed.The number of organic phases and the distribution of compounds between organic phases are then computed for each layer and each bin by assuming equilibrium between phases.To compute the concentrations with several organic phases at equilibrium we first study the conditions that have to be respected.The first condition is that the activities of each compound i are equal in each phase.For example, for two phases phase 1 and phase 2 : γ bin,layer,phase 1 i x bin,layer,phase 1 i = γ bin,layer,phase 2 i x bin,layer,phase 2 i . (70) The second condition, which is that each phase may be at equilibrium with the gas-phase (if condensation is too quick to be solved dynamically), can be written as To respect these two conditions, the Kelvin effect must be the same for each phase.
M bin surf , ρ organic and σ must be the same in the two phases.Therefore, if there are several organic phases, the partitioning coefficient must be computed with the following equation where M bin surf is the mean molar mass of all organic phases at the surface of particles.
The characteristic times are assumed to be the same for all the organic phases to prevent a compound from being absorbed dynamically into an organic phase when it is at equilibrium with another organic phase.The characteristic time

Absorption into the aqueous phase
For concentrations in the organic phases, the dynamic evolution follows Eq. ( 66) but the dynamic evolution of concentrations in the aqueous phase follows Eq. ( 76), because condensation-evaporation is assumed to not be limited by diffusion in the aqueous phase and a multi-layer representation is then not useful.
A bin aq,i K bin aq,i AQ bin (76)

Absorption into a particle with an aqueous phase and organic phases
Under most atmospheric conditions, particles are probably not entirely organic or entirely aqueous.The surface of particles is probably covered partially by both the OM and by the aqueous phase.Equations ( 76) and ( 66) are still valid but k absorption in Eqs. ( 76) and ( 45) must be corrected to take into account that there is a chance that a compound trying to condense into a phase encounters the wrong phase, i.e., a phase into which it may not condense.
The chance for a compound to encounter an aqueous phase fsurf aq is computed with Eq. ( 77) where S aq is the surface of particles that is aqueous and S tot the total surface of particles.
To evaluate this parameter, we assume that the surface of the particle is only covered by the aqueous phase and the organic phases and that the ratio of the aqueous surface over the organic surface is equal to the ratio of the volume of the aqueous phase over the volume of the organic phases: For the condensation in the aqueous phase, by taking into account the chance for a compound to encounter the aqueous phase, k absorption is For the condensation in the organic phases, by taking into account the chance for a compound to encounter the organic phases, k absorption is

Redistribution of compounds between layers
To use the approach described in this paper, the mass fraction (ratio of the mass of the layer over the mass of the particle) of layers must stay constant throughout the simulation and the mass of each layer must respect the condition given by Eq. ( 37), which specifies the mass of the layer with respect to the total mass of the particle.However, due to rapid condensation or evaporation of the layer near the interface, concentrations of organic compounds may need to be redistributed over layers to respect this condition.If compounds are not redistributed some layers may become larger (due to the differences in fluxes) and the layer near the interface may for example become larger than the layer at the center of the particle.Moreover, as the volume of the particle changes with condensation-evaporation, the concentration of a layer can be transferred to other layers.If a particle grows due to condensation into the layer at the interface, compounds that were previously in the interface layer are pushed into more internal layers and the newly condensed compounds remain at the www.geosci-model-dev.net/8/1111/2015/, the redistribution is done from the interface layer to the center of the particle: from ilayer = N layer to 2, ilayer2 = ilayer-1.interface layer.On the other hand, if a particle shrinks due to evaporation of the layer at the interface, the missing mass of the layer will be taken from more internal layers.For a case of evaporation, the mass of the layer at the interface may be too low (due to the more rapid evolution at the interface) and the missing mass of the layer is taken from layers at the inside of the particle, i.e., concentrations are redistributed from the outside to the inside.For a case of condensation, the mass of the layer at the interface may be too high, the exceeding mass of the layer is redistributed over the layers at the inside of the particle, i.e., concentrations are redistributed from the inside to the outside.The algorithm is detailed in Table 3.The redistribution algorithm that creates numerical diffusion as a small fraction of the mass of a layer is always transferred to other layers.This redistribution effect should decrease if the number of layers increases but it is necessary for a 3-D application to keep the number of layers low.However, the redistributed amount should be low compared to the absorbed amount of organic compounds.Similarly, in a 3-D air quality model, concentrations and number of particle in size sections have to be redistributed between sections so that the bounds diameters of sections are kept constant, which does also create some numerical issues.In SOAP, the diameters of sections evolve according to the mass that condenses or evaporates (without changing the number in sections) and there is no size redistribution between sections.However, a size-section redistribution algorithm should be added to the code if coagulation is added (for modeling purposes) or in 3-D models (in which SOAP would be implemented) after the call of SOAP.

Thermodynamic equilibrium
For numerical stability, some compounds in some bins and layers are assumed to be at equilibrium with the gas phase because equilibrium is reached very quickly (for example for very volatile compounds).To identify cases where equilibrium with the gas phase should be assumed, a criterion t equilibrium is used.If the characteristic time to reach equilibrium with the gas phase is lower than t equilibrium , the case is considered at equilibrium, whereas if it is higher than t equilibrium , the case is solved dynamically.
Concentrations of organic compounds in the organic phases are computed according to Eq. ( 81) (obtained by generalizing Eq. 2 for several bins, layers and phases): where conc i,eq is the total concentration of compound i at equilibrium with the gas phase (sum of the gas-phase concentration and of organic-phase concentrations in layers and bins at equilibrium) computed with Eq. ( 82) and ratio i,eq the ratio of the concentration of compound i at equilibrium with the gas phase in the particle phase to the concentration of i in the gas phase computed with Eq. ( 83) (similar to Eq. 15 for several bins, layers and phases).
conc i,eq =A tot,i − bin layer phase (1 − λ bin,layer,phase )A With λ bin,layer,phase defined such as λ bin,layer,phase = 1 if τ bin,layer,phase eq < t equilibrium (case at equilibrium) = 0 if τ bin,layer,phase eq >=t equilibrium (dynamic case). (84 The system is solved iteratively, as now detailed.The composition of the particles is first estimated using Eqs.( 81)-( 84).If the concentrations computed from this estimation are different from those obtained in the previous iteration, a new estimation of concentrations is computed.The algorithm is detailed in Table 1.Step 7 and step 2 prevent the non-convergence due to high variations of concentrations by reducing the variations between two iterations.Some numerical issues can arise from the equilibrium representation especially for low volatility compounds with high value of t equilibrium (for example 100 s).For these compounds, assuming equilibrium with the gas phase will give errors because these compounds will condense almost entirely on the bin with the higher organic mass instead of condensing on each bin according to the kinetics of condensation.To prevent this problem, low volatility cases (with K bin,layer,phase p,i > 10) are assumed dynamic.

Methodology used to compute the evolution of concentrations
The method used to solve the evolution of concentrations is shown in Fig. 5.As mentioned in the previous section, the model uses a hybrid method combining a dynamic representation where concentrations evolve as a function of time (for cases with characteristic times higher than a t equilibrium value specified by the user) with an equilibrium representation (for cases with characteristic times lower than the t equilibrium value).In the equilibrium representation, the distribution of organic compounds between phases and the gas/particle partitioning of cases with characteristic times lower than t equilibrium are computed.
In the dynamic representation, concentrations evolving dynamically (cases with characteristic times higher than t equilibrium ) are computed as a function of the time step with the second-order Rosenbrock scheme ROS2 (Verwer et al., 1999) for time integration.The time step can be rejected (the computation is redone with another time step) if the error between the second-order and the first-order concentrations is too high (higher than an EPSER parameter specified by the user).In that case, the time step decreases until the error is lower than EPSER (until the time step is accepted).If the time step is very low, the ROS2 scheme increases the time step to obtain the optimal time step.
To compute concentrations, the model first initializes required properties (activity coefficients, the number of phases in each layer and the characteristic times).For the first time step, the model calls for the first time the equilibrium representation.Once the equilibrium is reached, the model computes the dynamic evolution with the dynamic representation.Because concentrations changed, the gas/particle partitioning of cases at equilibrium changed, the equilibrium representation is then called again.For the next time steps, the dynamic and equilibrium representations are called once each, until the ending time of the simulation is reached.

Results
In this section, the dynamic implicit representation of SOAP is first compared to an explicit representation to check that the dynamic evolution of surrogates between the gas and particle phases is well represented by the implicit representation.Two test cases representative of European summer periods with high biogenic concentrations and high anthropogenic concentrations are then defined.Simulations are performed to study the impact of ideality, saturation, phase separation and thermodynamic equilibrium on the model results.Comparisons of the implicit and explicit representations for one of the test case are also performed.

Validation of the dynamic implicit representation
To validate the dynamic implicit representation of SOAP, we study the condensation of an hydrophobic organic surrogate of different volatilities (with a partitioning coefficient of 100, 1, or 0.01 m 3 µg −1 ) for different organic-phase diffusion coefficients (10 −19 , 10 −20 , 10 −21 , 10 −22 and 10 −24 m 2 s −1 ) into 5 µg m −3 of non-volatile organic compounds.The organic surrogate that condenses is initially only in the gas phase, with a mass of 5 µg m −3 .The size distribution of particles is divided into 5 bins.The number of particles inside each bin is assumed to be equal to 3.02 × 10 8 , 1.40 × 10 9 , 2.24 × 10 8 , 1.69 × 10 6 ; 12 224.0particles m −3 based on average simulation results over Europe (Couvidat et al., 2012).The initial mass concentrations of the surrogate are distributed between each section with the following fraction: 0.7, 26.6, 61.4,11.2 and 0.1 %.The initial diameter of each bin is computed from the initial mass and number concentrations in the bin.
The results of the implicit representation are compared with an explicit representation where diffusion between layers is treated explicitly.In the explicit representation, the diffusion of organic compounds in the particle is represented explicitly as in KM-GAP (with no equilibrium assumption).The particles are discretized with 100 layers, each of them having a volume representing 1 % of the total volume of the particles.The evolution of concentrations inside each layer i is computed by where F i,i+1 and F i,i−1 are the fluxes of diffusion between the layer i and the nearby layers (the computation of those fluxes is given in the Supplement) and V layer the volume of the layer.As in the implicit representation, the volume of each layer is kept constant throughout the simulation using the redistribution scheme presented in Sect.2.3.8. Figure 6 shows the comparison between the explicit representation with 100 layers and the implicit representation with four and five layers (one of these layer being the interface layer not affected by diffusion).The implicit representation reproduces the result of the explicit representation well.For non-volatile compounds (with a partitioning coefficient around 100 m 3 µg −1 ), there is no significant differences between the two methods.For such compounds, the condensation is not limited by diffusion and therefore they condense independently of the organic-phase diffusion coefficient.For low volatility compounds (with a partitioning coefficient around 1 m 3 µg −1 ), first, condensation into the first layer occurs without a limitation due to diffusion and then diffusion over all the particle occurs.The implicit representation with 5 layers slightly underestimates the diffusion of low volatility compounds for low organic-phase diffusion coefficients (10 −19 m 2 s −1 ) but reproduces very well cases with lower diffusion coefficients (lower than 10 −20 m 2 s −1 ), whereas the implicit representation with 4 layers slightly overestimates the diffusion of low volatility compounds.For more volatile compounds (with a partitioning coefficient around 0.01 m 3 µg −1 ), the implicit representation reproduces the explicit representation well with both four and five layers but underestimates low concentrations (lower than 10 % of the concentration at equilibrium) with four layers.
To fully validate the implicit representation, a comparison of the implicit and explicit representations is also performed for one of the test cases (see Sect. 3.5).

Setup of the test cases
The behavior of the model is tested using two test cases.The first test case corresponds to a summer period with high concentrations of biogenic SOA, and the second test case corresponds to a summer period with high concentrations of anthropogenic SOA.These two test cases hereafter are referred as the biogenic test case and the anthropogenic test case.For the biogenic and the anthropogenic test cases, pollutant mass The surrogates are the same as in the aerosol model H 2 O, which was used by Couvidat et al. (2012Couvidat et al. ( , 2013)).Surrogates are described in Table 4 and their total concentrations (gas + particle) are given in Table 5 for both test cases.For both test cases, the liquid water content of aerosol, the pH, the ionic strength and the concentrations of inorganic ions are computed with the ISORROPIA model (Nenes et al., 1998) at a specified RH.In H 2 O, the surrogates that are representative of primary and aged SVOCs (POAlP, POAmP, POAhP, SOAlP, SOAmP and SOAhP) do not have a molecular structure attached.Therefore, processes depending on the molecular structure (like absorption on the aqueous phase) are not estimated for those surrogates.They are assumed to be hydrophobic and their influence on the activity coefficients of other species is taken into account by assigning to them a default molecular structure representative of primary aerosol and lowly oxidized compounds.This default structure consists of 40 % of C 23 H 47 COOH, 5 % of C 8 H 17 CH = CHC 7 H 14 COOH, 15 % of 4-(2-propio)syringone, 12 % of C 29 H 60 and 28 % of 2-carboxybenzoic acid based on EPRI (1999).In our test cases, the same default structure as in H 2 O is used.
The size distribution and the number of particles are the same as those used for the validation.The volume of the solid core of each size section is computed from the mass concentrations of solid species (dust, black carbon and solid inorganic are given by ISORROPIA) to compute morphology factors.

Equilibrium representation
To study the impact of activity coefficients, hydrophilic properties of surrogates, saturation and phase separation, concentrations in the particle phase are assumed to be at equilibrium with the gas phase for all surrogates.

Influence of activity coefficients on organic aerosol formation (without phase separation of the organic phase)
The effect of non-ideality on aerosol concentrations is strong and complex.To determine the impact of non-ideality, Tables 6 and 7 show the concentrations of organic aerosol for both test cases formed from the various precursors as well as the concentrations of water, with and without the ideality assumption, respectively.The compounds are assumed to be both hydrophilic and hydrophobic except for the species POAlP, POAmP, POAhP, SOAlP, SOAmP and SOAhP, which are simply assumed to be hydrophobic, with the default molecular structure used in Couvidat et al. (2012).
Concentrations of compounds in the organic phase tend to decrease strongly when non-ideality is assumed (except for aromatics in the biogenic case due to non-linear effects), especially for the compounds formed from isoprene oxidation (most of them are very hydrophilic and therefore have low affinity with very hydrophobic compounds) and for some of the compounds formed from monoterpenes.Concentrations of hydrophilic compounds in the aqueous phase either increase or decrease depending on conditions.For the compounds formed from isoprene oxidation, their concentrations increase in the biogenic case from 0.70 µg m −3 for ideality to 1.24 µg m −3 for non-ideality with long-, medium-and short-ranges interactions and to 1.29 µg m −3 for non-ideality with only short-range interactions.It seems to indicate that in this case short-range interactions between organic compounds tend to stabilize hydrophilic organic compounds in the aqueous phase, whereas medium-range and long-range interactions between organic and inorganic compounds tend to destabilize hydrophilic organic compounds.These results are consistent with the result of Couvidat et al. (2012) where assuming ideality at infinite dilution can lead to an underestimation of the concentrations of hydrophilic species inside the aqueous phase of particles.Therefore, the concentrations of hydrophilic organic compounds in the aqueous phase probably depend strongly on other compounds and on the nonideality of aerosols.

Hydrophobic versus hydrophilic
Determining a priori if a compound is hydrophilic or hydrophobic is not straightforward.Table 8 shows that some compounds seem clearly hydrophobic (BiA2D, BiA1D, An-BlP, BiBlP, BiBmP, AnClP, BiNGA, BiNIT3, BiNIT) or hydrophilic (BiA0D) as they partition only into one phase.However, some compounds are present in both the organic and aqueous phases.AnBmP seems to be both hydrophilic and hydrophobic and can change phase depending on conditions.Moreover, BiA2D, BiA1D and BiMGA were assumed to be hydrophilic in H 2 O whereas it seems from these test cases that they are mainly hydrophobic.The fact that these compounds were assumed hydrophilic is probably due the choice of a criterion not representative of all atmospheric conditions.BiA2D and BiA1D were assumed hydrophilic based on their octanol/water coefficient (Pun et al., 2006), which is probably not representative of atmospheric conditions.BiMGA was assumed hydrophilic based on the results of Couvidat and Seigneur (2011).According to this study, BiMGA condenses mostly on the aqueous phase if the liquid water content of aerosols is high (superior to the concentration of the organic phase), which is not the case in the two test cases.Moreover, if medium-range and long-range interactions are not taken into account, the distribution of the compounds between phases changes significantly.For the biogenic test case, if medium-range and long-range interactions are not taken into account as in Couvidat and Seigneur (2011), the fraction in the organic phase decreases from 75 to 33 % for particulate BiA2D and from 89 to 46 % for particulate BiA1D.For the anthropogenic test case, the fraction in the organic decreases, from 99 to 76 % for particulate BiA2D, from 100 to 85 % for particulate BiA1D and from 59 to 46 % for particulate BiMGA.It is therefore possible for these compounds to be present in both phases depending on conditions.

Saturation and phase separation
Species having very different properties do not mix well together and phase separation can be computed by Gibbs energy minimization (see Sect. 2.2.4).Table 9 presents the concentrations with or without phase separation for the bio- genic test case at RH = 30 % without an aqueous phase.
In this case, without phase separation, both hydrophilic and hydrophobic compounds want to condense into the organic phase which is mainly constituted by hydrophobic compounds.Assuming phase separation in this case does not strongly influence concentrations of hydrophobic compounds, which decrease slightly.However, a second organic phase is created by phase separation which is composed of mainly very oxidized compounds (BiPER, BiDER and BiMGA).For the anthropogenic case, at RH = 30 %, phase separation does not occur because the concentrations of hydrophilic compounds are too low for the organic phase to be saturated.
Another organic phase may be created if there are compounds with a low oxidation state.For the biogenic case at low humidity, if the structure of nonadecane is used to represent the compounds POAlP, POAmP, POAhP, SOAlP, SOAmP and SOAhP (without a molecular structure; these compounds only condense into the less oxidized phase without impacting phase separation), a third organic phase may be created.Compounds with a low oxidation state may not readily mix with slightly oxidized compounds which in turn may not readily mix with more oxidized compounds.

Dynamic representation
For dynamic modeling, users must choose the organicphase diffusion coefficients because there is currently, to our knowledge, no method to estimate diffusion coefficients of the organic phase as a function of the composition of organic  aerosols.In this section, the same organic-phase diffusion coefficient is used for all organic compounds.For all the compounds, the gas-phase diffusion coefficient is assumed to be equal to 0.1 cm 2 s −1 which is the order of magnitude of this parameter (Seinfeld and Pandis, 1998) and the accommodation coefficient (value between 0 and 1) is assumed to be equal to 0.01 so that it corresponds to the order of magnitude of accommodation coefficients for organic compounds given in other studies, i.e., between 0.1 (Saleh et al., 2013) and 0.001 (Lee et al., 2011).The gas-phase diffusion coeffi-cient and the accommodation coefficient are used in Eqs. ( 79) and ( 80).In the regional-scale simulations from which initial concentrations are extracted, equilibrium is assumed between gas and particles.For the dynamic modeling, 5 µg m −3 of non-volatile hydrophobic organic compounds is added to have a mass into which compounds are not present and will have to diffuse.For each surrogate, the total mass gas + particles is kept the same as in the equilibrium simulations.However, the particle-phase concentrations are assumed to be equal to 10 % of the equilibrium particle-phase concentrations.
To fully validate the implicit representation of the dynamic evolution of surrogates, the results of SOAP are compared to the results of the explicit representation for the biogenic test case assuming ideality.RH is assumed to be equal to 30 %. Results of this comparison are shown in Fig. 7.This results shows good results with five layers and a slight overestimation with four layers of OM concentrations due to the overestimation of condensation for low volatility compounds.
The dynamic implicit representation is tested for various humidities and various diffusion coefficients for the biogenic and anthropogenic cases.Figures 8 and 9 show the concentrations in the organic phases for the biogenic and anthropogenic cases, respectively, whereas Figs. 10 and 11 show concentrations in the aqueous phase for the biogenic and anthropogenic cases, respectively.The influence of the number of layers is also tested in these figures.The two test cases show similar results.For the condensation into the organic phases, the growth of particles can be divided into two parts.The particle grows by condensation into the layer at the interface independently from the viscosity of the organic phases until the growth of particle become limited by the diffusion of the more volatile compounds into the organic phases.For the anthropogenic test case, humidity has a low effect between 30 and 70 % on the condensation of organic compounds into the organic phases because they are mainly constituted of hydrophobic compounds in the anthropogenic test case.However, at 70 % more hydrophilic compounds manage to condense into the aqueous phase.Between 70 to 90 %, the amount of organic compounds in the organic phases decreases whereas the amount of organic compounds increases in the aqueous phase.For the biogenic test case, the organic-phase concentration increases between 30 and 50 % of humidity because hydrophilic compounds (which are in higher concentrations than in the anthropogenic test case) partition more efficiently.Between 30 and 70 %, the organicphase concentration decreases but the aqueous-phase concentration increases because hydrophilic compounds condense more efficiently into the aqueous phase.However, diffusion coefficients are probably much higher at high humidity than at low humidity because the organic phases would be less viscous (less oligomerization due to esterification for example and more water in the organic phases which would decrease viscosity) as shown by Saukko et al. (2012) and Renbaum-Wolff et al. (2013).It may then be possible that at high humidity, organic aerosols reach equilibrium quicker than at low humidity.A method to estimate diffusion coefficients as a function of composition (or at least of as a function of humidity) is needed to properly represent this phenomenon.

Time analysis
Table 10 shows the computation time necessary to solve the biogenic case at RH = 70 % (relative to the time necessary to solve the case with an uncoupled and ideal system) for different configurations of the model.For the dynamic approach, 4 layers are used, the length of the computation corresponds to 600 s (corresponding to a time step of a simulation with the Polyphemus air quality platform over Europe; Sartelet et al., 2007), initial concentrations are assumed to be 80 % of concentrations at equilibrium and an organic-phase diffusion coefficient of 10 −21 m 2 s −1 is used.The times given here correspond to a specific case and therefore can greatly vary with initial conditions and with the chosen parameters for the numerical resolution of the system.They are provided here on an indicative basis.The dynamic approach is around 400 times slower than the equilibrium approach, making its applicability limited to short-term 3-D simulations.Table 10.Time necessary to solve the biogenic case at RH = 70 % (relative to the time necessary to solve the case with an uncoupled system with ideality) for different configurations of the model.

Phenomena taken into account
Time for Coupled Short-range Long-and medium-Organic phase Equilibrium Dynamic system interactions range interactions separation approach More reliable and complete information about the computation time cost will be provided in future studies about the implementation of the SOAP in an air quality model.

Conclusions
The Secondary Organic Aerosol Processor (SOAP) model is a modular model, which can compute the condensationevaporation of organic aerosol according to two different approaches: an equilibrium approach and a dynamic ap-  proach.In the equilibrium approach, concentrations in the particle phase are assumed at equilibrium with concentrations in the gas phase.In the dynamic approach, concentrations evolve according to the kinetics of condensationevaporation-diffusion of organic compounds.The dynamic approach uses a multi-layer representation of particles to represent the particle-phase diffusion of organic compounds.Future works will focus on improving the framework of this dynamic approach to take into account varying diffusion coefficients with layers, to allow for the entrapment of compounds in inner layers and to represent layer exchanges and transfers between the organic and the aqueous phases.Simulations with SOAP and comparisons to measurements should be performed to validate the model and to test the influence of each process and parameter on organic aerosol formation.
To improve the representation of aerosols, several processes should be added to the model.First, interactions between organic and inorganic compounds should be fully taken into account via activity coefficients.Currently, only the influence of inorganic compounds on organic compounds is taken into account.However, organic compounds can also impact the formation of inorganic compounds due to those interactions.This process could be taken into account by adding inorganic aerosol formation to SOAP.Second, a method to estimate diffusion coefficients in organic phases should be developed, as it is expected that composition of the organic phases greatly influences the viscosity and therefore diffusion coefficients of organic compounds.Finally, the model could be coupled to a solver for particle-phase chemistry and then represent processes such as oligomerization, which could affect the viscosity of the organic phases.

Appendix A: Notations
A eq the concentration in the organic phase at equilibrium A g,i the concentration of i in the gas phase (in µg m −3 ) A eq,bin g,i the concentration of iin the gas phase (in µg m −3 ) at equilibrium with the whole particle A eq,bin,interface g,i the concentration of i in the gas phase (in µg m −3 ) at equilibrium with the interface layer A eq,bin,layer g,i the concentration of i in the gas phase (in µg m −3 ) at equilibrium with the layer number of particles in a diameter bin (in m −3 ) P 0 i the saturation vapor pressure of i P i partial pressure of i P eq,i partial pressure of i at equilibrium (taking into account the Kelvin effect) R the ideal gas constant R p the rayon of the particle (in m) S aq surface of particles covered by an aqueous phase S tot surface of particles RH the relative humidity T the temperature T ref the temperature of reference at which P 0 i is determined V layer volume fraction of the layer V aq volume of the aqueous phase V org volume of the organic phase X i,aq molar fraction of i in the aqueous phase X i,org molar fraction of i in the organic phase α accommodation coefficient α layer ratio of the characteristic time for diffusion of the layer to the characteristic time for diffusion of the center of the particle γ i,aq activity coefficient of i in the aqueous phase γ ∞ i,aq activity coefficient of i at infinite dilution in water γ i,org activity coefficient of i in the organic phase γ water,aq activity coefficient of water in the aqueous phase γ water,org activity coefficient of water in the organic phase H i enthalpy of vaporization of i (in J mol −1 ) µ φ i chemical potential of i ξ i activity coefficient of i in the aqueous phase by reference to infinite dilution τ dif characteristic time for diffusion τ eq characteristic time to reach equilibrium The Supplement related to this article is available online at doi:10.5194/gmd-8-1111-2015-supplement.

Figure 1 .
Figure 1.Evolution of the ratio A p (t)/A eq as a function of the ratio t/τ dif .

Figure 2 .
Figure 2. Morphology factors as a function of the volume fraction of the solid phase (layer 1 corresponds to the layer near the solid core of the particle and layer 3 to the layer near the interface) for three diffusion layers (left) and four diffusion layers (right).
the system can not converge (the system is in a situation where relative errors on concentrations between two steps does not change or return to old values), n = n + 1. 8. If the system has not converged (relative errors on concentrations between two steps are too high), return to step 2.

Figure 3 .
Figure 3. Diagram of the first step of the algorithm to compute evaporation-condensation-diffusion fluxes.

Figure 4 .
Figure 4. Diagram of the second step of the algorithm to compute evaporation-condensation-diffusion fluxes.
that depends on the diffusion flux sign in each layer:If the flux of diffusion of compound i into the particle is positive (J bin diff,i > 0) when for the interface layer J of diffusion of compound i into the particle is negative (J bin diff,i < 0) when for the interface layer J bin,layer i = J bin,interface i if J bin,interface i < 0 else J bin,layer i If the new mass of organics M bin,new o = layer phase M bin,layer,phase,new o is lower than the mass at previous iteration or time step (before adding new compounds) M bin,old o , the redistribution is done from the center of the particle to the interface layer: from ilayer = 1 to N layer − 1, ilayer2 = ilayer+1: of the layer decreased: 3.1.1Mass to be redistributed: M = α layer × M bin,new o − M bin,layer,new o 3.1.2Fraction of concentrations to redistribute from the next layer (ilayer2): f = M M bin,layer2,new o 3.1.3Themissing mass of the layer is taken from the next layer (ilayer 2):

Figure 5 .
Figure 5. Diagram of the method used to compute the evolution of concentrations.

AF
p,i the concentration of i in the organic phase (in µg m −3 ) i in the organic phase (in µg m −3 ) in a phase A bin p,i the concentration of i in the organic phase (in µg m −3 ) in a diameter bin A bin,layer p,i the concentration of i in the organic phase (in µg m −3 ) in a diameter bin/layer A bin,layer,phase p,i the concentration of i in the organic phase (in µg m −3 ) in a bin/layer/phase A bin,interface p,i the concentration of i in the organic phase (in µg m −3 ) in a diameter bin at the interface layer A bin,interface,phase p,i the concentration of i in the organic phase (in µg m −3 ) in a bin/phase at the interface layer A p,water the concentration of water in the organic phase (in µg m −3 ) A aq,i the concentration of i in the aqueous phase (in µg m −3 ) A bin aq,i the concentration of i in the aqueous phase in a diameter bin (in µg m −3 ) A tot,i the concentration of i in all phases (in µg m −3 ) AQ total mass of the aqueous phase including organic compounds (in µg m −3 ) AQ inorg total mass of the aqueous phase excluding organic compounds (in µg m −3 ) AQ bin total mass of the aqueous phase including organic compounds in a bin (in µg m −3 ) C concentration (in M) C s concentration at the surface of the particle (in M) D p diameter of the particle (in m) D air diffusivity of the compound in air (in m 2 s −1 ) D org diffusivity of the compound in the organic phase (in m 2 s −1 ) a compound to encounter an aqueous phase at the surface of the particle f (Kn, α) transition regime formula f s volume of the solid phase in the particle over the volume of the particle G Gibbs energy H i Henry's law constant (in M atm −1 ) J bin cond/evap,i Flux of condensation-evaporation (in µg m −3 s −1 ) inside the whole particleJ bin diff,iFlux of diffusion (in µg m −3 s −1 ) inside the whole particle J bin,layer diff,i Flux of diffusion (in µg m −3 s −1 ) inside a bin/layer J bin,layer,phase diff,i Flux of diffusion (in µg m −3 s −1 ) inside a bin/layer/phase J bin,interface i Flux of condensation-evaporation (in µg m −3 s −1 ) at the interface layer J bin tot,i Flux of condensation-evaporation-diffusion (in µg m −3 s −1 ) inside the whole particle J bin,layer tot,i Flux of condensation-evaporation-diffusion (in µg m −3 s −1 ) inside a bin/layer J bin,layer,phase tot,i Flux of condensation-evaporation-diffusion (in µg m −3 s −1 ) inside a bin/layer/phase the aqueous-phase partitioning coefficient of a diameter bin K p,i the organic-phase partitioning coefficient K phase p,i the organic-phase partitioning coefficient for a phase www.geosci-model-dev.net/8phasepartitioning coefficient of a diameter bin/layer/phase K bin,interface p,i the organic-phase partitioning coefficient of a diameter bin at the interface layer K bin,interface,phase p,i the organic-phase partitioning coefficient of a diameter bin/phase at the interface layer k kinetic rate parameter of the absorption-diffusion equation k absorption kinetic rate parameter of the absorption equation k diffusion kinetic rate parameter of the diffusion equation K into which diffusion occurs M aq mean molar mass of the aqueous phase (in g mol −1 ) M i mean molar mass of i (in g mol −1 ) m i mass of i in a particle m bin i mass of i in a particle of the bin M o concentration of the organic phase (in µg m −3 ) M phase o concentration of the organic phase (in µg m −3 ) for a phase M bin o concentration of the organic phase in a diameter bin (in µg m −3 ) M bin,layer o concentration of the organic phase in a diameter bin/layer (in µg m −3 ) M bin,layer,phase o concentration of the organic phase in a diameter bin/layer/phase (in µg m −3 ) M bin,interface o concentration of the organic phase at the interface layer in a diameter bin (in µg m −3 ) M bin,interface,phase o concentration of the organic phase at the interface layer in a diameter bin/phase (in µg m −3 ) M ow mean molar mass of the organic phase (in g mol −1 ) M water molar mass of water n φ i number of moles of compound i in the phase φ N bin

Table 3 .
Algorithm to compute the redistribution of compounds between layers.

Table 5 .
Conditions of the test cases.

Table 6 .
Concentrations of organic aerosols formed for each precursor for both test cases at RH = 70 % if ideality is assumed.Biogenic test caseAnthropogenic test cases A p (µg m −3 ) A aq (µg m −3 ) A p (µg m −3 ) A aq (µg m −3 )

Table 7 .
Concentrations of organic aerosols formed for each precursor for both test cases at RH = 70 % if non-ideality is assumed (with short-range, medium-range and long-range interactions).−3 ) A aq (µg m −3 ) A p (µg m −3 ) A aq (µg m −3 )

Table 8 .
Distributions of surrogate organic compounds between the aqueous phase and the organic phase.

Table 9 .
Concentrations (in µg m −3 ) with or without phase separation for the biogenic test case at RH = 30 %.