1 Introduction

Surfactants, molecules containing both hydrophobic and hydrophilic parts, are important for a number of atmospheric processes (Brimblecombe and Latif 2004), including cloud formation (McFiggans et al. 2006; Lin et al. 2018), interfacial transport (Abbatt et al. 2012) and heterogeneous chemistry (Rossignol et al. 2016). Various surfactants have been found in cloud water and precipitation samples (Facchini et al. 1999; Orlović-Leko et al. 2010) and more recently also directly in aerosol particles with sizes relevant for cloud activation (Kroflič et al. 2018; Gérard et al. 2019). For example, fatty acids and their salts are common surfactants found in marine background aerosols (Mochida et al. 2002; Tervahattu et al. 2002; Cochran et al. 2017; Shaharom et al. 2018; Bikkina et al. 2019). The amphiphilic nature leads to a preferential adsorption of dissolved surfactant molecules at the solution–air interface. The ensuing partitioning surfactant solute between solution bulk and interface is particularly important in microscopic systems with relatively large surface areas, such as atmospheric droplets (cloud, fog, mist and rain) (Prisle et al. 2010; Malila and Prisle 2018; Lin et al. 2020; Bzdek et al. 2020).

The bulk/interface partitioning equilibrium is challenging to model, not only because of the scarcity of relevant experimental data, such as composition-dependent surface tension and density, for surfactant solutions, but also because of the tendency of surfactants to self–aggregate. Above a limit called the critical micelle concentration (CMC), surfactants form closed structures of spherical or tubular shape that coexist in solution with the remaining surfactant monomers (Langevin 1992). The self–aggregation process is driven by the balance of forces especially around the surfactant heads (Vlachy et al. 2008), leading to aggregates of different shapes and packing structures, including micelles and vesicles (Hashidzume and Harada 2015). Here, we refer to the process in general as micellization and the aggregate structures as micelles.

During micellization, the interaction forces between surfactant hydrophilic parts and other ions in solution are strongly affected by ion specificity and surface exclusion effects related to ion solvation (Warszyński et al. 2002). Surfactant heads develop different interactions with ions in solution that modify its surface-activity. Ions compete for the adsorption sites at the micellar surface depending on their ability to form contact–pair or solvent–separated pairs (Collins 1997; Vlachy et al. 2008; Kunz 2009; Salis and Ninham 2014). Aqueous droplets formed from atmospheric aerosols typically comprise multicomponent solutions of inorganic salts and surfactants (Gérard et al. 2016; Nozière et al. 2017; Gérard et al. 2019) and modeling of the micellization process in these systems is therefore not straightforward. For example, it has been experimentally determined by changing the relative humidity around levitated droplets comprising atmospheric fatty acids/fatty acid salts, that these surfactants can form complex tridimensional nanostructures including hexagonal and cubic close-packed arrangements of spherical and cylindrical micelles into bilayer stacks (Pfrang et al. 2017).

Even at very low concentrations, without micellization, aqueous surfactant solutions show significant deviations from ideal solution behaviour, especially in the case of ionic surfactants. The interaction of fatty acids with divalent cations of calcium and magnesium can lead to formation of stable insoluble complexes that precipitate even at very low concentrations (Rodriguez et al. 2001; Pereira et al. 2012). The greater the difference in size and shape between the hydrophobic and hydrophilic parts of the surfactant molecule, in general the more significant the nonidealities are. For example, aqueous solutions of the first three members of the fatty acid salts (formate, acetate, propionate) show very similar nonidealities to those observed in solutions of equivalent concentrations of a non-surface active 1:1 electrolyte, such as sodium chloride. With increasing number of carbon atoms in the aliphatic chain, the deviations from ideal solution behavior become increasingly pronounced, with a strong non-linear concentration dependence (Smith and Robinson 1942).

Many atmospheric processes, such as in–cloud and below–cloud gas scavenging, gas-to-particle conversion, deliquescence and crystallization of salts, involve non-ideal solution systems (Zuend et al. 2011) and therefore a proper description requires activity coefficients. For atmospheric applications, the activity coefficient model must typically be able to continuously describe solution properties from the early stages of droplet growth, when micellization is likely to occur, to highly diluted systems such as cloud and rain droplets. Several activity coefficient models exist for aqueous solutions of electrolytes mixed with organic compounds, but these are mainly devoted to short-chain carboxylic acids of weak surfactant strength. Raatikainen and Laaksonen (2005) and Tong et al. (2008) presented comparative analyses of these models as variations of the group contribution method UNIFAC or UNIQUAC functional group activity coefficients. Zuend et al. (2011) developed an extended parametrization combining the Pitzer and Kim (1974) electrolyte model with a UNIFAC-based approach for a variety of organic groups. However, to our knowledge, these parametrizations have not been used to estimate the properties of atmospheric surfactants such as the long–chain carboxylic acids or their salts. A different approach was taken by Kim and Frederick (1988a), who applied the Pitzer model to aqueous surfactant solutions and reported interaction parameters for C1–C8 fatty acid sodium salts (formate, propionate, butyrate, valerate, caprylate, pelargonate, caprate). These were obtained by fitting experimental data of osmotic coefficients, but disregarded micelle formation as an important source of non-idealities for especially the larger anions.

Until now, none of the approaches for atmospheric applications have specifically focused on multielectrolyte surfactant–containing systems, where phenomena such as self-aggregation, co-ion pairing, counterion binding and salting-out become particularly relevant. Activity coefficients for aqueous solutions with one ionic surfactant and one inorganic salt have been theoretically predicted by combining the mass-action model for the micellization reaction with classic Debye–Hückel theory (Tajima et al. 1970; Douhéret and Viallard 1982; Burchfield and Woolley 1984; Sharma et al. 2011) and also with more complex excess Gibbs free energy models such has the electrolyte non-random two-liquid model (Chen et al. 2001). Mass–action models require numerical solution for the concentrations of free monomers and aggregates in solution from the equations for the activity coefficients and the micellization reaction equilibrium constant. However, values for micellization equilibrium constants from experimental, modelling or theoretical studies are scarce. The micellization has also been described with the pseudo-phase separation method of Shinoda and Hutchinson (1962). Here, the properties of a solution above the CMC are averaged between those of the “saturated” solution and of a micellar pseudo–phase, using the micellization fraction of surfactant molecules as a weighting factor. Instead of solving the equilibrium constant for the reaction, the micellization fraction is directly obtained from CMC values assuming that all micelles contain the same number of surfactant monomers (molecules or ions), i.e. they have a uniform aggregation number.

In this work, we present an activity coefficient model that uses interaction parameters from binary electrolyte–water systems to predict the properties of aqueous multielectrolyte ionic surfactant solutions including micellization. Binary data is much more abundant in the literature and easy to obtain from experiments, than data for ternary and higher order solutions. The properties of the monomer–micelle equilibria are calculated with the phase-separation method using a composition-dependent parametrization for the CMC as a function of salt molality. From here on, we will refer our model as the CMC based Ionic Surfactant Activity, or CISA, model.

The CISA model represents departures from ideal solution behaviour through the Pitzer–Debye–Hückel (PDH) thermodynamic framework for mixed electrolytes (Pitzer and Kim 1974). The model is used to predict water activity and activity coefficients of ionic solute species in binary and ternary solutions comprising water, surfactant and inorganic salt. We explore the performance of the model from dilute systems to surfactant concentrations as high as ten times the expected CMC. Results from the CISA model are validated against experimental data for binary aqueous solutions of three surfactants of atmospheric relevance (Prisle et al. 2011; Petters and Petters 2016; Forestieri et al. 2018), sodium dodecylsulfate, sodium octanoate and sodium decanoate, and for ternary solutions of sodium decanoate mixed with sodium chloride, which is abundantly present in aerosols and droplets from marine environments. The CISA model is also compared with results from the model by Burchfield and Woolley (1984), which is based on the Guggenheim–Scatchard theory to describe deviation from ideal solution behaviour and numerically simpler, but requires specific model parameters for higher order solutions.

2 Description of the CISA model

2.1 Dissociation and micellization of an ionic surfactant

Here, we derive the equations for the case of aqueous solutions comprising an organic salt (MX) with a surfactant ion (X) and an inorganic salt (MY) sharing the surfactant counterion (M) as a common ion. This type of system is a common surrogate for marine background cloud condensation nuclei (Cochran et al. 2017; Shaharom et al. 2018; Cravigan et al. 2020). The surfactant is here assumed to be anionic, but may also be cationic without changing the overall derivation. The 1:1 organic and inorganic salts dissociate as

$$ \text{MX}\quad \xrightarrow{} \quad \mathrm{M^{+}} \quad + \quad \mathrm{X^{-}} $$
(1a)

and

$$ \text{MY}\quad \xrightarrow{} \quad \mathrm{M^{+}} \quad + \quad \mathrm{Y^{-}}\textrm, $$
(1b)

where X− is the surfactant anion, M+ is the counterion, and Y− is the co-ion.

If the molecular concentration of the surfactant mMX is above the CMC, surfactant anions and counterions associate to form a new ionic species in solution, micelles (labeled Mic), according to the reaction

$$ q\mathrm{M^{+}}\quad + \quad n\mathrm{X^{-}} \quad = \quad \left( \mathrm{M}_{q}\mathrm{X}_{n}\right)^{q-n}\quad = \quad \text{Mic}\textrm, $$
(2)

where n is the aggregation number, or number of surfactant ions per micelle, and q is the number of counterions bound to each micelle, which can be expressed in terms of the counterion binding coefficient β as q = βn (Burchfield and Woolley 1984). Here, β represents the fraction of associated counterions per micelle and for a series of similar compounds vary according to the number of carbon atoms in the hydrophobic part of the surfactant ion. For example, sodium n-alkylcarboxylate surfactants between C6H11O2–C14H25O2 have β values increasing from 0.57 to 0.74 (Vikingstad et al. 1978). Due to the surfactant aggregation, a mole of micelles has an equivalent electrical charge that can be significantly higher than the electrical charge for a mole of the individual surfactant ions. This electrical charge also depends on the degree of counterion binding to the micelle surface and therefore does not have a unique value, but decreases as the extent of counterion binding increases with increasing MY concentrations in the solution (Corrin and Harkins 1947; Vikingstad et al. 1978; Vikingstad and Høiland 1978; Høiland and Vikingstad 1978). With a larger number of M–X contact ion pairs at the micelle surface, the more effective screening by the surface electrical charge lead to lower surfactant areas and higher aggregation numbers (Vlachy et al. 2008; Kunz 2009; Salis and Ninham 2014).

Here, we assume that all micelles in a solution have uniform properties which do not change with increasing surfactant concentration above the CMC. Specificaly, values for β and n are assumed to be constant and all micelles have a uniform electric valence in calculations of the solution ionic strength. Deviations from these assumptions are implicitly accounted for by the composition-dependent terms of the solution excess Gibbs free energy used to derive equations for the activity coefficients.

To estimate the extent of micellization, we define the fraction of surfactant molecules ξ that reside in micelles at given conditions. In the mass–action models, the fraction of micellization can be obtained by solving the activity coefficient model together with the equilibrium constant of the micellization reaction (Eq. 2) defined as

$$ K=\frac{a_{\text{Mic}}}{a_{\mathrm{M}}^{q} a_{\mathrm{X}}^{n}}=\frac{\gamma_{\text{Mic}}}{\gamma_{\mathrm{M}}^{q} \gamma_{\mathrm{X}}^{n}}\frac{m_{\text{Mic}}/m_{\text{Mic}}^{o}}{\left( m_{\mathrm{M}}/m_{\mathrm{M}}^{o}\right)^{q} \left( m_{\mathrm{X}}/m_{\mathrm{X}}^{o}\right)^{n}}\textrm, $$
(3)

where \({m_{i}^{o}}\) is the molality of ionic species i at the reference state of 1 mol kg− 1 where its molality-based activity coefficient γi = 1 and activity ai = 1.

Solving Eq. 3 is a resource consuming computational task due to the highly nonlinear nature of the equations involved. For atmospheric applications, the implementation of an activity coefficient model that includes micellization or ion–surfactant specific interactions has therefore so far not been done. For example, this is a significant limitation when activity coefficients are required for modelling cloud droplet activation with considerations of bulk/surface partitioning, as calculation of the Köhler curve in such cases requires iterative solution of droplet bulk and surface concentrations for all droplet sizes (Prisle et al. 2008; Prisle et al. 2010), unless an approximate parametrization is used (e.g. Raatikainen and Laaksonen, 2011. Ovadnevaite et al. (2017) used the AIOMFAC model to calculate cloud activation of surfactant–salt droplets, but only under the assumption of complete phase separation of water- and organic-rich phases.

Instead of full equilibrium calculations, we here use a pseudo-phase separation approach, where the fraction of micellization ξ is estimated from a concentration-dependent parametrization of the CMC as a function of the salt molality (De Lisi et al. 1981) as

$$ \xi=\left( \frac{\left( \text{sgn} \left( \text{CMC}-m_{\text{MX}}\right)+1\right)}{2}\right)\frac{m_{\text{MX}}-\text{CMC}}{m_{\text{MX}}}, $$
(4)

Here, CMC is the surfactant critical micelle concentration in molal units and sgn is the sign function. Comparing the calculated values of ξ for the aqueous surfactant systems considered here (as shown in Figs. S1 and S2 in the Electronic Supplementary Material), we have not found significant differences between the approaches of Eqs. 3 and 4.

We can simplify the system representation by treating it as a mixture of 1 kg of water were the molality mi is equal to the moles of the different dissociating species \(i=2{\dots } d\) (e.g. i = 2 for the surfactant, i = 3 for the inorganic salt) and \(m_{1}\equiv s= 1/\left (18.0153\times 10^{-3}\right )\) moles is the water molality. Once ξ is found from Eq. 4, the concentrations of molecular and ionic forms of all species in solution can be related as

$$ m_{2}=m_{\mathrm{X}}+nm_{\text{Mic}}=m_{\mathrm{X}}+\xi m_{2}\textrm, $$
(5a)
$$ m_{2}+m_{3}=m_{\mathrm{M}}+qm_{\text{Mic}} =m_{\mathrm{M}}+\xi m_{2}/n\textrm, $$
(5b)
$$ m_{3}=m_{\mathrm{Y}}\textrm, $$
(5c)

and

$$ m_{\text{Mic}}=\xi\frac{m_{2}}{n}\textrm. $$
(5d)

The activity of MX can also be defined via the activities of individual ionic species as

$$ a_{\text{MX}}=\left( a_{\mathrm{M^{+}}}a_{\mathrm{X^{-}}}\right)=\left( \gamma_{\mathrm{M^{+}}}m_{\mathrm{M^{+}}}\right)\left( \gamma_{\mathrm{X^{-}}}m_{\mathrm{X^{-}}}\right)=\left( \gamma_{\mathrm{M^{+}}}\gamma_{\mathrm{X^{-}}}\right)m_{\text{MX}}^{2}\textrm, $$
(6)

where ai and mi are the activity and the molality of species i in the solution, respectively. To simplify Eq. 6 the activity of MX are expressed in terms of the mean ionic activity coefficient \(\gamma ^{\pm }_{\text {MX}}=\left (\gamma _{\mathrm {M^{+}}}\gamma _{\mathrm {X^{-}}}\right )^{\frac {1}{2}}\) (Soustelle 2015) as

$$ a_{\text{MX}}={\gamma^{\pm}_{\text{MX}}}^{2} m_{\text{MX}}^{2} \textrm. $$
(7)

2.2 Solute activity coefficients

The equations for activity coefficients in the CISA model are derived from Pitzer–Debye–Hückel (PDH) theory. Debye–Hückel theory is used to estimate the nonidealities due to long-range electrostatic forces between ions. The effect of short-range forces is represented by a virial expansion in terms of the solution ionic strength. Virial terms represent the pair and triple interactions between ions of different sign. Here, triple interactions of ions with the same charge are disregarded. The interactions between neutral and ionic species can also be included in the virial series (Pitzer 2017). Our current framework does not involve neutral species and the present formulation therefore does not include these interaction terms.

The CISA model was formulated to describe properties of ternary solutions with multiple electrolytes using only model parameters of the corresponding binary systems. It can be extended to calculations of other thermodynamic and volumetric properties of the solution in partial and total forms, such as density, enthalpy and entropy. Unlike other approaches, including that of Burchfield and Woolley (1984), the effect of the short-range forces is not linear with respect to the solution ionic strength. This allows representing the solution nonidealities at high concentrations, nearly up to the solubility limits (Pitzer and Kim 1974; Kim and Frederick 1988a, b).

Model parameters can be fitted to experimental data, but are already available in the literature for more than 300 binary systems of inorganic and organic salts in water (Kim and Frederick 1988a, b). The γi values in aqueous solutions with three or more electrolytes can be predicted using binary parameters with uncertainties under 5% for solutions at moderate ionic strength (I = 2–6 mol kg− 1). The effect of cation–cation–anion or anion–anion–cation triple interactions is directly proportional to the concentration of the ions involved and can be disregarded if the individual concentrations are low (Pitzer and Kim 1974).

The expressions to calculate activity coefficients include the Debye–Hückel limiting law and the second and third virial coefficients BMX and CMX to account for pair interactions of the type MM, MX and XX and triple interactions of the type MMX and XXM, respectively (Pitzer 1973). After neglecting the virial terms for cation–cation–cation or anion–anion–anion interactions, the expressions for the activity coefficients of the ions from the surfactant MX and the inorganic salt MY in a ternary water–MX–MY system are

$$ \begin{array}{@{}rcl@{}} \ln\gamma_{\mathrm{M}}=&F+m_{\mathrm{X}}\left( 2B_{\text{MX}}^{\gamma}+EC_{\text{MX}}^{\gamma}\right)+m_{\mathrm{Y}}\left( 2B_{\text{MY}}^{\gamma}+EC_{\text{MY}}^{\gamma}\right)\\ &+m_{\text{Mic}}\left( 2B_{\text{MMic}}^{\gamma}+EC_{\text{MMic}}^{\gamma}\right)\textrm, \end{array} $$
(8a)
$$ \ln\gamma_{\mathrm{X}}=F+m_{\mathrm{M}}\left( 2B_{\text{MX}}^{\gamma}+EC_{\text{MX}}^{\gamma}\right) \textrm, $$
(8b)
$$ \ln\gamma_{\mathrm{Y^{-}}}=F+m_{\mathrm{M}}\left( 2B_{\text{MY}}^{\gamma}+EC_{\text{MY}}^{\gamma}\right) $$
(8c)

and for micelles

$$ \ln\gamma_{\text{Mic}}^{-}=\delta^{2}\left( q-n\right)^{2}F+m_{\mathrm{M}}\left( 2B_{\text{MMic}}^{\gamma}+EC_{\text{MMic}}^{\gamma}\right)\textrm, $$
(8d)

with

$$ F=f^{\gamma}+m_{\mathrm{M}}\left( m_{\mathrm{X}}B^{\prime}_{\text{MX}}+m_{\mathrm{Y}}B^{\prime}_{\text{MY}}+m_{\text{Mic}}B^{\prime}_{\text{MMic}}\right)\textrm, $$
(8e)
$$ f^{\gamma}=-\frac{A_{\phi}\sqrt{I}}{1+b\sqrt{I}}+\frac{2}{b}\ln\left( 1+b\sqrt{I}\right) $$
(8f)

and

$$ E=m_{\mathrm{X}}+m_{\mathrm{Y}}+\delta\left( n-q\right)m_{\text{Mic}}, $$
(8g)

where δ is the composition-independent shielding factor, a corrective factor introduced by Burchfield and Woolley (1984) to reduce the electric charge of the micelles in solution and modulate its contribution to the total ionic strength of the solutions. The value for the shielding factor can range between 0.4 and 0.7 and is found by fitting the model results to experimental data of osmotic and activity coefficients in aqueous surfactant solutions (Burchfield and Woolley 1984). Without the shielding factor, the charge of a single NaDec-micelle would be -12.48 (in elementary charge units) for an aggregation number of 39 and a counterion binding coefficient of 0.68 (Burchfield and Woolley 1984) at 298 K (Vikingstad et al. 1978). However, better predictions of the osmotic coefficient in aqueous solutions of NaDec have been obtained using a shielding factor of 0.50 that corresponds to a micelle charge of -6.24 (Burchfield and Woolley 1984; Sharma et al. 2011).

The third virial coefficient of the form Cca (“c” and “a” stand for cation and anion, respectively) is composition-independent but the second virial coefficient Bca depends on the ionic strength of the solution as

$$ B_{\text{ca}}=2\beta_{\text{ca}}^{\left( 0\right)}+ \frac{2\beta_{\text{ca}}^{\left( 1\right)}}{{\alpha_{1}^{2}} I}g\left( \alpha_{1} \sqrt{I}\right), $$
(9a)
$$ B^{\prime}_{\text{ca}} =\frac{2\beta_{\text{ca}}^{\left( 1\right)}}{{\alpha_{1}^{2}} I}\left( -1+\left( 1+\alpha_{1}\sqrt{I}\right)\exp{\left( -\alpha_{1}\sqrt{I}\right)}\right), $$
(9b)
$$ C_{\text{ca}}^{\gamma} =\frac{3}{2}C_{\text{ca}}^{\phi} $$
(9c)

with

$$ g\left( x\right)=1-\exp\left( -x\right)\left( 1+x+\frac{x^{2}}{2}\right), $$
(9d)

where α1 is a model parameter equal to 2.0 kg0.5 mol− 0.5 for all electrolyte systems. The model parameters \(\beta _{\text {ca}}^{\left (0\right )}\), \(\beta _{\text {ca}}^{\left (1\right )}\) and \(C_{\text {ca}}^{\phi }\) can be found by fitting experimental data for the osmotic coefficient or mean ionic activity coefficients.

For all solutes, the reference state is the hypothetical one molal solution of the ion, for which the ratio ai(mi)− 1 tends to unity as the molality approaches zero (Pitzer 2017). In this hypothetical solution, the solvent behaves as in its pure form and the electrolyte behaves ideally with an activity coefficient equal to unity. Because the density of electrolyte solutions change with concentration, molality is the preferred concentration unit. This corresponds to the asymmetric convention III by Levine (2008).

In addition to the micellization and counterion binding processes, the CISA model also considers the increasing effect of intermicellar forces on solution non-ideality with increasing surfactant concentration in solution. Via the third virial coefficient \(C_{\text {MMic}}^{\gamma }\) and the equivalent electrical molality E of the solution, the CISA model can capture how interactions of the type Mic–Mic–M and M–M–Mic become dominant because the effective electrical charge of micelles \(\delta \left (n-q\right )\) is much larger than for individual surfactant ions.

As the CISA model is based on PDH theory, it can be fully extended to calculate the properties of aqueous surfactant solutions mixed with different types of inorganic salts. This implementation would require only the binary water–surfactant interaction parameters and a parametrization of the CMC as function of composition in the ternary system.

2.3 Water activity

The water activity in solution can be related to the ion activity coefficients via the osmotic coefficient according to the Gibbs–Duhem equation for thermodynamic consistency (Pitzer 2017). The osmotic coefficient ϕ is defined as the ratio between the observed osmotic pressure πobs for a solution and the expected osmotic pressure πid for an ideal solution of the same composition and temperature. The activity of water in the solution is then calculated as

$$ \phi\equiv \frac{\pi_{\mathrm{{obs}}}}{\pi_{\mathrm{{id}}}}=\frac{\ln a_{1}}{\ln x_{1}}=-\frac{s}{{\sum}_{i>1}m_{i}} \ln a_{1}\textrm. $$
(10)

If the molality mi is used based on the concentration of dissolved solute in molecular form (disregarding dissociation and aggregation), the osmotic coefficient is designated ϕmol, whereas it is designated ϕion if the molality is based on the ionic species in solution (after dissociation and/or aggregation). The values of ϕmol and ϕion are related as

$$ \phi^{\text{mol}}=\frac{m_{\mathrm{X}}+m_{\mathrm{M}}+m_{\mathrm{Y}}+m_{\text{Mic}}}{2m_{2}+2m_{3}}\phi^{\text{ion}}. $$
(11)

The definition of the osmotic coefficient (Eq. 10) can be combined with the Gibbs-Duhem equation (Eq. S4) to obtain an expression for the water activity (in terms of ϕ) based on the nonidealities of ionic solute species (in terms of γi) as

$$ d\left( \left( \phi-1\right)\sum m_{i}\right)=\sum m_{i}d\ln\gamma_{i}. $$
(12)

A full derivation is provided in Eqs. S1–S5 in the Electronic Supplementary Material.

The ionic osmotic coefficient is then calculated as

$$ \begin{array}{@{}rcl@{}} \phi^{\text{ion}}&=&1+\frac{2}{\sum m_{i}}\left( -\frac{A_{\phi}I^{3/2}}{1+b\sqrt{I}}+m_{\mathrm{M}}\left( m_{\mathrm{X}}\left( B_{\text{MX}}^{\phi}+EC_{\text{MX}}^{\phi}\right)+m_{\mathrm{Y}}\left( B_{\text{MY}}^{\phi}+EC_{\text{MY}}^{\phi}\right)\right.\right.\\ &&+\left.\left.m_{\text{Mic}}\left( B_{\text{MMic}}^{\phi}+EC_{\text{MMic}}^{\phi}\right)\right)\right)\textrm. \end{array} $$
(13)

In all expressions, the Debye–Hückel limiting slope Aϕ = Aγ/3 is given as

$$ A_{\gamma}=\left( 2\pi N_{A}\rho_{1}\right)^{\frac{1}{2}}\left( \frac{e^{2}}{4\pi\epsilon_{0}Dk_{B}T}\right)^{\frac{3}{2}}\\ = 1.17374046 \text{kg}^{0.5} \text{mol}^{-0.5} $$
(14)

at 298 K. In Eq. 14, kB is the Boltzmann constant, NA the Avogadro’s number, T the temperature, ρ1 the density of water (Pátek et al. 2009), e the electron charge, D the dielectric constant of water (Fernández et al. 1997) and 𝜖0 the permittivity of vacuum.

2.4 Critical micelle concentration

The CMC of the aqueous surfactant solution is a key variable in the CISA calculation framework. CMC values for binary aqueous systems comprising the selected model atmospheric surfactants are listed in Table 1. These values were calculated with the semi-empirical model of Muller (1993) and parameters fitted to experimental data presented by Blanco et al. (2005) for binary NaOct and NaDec aqueous solutions. The CMC values obtained in this way for NaDec solutions did not match the experimental values provided by Sharma et al. (2011) at the same temperature. Because values of the CMC proved to strongly affect the quality of closure between experimental and modeled activity coefficient values for ternary water–NaDec–NaCl solutions, we opted to use the experimentally obtained binary CMC values of Sharma et al. (2011). CMC values for NaDS in binary aqueous solutions are estimated from curve fitting of the experimental data presented by Mukerjee (1971) using the model of Onori et al. (1994).

Table 1 Physico–chemical properties relevant for this study
Full size table

For ternary solutions comprising water, NaDec and NaCl, we fitted a rational function to the experimental CMC values measured by Sharma et al. (2011). At 298 K, the NaDec CMC values were 0.113, 0.110, 0.107, 0.098, 0.092 and 0.089 mol kg− 1 at NaCl concentrations of 0, 0.0122, 0.0268, 0.042, 0.0613 and 0.089 mol kg− 1, respectively. The composition-dependent function for the CMC of ternary aqueous NaDec–NaCl solutions is thenFootnote 1

$$ \text{CMC}=\frac{0.0335}{0.2941+m_{3}}\textrm, $$
(15)

where m3 is the molality of NaCl. This model tends to a limiting value for the ternary CMC of 0.0053 mol kg− 1 at the water solubility limit of NaCl.

3 Comparison to the activity coefficient model of Wolley (1984)

We compare results from the presented CISA model with those obtained from the Burchfield and Woolley (1984) thermodynamic model, from hereon referred to as the BW model. The calculation framework of the BW model combines representation of long-range interactions with Debye–Hückel theory and short-range interactions with Guggehheim–Scatchard theory. The main difference between the CISA and BW models is that interaction parameters Bca in the BW model are composition-independent as opposed to changing with the solution ionic strength as in CISA. The advantage of the BW model is the numerical simplicity, compared to the CISA model. However, it gives negative deviations from experimental osmotic coefficients in solutions already at moderate solute concentrations above 0.1 M (Guggenheim 1967; Pitzer 1973), where triple interactions contribute to the nonideal solution behaviour.

The expressions from the BW model for the activity coefficients of each ionic species in a ternary system of water, surfactant and an inorganic salt were adapted by Sharma et al. (2011) as

$$ \ln\gamma_{\mathrm{X}}=-\frac{A_{\gamma}I^{\frac{1}{2}}}{1+b_{\mathrm{X}}I^{\frac{1}{2}}}+B_{\text{MX}}m_{\mathrm{M}}\textrm, $$
(16a)
$$ \ln\gamma_{\mathrm{M}}=-\frac{A_{\gamma}I^{\frac{1}{2}}}{1+b_{\mathrm{M}}I^{\frac{1}{2}}}+B_{\text{MX}}m_{\mathrm{X}}+B_{\text{MY}}m_{\mathrm{Y}}+B_{\text{MMic}}m_{\text{Mic}}\textrm, $$
(16b)
$$ \ln\gamma_{\text{Mic}}=-\delta^{2}\left( q-n\right)^{2}\frac{A_{\gamma}I^{\frac{1}{2}}}{1+b_{\text{Mic}}I^{\frac{1}{2}}}+B_{\text{MMic}}m_{\mathrm{M}} $$
(16c)

and

$$ \ln\gamma_{\mathrm{Y}}=-\frac{A_{\gamma}I^{\frac{1}{2}}}{1+b_{\mathrm{Y}}I^{\frac{1}{2}}}+B_{\text{MY}}m_{\mathrm{M}}\textrm. $$
(16d)

Ion size parameters were set by Burchfield and Woolley (1984) to be bX = bM = bY = bMic = b, all equal to 1 mol0.5 kg− 0.5, to assure thermodynamic consistency via the required equality of the mixed partial derivatives of the activity coefficients in terms of concentration.

The expression for the osmotic coefficient in terms of ionic concentrations is

$$ \begin{array}{@{}rcl@{}} \phi^{\text{ion}}&=&1+\frac{-\frac{2}{3}A_{\gamma}I^{\frac{3}{2}}\tau\left( \sqrt{I}\right)}{m_{\mathrm{X}}+m_{\mathrm{M}}+m_{\text{Mic}}+m_{\mathrm{Y}}}\\ &&+\frac{B_{\text{MX}}m_{\mathrm{X}}m_{\mathrm{M}}+B_{\text{MMic}}m_{\mathrm{M}}m_{\text{Mic}}+B_{\text{MY}}m_{\mathrm{M}}m_{\mathrm{Y}}}{m_{\mathrm{X}}+m_{\mathrm{M}}+m_{\text{Mic}}+m_{\mathrm{Y}}}\textrm. \end{array} $$
(17)

The composition-independent model parameters Bca can be found by fitting experimental data for osmotic coefficients in aqueous solutions of MX and MY. Sharma et al. (2011) reported values of BMX = 0 kg mol− 1, BMY = 0 kg mol− 1 and BMMic = 32 kg mol− 1 for aqueous NaDec–NaCl solutions at 298.15 K. The water activity and activity and osmotic coefficients are calculated for the same concentrations in both the BW and CISA models, by keeping the framework for representation of micelle formation as a common point. Concentrations are calculated for aqueous binary solutions using Eq. 4 for the fraction of micellization ξ together with the CMC values given in Table 1, and for aqueous ternary NaDec–NaCl solutions using the parametrization for the CMC given in Eq. 15.

4 Fitting of parameters for the activity coefficient models

To compare the performance of the CISA and BW models for describing activity coefficients in aqueous surfactant solutions, the model parameters were fitted to the same experimental data sets. We used activity coefficients and water activities derived from measurements of the osmotic coefficient. The osmotic coefficient can be measured with good precision at very low concentrations using relatively simple techniques such as freezing point depression (De Lisi et al. 1981; Sharma et al. 2011), the isopiestic method and vapor-pressure osmometry (Smith and Robinson 1942; Widera et al. 2003). Experimental activity coefficients can also be derived from electromotive force measurements of aqueous solutions (Sasaki et al. 1975; Cutler et al. 1978; Kale et al. 1980). The procedure to obtain the activity coefficients from these data is described in detail elsewhere (e.g. Pitzer 2017). A list of the experimental data used to fit model parameters for each of the atmospherically relevant surfactants studied here is given in Table 2.

Table 2 Sources of experimental data for the activity of water and/or the mean activity ionic coefficients of solutes in aqueous surfactant solutions as a function of solute molality
Full size table

Model parameters were found by minimizing the objective function

$$ \frac{{{\sum}_{j}^{N}}\left( \left( \gamma^{\pm}_{\exp}\left( \boldsymbol{x}_{j}\right)-\gamma^{\pm}_{\text{calc}}\left( \boldsymbol{x}_{j}\right)\right)w_{j}\right)^{2}}{N-N_{p}-1}=0\textrm, $$
(18)

where N is the number of experimental values, Np is the number of model parameters, and \(\gamma ^{\pm }_{\exp }\) and \(\gamma ^{\pm }_{\text {calc}}\) are the mean ionic activity coefficients for an electrolyte in an aqueous solution with mole fractions xj, obtained from experimental measurements and calculated from the models, respectively. A weighting factor of wj = 1.5 was used for all data points corresponding to concentrations below the CMC, which improved the quality of the fitting. We tried also with a similar objective function but based on the experimental values for the osmotic coefficients. Although the quality of the fitting was better, these model parameters did not represent the activity coefficient data well. This is likely due to nonlinearities between these variables, where a small change in the value of ϕ is magnified in the values obtained for γ. The model parameters obtained for the CISA and the BW models are reported in Tables 3 and 4, respectively, together with descriptive parameters of the fitting quality.

Table 3 Model parameters used to calculate the activity coefficients of ionic species in aqueous surfactant solutions
Full size table
Table 4 Model parameters BMX, BMMic and \(m_{2}^{\max \limits }\) (in kg mol− 1) for the BW model
Full size table

For binary aqueous surfactant solutions (no inorganic salt added), when activity coefficients were not directly available, experimental values of ϕmol were fitted to rational or polynomial functions describing the ratio of (ϕmol-1) / \(\sqrt {m}\). The fitted functions are shown in Eqs. S13–S16 in the Electronic Supplementary Material and were used to find the mean ionic activity coefficient of the surfactant as

$$ \ln\gamma^{\pm}_{2}=\phi^{\text{mol}}-1+2{\int}_{0}^{m_{2}}\frac{\phi^{\text{mol}}-1}{\sqrt{m}}d\sqrt{m}. $$
(19)

The integral was solved numerically and the mean activity coefficients are reported in Tables S1–S7 in the Electronic Supplementary Material. Because it is very difficult to obtain reliable experimental data for solutions approaching the infinite dilution limmit \(m_{2}\rightarrow 0\), we instead impose the condition, proved by Debye–Hückel theory at this limiting state, \(\phi \rightarrow 1-A_{\phi }\) as \(m\rightarrow 0\) (Pitzer 2017).

5 Results and discussion

5.1 Binary aqueous surfactant solutions

The performance of the CISA model for describing binary aqueous surfactant systems is seen in Fig. 1. Calculated values for the mean ionic activity coefficients of NaOct, NaDec and NaDS are shown in panels a), c) and e), respectively, while the corresponding activity of water as a function of the surfactant mole fraction in their aqueous solutions is shown in panels b), d) and f). For each binary system, two sets of experimental data are shown alongside model calculations, the first is used for validation (yellow triangles) and the second to fit the model parameters (blue circles). CMC values for binary solutions of each surfactant as given in Table 1 are also shown in the respective panels, to highlight the strong effect of micellization on solution nonideality. We have added the superscript “bin” to the binary CMC, to distinguish values from those observed in ternary water–surfactant–salt solutions. For each of the curves describing the surfactant mean ionic activity coefficient and the water activity, a sharp change in the slope with respect to composition is seen at the CMC. This is caused by the onset of micellization, which induces a drastic reduction in the molalities of both free surfactant ions and counterions, due to their association.

Fig. 1
figure 1

Mean ionic activity coefficients of ionic surfactants and water activity in aqueous solutions (a,b) NaOct (c,d) NaDec (e,f) NaDS. Experimental values and calculated values from the CISA model and the BW model. Experimental data for fitting model parameters is shown with blue circles. The CMCbin indicates the critical micelle concentration in the binary water–surfactant system

Full size image

For aqueous solutions of sodium octanoate, Fig. 1a) shows a good agreement between the model predictions and experimental values for the surfactant mean ionic activity coefficient at concentrations both below and above the CMC. Comparing the calculated values using each of the two activity models, there are no significant differences for the surfactant mean ionic activity coefficients, whereas the CISA model shows a better performance predicting the water activity in Fig. 1b). The CISA model closely follows the trend in experimental values of a1 at surfactant concentrations as high as three times the CMC, x2 < 3CMC. In contrast, the BW model is overestimating the water activity at surfactant concentrations immediately above the CMC.

In the case of sodium decanoate solutions, we see in Figs. 1c) and 1d) that the calculated values of \(\gamma ^{\pm }_{2}\) and a1 from both models follow closely the experimental values across the full concentration range. Predictions deviate between the models at very high surfactant concentrations (x2 > 12CMC), where transitions to complex tridimensional structures may occur (Pfrang et al. 2017), which are not described by our current framework. As it was observed also for NaOct solutions, the CISA model predicts the reduction of water activity with increasing surfactant concentrations, while the BW model tends to overestimate the water activity when surfactant concentrations approach the CMC and with larger positive biases above the CMC.

For aqueous solutions of sodium dodecyl sulfate, the situation is very different. As can be seen in Fig.1e), both models overestimate the surfactant mean ionic activity coefficient at concentrations below the CMC, especially in the limit of highly diluted solutions. Unfortunately, we could not find additional experimental values for the water activity in this concentration range in order to verify potential biases. When the surfactant concentration exceeds the CMC, both models show good performance in predicting the surfactant mean ionic coefficient, with no significant mutual differences in the calculated values. However, as seen in Fig. 1f), the CISA model also performs well for calculations of the water activity, while the BW model underestimates these values, when x2 > CMC.

To assure thermodynamic consistency in terms of the Gibbs–Duhem relation, ∂γi/∂mj = ∂γj/∂mi, the BW model assumes a unique ion-size parameter b equal to unity for all ions in the solution. We tried to improve performance of the BW model by fitting different values for b to describe the ions in aqueous SDS solutions. However, there was no significant improvement of either the quality of fitting or on the model predictions. As we could not find an alternative set of experimental data for activity coefficients in dilute SDS solutions, we could not confirm the experimental values of Tajima et al. (1970) or verify potential model biases from them.

In general, the CISA model gives better predictions than the BW model of the binary surfactant solution water activity at any concentration, both below and above the CMC, while there are not significant differences between the models for calculated surfactant mean ionic activity coefficients. Even when both CISA and BW models use the same parameters to represent the surfactant aggregation and counterion binding, the CISA model better captures the effect of triple solute interactions at high surfactant concentrations. This model uses the third virial coefficient Cca term to account for triple interactions of the type MMX and XXM that lead to the surfactant aggregation. It also accounts for the MMMic and MMicMic triple interactions (Eq. 8d) above the CMC, which have a profound effect on the solution behaviour. Inter-micelle forces are strong because micelles are multiply charged ionic species and not only raise the equivalent electrical molality and the ionic strength of the solution, but also alter the mobility of all ions in solution. These effects combined result in strong negative deviations from ideal solution behaviour and surfactant mean ionic activity coefficients that are consistently below unity.

It is clear that the surfactant aggregation cannot be disregarded in this type of calculations, especially in atmospheric applications where the water activity is used to estimate solution equilibrium with water in the gas-phase, or surfactant activity coefficients are used to estimate the surface tension and composition of aqueous solution–air interfaces. Results for calculations of mean activity coefficient and water activity in the studied systems without consideration of micellization are shown in Fig. S3 of the Electronic Supplementary Material.

5.2 Ternary aqueous solutions

We then applied the CISA model to predict nonideal behavior of ternary NaDec–NaCl aqueous solutions using only the ion–ion interaction parameters βca and Cca obtained for the corresponding binary systems. As described above, water–surfactant parameters were obtained from our fittings (as listed in Table 3) and water–inorganic salt parameters from literature. Results from the CISA model are compared with those from the BW model as presented by Sharma et al. (2011) and using model parameters fitted to experimental activity coefficient data for ternary aqueous NaDec–NaCl solutions.

Figure 2 shows the mean ionic activity coefficients \(\gamma ^{\pm }_{2}\) for NaDec and \(\gamma ^{\pm }_{3}\) for NaCl in aqueous solutions with different relative mixing ratios of surfactant and inorganic salt solutes, here expressed as \(x_{3}^{\text {dry}}=m_{3}/(m_{2}+m_{3})\). Values of activity coefficients calculated with the CISA model are compared to those of the BW model and to the experimental values of (Sharma et al. 2011). In panels a), c), e) and g), we also include the NaDec mean ionic activity coefficient in a binary solution with equivalent concentration and similarly, in the panels b), d), f) and h), the binary NaCl mean ionic activity coefficient. This allows us to highlight the changes in the mean ionic activity coefficients caused by solute–solute interactions specific to the ternary solutions. There are no significant differences between the calculated \(\gamma ^{\pm }_{2}\) and \(\gamma ^{\pm }_{3}\) obtained from the CISA model, the BW model and experimentally derived values. This is highly encouraging, because calculations with the CISA model are based on much the more readily available parameters from corresponding binary aqueous systems, while the BW model requires parameters from ternary solution data to reach the same degree of agreement with experimental values.

Fig. 2
figure 2

Calculated mean ionic activity coefficients of NaDec and NaCl in ternary aqueous solutions at different NaDec–NaCl relative mixing ratios, compared to experimental values from Sharma et al. (2011), at 298 K. The CMCter indicates the critical micelle concentration in the ternary water–surfactant–salt system given by Eq. 15

Full size image

In each panel of Fig. 2, the critical micelle concentration for the corresponding \(x_{3}^{\text {dry}}\) is indicated, with the superscript “ter” added to distinguish values from the CMC observed in binary water–NaDec solutions. The CMCter decreases with increasing inorganic salt molality according to Eq. 15, as discussed in more detail in connection with Fig. 4 below. The indicated CMCter therefore corresponds to the first time where the condition x2 > x2,CMC is fulfilled for each \(x_{3}^{\text {dry}}\). As it was also observed for the binary systems, the slope of the mean activity ionic coefficient curves changes significantly at the indicated CMCter.

In all panels of Fig. 2 we notice how both models follow the trend of the experimental data both below and above the CMC. We also notice that when the surfactant–salt mixing ratio decreases (increasing \(x_{3}^{\text {dry}}\)), the calculated \(\gamma ^{\pm }_{3}\) in the ternary system (panels b, d, f and h of Fig. 2) deviate more from those of the binary water–NaCl system of equivalent molality. A similar tendency is also observed for \(\gamma ^{\pm }_{2}\) (panels a, c, e and g of Fig. 2), but most significantly at the lowest surfactant–salt mixing ratios (highest \(x_{3}^{\text {dry}}\)) studied. This suggests that the strongest departures from ideal solution behavior for both species occur when the organic and inorganic electrolytes are in approximately equal molar amounts in solution.

Figure 3 shows how the calculated water activity (panels a, c, e and g) and the osmotic coefficient (panels b, d, f and h) from the CISA model follow the experimental values with a good performance across the whole concentration range. Again, the CMCter for each \(x_{3}^{\text {dry}}\) is indicated and it is clear how especially the osmotic coefficient changes sharply as function of NaDec concentration at the onset of micellization.

Fig. 3
figure 3

Calculated activity and osmotic coefficients of water in ternary NaDec–NaCl aqueous solutions at different NaDec–NaCl relative mixing ratios, compared to experimental values from Sharma et al. (2011), at 298 K. The CMCter indicates the critical micelle concentration in the ternary water–surfactant–salt system given by Eq. 15

Full size image

In general, the CISA model describes nonidealities in ternary aqueous NaDec–NaCl solutions well, using just the parameters for the corresponding binary aqueous systems. This represents a great advantage for applications to more complex systems, such as atmospheric solution droplets, as the calculation framework of the CISA model can be extended to predict activity coefficients of more general multiple-electrolyte systems. Many binary solution interaction parameters, in particular for well-known inorganic salts and atmospheric surfactants, are already available in the literature or may alternatively be obtained experimentally from relatively simple measurements of the solution osmotic coefficient and electrical conductivity. The main limitation of the CISA framework lies in the construction of robust composition-dependent functions for the CMC in solutions of multiple surfactants and electrolytes, because surfactant aggregation can here occur in multiple parallel reactions with mixed-surfactant micelle formation (MacNeil et al. 2014) and because experimental CMC values can vary significantly, depending on the property and the technique used to determine them (Álvarez-Silva et al. 2010).

We have shown that with a parametrization of the surfactant CMC as function of the inorganic salt content in the solution, the CISA model can be used to estimate the water activity and mean ionic coefficients of both surfactant and inorganic salts in solutions spanning concentrations from highly dilute conditions up to several times the surfactant CMC. Both water activity and surface tension (as a function of surface or bulk surfactant activity, depending on the equation of state) play important roles in controlling the solution–air interface equilibrium of atmospheric droplets (Malila and Prisle 2018; Lin et al. 2018; Bzdek et al. 2020). Their proper estimation is therefore essential for predicting the formation of clouds and their radiative effects in the atmosphere (Prisle et al. 2012). When droplet solutions contain surfactants, their adsorption at the droplet–air interface must be explicitly evaluated from adsorption models that require the activity coefficients as inputs, especially if surfactants are mixed with salts where ion–ion interactions modify the surfactant aggregation and counterion binding and with that the adsorption isotherm (Sorjamaa et al. 2004; Prisle et al. 2010; Ovadnevaite et al. 2017; Prisle and Molgaard 2018).

5.3 Activity coefficients in ternary solutions of high ionic strength

Figure 4 shows the variation of the mean ionic activity coefficients of NaDec and NaCl in their ternary aqueous mixtures, for higher relative amounts of the inorganic salt (\(x_{3}^{\text {dry}}>0.4\)) than currently reported by experiments (Smith and Robinson 1942; Sharma et al. 2011). In each case, the variation of the surfactant CMC with the mole fraction concentration of the inorganic salt x3 in the ternary solutions is indicated with black curves. The CMC decreases with increasing NaCl molalities in ternary solutions and defines a local change in the derivative of the surfactant activity coefficient with respect to composition (x2 or x3).

Fig. 4
figure 4

Mean ionic activity coefficient of the surfactant and the salt in NaDec–NaCl solutions at 298 K at high salt–surfactant mixing ratios. CMCter values are depicted by the black curve. Note that the color of the surface lines changes with increasing NaCl mole fractions in the solution

Full size image

The variation of both the surfactant and inorganic salt activity coefficients, as well as the CMC, is more pronounced at the higher value of \(x_{3}^{\text {dry}}\) (Figs. 4c and 4d), due to stronger ion–ion interactions in solutions of higher ionic strength. However, while the mean ionic activity coefficient of the inorganic salt varies smoothly around unity, the surfactant mean ionic activity coefficient increases steeply by several orders of magnitude for relatively small increments in NaCl concentration. Such asymmetric behaviour and strong non-idealities are typically associated with miscibility gaps (Van Ness et al. 2000), which may manifest as salt-induced liquid–liquid phase separation (Imae et al. 1988) or in a form of surfactant precipitation, depending on the total concentration of the ternary water–surfactant–salt system. For example, when the concentration of the surfactant counterion is increased by addition of the inorganic salt sharing the common counterion, the ionic surfactant can precipitate if the activity product of anion and its counterion exceeds the solubility product (Scamehorn and Harwell 2005). The CISA model could be extended to represent phase separation in ternary systems at specific surfactant–salt mixing ratios, by establishing the required conditions for spinodal decomposition in terms of the second derivatives of the Gibbs free energy (Zuend et al. 2010).

In all panels of Fig. 4, predictions of activity coefficients with the CISA model are continuous across a concentration range above the CMC up to surfactant molalities of 0.5–0.7 mol kg− 1. This is a crucial feature of the CISA model for atmospheric applications, in order to continuously describe aqueous solutions spanning a wide concentration range representative of growing droplets in the atmosphere (Li et al. 1998; Prisle et al. 2008; Prisle et al. 2010; Lin et al. 2018; Lin et al. 2020; Bzdek et al. 2020). From the early stages of hygroscopic growth until the point of cloud activation, activity coefficients are essential to fully describe the dilution of solid phases, micellization and other changes, such as the packing factor and shape of the surfactant aggregates that have been observed in experiments for surfactant-containing atmospheric droplet model systems (Pfrang et al. 2017). In an extension of the CISA model to account for such phenomena, the steep increase in the surfactant activity coefficients at certain compositions, as seen in Fig. 4, could be used as criterion to switch between models that include/exclude different surfactant aggregates from equilibrium surface tension equations. For example, (molality based) activity coefficient values higher than 100 could be used to indicate the possible existence of micelles that are not accounted for by typical adsorption isotherms (cf. Table 1 in Malila and Prisle, 2018).

5.4 The CISA model at different temperatures

To perform calculations with the CISA model at different solution temperatures, parametrizations of the model parameters for binary water–surfactant and water–inorganic salt systems, as well as the CMC in binary and ternary surfactant solutions at the corresponding temperatures are needed. While osmotic coefficients for aqueous solutions at different temperatures are available for many inorganic electrolytes, this type of experimental data for aqueous surfactant solutions of atmospheric relevance is scarce. We found a single set of temperature-dependent data for osmotic coefficients in binary aqueous solutions of NaDec presented by De Lisi et al. (1981). The CISA model interaction parameters \(\beta ^{0}_{\text {NaDec}}\), \(\beta ^{1}_{\text {NaDec}}\) and \(C^{\phi }_{\text {NaDec}}\), \(\beta ^{0}_{\text {NaMic}}\), \(\beta ^{1}_{\text {NaMic}}\) and \(C^{\phi }_{\text {NaMic}}\) were fitted at temperatures 278.15 K and 288.15 K, in addition to those at 298.15 K already used in the main part of this work. Results for the fits and implementation of the CISA model to calculate activity coefficients and water activities in binary aqueous NaDec solutions at these temperatures can be found in the Electronic Supplementary Material (Figs. S4–S6).

Binary solution NaDec and NaMic pair interaction parameters were then extended to 273.15 K by linear extrapolation and used in the CISA model to calculate NaDec and NaCl mean ionic activity coefficients and water activity and osmotic coefficient in ternary aqueous solutions with different NaDec–NaCl mixing ratios. Results are shown in Supplementary Figs. S7 and S8, together with experimental results for ternary NaDec–NaCl solutions at 273.15 K from Sharma et al. (2011). The corresponding binary water–NaCl interaction parameters \(\beta ^{0}_{\text {NaCl}}\), \(\beta ^{1}_{\text {NaCl}}\) and \(C^{\phi }_{\text {NaCl}}\) at 273.15 K used for the calculations were obtained from temperature-dependent functions built using the values reported by Pitzer et al. (1984) for aqueous NaCl solutions at 1 bar between 273 K and 313 K.

The CISA model closely follows the experimental NaDec activity coefficients and ternary solution water activity at all the surfactant–inorganic salt mixing ratios, but tends to overestimate the salt mean ionic activity coefficient at concentrations above the CMC. Without additional experimental data, it is difficult to assess whether these discrepancies at lower temperature indicate the existence of stronger ion–micelle interactions in ternary water–surfactant–salt solutions, or reflect an underestimation of the NaMic pair interactions fitted from binary water–NaDec solutions. We also note that the experimental data of Sharma et al. (2011) for osmotic coefficients was measured at the solution freezing temperatures (which deviate from the water freezing due to freezing point depression) and then converted to values at 273 K, which may also introduce a source of variance.

5.5 Extensions of the CISA model

Application of the CISA model to ternary solutions of an ionic surfactant (MX) mixed with an inorganic salt (NY), which do not share a common ion (M+≠N+), requires a description of the CMC and micelle composition, aggregation number, and degree of counterion binding. These variables can however not be readily predicted using ideal mixing rules for the binary aqueous systems.

The main challenge lies in the ion specificity of ion pairing and counterion binding affecting the surfactant adsorption and micellization processes. In ternary aqueous solutions, dissimilar surfactant and inorganic salt counterions compete for adsorption sites at the solution–air interface or micellar surfaces, depending on the degree of coion pairing with the surfactant hydrophilic part. The micellization occurs as a competitive process between mixed counterion adsorption at the micellar surface and micelle compositions and sizes spanning the range from those observed in binary MX and NX aqueous solutions, respectively. The law of matching water affinities establishes that “inner sphere ion pairs are preferentially formed between oppositely charged ions with matching absolute enthalpies of hydration” (Collins 1997). In ternary systems, ions that forms the closest ion-pair with the surfactant hydrophilic part will dominate the counterion binding to micelles (Kim et al. 2001; Moreira and Firoozabadi 2010). Between M–X and N–X, the closest ion-pair will also lead to the largest reduction in the surfactant partial molar area, increasing micelle aggregation numbers and decreasing the CMC in the ternary solution (Vlachy et al. 2008). However, the extent of such changes ultimately depends on the counterion mixing ratio and total electrolyte content in the ternary solution.

Effects of ion specificity on surface activity in ternary MX–NY solutions poses a significant challenge to atmospheric calculations, where droplet solutions can contain dozens of inorganic and organic ions. To our knowledge, these effects have so far not been taken into account in models for atmospheric droplet systems.

There are several molecular thermodynamic theories to account for ion-specific effects on micellization (Warszyński et al. 2002; Moreira and Firoozabadi 2010; Danov et al. 2014; Lukanov and Firoozabadi 2014; Zhu and Free 2015; Hao et al. 2016; Khoshnood et al. 2016; Koroleva et al. 2019), but their numerical implementation is complex and computationally prohibitive to atmospheric models. As a way forward, a simple approach to estimate the aggregation number n, thermodynamic potentials and equilibrium constant of the micellization reaction could be based e.g. on microcalorimetric (Olesen et al. 2015; Durand-Vidal et al. 2020) and conductance (Jocić 1995) measurements. These would enable estimates of the CMC and average degree of counterion binding in series of ternary solutions with different M+ −N+ ion mixing ratios. Together with the mass–action model of micellization, the CMC model of Karakashev and Smoukov (2017) and the activity coefficient equations of the CISA model, these parameters would constrain the extended framework. We further anticipate that observed deviations from mixing rules based on micellization parameters for binary solutions can be used to identify ion–ion synergistic effects and guide an intelligent inference system (fuzzy logic) for extending the CISA framework.

6 Conclusions

We present the CISA model to estimate ionic activity coefficients and water activity in aqueous surfactant solutions with or without added inorganic salt. The model can be implemented to describe the properties of ternary ionic surfactant–inorganic salt solutions using only model parameters obtained from the corresponding binary water–surfactant and water–salt systems. The concentration of all solute species can be solved without iteration using composition-dependent functions for the surfactant CMC as a function of the inorganic salt molality. This flexibility adds significant potential for application to complex aqueous surfactant mixtures of atmospheric relevance. For the investigated ternary sodium decanoate–sodium chloride aqueous solutions, a common model system for marine atmospheric droplets, the CISA model is able to capture solution non-idealities induced by micellization at surfactant concentrations up to 3 times the CMC.

When compared to the model of Burchfield and Woolley (1984), the CISA model gives an improved description of the activity of water above the CMC, for both binary and ternary surfactant systems. This is particularly important when describing early stages of droplet growth for surfactant–enriched aerosols, as observed in both laboratory and ambient atmospheric settings. Because the CISA model is based on Pitzer–Debye–Hückel theory, it can be extended to estimate other thermodynamic and volumetric properties of aqueous electrolyte solutions, including one or more inorganic salts. The PDH theory is very well documented and interaction parameters are available for more than three hundred different binary aqueous electrolyte solutions. PDH equations have been used in other frameworks for atmospheric modeling purposes, including E–AIM (Clegg et al. 2001; Clegg et al. 2003) and AIOMFAC (Zuend et al. 2008; Zuend and Seinfeld 2012; 2013), to describe multi-component mixtures of both symmetrical and unsymmetrical electrolytes with organic compounds. The CISA model can be used to extend the capabilities of these frameworks to describe properties of atmospheric surfactant-enriched aerosols and droplets.