A thermodynamic framework for bulk–surface partitioning in finite-volume mixed organic–inorganic aerosol particles and cloud droplets

Atmospheric aerosol particles and their interactions with clouds are among the largest sources of uncertainty in global climate modeling. Aerosol particles in the ultrafine size range with diameters less than 100 nm have very high surface-area-to-volume ratios, with a substantial fraction of molecules occupying the air–droplet interface. The partitioning of surface-active species between the interior bulk of a droplet and the interface with the surrounding air plays a large role in the physicochemical properties of a particle and in the activation of ultrafine particles, especially those of less than 50 nm diameter, into cloud droplets. In this work, a novel and thermodynamically rigorous treatment of bulk–surface equilibrium partitioning is developed through the use of a framework based on the Aerosol Inorganic–Organic Mixtures Functional groups Activity Coefficients (AIOMFAC) model in combination with a finite-depth Guggenheim interface region on spherical, finite-volume droplets. We outline our numerical implementation of the resulting modified Butler equation, including accounting for challenging extreme cases when certain compounds have very limited solubility in either the surface or the bulk phase. This model, which uses a single, physically constrained interface thickness parameter, is capable of predicting the size-dependent surface tension of complex multicomponent solutions containing organic and inorganic species. We explore the impacts of coupled surface tension changes and changes in bulk–surface partitioning coefficients for aerosol particles ranging in diameters from several micrometers to as small as 10 nm and across atmospherically relevant relative humidity ranges. The treatment of bulk–surface equilibrium leads to deviations from classical cloud droplet activation behavior as modeled by simplified treatments of the Köhler equation that do not account for bulk–surface partitioning. The treatments for bulk–surface partitioning laid out in this work, when applied to the Köhler equation, are in agreement with measured critical supersaturations of a range of different systems. However, we also find that challenges remain in accurately modeling the growth behavior of certain systems containing small dicarboxylic acids, especially in a predictive manner. Furthermore, it was determined that the thickness of the interfacial phase is a sensitive parameter in this treatment; however, constraining it to a meaningful range allows for predictive modeling of aerosol particle activation into cloud droplets, including cases with consideration of co-condensation of semivolatile organics.

1 Introduction

Atmospheric aerosols are suspensions of condensed particles and the gas which surrounds them. They can be a major contributor to poor indoor and outdoor air quality. Exposure to poor air quality is the leading environmental risk factor for premature mortality globally (Cohen et al., 2017; Burnett et al., 2018). Beyond their public health effects, the role of atmospheric aerosols and aerosol–cloud interactions in the global climate system remains one of the least understood processes in climate models (Boucher et al., 2013). Some light-absorbing aerosols exhibit a positive direct radiative forcing, while others mostly scatter solar radiation and exhibit negative direct radiative forcing. The magnitude and sign of the local or regional direct radiative forcing depend on numerous factors, including an aerosol's composition, size distribution, geographic location, season, and altitude in the atmosphere (Bellouin et al., 2020). Furthermore, atmospheric aerosols interact with clouds and precipitation, thereby contributing to important indirect effects on Earth's climate that are poorly constrained (Seinfeld et al., 2016). Therefore, it is of the utmost importance to understand the physicochemical properties and microphysical processes of aerosol particles, including aerosol–cloud–radiation interactions, as these drive their impacts on the climate system.

Aerosol particles exist over a broad range of sizes. Differently sized particles frequently originate from distinct sources or form through unique processes. Particles in the nucleation mode are primarily formed via the spontaneous condensation of gaseous compounds to create or extend a tiny condensed phase. Similarly, Aitken mode particles may form by the continued growth of nucleation mode particles through condensation of gaseous species or coagulation of condensed particles. Particles in the nucleation mode or Aitken mode are collectively referred to as ultrafine particles. Accumulation mode particles form through coagulation and agglomeration of smaller Aitken mode particles in the air. The largest aerosols, such as mineral dust and sea spray particles, are primarily composed of inorganic compounds and typically form via mechanical processes.

As aerosol particles decrease in size, they are more likely to originate from and grow through complex multiphase chemical reactions. Many fine and ultrafine aerosol particles are composed of organic compounds that form via oxidation reactions of volatile organic compounds (VOCs) or intermediate-volatility organic compounds (IVOCs) in the gas phase. Many of the products of these reactions are lower in volatility and partition into condensed particles. Semivolatile and low-volatility organic compounds (SVOCs and LVOCs) may also be emitted and partition into particles with minimal processing in the gas phase. Organic aerosol particles, or components thereof, formed through gas-phase or multiphase chemistry are known as secondary organic aerosols (SOAs) (Jimenez et al., 2009; Hallquist et al., 2009); for a complete list of abbreviations and symbols, please refer to Tables S1 and S2 of the Supplement. Within a condensed phase additional chemical reactions may proceed. This is particularly notable in low-viscosity aqueous aerosol particles and cloud droplets; those droplets provide aqueous media wherein numerous laboratory and field observations have observed the formation or transformation of SOAs through chemical processing (e.g., Ervens et al., 2011). For certain species these reactions can further reduce their volatility, effectively “trapping” said species in the particle phase. For example, isoprene epoxydiols (IEPOXs) are known to undergo acid-catalyzed reactions in aqueous ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$-containing particle phases (Surratt et al., 2010). The IVOCs and VOCs which lead to the formation of SOAs can be biogenic (e.g., isoprene or monoterpenes) or anthropogenic (e.g., toluene, xylene, naphthalene) (Ng et al., 2007; Chan et al., 2009). Observational studies have concluded that the majority of SOA material is derived from biogenic precursors (Zhang et al., 2007).

The internal mixing state and geometries of aerosol particles are varied. Aerosol particles that are emitted by primary emission from biomass burning often take the shape of long agglomerations of individual spherical particles (Pokhrel et al., 2021), while aqueous inorganic-rich particles and soluble SOAs are thought to be approximately spherical. The frequent presence of liquid–liquid-phase separation (LLPS) can complicate the shape aspect. LLPS particles may form radially symmetric structures with an organic-rich “shell” (phase β) covering an inorganic-rich “core” (phase α) or other geometrically more complex structures, such as a partially engulfed morphology, wherein the particle is no longer radially symmetric nor perfectly spherical. In a partially engulfed particle, phase β does not spread completely over phase α and leaves a portion of phase α exposed (Song et al., 2013, 2012; Ciobanu et al., 2009; Shiraiwa et al., 2013). These particles are also called “Russian dolls” or “janus particles” and interact with radiation differently than their core–shell counterparts (Lang-Yona et al., 2010). Furthermore, should phase β of core–shell-phase-separated particles become a highly viscous liquid or even glassy, it can limit the reactive uptake of different species (Schmedding et al., 2020, 2019; Zhou et al., 2019; Kuwata and Martin, 2012).

Since the surface-area-to-volume ratio of a sphere scales with the inverse of its radius, said ratio becomes of great importance for ultrafine particles. If the assumption is made that the surface of a particle is represented by a monolayer of molecules in contact with the gas phase, an assumption further explored in this study, then spherical particles with a diameter less than 100 nm will have a non-negligible quantity of molecules present at their surface. For example, in a pure water droplet with a diameter of 50 nm, ∼4 % of the water molecules will be present in the surface monolayer, while more than half of the water molecules in a pure droplet with a diameter of 3 nm will be present in the surface monolayer.

Many ultrafine aerosol particles are complex multicomponent systems and contain numerous surfactant and tensoionic species in differing concentrations. These species modify the surface tension of a particle (Sorjamaa et al., 2004; Nozière et al., 2014; Gérard et al., 2016; Petters and Petters, 2016; Ruehl et al., 2016; Ovadnevaite et al., 2017; Kroflič et al., 2018; Malila and Prisle, 2018; Gérard et al., 2019). In finite-volume systems, surfactants may become depleted out of the particle interior; likewise tensoionic species may become increasingly concentrated in the particle interior. This may influence the equilibrium morphology of an aerosol particle for a given size. Even in macroscopic systems, it is possible for LLPS particles to adopt non-spherical partially engulfed morphologies based on their surface properties (Binyaminov et al., 2021). The surface enrichment and depletion of species may affect the conditions under which an aerosol particle will activate and quickly grow into a cloud droplet.

The impact of surface tension modification by surfactants on cloud droplet activation has been known and studied for decades (Facchini et al., 1999, 2000; Topping et al., 2007; Ovadnevaite et al., 2017). The conditions under which a hygroscopically growing aerosol particle will activate into a cloud droplet were first stated by Köhler (1936), recognizing the importance of the global maximum of a particle's equilibrium saturation ratio, S_crit, as expressed by the following equation:

$$\begin{array}{}\text{(1)}& S=a_{\mathrm{w}}\exp \left({\displaystyle \frac{\mathrm{4}\mathit{\sigma }{M}_{\mathrm{w}}}{RT{\mathit{\rho }}_{\mathrm{w}}{D}_{\mathrm{p}}}}\right).\end{array}$$

Here, σ denotes the effective surface tension of the particle at the air–liquid interface. M_w and ρ_w are the molar mass and density of water, respectively. R is the gas constant, T is the temperature, and D_p is the diameter of the particle. The so-called solute or Raoult effect is captured by the water activity, a_w, whereas the exponential factor captures the Kelvin effect. These two effects are often considered to be in competition with one another, since the Raoult effect decreases the supersaturation necessary for an aerosol (or cloud condensation nucleus, CCN) to reach activation into a cloud droplet, while the Kelvin effect increases the critical supersaturation. Note that values of S are often reported as percentages, in which case the value of Eq. (1) is multiplied by 100 %. A common reason for considering this a competition is the following: as surface-active species are depleted from a droplet's interior “bulk” to populate the surface during hygroscopic droplet growth, their lowered concentration in the bulk will typically raise the value of a_w while simultaneously decreasing the value of σ by lowering the effective surface tension. However, we note that there is not necessarily a competition. In particular, in multicomponent organic–inorganic particles with substantial nonideal mixing, such as in cases with LLPS, the bulk–surface partitioning of low-polarity surfactants can lead to a relative lowering of both surface tension and a_w (e.g., Ovadnevaite et al., 2017).

In an effort to simplify calculations of aerosol water uptake and CCN activation, Petters and Kreidenweis (2007) developed a single-parameter model of hygroscopic growth, commonly known as κ-Köhler theory, wherein a_w can be related to the current volume of water in the particle V_w and the starting dry particle (solute) volume, V_dry, through the following equation:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{1}}{a_{\mathrm{w}}}}=\mathrm{1}+\mathit{\kappa }{\displaystyle \frac{V_{\mathrm{dry}}}{V_{\mathrm{w}}}}.\end{array}$$

The value of κ for a multicomponent solution can be found from a volume-fraction-based linear weighting of the κ values of individual components. The combination of Eqs. (1) and (2), along with an assumption about σ, allows for the calculation of critical supersaturations of aerosol particles given only their composition and dry sizes (Petters and Kreidenweis, 2007); however, the utility of κ-Köhler theory is thought to be more limited for ultrafine particles (Topping et al., 2016).

While κ-Köhler theory may provide reasonable predictions for many systems, including particles in the accumulation mode or of even larger sizes that are composed of highly soluble components, its neglect of bulk–surface partitioning treatments in combination with the high surface-area-to-volume ratios of ultrafine particles is worth exploring in greater detail. There are numerous methods for determining the surface composition at a gas–liquid interface, all of which rely on various assumptions, such as the location and dimensionality of said interface and the inclusion of various system-dependent fit parameters. One of the classical methods is the Gibbs 2-dimensional (2D) dividing surface (or plane) approach, which relates the change in surface tension of a solution (dσ) to the surface excess concentration (Γ_i) and change in chemical potential dμ_i of component i in the bulk solution via (Gibbs, 1874)

$$\begin{array}{}\text{(3)}& \mathrm{d}\mathit{\sigma }=-\sum _i{\mathrm{\Gamma }}_i\mathrm{d}{\mathit{\mu }}_i.\end{array}$$

The term surface excess concentration may be misleading in this case, since it is defined as the difference between the concentration of i on the Gibbs dividing plane and that in the interior volume of a bulk phase adjacent to the plane. Thus negative values of Γ_i are possible if i partitions preferentially into the bulk of a solution and consequently has a lower concentration at the surface. The definition of the location of the Gibbs surface is of the utmost importance when utilizing this approach. Its location is typically selected such that the surface excess concentration of the first species (usually the main solvent, e.g., water) is precisely 0, thus making that 2D plane location a system-specific as well as composition- and size-dependent parameter. Determining the location of the Gibbs dividing plane presents additional challenges, since the location of the surface is defined to be located within the interfacial region between two phases but may not necessarily be found at the same radial position as that of the “outside” edge of the monolayer of molecules found at the boundary of a phase.

From the theoretical framework laid out by Gibbs, the semi-empirical Szyszkowski–Langmuir isotherm for bulk–surface partitioning was developed (Szyszkowski, 1908):

$$\begin{array}{}\text{(4)}& \mathit{\sigma }={\mathit{\sigma }}_{\mathrm{w}}\left(T\right)-A^{\mathrm{SL}}\ln \left(\mathrm{1}+{\displaystyle \frac{B^{\mathrm{SL}}}{C_i^{\mathrm{SL}}}}\right).\end{array}$$

Here, σ_w is the surface tension of pure water at the temperature of interest, A^SL and B^SL are system-dependent fit parameters, and $C_i^{\mathrm{SL}}$ is the concentration of solute i in the bulk liquid of the system. While simpler as an approach and adequate in many system-specific cases (such as the interpretation of laboratory studies), the predictive power of this equation is limited due to the inclusion of two fit parameters for each chemical system.

A third treatment for modeling bulk–surface partitioning of different species was introduced by Jura and Harkins (1946). It combined semi-empirical fitted models, relating surface concentrations to surface tension as follows:

$$\begin{array}{}\text{(5)}& \mathit{\sigma }={\mathit{\sigma }}_{\mathrm{w}}\left(T\right)-(A_{\mathrm{0}}^{\mathrm{CF}}-A_i^{\mathrm{CF}})m_{\mathit{\sigma }}.\end{array}$$

This equation describes an insoluble compressed film (CF) at the interface between a liquid and a gas with a 2D equation of state for bulk–surface partitioning. $A_{\mathrm{0}}^{\mathrm{CF}}$ is the maximum surface adsorption possible, $A_i^{\mathrm{CF}}$ is the current surface adsorption of i, and m_σ is a term relating the change in surface tension to the change in surface concentration. This equation, when coupled with an isotherm relating bulk and surface concentrations, is capable of describing the surface tension as aerosol particles grow hygroscopically when exposed to increasing relative humidity (RH) (Ruehl et al., 2016). In macroscopic systems, the approaches of Szyszkowski (1908) and Jura and Harkins (1946) assume that the enrichment or depletion of species at the particle surface has a negligible effect on the bulk particle composition; however, in ultrafine particles, depletion of surface-active species from the bulk phase may become important. This bulk depletion effect has recently come under scrutiny (Prisle et al., 2010; Bzdek et al., 2020; Lin et al., 2020) and has already been considered previously by Sorjamaa et al. (2004), who considered binary systems of water and the surfactant sodium dodecyl sulfate (SDS). They noted that Köhler curve calculations should include the effect of bulk depletion, since neglecting this effect led to unrealistically low critical supersaturation conditions. These results were supported by laboratory measurements of ternary water–SDS–sodium-chloride systems taken by Prisle et al. (2008).

While many organic species present in the atmosphere are expected to be strongly surface active, there are others classified as weak surfactants, including organosulfates (Hansen et al., 2015), certain components in mixtures of marine SOA and POA (Ovadnevaite et al., 2017), and aliphatic dicarboxylic acids (Ruehl et al., 2016). The surface enrichment of these compounds is difficult to predict, currently inaccessible experimentally (for airborne particles), and potentially showing complex interactions with particle size (Sorjamaa et al., 2004; Sorjamaa and Laaksonen, 2007; Davies et al., 2019; Ovadnevaite et al., 2017).

Prior treatments of bulk–surface partitioning in aerosol systems

There have been numerous methods developed for predicting the bulk–surface partitioning of aerosol chemical species; the following is a summary of recent works. Briefly, Sorjamaa et al. (2004) used an approach based on the 2D Gibbs dividing surface theory and found that the depletion of water-soluble surfactants may have a substantial impact on Köhler curves through both surface tension depression and modification of the Raoult effect as surfactants are depleted from the bulk. They report that failing to account for both surface and bulk effects in growing aerosol particles may lead to under-predictions of the critical supersaturation, particularly at higher organic mass fractions. Prisle et al. (2008) used the semi-empirical Szyszkowski equation (Szyszkowski, 1908) to model the bulk–surface partitioning of various fatty acids with increasing carbon chain length. Importantly, they noted that using the surface tension of pure water while accounting for changes to the Raoult effect from bulk–surface partitioning led to good agreement with experimental data and less costly calculations; however, the use of accurate values of a_w that accounted for bulk–surface partitioning was critical. These assumptions were more thoroughly explored by Prisle et al. (2010), wherein they found that particles with at least 50 % of their mass composed of surfactant species required more accurate treatments of surface tension depression than considering σ to be the same as that of pure water. In systems with lower concentrations of surfactants it was noted that the surface tension may be similar to that of pure water at the point of CCN activation. Despite this, an evolving surface tension depression in growing aerosol particles that have not yet activated may play a role in their growth and activation behavior (Ovadnevaite et al., 2017; Davies et al., 2019). Further complicating the issue is a lack of measurements of surface tension on sub-500 nm diameter particles.

Romakkaniemi et al. (2011) included a Szyszkowski–Langmuir treatment of bulk–surface partitioning to droplets containing methylglyoxal, a semivolatile species and moderate surfactant. They found that treating the surface as a hypothetical 2D plane that mixes ideally led to an over 1 order of magnitude increase in the total particle-phase concentration of methylglyoxal, further confirming the importance of bulk–surface exchange treatments. Beyond Szyszkowski–Langmuir isotherm-based treatments of bulk–surface partitioning, other equations of state have been employed, including 2D Van der Waals models, compressed film models, and LLPS-based models of bulk–surface partitioning (Ruehl and Wilson, 2014; Ruehl et al., 2016; Ovadnevaite et al., 2017). The compressed film model of Jura and Harkins (1946) (Ruehl et al., 2016) was also utilized by Forestieri et al. (2018) to examine long-chain fatty acid coatings on sodium chloride particles to mimic sea spray aerosols. It was found that different surfactant species may have large effects on CCN activation through both their impacts on σ and the effective hygroscopic growth parameter κ under high-RH conditions. They note that compounds traditionally thought of as highly surface active, like the fatty acids in their study, may not have as large of an impact on CCN activation as others (Forestieri et al., 2018).

The treatment of bulk–surface partitioning was also studied using an AIOMFAC-based coupled liquid–liquid equilibrium and gas–particle partitioning calculation (Ovadnevaite et al., 2017; Davies et al., 2019). The effective surface tension was estimated based on the predicted LLPS phase compositions and the surface coverage by an organic-rich shell phase, constrained to be of a defined minimum thickness δ_β,min. This approach found that LLPS aerosols can be strongly affected by surface tension reductions, yet weakly by changes to the Raoult effect following bulk–surface partitioning (Ovadnevaite et al., 2017; Davies et al., 2019). The treatment of bulk–surface partitioning in the LLPS-based approach was primarily due to bulk equilibrium LLPS, organic surface coverage, and gas–particle partitioning, which directly account for substantial nonideal mixing. However, other more detailed effects, such as a size-dependent feedback from bulk–surface partitioning on LLPS phase compositions as well as liquid–liquid–interfacial-phase partitioning and energy effects were not accounted for. Furthermore, this approach did not depend on any assumptions about the maximum thickness of the interfacial region, only a prescribed minimum thickness (Ovadnevaite et al., 2017). Likewise, a subsequent study by Davies et al. (2019) compared a compressed film model and three versions of AIOMFAC-based bulk–surface partitioning models. The first AIOMFAC-based approach involved using a full (unconstrained) liquid–liquid equilibrium calculation. This methodology assumed that the organic-rich phase β would form a spherical shell around phase α if there was sufficient material; otherwise phase β formed a partial spherical shell of thickness δ_β,min over phase α. An area-weighted mean of pure-component surface tensions based on the areas of each phase exposed to the gas phase was used to determine the effective droplet surface tension. The second approach involved assuming a complete phase separation among organics and aqueous inorganic electrolytes, with only water allowed to partition between both phases and the assumption that the organic species formed a film, i.e., a partial monolayer or up to multiple molecular layers at the droplet surface. In the case of insufficient organic material for forming a complete monolayer over the droplet, the effective surface tension of the droplet is computed as the surface coverage area-weighted average of the pure organic species surface tension and the surface tension of water. These two treatments of surface tension provided a good and predictive estimate of the upper and lower bounds of the measured critical supersaturation of a CCN, given the initial dry composition of the particle.

Malila and Prisle (2018) developed a semi-empirical monolayer-based bulk–surface partitioning model based on an extension of an earlier method developed by Laaksonen and Kulmala (1991) wherein the surface tension of a droplet was related to a surface-composition-weighted average of pure-component surface tensions (σ_i^∘):

$$\begin{array}{}\text{(6)}& \mathit{\sigma }(x^{\mathrm{b}},T)={\displaystyle \frac{{\sum }_i{\mathit{\sigma }}_i{}^{\circ }{\mathcal{V}}_i{x}_i^{\mathrm{s}}}{{\sum }_i{\mathcal{V}}_i{x}_i^{\mathrm{s}}}}.\end{array}$$

Here, 𝒱_i and $x_i^{\mathrm{s}}$ are the molar volumes and surface mole fractions of i, respectively. Coupled with mass conservation, this equation must be solved iteratively from a given droplet size and overall composition. It is important to note that pseudo-binary approximations must be made for Eq. (6) for systems with three or more components and that the iterative solution must be optimized for each system tested, which somewhat limits the utility of this approach as a predictive tool.

Vepsäläinen et al. (2022) conducted a comparative study examining the differences between the 2D Gibbs adsorption model of Prisle et al. (2011), the simplified complete partitioning model of Prisle et al. (2011), the compressed film model of Ruehl et al. (2016), a partial monolayer model based on Ovadnevaite et al. (2017), the monolayer model of Malila and Prisle (2018), and a simple bulk-composition-based model that did not allow for bulk–surface partitioning. In their study, the hygroscopic growth of 50 nm particles of different moderately surface-active dicarboxylic acid species was modeled. It was noted that the more complex models for bulk–surface partitioning, those which allowed for the partial partitioning of species between the bulk and surface, had the best agreement with measured critical supersaturations. Despite the agreement between the more complex models with measured critical supersaturations as a function of dry particle size, they predicted varying degrees of bulk–surface partitioning for different species and, thus, different equilibrium compositions for the bulk and surface of the droplet. A related model comparison by Vepsäläinen et al. (2023) was also extended to systems including stronger surfactant species, such as myristic acid. In that recent work, it was also noted that the various bulk–surface partitioning treatments from Vepsäläinen et al. (2022) were in agreement at low surfactant mass fractions; however, there was disagreement at higher surfactant mass fractions.

Other recent work with thermodynamics-based models has attempted to predict the degree to which surfactants may cover an aerosol particle (surface) at equilibrium, with many surfactant species expected to have a surface coverage on the order of 60 %–85 %, while some surfactant remains dissolved in the particle bulk, for compositions both below and above the critical micelle concentration (McGraw and Wang, 2021).

While simplifying assumptions about bulk–surface partitioning are reasonable for macroscopic systems, in the case of ultrafine particles, the competing effects of interior bulk-phase depletion and surface accumulation are complex and must be considered. Recent experimental findings have highlighted that even for larger droplets on the order of several micrometers, which contain non-ionic surfactant species similar to those found in atmospheric aerosols, varying surface tensions may be exhibited across particle sizes (Bzdek et al., 2020). Thus, a rigorous framework for determining the equilibrium bulk–surface partitioning of a droplet is necessary to better understand the size dependency of bulk–surface partitioning and related surface tension impacts in aerosol systems. Ideally, such a framework is predictive in design; i.e., it does not need to be fine-tuned or fitted specifically for every system or size range of interest. In the following, we introduce and evaluate such a framework by extending the AIOMFAC-based gas–particle partitioning calculations to include a distinct surface phase present between aerosol particles and the surrounding gas phase; see Fig. 1.

Figure 1 Conceptual diagram of bulk–surface partitioning in a single-bulk-phase spherical aerosol particle or cloud droplet. The thickness of the finite-volume surface phase is represented by δ.

2 Theory and methods

In his seminal works, Gibbs (1874) defined the following equation capable of expressing the so-called “free” energy of a thermodynamic system (G) with negligible energetic contributions from surfaces and interfaces as

$$\begin{array}{}\text{(7)}& G=U-TS-PV.\end{array}$$

Here, U is the internal energy of the system, T is temperature, S is entropy, P is pressure, and V is the volume of the system. Taking the derivative of Eq. (7) and then integrating with respect to the number of moles of i in the system (n_i) at constant T and P gives the following equation for a system with k components:

$$\begin{array}{}\text{(8)}& G=\sum _k{\mathit{\mu }}_k{n}_k.\end{array}$$

Using the information given in Eq. (8), the equilibrium composition of a closed multiphase system can be found when its energy is at its global minimum. However, in order to accurately predict the chemical potential of each component in a system, physicochemical mixing models must be employed to account for non-ideal mixing effects on the chemical potential of i in a liquid-state bulk phase b as given by

$$\begin{array}{}\text{(9)}& {\mathit{\mu }}_i^{\mathrm{b}}={\mathit{\mu }}_i^{\circ ,\mathrm{b}}+RT\ln \left(a_i^{\mathrm{b}}\right),\end{array}$$

where ${\mathit{\mu }}_i^{\circ ,\mathrm{b}}$ is the standard chemical potential, R is the gas constant, and $a_i^{\mathrm{b}}$ is the (chemical) activity of i. The activity of a species is a unitless value which represents the “effective concentration” of component i on a chosen composition scale (e.g., mole fraction or molality) as determined by non-ideal mixing in a given phase of a system. Under equilibrium conditions, the chemical potentials of a component must be equivalent across all phases (e.g., Zuend et al., 2010). A useful choice, often employed when calculating μ_i for multiple phases of a system, is that μ_i^∘ is defined to be identical for all phases. With that definition, the activities of a species must also be equal in all coexisting phases for a system at equilibrium. In order to accurately predict the various $a_i^{\mathrm{b}}$ values over a wide range of solution compositions, complex mixing models, such as the AIOMFAC model, must be employed to account for molecular interactions and shape and size differences.

The AIOMFAC model is a state-of-the-art thermodynamic mixing model capable of accurately predicting the chemical activity of species in atmospheric aerosol systems in a thermodynamically rigorous and consistent manner. AIOMFAC can be considered to be a major extension of the UNIversal quasichemical Functional Group Activity Coefficients (UNIFAC) model by Fredenslund et al. (1975), since AIOMFAC includes the treatment of aqueous electrolytes and interactions between ions and organic molecules. UNIFAC itself is a group-contribution model derived from the UNIversal QUAsi-Chemical (UNIQUAC) theory of liquid mixtures developed by Abrams and Prausnitz (1975). The UNIQUAC theory and model represent a local composition model for mixtures of non-electrolyte components that generalizes the original quasi-chemical theory of Guggenheim (1952) to mixtures consisting of molecules of various shapes and sizes (while Guggenheim's quasi-chemical lattice model was restricted to spherical molecules of approximately equal sizes). The UNIQUAC model, the UNIFAC group-contribution version, and subsequently extended variants like AIOMFAC are therefore all deeply rooted in an advanced statistical mechanics treatment of mixing and interactions among different molecules in solution based on the local composition principle and are more rigorous than Guggenheim's two-component lattice gas model. AIOMFAC has been updated numerous times to include many atmospherically relevant inorganic cations, anions, and organic functional groups so that this model can represent a large number of mixed organic–inorganic aerosol systems, including systems with tens to thousands of components (Zuend et al., 2008, 2011; Yin et al., 2022).

If two or more distinct phases (phases are here indexed by superscript ϕ) are present in a system at constant temperature and pressure, Eq. (8) may need to be expanded to account for the energy associated with the boundaries between the phases. This can be achieved via the addition of a term accounting for the interfacial energy per unit area (σ) scaled by the area of the interface (A) (e.g., Aston and Herrington, 1991):

$$\begin{array}{}\text{(10)}& G=\sum _{\mathit{\varphi }}\sum _j{\mathit{\mu }}_j^{\mathit{\varphi }}{N}_j^{\mathit{\varphi }}+\mathit{\sigma }A.\end{array}$$

The interfacial Gibbs energy contribution can be important for microscopic systems, in which the contribution gains in magnitude relative to the collective energy content of the bulk volume of the phases. This interfacial Gibbs energy is therefore of interest in the context of small aerosol particles and cloud droplets. In the case of macroscopic systems, the interfacial energy term can usually be neglected due to its comparably tiny contribution.

Gibbs defined interfacial energy (or tension) as the excess energy attributed to the 2D boundary present between two phases. However, the choice to describe the boundary between two phases as a 2D surface and deciding its exact location presented several problems because real systems may exhibit concentration gradients near an interface. To account for this, Guggenheim introduced the concept of a 3-dimensional (3D) interfacial phase, which includes the boundary between two distinct phases and the concentration gradients of species therein (Guggenheim, 1940). While often on the order of a molecular monolayer or bilayer in depth, such an interfacial compartment can be treated as a distinct phase and must abide by the thermodynamic equilibrium conditions. Given the abovementioned challenges when using Gibbs' 2D surface treatment and in order to more easily connect geometrically confined interfacial and bulk phases in a mass-conserving finite-volume system, the approach of Guggenheim's 3D interface is adopted in our framework. Thus, for application to fine and ultrafine particles, we introduce an interface with a small but finite thickness (δ). We will focus in the following derivations on particles of spherical shape (droplets).

Following the approach of Cai and Griffin (2005), species adsorbed onto or absorbed into the surface of an aerosol particle contribute to the total volume of the particle. Thus, the interface only extends from the limit of the particle's radius inward, and the total number of condensed-phase molecules of i ($n_i^{\mathrm{tot}}$) is defined as the sum of the number of molecules of i in the bulk ($n_i^{\mathrm{b}}$) and in the surface ($n_i^{\mathrm{s}}$):

$$\begin{array}{}\text{(11)}& n_i^{\mathrm{s}}+n_i^{\mathrm{b}}=n_i^{\mathrm{tot}}.\end{array}$$

If the surface of an aerosol is chemically distinct from the underlying bulk phase, then a particle with one condensed (here liquid) bulk phase, a distinct surface phase, and the surrounding gas phase will satisfy the criterion of equivalency of chemical potentials at equilibrium, as mentioned previously. As described by Aston and Herrington (1991), the chemical potential of a surface can be defined by differentiating Eq. (10) for a surface phase with respect to the molar amount of component i in the surface $n_i^{\mathrm{s}}$:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\left({\displaystyle \frac{\partial G^{\mathrm{s}}}{\partial n_i^{\mathrm{s}}}}\right)}_{T,P,n_j^{\mathrm{s}},A}& ={\left({\displaystyle \frac{\partial G^{\mathrm{s}}}{\partial n_i^{\mathrm{s}}}}\right)}_{T,P,n_j^{\mathrm{s}},\mathit{\sigma }}\\ & -{\left({\displaystyle \frac{\partial G^{\mathrm{s}}}{\partial A}}\right)}_{T,P,n_i^{\mathrm{s}}}\times {\left({\displaystyle \frac{\partial A}{\partial n_i^{\mathrm{s}}}}\right)}_{T,P,n_j^{\mathrm{s}},\mathit{\sigma }}.\end{array}\end{array}$$

This leads to the following expression for the chemical potential of a surface phase:

$$\begin{array}{}\text{(13)}& {\mathit{\mu }}_i^{\mathrm{s}}={\mathit{\xi }}_i^{\mathrm{s}}-\mathit{\sigma }{\mathcal{A}}_i,\end{array}$$

where ${\mathit{\xi }}_i^{\mathrm{s}}={\left(\frac{\partial G^{\mathrm{s}}}{\partial n_i^{\mathrm{s}}}\right)}_{T,P,n_j^{\mathrm{s}},\mathit{\sigma }}$ is the surface energy per unit surface area of the solution (i.e., the interfacial energy per unit area of the gas–liquid interface), and ${\mathcal{A}}_i={\left(\frac{\partial A}{\partial n_i^{\mathrm{s}}}\right)}_{T,P,n_j^{\mathrm{s}},\mathit{\sigma }}$ is the partial molar area of i. 𝒜_i can be found for a finite-depth Guggenheim surface phase of a spherical particle as follows:

$$\begin{array}{}\text{(14)}& {\left({\displaystyle \frac{\partial A}{\partial n_i^{\mathrm{s}}}}\right)}_{T,P,n_i^{\mathrm{s}},\mathit{\delta }}={\displaystyle \frac{{\left(\frac{\partial A}{\partial r}\right)}_{T,P,n_i^{\mathrm{s}},\mathit{\delta }}}{\frac{\partial n_i^{\mathrm{s}}}{\partial r}}}.\end{array}$$

Here, the change in surface area of a sphere with respect to its radius is

$$\begin{array}{}\text{(15)}& {\left({\displaystyle \frac{\partial A}{\partial r}}\right)}_{T,P,n_i^{\mathrm{s}},\mathit{\delta }}=\mathrm{8}\mathit{\pi }r,\end{array}$$

and the volume of the finite-depth surface phase on a spherical droplet is

$$\begin{array}{}\text{(16)}& v^{\mathrm{s}}={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }(r^{\mathrm{3}}-(r-\mathit{\delta }{)}^{\mathrm{3}}).\end{array}$$

Differentiating Eq. (16) with respect to the particle radius assuming a constant surface depth leads to

$$\begin{array}{}\text{(17)}& {\left({\displaystyle \frac{\partial v^{\mathrm{s}}}{\partial r}}\right)}_{\mathit{\delta }}=\mathrm{4}\mathit{\pi }(\mathrm{2}\mathit{\delta }r-{\mathit{\delta }}^{\mathrm{2}}).\end{array}$$

Equation (17) can be converted into an expression for $\frac{\partial n_i^{\mathrm{s}}}{\partial r}$ with the inclusion of the molar volume of i in the surface (${\mathcal{V}}_i^{\mathrm{s}}$). This value is assumed to be the same as the bulk molar volume of pure i, which can be calculated by dividing the molar mass (M_i) by the liquid-state density (ρ_i). This yields

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\partial n_i^{\mathrm{s}}}{\partial r}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }(\mathrm{2}\mathit{\delta }r-{\mathit{\delta }}^{\mathrm{2}})}{{\mathcal{V}}_i^{\mathrm{s}}}}.\end{array}$$

Combining Eqs. (15) and (18) along with the assumption that the molar volumes are additive in order to determine $n_i^{\mathrm{s}}$ leads to Eq. (19) for the partial molar area 𝒜_i:

$$\begin{array}{}\text{(19)}& {\mathcal{A}}_i={\left({\displaystyle \frac{\partial A}{\partial n_i^{\mathrm{s}}}}\right)}_{T,P,n_j^{\mathrm{s}},\mathit{\sigma }}={\mathcal{V}}_i^{\mathrm{s}}{\displaystyle \frac{\mathrm{2}r}{\mathrm{2}\mathit{\delta }r-{\mathit{\delta }}^{\mathrm{2}}}}.\end{array}$$

In the extreme case of the depth of the interface approaching the value of a droplet's radius, ${lim}_{\mathit{\delta }\to r}$ of Eq. (19) reduces to $\frac{\mathrm{2}}{r}{\mathcal{V}}_i^{\mathrm{s}}$. In the opposite limiting case of a macroscopic system, in which δ≪r and δ²≈0, the ${lim}_{r\to \mathrm{\infty }}$ of Eq. (19) is $\frac{\mathrm{1}}{\mathit{\delta }}\times {\mathcal{V}}_i^{\mathrm{s}}$. Integrating partial molar areas by Euler's theorem leads to

$$\begin{array}{}\text{(20)}& A=\sum _j{n}_j^{\mathrm{s}}{\mathcal{A}}_j\end{array}$$

for a multicomponent system with $j=\mathrm{1},\mathrm{\dots },k$ different species present. Note that for cases in which the surface phase extends beyond monolayer thickness, Eq. (20) indicates that the 𝒜_i values will not necessarily correspond to the area per molecule located directly at the 2D gas–liquid surface.

Returning to Eq. (13), the intrinsic chemical potential, ${\mathit{\xi }}_i^{\mathrm{s}}$, of components in the surface phase can be formulated analogously to that of a liquid phase (Aston and Herrington, 1991):

$$\begin{array}{}\text{(21)}& {\mathit{\xi }}_i^{\mathrm{s}}={\mathit{\xi }}_i^{\circ ,\mathrm{s}}+RT\ln \left(a_i^{\mathrm{s}}\right).\end{array}$$

Here, ${\mathit{\xi }}_i^{\circ ,\mathrm{s}}$ is the intrinsic standard chemical potential, and $a_i^{\mathrm{s}}$ is the activity of i in the surface phase. In the case of pure component i, Eq. (13) yields

$$\begin{array}{}\text{(22)}& {\mathit{\mu }}_i^{\circ ,\mathrm{s}}={\mathit{\xi }}_i^{\circ ,\mathrm{s}}-{\mathit{\sigma }}_i{}^{\circ }{\mathcal{A}}_i{}^{\circ },\end{array}$$

where σ_i^∘ is the surface energy per unit area of pure component i. For a more general multicomponent surface phase, Eq. (13) becomes

$$\begin{array}{}\text{(23)}& {\mathit{\mu }}_i^{\mathrm{s}}={\mathit{\mu }}_i^{\circ ,\mathrm{s}}+RT\ln \left(a_i^{\mathrm{s}}\right)+{\mathit{\sigma }}_i{}^{\circ }{\mathcal{A}}_i{}^{\circ }-{\mathit{\sigma }}_i{\mathcal{A}}_i.\end{array}$$

Assuming that ${\mathcal{A}}_i{}^{\circ }\approx {\mathcal{A}}_i$ under all conditions, defining ${\mathit{\mu }}_i^{\circ ,\mathrm{s}}={\mathit{\mu }}_i^{\circ ,\mathrm{b}}$, and further requiring that at equilibrium the chemical potentials of species across phases are equal, the effective surface tension of a multicomponent solution can be isolated to form the following equation:

$$\begin{array}{}\text{(24)}& {\mathit{\sigma }}_i={\mathit{\sigma }}_i{}^{\circ }+{\displaystyle \frac{RT}{{\mathcal{A}}_i}}\ln \left({\displaystyle \frac{a_i^{\mathrm{s}}}{a_i^{\mathrm{b}}}}\right).\end{array}$$

At equilibrium all σ_i values must be equivalent to a common value, the effective solution surface tension (σ) . In this case, Eq. (24) can be recognized as the Butler equation (Butler and Kendall, 1932; Sprow and Prausnitz, 1966). Employing Eq. (24) for each component in a system, along with Eq. (11), leads to a system of equations for the equilibrium composition of a 3D Guggenheim interface. These equations can be solved iteratively in a nested manner in combination with the AIOMFAC-based gas–particle partitioning and liquid–liquid-phase-separation algorithms laid out in Zuend et al. (2010) and Zuend and Seinfeld (2013).

2.1 Calculation of surface composition

In order to determine the composition of the surface phase in practice, a system of equations which allows maximum freedom of variable ranges, yet simultaneously satisfies the volume balance constraints of the 3D interface, is required. Furthermore those equations must not modify variables that should be held constant during the calculation of partial derivatives. For this purpose, we have developed a conforming implementation within the extended AIOMFAC equilibrium model. We begin by introducing a term representing the fractional amount of component i relative to the total amount of i in the particle phase:

$$\begin{array}{}\text{(25)}& {\mathit{\epsilon }}_i={\displaystyle \frac{n_i^{\mathrm{s}}}{n_i^{\mathrm{tot}}}},\end{array}$$

where $n_i^{\mathrm{s}}$ denotes the molar amount in the surface phase, and $n_i^{\mathrm{tot}}$ is the total molar amount of i available for partitioning (a constant during partitioning calculations). From the values of r_p and δ, the total volume of the surface phase is calculated according to Eq. (16). Likewise $n_i^{\mathrm{s}}$, ρ_i, and M_i can be used to determine the volume of i in the surface phase. The fraction of the total surface volume occupied by i is given by

$$\begin{array}{}\text{(26)}& f_i={\displaystyle \frac{v_i^{\mathrm{s}}}{v^{\mathrm{s}}}}={\displaystyle \frac{v_i^{\mathrm{s}}}{{\sum }_j{v}_j^{\mathrm{s}}}}.\end{array}$$

We now define a new variable, ζ_i, expressing the fraction of i in the surface relative to its assigned volume range such that

$$\begin{array}{}\text{(27)}& {\mathit{\zeta }}_i{v}_i^{\mathrm{rg}}+v_i^{\min }=v_i^{\mathrm{s}\star },\end{array}$$

where $v_i^{\mathrm{rg}}=v_i^{\max }-v_i^{\min }$, and $v_i^{\max }$ and $v_i^{\min }$ are the respective maximum and minimum possible surface volume contributions. $v_i^{\mathrm{s}\star }$ is the unnormalized surface volume contribution of i. The values of ζ_i can vary in the range from 0.0 to 1.0; however, there can be non-zero $v_i^{\min }$ values for some components (e.g., for water at very high RH) in order to achieve volume closure between the targeted geometric surface shell volume and the (unnormalized) volume as calculated by summing up the surface quantities such that $v_i^{\min }\ge \mathrm{0}$ and $v_i^{\max }\ge v_i^{\min }$. For a given particle composition, r_p, and δ, both $v_i^{\min }$ and $v_i^{\max }$ are species-specific constants. As such, $v_i^{\min }$ can be determined depending on whether all other components j≠i can occupy all of the surface volume if they are at their maximum abundance in the surface. If this is the case, then $v_i^{\min }=\mathrm{0}$. Otherwise, $v_i^{\min }>\mathrm{0}$. This leads to

$$\begin{array}{}\text{(28)}& v_i^{\min }=max\left[v^{\mathrm{s}}-\sum _{j,j\ne i}{v}_j^{\max },\mathrm{0}\right].\end{array}$$

Here, $v_i^{\max }$ is the maximum possible volume in the surface phase, less than or equal to v^s, such that

$$\begin{array}{}\text{(29)}& v_i^{\max }=min\left[n_i^{\mathrm{tot}}\times {\mathcal{V}}_i,v^{\mathrm{s}}\right].\end{array}$$

Returning to ζ_i, Eq. (27) can be rearranged as

$$\begin{array}{}\text{(30)}& {\mathit{\zeta }}_i={\displaystyle \frac{v_i^{\mathrm{s}\star }-v_i^{\min }}{v_i^{\max }-v_i^{\min }}}.\end{array}$$

For a given value of ζ_i, $v_i^{\mathrm{s}\star }$ can be computed and then normalized to find f_i via

$$\begin{array}{}\text{(31)}& f_i={\displaystyle \frac{v_i^{\mathrm{s}\star }}{{\sum }_j{v}_j^{\mathrm{s}\star }}},\end{array}$$

and a value for $v_i^{\mathrm{s}}$ can be determined via

$$\begin{array}{}\text{(32)}& v_i^{\mathrm{s}}=f_i\times v^{\mathrm{s}}.\end{array}$$

We note that such calculated f_i and $v_i^{\mathrm{s}}$ values may violate the condition that $v_i^{\mathrm{s}}\le v_i^{\max }$. This issue is remedied by introducing the relative deviation of actual to targeted surface-phase volume, an additional equation (constraint) to be solved alongside the equations describing the partitioning of k−1 components. From $v_i^{\mathrm{s}}$, variables $n_i^{\mathrm{s}}$ and ε_i can be computed via $n_i^{\mathrm{s}}=\frac{v_i^{\mathrm{s}}}{{\mathcal{V}}_i}$ and then applied in Eq. (25) to obtain ε_i.

2.2 Initial guess generation

An initial guess for the surface composition of a given particle can be derived from a first calculation for a non-partitioning (superscript “np”) case. In this trivial case, the relative compositions of the surface and bulk phases are set to be identical such that $x_i^{\mathrm{s}}=x_i^{\mathrm{b}}$ and $f_i^{\mathrm{s}}=f_i^{\mathrm{b}}$. From this case, an estimation of ${\mathit{\zeta }}_i^{\mathrm{np}}$ can be computed given that

$$\begin{array}{}\text{(33)}& f_i^{\mathrm{np}}={\displaystyle \frac{n_i^{\mathrm{tot}}{\mathcal{V}}_i}{{\sum }_j{n}_j^{\mathrm{tot}}{\mathcal{V}}_j}}\end{array}$$

and

$$\begin{array}{}\text{(34)}& v_i^{\mathrm{s}\star ,\mathrm{np}}=min\left(f_i^{\mathrm{np}}{v}^{\mathrm{s}},v_i^{\max }\right)\end{array}$$

to yield

$$\begin{array}{}\text{(35)}& {\mathit{\zeta }}_i^{\mathrm{np}}={\displaystyle \frac{v_i^{\mathrm{s}\star ,\mathrm{np}}-v_i^{\min }}{v_i^{\max }-v_i^{\min }}}.\end{array}$$

Given the definition that a_i=x_iγ_i for neutral components (and $a_{\pm ,i}=\frac{m_{\pm ,i}}{m^{\circ }}{\mathit{\gamma }}_{\pm ,i}$ for electrolytes), Eq. (24) can be rearranged to the following form:

$$\begin{array}{}\text{(36)}& \ln \left({\displaystyle \frac{x_i^{\mathrm{s}}}{x_i^{\mathrm{b}}}}\right)=(\mathit{\sigma }-{\mathit{\sigma }}_i{}^{\circ }){\displaystyle \frac{{\mathcal{A}}_i}{RT}}-\ln \left({\displaystyle \frac{{\mathit{\gamma }}_i^{\mathrm{s}}}{{\mathit{\gamma }}_i^{\mathrm{b}}}}\right).\end{array}$$

As long as the same nonideal mixing model is used for bulk and surface activity coefficients (in this case, AIOMFAC), the activity coefficient ratio on the right-hand side of Eq. (36) is equal to 1 (only in this non-partitioning case). For a spherical particle of known surface volume v^s as well as $f_i^{\mathrm{s}}$ and 𝒱_i values, the non-partitioning assumption enables the calculation of the molar-phase amounts and the surface-to-bulk-molar ratios, $\frac{n_i^{\mathrm{s}}}{n_i^{\mathrm{b}}}$. Using a composition-weighted mean of the pure-component surface tensions for σ allows for the evaluation of Eq. (36), which yields $\frac{x_i^{\mathrm{s}}}{x_i^{\mathrm{b}}}$. For neutral components (mole fraction scale), the obtained mole fraction ratios can then be converted into a (new) guess for the set of ε_i values as follows:

$$\begin{array}{}\text{(37)}& {\mathit{\epsilon }}_i^{\mathrm{guess}}={\displaystyle \frac{n_i^{\mathrm{s}}/n_i^{\mathrm{b}}}{\mathrm{1}+(n_i^{\mathrm{s}}/n_i^{\mathrm{b}})}}={\displaystyle \frac{x_i^{\mathrm{s}}}{x_i^{\mathrm{b}}}}{\displaystyle \frac{{\sum }_j{n}_j^{\mathrm{s}}}{{\sum }_j{n}_j^{\mathrm{b}}}}.\end{array}$$

For electrolyte components with activities defined on molality scale, the corresponding bulk–surface partitioning guess is generated via scaling by the surface-to-bulk-phase mass ratio,

$$\begin{array}{}\text{(38)}& {\mathit{\epsilon }}_i^{\mathrm{guess}}={\displaystyle \frac{x_i^{\mathrm{s}}}{x_i^{\mathrm{b}}}}{\displaystyle \frac{{\sum }_l{n}_l^{\mathrm{s}}{M}_l}{{\sum }_j{n}_l^{\mathrm{b}}{M}_l}},\end{array}$$

where the summation index l covers non-electrolyte components (i.e., solvents) only. Using the determined set of ${\mathit{\epsilon }}_i^{\mathrm{guess}}$ values, Eq. (24) can be evaluated to obtain updated activity coefficient ratios, $\frac{{\mathit{\gamma }}_i^{\mathrm{s}}}{{\mathit{\gamma }}_i^{\mathrm{b}}}$, as well as an updated weighted mean estimate of the surface tension, which can then once more be evaluated with Eq. (36) and processed to obtain an updated ε^guess vector. If desired, one can expand on this approach by using the determined activity coefficient ratios (in this case kept fixed) together with a set of distinct guesses for potential equilibrium σ values in Eq. (36), yielding a set of initial guesses for the ε_i values. Systematically generating more than one initial guess is useful when the subsequent numerical solution of the system of nonlinear equations, given by Eqs. (24) or (40), is unsuccessful in the case of the first initial guess evaluated – or to further explore whether more than one solution may exist. In our modern Fortran implementation, the system of equations is solved by a modified, bound-constrained version of Powell's hybrid method (Moré et al., 1980, 1984). Our extensive numerical testing suggests that this approach results in a fast and robust method for finding the equilibrium bulk–surface partitioning state for a given overall particle composition, radius, and interfacial thickness.

2.3 Model assumptions

A key piece of information that is necessary to solve Eq. (24) is the liquid-state pure-component surface tension, σ_i^∘, at given temperature. In this work, σ_i^∘ values of organic components were taken from published data (Hyvärinen et al., 2006; Riipinen et al., 2007; Booth et al., 2009) or, in the case of glutaric acid, extrapolated from high-concentration data of binary aqueous solutions (Booth et al., 2009). For inorganic electrolyte components, σ_i^∘ was calculated using the approach for estimating pure molten salt surface tensions as described in Dutcher et al. (2010). Organic compounds with poorly constrained σ_i^∘ values were assumed to have values of 35 mJ m⁻². In addition, there are physical constraints applied to the surface-phase thickness δ in this study. The lower limit of δ was selected to be 0.15 nm or the approximate length of a single alkane C−C bond. The upper bound for δ was selected to be 1.0 nm for all systems tested, which is the approximate length of three water molecule diameters. An alternative assumption was also tested: a treatment where δ is a function of particle- or surface-phase composition; however, there is limited information on how exactly δ should change as a function of composition. Therefore, in this test case it was assumed that a simple weighted average of molecular lengths based on the surface mole fraction in the particle phase could be used. These molecular lengths were computed based on 𝒱_i for each species and the assumption that the molecular length scale can be approximated by the side length of a volume-equivalent cube. Following the calculation of the surface composition, δ was updated, and the surface composition was recalculated. This process was repeated iteratively until convergence (within a set tolerance) to a stable δ value. It is also important to note that Eqs. (24) and (40) cannot be solved directly for a species which is completely insoluble in either the surface or the bulk phase. As such, the relative surface tension deviations, to be solved for from these equations, are scaled by a smooth (∼ rectangular step) weighting factor expressed by the following function:

$$\begin{array}{}\text{(39)}& \mathrm{\Delta }{\mathit{\sigma }}_i={\displaystyle \frac{{\mathit{\sigma }}_i-\overline{\mathit{\sigma }}}{\left|\overline{\mathit{\sigma }}\right|+{\mathit{\tau }}_{\mathit{\sigma }}}}\times {\displaystyle \frac{{\mathit{\epsilon }}_i(\mathrm{1}-{\mathit{\epsilon }}_i)}{{\mathit{\epsilon }}_i(\mathrm{1}-{\mathit{\epsilon }}_i)+\mathrm{100}\sqrt{\mathit{ϵ}}\cdot \exp \left[-\sqrt{{\mathit{\epsilon }}_i(\mathrm{1}-{\mathit{\epsilon }}_i)}\right]}}.\end{array}$$

Here, ϵ is the floating point machine precision employed, $\overline{\mathit{\sigma }}$ the weighted mean surface tension, and τ_σ a tolerance value, typically set to 0.1 J m⁻². The weighting factor (expression after ×) on the right-hand side of Eq. (39) is evaluated to be near 1.0 in most cases and smoothly transitions to substantially smaller values only as i becomes very close to insoluble in either the surface phase or the bulk phase, in which case ε_i approaches either 1 or 0. Therefore, in cases of extremely limited solubility of i in the surface or bulk, the contribution of i to the system of equations used to solve for bulk–surface equilibrium is diminished (a desired property, since the numerical uncertainty grows near the limits of solubility, and numerical precision limitations become substantial). We note that the weighting factor is computed prior to (but not during) numerically solving the system of equations and is only updated if deemed necessary afterwards, such as when the solver was unsuccessful for set numerical tolerances and the equation solving needed to be repeated. Furthermore, a closely related weighting factor, which is normalized by the sum of weightings such that the resulting fractional weights sum to 1.0, is used in the calculation of the weighted mean $\overline{\mathit{\sigma }}$ value during the process of iteratively solving the system of equations (i.e., simultaneously solving Eq. 39 via Eq. 24 for all components).

An additional assumption made in this study concerns the computation of activity coefficients. It is assumed that there are no modifications to the calculation (by AIOMFAC) of activity coefficients in the surface phase compared to the bulk phase. Lane (1983) introduced a common exponential correction factor, t, which is applied to the calculated activity coefficients of a surface phase such that ${\mathit{\gamma }}_i^{\mathrm{s}}={\left({\mathit{\gamma }}_i^{\mathrm{s},\mathrm{calc}}\right)}^t$. The introduction of such an exponent is motivated by the idea that activity coefficients in a surface, affected by some limitations in the molecular packing options, may deviate slightly from those calculated for a bulk phase of identical composition and temperature. However, estimated values of t are system-dependent, yet often close to unity (Lane, 1983). Because the inclusion of additional semi-empirical terms limits the flexibility of the targeted predictive capability of the model developed in this work, it is assumed that t=1 for all systems.

Returning to Eq. (24), if an alternate assumption is made about the partial molar areas, specifically that ${\mathcal{A}}_i{}^{\circ }\ne {\mathcal{A}}_i$, then the expressions previously leading to Eq. (24) result in

$$\begin{array}{}\text{(40)}& {\mathit{\sigma }}_i={\mathit{\sigma }}_i{}^{\circ }{\displaystyle \frac{{\mathcal{A}}_i{}^{\circ }}{{\mathcal{A}}_i}}+{\displaystyle \frac{\mathrm{1}}{{\mathcal{A}}_i}}RT\ln \left({\displaystyle \frac{a_i^{\mathrm{s}}}{a_i^{\mathrm{b}}}}\right),\end{array}$$

and 𝒜_i^∘ can be found for a pure droplet of i analogous to Eq. (19), with the only modification being that the interfacial thickness used is that of the pure component, i.e., set $\mathit{\delta }={\mathit{\delta }}_i{}^{\circ }$. If it is assumed that only a monolayer of molecules forms the surface of a pure-component droplet, then δ_i^∘ can be estimated based on the molecular size of i. With this information, Eq. (40) can be employed with Eqs. (11) and (39) in the same manner as Eq. (24) to form a system of equations for numerically solving the equilibrium bulk–surface partitioning problem.

Figure 2 Predicted and measured surface tensions for (a) a binary water–glutaric acid system and (b) a binary water–sodium-chloride system as a function of solute concentration in water at 298 K. The area shaded in grey represents the predicted surface tension bounded by the proposed limits for δ of 0.1 and 1.0 nm. The solid black line is the predicted surface tension for δ=0.3 nm. Also shown in (c) and (d) are the same systems as (a) and (b), respectively, but with the assumption that ${\mathcal{A}}_i\ne {\mathcal{A}}_i{}^{\circ }$, leading to a modified form of the Butler equation (Eq. 40) and the testing of the surface-composition-based approach for determining δ.

3 Results and discussion

3.1 Comparison of measured and predicted surface tension

In order to determine the validity of the described predictive model, comparisons were made to measurements for a selection of atmospherically relevant binary systems. Figure 2a shows the predicted surface tension (utilizing AIOMFAC with Eq. 24) as a function of the total particle-phase concentration of glutaric acid in a binary water–glutaric acid droplet that was allowed to grow hygroscopically from a starting dry diameter of 5 µm to the point of cloud droplet activation, corresponding to a diameter of approximately 10 µm. This scenario allowed us to compare predicted surface tensions to the bulk tensiometry and optical tweezer measurements taken by Bzdek et al. (2016). The three curves shown in Fig. 2a correspond to three different values of δ: 0.1 nm, which is the approximate length of a single carbon–carbon bond; 0.3 nm, which is the approximate length scale of a single water molecule (the same value was also used as a minimum thickness in models by Davies et al., 2019 and Ovadnevaite et al., 2017); and 1.0 nm, which corresponds to the approximate size of a single glutaric acid molecule along its longest axis. It is shown that assuming a thinner interfacial thickness value leads to a surface tension curve which is highly sensitive to the overall concentration of glutaric acid in the droplet. Analogously, an interfacial thickness that is larger will require greater changes in the total particle-phase concentration of an organic species in order to observe a similar decrease in surface tension. In addition to the binary water–glutaric acid system, Bzdek et al. (2016) also analyzed a binary water–sodium-chloride system of the same size using the bulk tensiometry and optical tweezer approaches. Shown along with these data in Fig. 2b are bulk measurements by Ozdemir et al. (2009). While the surface tension increases with salt concentration in this case, a similar behavior can be seen in the surface tension vs. concentration curves for the inorganic electrolyte system, wherein the curve is sensitive to the selection of δ. However, for the system shown in Fig. 2b, the sensitivity of the modeled surface tension to the value of δ is opposite of that for surface-active organic species. This is likely due to the fact that electrolytes preferentially partition into the bulk phase; therefore, a thinner interface will contain a higher mole fraction of water (at a specific wet diameter) and requires greater concentrations of electrolytes to increase a droplet's surface tension. Additional sensitivity comparisons were performed to determine the effect of modifying σ_i^∘ by ±10 %. It was found that glutaric acid was more sensitive to modifications in the value of σ_glutaric^∘ than NaCl, which was very weakly sensitive to increases in σ_NaCl^∘ and somewhat sensitive to reductions in σ_NaCl^∘. Overall, using δ=0.3 nm leads to surface tension predictions in better agreement with the measurements shown for large droplets or macroscopic solutions. This suggests that a molecular monolayer assumption for representing the surface phase is a relatively good model, at least for the systems shown in Fig. 2a and b.

Figure 2c and d show the effects of exploring different assumptions of Eq. (24) along with the same measurements from Fig. 2a and b. The modified Butler equation, which assumes that ${\mathcal{A}}_i{}^{\circ }\ne {\mathcal{A}}_i$ (Eq. 40 in this work), consistently overestimates the surface tension of the water–glutaric-acid system and only agrees with the water–NaCl system at very low concentrations. Assuming that δ is surface composition-dependent also leads to poor agreement with experimental data at high concentrations for the water–glutaric-acid system and better agreement at low concentrations. For the water–NaCl system, it can be seen that the opposite is true; the composition-dependent δ leads to good agreement at higher concentrations than at lower concentrations. The combination of Eq. (40) and a surface-composition-dependent δ leads to poor agreement across all concentrations shown for the water–glutaric-acid system and only gives good agreement for the water–NaCl system at low concentrations.

It is important to note that the sizes of the droplets analyzed in Bzdek et al. (2016) are large relative to the more numerous but substantially smaller atmospheric aerosol particles of importance in cloud formation. However, there is a paucity of surface tension data available for droplets in the sub-500 nm size range, since it is extremely difficult to measure the surface tension of atmospherically relevant aerosol particles that are freely suspended while at sizes near or below the wavelengths of visible light. Recent developments in atomic force microscopy have allowed surface tension measurements of sub-500 nm particles to be taken (Lee et al., 2017) directly and are indeed promising as a source of measurement data at that size scale. In the present study, such measurements were not analyzed because placing the droplet on a glass substrate introduces a second interface, the substrate–droplet interface, which may modify the partitioning behavior of different species as well as affect the geometry of the droplet from spherical to approximately semi-spherical. At the moment, the framework laid out in this work is only capable of handling spherical geometries with a single gas–droplet surface; the inclusion of interfaces between two condensed phases and the treatment of non-spherical geometries are the subject of future work.

Figure 3 Predicted surface tensions of a binary water–adipic acid system across 8 orders of magnitude in dry particle diameter (see legend in panel c) for interfacial thickness values of (a) 0.1 nm, (b) 0.3 nm, and (c) 0.5 nm. (d) The specific concentrations of adipic acid for which the surface tension is the most sensitive to change in the overall mole fraction of adipic acid as a function of dry diameter and δ. Also shown in (a), (b), and (c) are measurements taken by Riipinen et al. (2007) and Booth et al. (2009).

3.2 Bulk-phase depletion

While the measured particles associated with the data from Fig. 2 are large, the effects of bulk-phase depletion cannot be neglected entirely. Figure 3a, b, and c show surface tension as a function of particle size and δ for binary water–adipic acid systems of particles with dry diameters of 1.0 nm, 10 nm, 100 nm, 1 µm, 10 µm, 100 µm, 1 mm, and 1 cm (ranging over 8 orders of magnitude). It is shown that even highly diluted particles with dry diameters below 100 µm have different surface tensions at the same total mole fraction of adipic acid in the condensed phase. If the value of δ decreases, the surface tension curves for larger dry diameters converge, and bulk-phase depletion is only noticeable for the smallest particles. If δ is increased to be more similar to the molecular length scale of the solute species (∼0.45 nm), there is better agreement between modeled and measured surface tension values in the larger-sized droplets. Increases in δ also increase the minimum concentration of solute necessary to decrease the surface tension from that of pure water for all particle sizes; however, smaller particle sizes are more responsive to this change than larger ones. A method for determining this dependence is by taking the concentration of solute ($c_i^{\mathrm{d}\mathit{\sigma }}$) at which the change in surface tension with each additional molecule of solute added is the greatest, in other words, the global extrema (minimum for surfactant species and maximum for tensoionic species) of $\frac{\mathrm{d}\mathit{\sigma }}{\mathrm{d}{X}_i^{\mathrm{total}}}$. Figure 3d shows how $X_i^{\mathrm{d}\mathit{\sigma }}$ varies as a function of size for the binary water–adipic acid system. Droplets with initial dry diameters on the order of 1 µm may still experience mild bulk depletion effects depending on the selection of the interfacial thickness value, with larger δ values leading to more pronounced bulk-phase depletion effects at larger sizes. Moreover, we note that the value of $X_i^{{\mathrm{d}}^{\mathrm{2}}\mathit{\sigma }}$ is strongly dependent on δ at larger sizes compared to smaller ones, with a doubling of the interfacial thickness from 0.15 to 0.3 nm leading to approximately 3 orders of magnitude increase in the value of $X_i^{\mathrm{d}\mathit{\sigma }}$ for particles with dry diameters larger than 1 µm. This same change in δ for particles with dry sizes below 100 nm leads to differences of about 2 orders of magnitude or less. Thus, we demonstrate that bulk-phase depletion may modify surface tension as a function of particle composition on particles up to the micrometer size scale.

Figure 4 Predicted bulk–surface partitioning coefficient ($\frac{x_i^{\mathrm{surf}}}{x_i^{\mathrm{bulk}}}$) of (a) water, (b) glutaric acid, and (d) sodium chloride present in a forced single-bulk-phase particle at T=298 K. A molar dry solute ratio of 1:1 was used in all cases. (d) Predicted effective surface tension for several particle dry diameters as indicated by color. Right column (composition bar graphs): shown are the mole fractions of each species in the surface and the bulk phase (α) for a particle of 25 nm dry diameter.

The bulk–surface concentration ratio $\frac{x_i^{\mathrm{surf}}}{x_i^{\mathrm{bulk}}}$ is independent of the volume ratio and represents the physicochemical partitioning coefficient. Figure 4 shows these mole fraction ratios for each species in a ternary water–glutaric-acid–sodium-chloride system as a function of equilibrium saturation ratio (S). Monodisperse particles consisting of a 1:1 molar ratio of glutaric acid to sodium chloride, with dry diameter values ranging from 10 to 500 nm, were allowed to grow hygroscopically until the point of cloud droplet activation (while maintaining the same solute masses, i.e., no gas–particle partitioning of glutaric acid considered). In this test, a forced one-phase calculation was performed (only allowing a single bulk liquid phase to exist), which prevented the particles from undergoing LLPS at compositions where that would be favorable. The value of δ was held constant at 0.3 nm, and Eq. (24) was employed. It can be seen that glutaric acid, while only a weakly surface-active compound compared to lower-polarity organics, is nevertheless strongly enriched in the surface of the particle across all particle sizes, especially at higher values of S. The predicted surface tension for each particle is also shown. The difference in σ of particles is the largest at both very low and very high values of S, with all of the σ values being the most similar for saturation ratios between 0.675 and 0.725. However, particles with diameters below 25 nm exhibit greater deviations from the behavior of their larger counterparts. These smallest particles achieve both a slightly lower minimum surface tension and a slightly higher surface tension value under low-a_w conditions. This may be driven by the fact that under these conditions, more tensoionic species must be present in the surface due to limited amounts of both water and organic species. However, this effect is still quite weak even at such small particle sizes.

Figure 5 The effect of implementing a modified form of the Butler equation (Eq. 40) with the assumption that ${\mathcal{A}}_i{}^{\circ }\ne {\mathcal{A}}_i$ for a water–glutaric-acid–sodium-chloride system with a 1:1 water-free molar ratio of glutaric acid to sodium chloride on (a) the bulk–surface partitioning coefficient, $\frac{x_i^{\mathrm{surf}}}{x_i^{\mathrm{bulk}}}$, and (b) predicted surface tension, σ. All calculations were performed at a temperature of 298 K.

If Eq. (40) is employed for bulk–surface equilibrium predictions, the δ_i^∘ values for water, glutaric acid, and sodium chloride are calculated based on their molecular sizes. This leads to modified bulk–surface partitioning behavior as shown in Fig. 5a and b. The assumption that ${\mathcal{A}}_i{}^{\circ }\ne {\mathcal{A}}_i$ effectively modifies the value of σ_i^∘ by multiplying it by the ratio $\frac{{\mathcal{A}}_i{}^{\circ }}{{\mathcal{A}}_i}$. In the macroscopic case, where r≫δ, then ${lim}_{r\to \mathrm{\infty }}\frac{{\mathcal{A}}_i{}^{\circ }}{{\mathcal{A}}_i}=\frac{{\mathit{\delta }}_i{}^{\circ }}{\mathit{\delta }}$, again showing an important effect of the calculated or assumed value of δ. In the case of species for which ${\mathit{\delta }}_i{}^{\circ }<\mathit{\delta }$, a stronger affinity for the surface results, while for ${\mathit{\delta }}_i{}^{\circ }>\mathit{\delta }$, it is expected that species i would have a weaker affinity for the surface phase. Notably, the feedbacks on the partitioning behavior of one species can still modify the bulk–surface partitioning of other species, since modifications to the composition of the surface phase will modify the activity coefficients for all species in the phase. As such, in Fig. 5, it can be seen that the assumption that $\mathit{\delta }={\mathit{\delta }}_{\mathrm{water}}{}^{\circ }$ leads to decreased partitioning of water, as both glutaric acid and sodium chloride show an increased affinity for the surface phase. This also modifies the effective solution surface tension. In the case where $\mathit{\delta }={\mathit{\delta }}_{\mathrm{water}}{}^{\circ }$, σ is consistently lower than the values predicted by Eq. (24). In the case where δ lies between all of the values of δ_i^∘, σ predicted by Eq. (40) is lower than the value predicted by Eq. (24) when the solution is highly concentrated in both organic and inorganic solutes (low a_w) and higher than Eq. (24) in more dilute cases (high a_w). These findings suggest that it is best to assume that ${\mathcal{A}}_i{}^{\circ }\approx {\mathcal{A}}_i$ for most systems.

Figure 6 Köhler curves for ternary water–suberic-acid–ammonium-sulfate particles at 298 K at an organic volume fraction of 0.88 and with a water-free diameter of (a) 150 nm and (b) 40 nm. Also shown in (a) are the measured SS_crit and the critical wet diameter measured by Ruehl et al. (2016). Panels (c) and (d) show the calculated Köhler curves if the assumption is made that ${\mathcal{A}}_i{}^{\circ }\ne {\mathcal{A}}_i$, leading to the modified Butler equation (Eq. 40), as well as the effect of modifying δ as a function of overall particle composition. Panels (e) and (f) show the calculated σ values for the systems shown in (c) and (d). A complete list of SS_crit and D_crit values for this system and the system shown in Fig. 7 can be found in Table 2.

3.3 Köhler curves

Numerous different test systems have been used in laboratory experiments and modeling studies to better understand the role of bulk–surface partitioning in the behavior of CCN both before and after activation. One such system that will now be considered has been studied in experiments and theory by Ruehl et al. (2016). They generated 150 nm (dry diameter) mixed suberic-acid–ammonium-sulfate particles corresponding to a spherical 50 nm diameter ammonium sulfate core coated with a 50 nm layer of suberic acid. Figure 6a shows the Köhler curve calculated via Eq. (1) with the assumption that a_w in that equation is determined from the bulk-phase composition of the particle. Figure 6b also shows a particle with the same organic volume fraction as Fig. 6a but of a water-free diameter of 40 nm instead of 150 nm. Different interfacial thicknesses (δ) are shown, including the values determined by the compressed film model used in Ruehl et al. (2016). The resulting behavior of the Köhler curve is highly sensitive to the value of δ, particularly as δ approaches the lower limit of physically realistic values. The shapes of these Köhler curves are determined by the point at which the surface tension of the particle becomes similar to that of pure suberic acid relative to a given saturation ratio. In systems with δ values that are smaller, the surface remains are enriched in suberic acid at higher S values than in systems with larger δ values. This lowered surface tension at high S values may lead to modifications of the shape of the Köhler curve, including the branch at sizes greater than the critical wet diameter for CCN activation. Such behavior can be seen in the blue curve of Fig. 6a. Also shown in panels (a) and (b) of Fig. 6 are AIOMFAC-based predictive treatments of bulk–surface partitioning and CCN activation for such systems, as discussed in prior work (Ovadnevaite et al., 2017; Davies et al., 2019). These prior treatments both assume a surface with adjustable depth; consideration of LLPS; and, in the presence or absence of bulk LLPS at higher a_w values, the calculation of an effective surface tension using a volume-fraction-weighted mixing rule of the pure-component surface tension values based on the bulk-liquid-phase compositions (and area fractions in the case of LLPS) in contact with the droplet surface. In this case, both simulations used δ=0.3 nm and the same pure-component surface tension values as used in the more detailed approach developed in this work. Briefly, these are two thermodynamic model variants for the computation of droplet surface tension as a function of composition, size, and temperature, introduced by Ovadnevaite et al. (2017) (see their supplementary information) and also discussed and applied by Davies et al. (2019). AIOMFAC-Equil is a full bulk equilibrium model, including gas–particle partitioning and liquid–liquid-phase separation (LLPS), but without consideration of bulk–surface partitioning. In the context of surface tension predictions, a procedure for the post-processing of model outputs has been introduced by Ovadnevaite et al. (2017), which assumes a core–shell droplet morphology in the case of LLPS. The AIOMFAC-Equil model treats the droplet surface tension as the surface-area-fraction-weighted average of the surface tensions of the present liquid phases (when no complete shell is formed by an organic-rich phase). The initial surface tensions of those phases are computed based on a volume-fraction-weighted mean of the pure-component surface tensions. The AIOMFAC-CLLPS with organic film model variant assumes complete LLPS among organics and inorganics (except for water) at all RH levels. It further assumes that all organic species in the droplet are present in a water-free layer at the surface of the droplet, while all electrolytes are present in an (aqueous) core phase of the droplet. The surface tension of the droplet in this model is equal to the surface tension of the organic film (phase), assuming complete coverage by the organic phase, or the surface-area-fraction-weighted average of the organic phase and the electrolyte-rich aqueous phase, should there be insufficient organic material to completely cover the droplet (Davies et al., 2019). The AIOMFAC-Equil prediction leads to substantially higher critical supersaturations than observed because this model variant ignores bulk–surface partitioning in the case of a single bulk phase present at higher RH, as in the system of Fig. 6. In contrast, the AIOMFAC-CLLPS variant with an imposed organic film assumption agrees reasonably well with the measured critical supersaturations. Indeed, if δ is assumed to be 0.3 nm for both prior AIOMFAC-based treatments (in those cases setting the minimum thickness of the surface phase), then the calculation by the AIOMFAC-CLLPS variant with an organic film represents the measured data better than the approach laid out in this work. However, if the δ value is lowered, then the bulk–surface equilibrium approach from this study shows excellent agreement with the measured peak supersaturation and, importantly, does not rely on the same simplifying assumptions as the organic-film-based calculation. Figure 6b demonstrates the effect of particle size (dry diameter of 40 nm) and bulk–surface partitioning on CCN activation, since the deviations from classical Köhler curve behavior (using a fixed surface tension, that of pure water) are more pronounced in both the δ=0.3 nm and δ=0.15 nm cases. In the example of Fig. 6, the AIOMFAC-Equil prediction is representative of a classical Köhler curve, since that model does not predict LLPS in the high-a_w range close to the CCN activation point for this system. Clearly, in the case of ultrafine particles, a more detailed treatment of bulk–surface partitioning and associated surface tension evolution leads to notable deviations from classical behavior. Indeed a more detailed and thermodynamically rigorous treatment of bulk–surface partitioning leads to better agreement of the activation conditions in comparison to laboratory studies of water–suberic-acid–ammonium-sulfate particles than previous AIOMFAC-based treatments.

To further demonstrate the predictive power of the model developed in this work, an isoprene SOA system was considered as well. This system consists of water and 21 semivolatile isoprene photo-oxidation products, as proposed for a simplified isoprene-derived SOA representation in previous modeling work (Rastak et al., 2017; Gervasi et al., 2020). Because of the higher volatility of some of the isoprene SOA species, the effects of organic co-condensation (or more generally gas–particle partitioning) during hygroscopic growth should also be analyzed to better understand atmospheric implications of such aerosol systems. Concentrations simulated by the Master Chemical Mechanism (Jenkin et al., 2015, 1997, 2012) were used as inputs for the “co-condensation-enabled” case (Rastak et al., 2017) (see Table S3 in the Supplement). A distinct case, in which the organic composition of the particle was fixed, termed the “co-condensation-disabled” case, was also used for comparison. In that case, the equilibrium composition of the particle was taken from a bulk equilibrium gas–particle partitioning calculation at 0.1 % RH and then fixed for the organics (essentially rendering them nonvolatile) such that only water could partition between the gas and particle phases in subsequent computations. We note that for this isoprene SOA system, the values of the various σ_i^∘ values have not been measured. We therefore assumed that ${\mathit{\sigma }}_i{}^{\circ }=\mathrm{35}\mathrm{mJ}{\mathrm{m}}^{-\mathrm{2}}$ for all species, which is in line with the assumptions made in Davies et al. (2019) and Ovadnevaite et al. (2017) (see Fig. S7 in the Supplement for an analysis of the framework sensitivity to σ_i^∘ values). Based on experiments and AIOMFAC LLPS equilibrium computations for bulk solutions, this system is not expected to undergo LLPS at any RH (Rastak et al., 2017).

Figure 7 (a) The effect of including co-condensation of SVOC species for an aqueous 21-component system of isoprene-derived SOA (Rastak et al., 2017; Gervasi et al., 2020) with D_dry=23 nm. (b) The bulk-phase water activity, corresponding to the contribution of the Raoult effect to the Köhler curves shown in (a). Panel (c) shows the Kelvin effect term's contribution to the curves in (a). For the co-condensation-disabled case, the particle's water-free composition was taken at RH=0.1 %, and then only water was allowed to partition to and from particles with this dry composition. For a detailed description of the system components, see Table S1.

Figure 7a shows the Köhler curve for a particle of 25 nm in dry diameter with co-condensation enabled and disabled. Figure 7b shows the contribution of the Raoult effect for both systems shown in Fig. 7a. Likewise Fig. 7c shows the contribution of the Kelvin effect for the same system. The inclusion of co-condensation of organic species leads to substantial reductions in S_crit for this system through modifications to both the Raoult effect and the Kelvin effect.

Goldsack and White (1983)

Table 1 Partial molar areas (m² mol⁻¹) of different organic species calculated using a simplified geometric approach based on the values given in Topping et al. (2007) and the model developed by Goldsack and White (1983), as well as by Eq. (19) of this study for particle diameters of 5, 50, and 5000 nm and two distinct δ values (0.15 or 1.0 nm).

Table 2 Critical wet diameter and supersaturations for the Köhler curves shown in Fig. 6.

4 Theoretical and atmospheric implications

4.1 Theoretical implications

One source of uncertainty in the approaches to bulk–surface partitioning described in Sect. 1 is determining the effective partial molar area of a given species, 𝒜_i, in a mixed-surface phase and how that value may differ from the molar area of pure i, 𝒜_i^∘. A common assumption is that the apparent molar area can be calculated from the molar volume (𝒱_i) of a species as ${\mathcal{A}}_i={\mathcal{V}}_i^{\frac{\mathrm{2}}{\mathrm{3}}}$ and that the molar area of a species in solution is the same as its pure-component value ${\mathcal{A}}_i\approx {\mathcal{A}}_i{}^{\circ }$. Table 1 lists the partial molar areas of numerous organic species calculated using Eq. (19) for particles of three distinct diameters and δ=0.15 nm or δ=1.0 nm. Those diameter choices serve to demonstrate a size dependence in this parameter when calculated via Eq. (19). Also shown are the size-independent partial molar areas when computed under the ${\mathcal{A}}_i=n_A^{\frac{\mathrm{1}}{\mathrm{3}}}{\mathcal{V}}_{b,i}^{\frac{\mathrm{2}}{\mathrm{3}}}$ assumption or with the empirical approach developed by Goldsack and White (1983):

$$\begin{array}{}\text{(41)}& {\mathcal{A}}_i=\mathrm{1.021}\times {\mathrm{10}}^{\mathrm{8}}\times {\mathcal{V}}_{c,i}^{\frac{\mathrm{6}}{\mathrm{15}}}{\mathcal{V}}_{b,i}^{\frac{\mathrm{4}}{\mathrm{15}}},\end{array}$$

where 𝒱_c,i and 𝒱_b,i are the critical and bulk molar volumes. It should be noted that the scaling factor in Eq. (41) requires that the values of 𝒱_c,i and 𝒱_b,i are input in units of cm³ mol⁻¹; the equation then returns 𝒜_i in units of cm² mol⁻¹. In the mathematically sound framework developed in this work, the values of 𝒜_i are weak functions of particle radius r_p; however, they are stronger functions of the value of δ, with smaller δ values leading to smaller 𝒜_i values.

Another important assumption made regarding the treatment of 𝒜_i is the assumption that the density of i in the surface phase (${\mathit{\rho }}_i^{\mathrm{s}}$) is equivalent to that of the pure-component value of a (bulk) liquid state. Currently, all AIOMFAC-based bulk solution calculations make the assumption that the 𝒱_i and related ρ_i of a species do not change as a function of solution composition and that the total volume of a phase is a linearly additive function of the individual component molar volumes times their molar abundance. Other studies have noted that ${\mathit{\rho }}_i^{\mathrm{s}}$ may differ from ${\mathit{\rho }}_i^{\mathrm{b}}$ and that lower values of ${\mathit{\rho }}_i^{\mathrm{s}}$ may lead to better agreement between surface tension models and experimental data (Defay et al., 1966).

It should also be noted that deviations in activity coefficients are possible when comparing surface versus bulk phases of the same molar composition. As mentioned in Sect. 2.3, the introduction of a single exponential scaling factor (t) for all component activity coefficients in the surface phase has been used in the past as a fit parameter in order for binary solution surface tension curves to better match experimental data without violating the Gibbs–Duhem relation. This t value may be thought of as treating the surface-phase nonideality as taking place at a different temperature than that of the bulk phase, since $RT\ln \left[({\mathit{\gamma }}_i^{\mathrm{s}}{)}^t\right]=R\left(tT\right)\ln \left[\left({\mathit{\gamma }}_i^{\mathrm{s}}\right)\right]$. Therefore, one may argue that the value of t should be constrained such that the temperature change remains physically realistic. Thus, a value of t=25 is physically unrealistic for a system at 298 K which shows substantial nonideal mixing, since it would mean that the surface-phase nonideality would be behaving as if it were at 7450 K. The inclusion of t also introduces an additional fit parameter that is likely unique for each system, thus limiting the predictive power of the model introduced in this work. Use of exponent t in combination with an unconstrained fit leads to better agreement with points below cloud droplet activation for some of the Köhler curves presented in Ruehl et al. (2016), such as ternary water–succinic-acid–ammonium-sulfate, water–pimelic-acid–ammonium-sulfate, and water–glutaric-acid–ammonium-sulfate particles. However, those fitted t values must be combined with rather extreme values of δ and σ_org for good agreement with the experimental data both at the CCN activation point and at points below activation (see Figs. S4–S6 for examples). In addition, the framework laid out in this work is incapable of simultaneously matching the growth data points and critical supersaturation point reported by Ruehl et al. (2016) for water–malonic-acid–ammonium-sulfate particles (see Fig. S3). Further explorations of variations in the activity coefficients of species in the surface phase are warranted to better understand how these activity coefficients may differ in value from those of a bulk solution with the same molar composition. Likewise, the explicit treatment of the dissociation of organic acids under dilute aqueous conditions is not considered in regard to its role in bulk–surface partitioning in this study. Under highly diluted conditions, such as those found in activating CCN, many dicarboxylic acids may partially or fully dissociate. It is possible that consideration of such acid dissociation may lead to modifications of both the surface enrichment and the bulk depletion of different species as well as to enhancements of the solute effect via an increase in dissolved ionic species. Explicit treatments of organic acid dissociation and resulting interactions among various additional ions in an AIOMFAC-based model framework are thus a direction to be explored in future work.

An alternative approach to using the framework laid out in this work is to use other statistical mechanics models to predict the surface tension as a function of bulk solution composition. A simplified statistical-mechanics-based approach for surface tension predictions was developed by Wexler and Dutcher (2013), which relies on a single physically constrained fitted parameter, r, which represents the average number of water molecules displaced by a solute molecule at the surface. As shown in Fig. S6, this model had a root mean square error of 2.90 mJ m⁻² in comparison to measurements of surface tension for a (macroscopic) binary water–ethanol system (Ernst et al., 1935). In comparison, the model developed in this work has a root mean square error of 2.930 mJ m⁻², when using our default assumption that the thickness of the interface is δ=0.3 nm. The fitted value of r as reported by Wexler and Dutcher (2013) is 3.00. If the number of water molecules displaced by an ethanol molecule at the surface of a droplet is assumed to be determined based on the respective values of 𝒜_water and 𝒜_ethanol, then the number of water molecules displaced by an ethanol molecule in the surface phase is 3.02. Despite both models being in good agreement in the macroscopic case, it is important to note that their statistical mechanical model does not directly account for bulk-phase depletion in volume-constrained systems. For a comparison of various other frameworks for estimating the surface tension of liquid solutions and/or atmospheric aerosol particles, we refer the reader to the recent work by Kleinheins et al. (2023) and Vepsäläinen et al. (2022, 2023).

Figure 8 Critical supersaturations vs. dry diameter as predicted by the AIOMFAC-Equil model, the AIOMFAC-CLLPS with organic film model, and a model from this work (Eq. 24) for (a) water–suberic-acid–ammonium-sulfate particles and (b) water–isoprene SOA particles. The horizontal shaded bands represent distinct regimes of maximum supersaturations encountered by aerosols in either marine (blue), clean continental (green), background (orange-yellow), or urban polluted (brown) conditions at cloud base conditions, as reported by Pinsky et al. (2014).

4.2 Atmospheric implications

The effect of bulk–surface partitioning on Köhler curve shapes is evident for submicron-sized aerosol particles. Even the use of relatively “thick”, yet reasonable, δ values (>0.5 nm) still exhibits substantial suppression of the critical supersaturation and modifications to the shape of the Köhler curve for the particle size range prior to reaching the CCN activation point under growth conditions. The inclusion of the thermodynamic treatment of equilibrium bulk–surface partitioning outlined in this study leads to simulated droplets that will grow to larger diameters at lower relative humidities than classical Köhler theory would otherwise suggest. If the value of δ is lowered, Köhler curves may exhibit a second local maximum as the CCN surface tension approaches that of pure water after the point of droplet activation. This behavior suggests that particles with very thin surface phases are more likely to activate into cloud droplets (for a given dry diameter). Clouds which form from rising air parcels populated by surfactant-containing particles may exhibit substantially higher cloud droplet number concentrations than those forming from air parcels of comparable particle size distributions but lacking aerosol particles of lowered (yet evolving) surface tensions. Figure 8a shows the critical supersaturation for CCN activation of particles with the same condensed-phase composition as those of Fig. 6, with D_dry values ranging from 25 to 130 nm. Similarly, Fig. 8b shows particles that grow under the same input parameters as used for Fig. 7, with dry sizes from 25 to 130 nm. The colored horizontal bands shown in Fig. 8 correspond to the typical supersaturation values experienced by aerosol particles in marine (blue), clean continental (green), background (orange), and urban polluted (brown) cloud base conditions, according to the classification by Pinsky et al. (2014). If the AIOMFAC-Equil model is used to determine the CCN activation conditions, particles with D_dry between approximately 60 and 75 nm would not be predicted to activate in continental clouds for both systems shown. For the isoprene-derived SOA system, in particular, the inclusion of both equilibrium co-condensation of SVOCs and bulk–surface partitioning leads to larger modifications in predicted S_crit, especially for D_dry>50 nm. This can have important implications for both the radiative forcing effects of resulting clouds and for the precipitation formation microphysics in these clouds, since lowering the critical dry diameter for typical peak supersaturation experienced at cloud base conditions may lead to substantially increased cloud droplet number concentrations (depending on the present aerosol number–size distribution) (Ovadnevaite et al., 2017). It is also important to note that the new framework introduced in this study to date has only considered impacts of organic components on CCN activation, yet does not consider the role that bulk–surface partitioning may play (if any) in ice-nucleating particles in cirrus clouds or mixed-phase clouds. Likewise, the interplay of bulk–surface partitioning and co-condensation with liquid–liquid equilibria is not yet considered in a fully coupled manner. For example, in the case of phase separation, such interactions may influence particle morphology and the RH range within which LLPS occurs. Consequently, coupled size and LLPS effects may also change the interactions between aerosol particles and light.

5 Conclusions

The surface-area-to-volume ratio of atmospheric aerosol particles increases substantially as particle diameter decreases in the fine and ultrafine size ranges. Any unique properties of the exterior surface of an aerosol particle must be accounted for in order to accurately model the behavior of the smallest (sub-100 nm sized) particles. This study builds on the finite-depth Guggenheim surface-phase treatment in combination with variants of the Butler equation and AIOMFAC-based vapor–liquid equilibrium computations to create a thermodynamically rigorous treatment of bulk–surface partitioning in spherical aerosol particles with diameters as low as 10 nm. This model relies on one adjustable, loosely constrained parameter, the surface-phase thickness δ, and applies it consistently to any number of species in multicomponent organic–inorganic aerosol systems. The approach is capable of representing experimentally measured surface tension data for atmospherically relevant systems across a range of relative humidities. The inclusion of a thermodynamically sound treatment of interfacial regions leads to modified Köhler curve predictions that are in agreement with measured data, including for cases for which simpler approaches with fixed surface tension fail. For particles with diameters larger than ∼ 100 nm, the simpler AIOMFAC-CLLPS model variant with an organic film assumption agrees reasonably well with the more thermodynamically sound model in terms of predicted SS_crit values and may serve as a good approximation when computational efficiency is a key concern. For smaller particles, where bulk-phase depletion may play a larger role, larger disagreement arises between the AIOMFAC-CLLPS model with organic film and the approaches laid out in this work.

While measurements of physicochemical properties of particles, such as surface tension and chemical composition, in the size range below 1 µm and especially below 100 nm are rare or nonexistent, there have been numerous measurements made on larger particles. Models trained on data from bulk measurements and large microscopic droplets have been used to study sub-100 nm particles. Frequently, those measurements were done on droplets exhibiting a single liquid (bulk) phase and spherical shape when freely suspended, but phase separation and associated phase boundaries can affect the particle shape. Indeed, many systems have been modeled or observed to adopt more complex, non-spherical morphologies, in some cases involving multiple liquid phases (Huang et al., 2021; Kwamena et al., 2010; Reid et al., 2011). The basic thermodynamic theory introduced in Sect. 2 is generally applicable; however, we show that when applied to finite-volume droplets, geometric considerations introduce shape- and mass-balance constraints which impact the bulk–surface partitioning, particularly in submicron-sized particles. In this study, we have outlined the detailed expressions for spherical single-bulk-phase particles, which were implemented in our AIOMFAC-based bulk–surface partitioning model. To date, this model has not accounted for non-spherical shape or feedback effects from energy stored in liquid–liquid interfaces. However, our model provides a basis for future extensions to account for size-dependent feedback between droplet size, liquid–liquid interfaces, non-sphericity, and size effects on LLPS onset. Furthermore, the use of Eq. (24) or (40) for gas–liquid interfaces requires accurate measurements or predictions of a reference state value of σ, usually pure-component surface tension, for many atmospherically relevant species. These data do not exist or, in the case of inorganic electrolytes, are disputed as to what the correct value should be. This may limit the systems for which our approach can be used – or requires assumptions to be made. Therefore, this study highlights the need for, and benefit of, reliable data for pure-component surface tension. A key goal for applications of thermodynamic multiphase aerosol models in atmospheric chemistry is achieving predictive capability, unrestricted by system- and size-dependent fit parameters. The model introduced here marks a major step toward this goal. It enables us to better quantify the role of interfacial properties in environmental systems on the nanometer and micrometer size scales. Models like those introduced in this study can serve as a bridge between the measurable particle size range and the presently experimentally inaccessible ultrafine size range of interest for cloud formation.