Classical nucleation theory of immersion freezing: Sensitivity of contact angle schemes to thermodynamic and kinetic parameters

Heterogeneous ice formation by immersion freezing in mixed-phase clouds can be parameterized in general circulation models (GCMs) by Classical Nucleation Theory (CNT). CNT parameterization schemes describe ::::::::: immersion freezing as a stochastic process including the properties of insoluble aerosol particles , so called ice nuclei, in the droplets. There are different ways how to describe :::::::::: parameterize the properties of aerosol particles (i.e. contact angle schemes), which are compiled and tested in this paper. The goal of this study is to find a parameterization scheme for GCMs to describe immersion freezing 5 with the ability to shift and adjust the slope of the freezing curve compared to homogeneous freezing to match experimental data. The results of using :: We ::::::: showed :: in : a :::::::: previous ::::::::: publication :::: that ::: the ::::::: resulting ::::::: freezing :::::: curves ::::: from CNT are very sensitive to unconstrained kinetic and thermodynamic parameters in the case of homogeneous freezingleading to uncertainties in calculated nucleation rates Jhom of several orders of magnitude. Here we investigate how sensitive the outcome of a parameter estimation 10 for contact angle schemes from experimental data is to ::::::::::: unconstrained : kinetic and thermodynamic parameters. We show that additional free parameter :::::::::: demonstrate ::: that ::: the :::::::::: parameters ::::::::: describing ::: the :::::: contact ::::: angle :::::::: schemes can mask the uncertainty of Jimm due to:: in thermodynamic and kinetic parameters. Different CNT formulations are fitted to an extensive immersion freezing dataset :::::::: consisting :: of :::: size ::::::: selected :::::::::::: measurements as a function of particle diameter (d), temperature T and time t :::::::::: temperature ::: and :::: time : for different mineral dust types, namely 15 kaolinite, illite, montmorillonite, microcline (K-feldspar) and Arizona test dust. It is investigated how accurate different CNT formulations (with the estimated fit parameters :: for :::::::: different :::::: contact ::::: angle :::::::: schemes) reproduce the measured freezing curves :::: data, especially the time and particle size dependence of the freezing process. The results are compared to a simplified deterministic freezing scheme. It is evaluated in this context : In :::: this :::::: context :: it :: is :::::::: evaluated which CNT based parameterization scheme :::: able to represent particle properties is a good :: the :::: best : choice to describe immersion freezing in a GCM. 20


Introduction
In mixed-phase clouds freezing of cloud droplets occurs by different pathways of heterogeneous freezing/nucleation. The 20 nucleation process is initiated on the surface of an aerosol particle, called ice nucleus (IN), which either collides with a supercooled droplet (contact freezing), acts as cloud condensation nucleus (CCN) and causes freezing when the droplet is increasingly supercooled (immersion freezing), freezes immediately after CCN activation at supercooled conditions (condensation freezing), or provides a site where water vapor deposits as ice (deposition nucleation) [Vali (1985)].
In mid latitudes, where supercooled clouds are common, IN and their effect on precipitation formation through immersion freezing influence the hydrological cycle (Lohmann, 2002;Zeng et al., 2009;DeMott et al., 2010) and thereby e.g. the biosphere and agriculture. Aerosol particles determine the formation and ice-water ratio of mixed-phase clouds, thereby the cloud radiative properties and indirectly the radiation budget, which affects earth's climate. Therefore results of climate simulations 5 in regional and global models are sensitive to the parameterization scheme used for heterogeneous ice formation and in particular immersion freezing as it is the most abundant freezing pathway (Ansmann et al., 2009;Wiacek et al., 2010). One approach to parameterize immersion freezing in global and regional climate models is by Classical Nucleation theory (CNT), which requires approximations of the thermodynamics and kinetics of nucleation. Although computationally more expensive compared to empirical parameterization schemes, it allows a physical treatment of ice nucleation as function of temperature T , ice 10 supersaturation S i , time t and IN type (e.g. size, surface properties). Using a theoretical scheme has the advantage that the scheme is valid over the whole T -S i -space, which is mandatory for the use in a GCM, where all kinds of conditions occur (especially in certain regions like the Arctic, but also in simulations of future climate, where atmospheric conditions can be different from the present day or the pre-industrial epoch). Empirical schemes are often limited to the narrow conditions from which the scheme was derived and can lead to unphysical results when extrapolated. Therefore empirical schemes might not be 15 representative for the future atmosphere or untypical atmospheric conditions. One example is the empirical Meyers et al. (1992) scheme, which was developed using measurements in mid latitudes and has been found to be inaccurate when extrapolated to Arctic conditions (Prenni et al., 2007).
Some parameters in the framework of CNT are so far unconstrained. At the same time results from the CNT are very sensitive to the choice of thermodynamic and kinetic parameters, in particular interfacial tension between ice and water σ iw 20 and activation energy ∆g # . Sensitivity of CNT on σ iw and ∆g # in the case of homogeneous freezing has been discussed in Ickes et al. (2015). Using CNT as an approach to parameterize immersion freezing in aerosol-climate models raises the question of the sensitivity of the parameterization scheme to σ iw and ∆g # in the case of heterogeneous freezing. Additionally there is a need to include and represent IN properties. We test three different schemes to describe the effect of an IN population on immersion freezing conditions and investigate the impact of the chosen scheme on the parameterization of immersion freezing. 25 We also discuss strategies how to evaluate the applicability of different CNT formulations.
The formalism of CNT for immersion freezing is summarized in section 2. Advantages and disadvantages of certain formulations for the use in GCMs are discussed. In section 3 the sensitivity of the immersion freezing nucleation rate J imm [s −1 m −2 ] and the fit of the geometric term f to thermodynamic and kinetic parameters is investigated by fitting and comparing the results to an ice nucleation measurement dataset of kaolinite (Welti et al., 2012). The section is followed by suggestions for criteria 30 how to evaluate the quality of a CNT parameterization scheme (section 4). Section 5 presents CNT parameters estimated from experimental data for five different mineral dust types and in section 5.1 the criteria from section 4 are applied for three CNT formulations and compared to the performance of an empirical parameterization. Throughout this paper we refer to "CNT formulations" indicating a specific CNT framework for immersion freezing based on Eq. 4, including parameterization schemes for σ iw , ∆g # and the geometric term f (contact angle scheme). CNT formulation #1 as an example is later on abbreviated by CNT #1. We refer to "schemes" to discuss the different parameterization schemes for the geometric term f or the contact angle α to express the ice nucleating surface properties of aerosol particles.

Classical Nucleation Theory for immersion freezing
In CNT freezing is described as a stochastic process by a temperature dependent nucleation rate. In the case of homogeneous freezing (pure water droplets) statistical fluctuation of water molecules can lead to the formation of small ice-like structures 5 (ice embryos) that lead to freezing of the supercooled droplet if they reach a certain critical size (ice germ). The nucleation rate describes the formation of ice germs leading to freezing over time. It consists of a thermodynamic and a kinetic component.
The thermodynamic component describes the formation of ice embryos (determined by the thermodynamic energy barrier ∆G), the kinetic component describes the number of molecules, which can be incorporated into the ice embryo (determined by the activation energy barrier ∆g # ). 10 The presence of IN immersed in supercooled droplets facilitates ice nucleation compared to homogeneous nucleation by providing a catalytic surface. The IN surface reduces the thermodynamic energy barrier ∆G determined by T , S i and σ iw .
The difference in nucleation with and without an IN i.e. homogeneous or heterogeneous nucleation, is accounted for by the geometric term f , also called wettening factor, compatibility factor or contact parameter. This term indicates the increased probability to nucleate a stable ice germ due to the presence of the IN surface because of the reduced number of water molecules 15 necessary to form an ice germ. f describes by how much the IN properties (of unknown nature) reduce the energy barrier for the formation of ice embryos on its surface compared to homogeneous freezing: where v ice is the volume of a water molecule in the ice embryo, k B the Boltzmann constant, S i the saturation ratio with respect 20 to ice, r IN the radius of the catalytic IN surface and r germ the critical radius of an stable ice cluster which initiates freezing of the droplet.
The contact angle α has a value between 0 • and 180 • , where the latter is equal to the case of homogeneous freezing (f =1 → ∆G = ∆G hom ).
If the radius of the IN is significantly larger than the radius of the ice germ, curvature of the IN surface can be neglected leading 30 to a simplified form of f (Volmer, 1939): f (α) = (2 + cos α)(1 − cos α) 2 4 . (3) Whereas the thermodynamic term in the nucleation rate J imm (thermodynamic exponent determined by the energy barrier ∆G, see above) changes from homogeneous to heterogeneous freezing, the kinetic term is assumed to be the same for homogeneous and immersion freezing. The prefactor of the nucleation rate changes as well. It is defined as: where n s is the number of water molecules in contact with the unit area of the ice cluster, Z the non-equilibrium Zeldovich factor, h the Planck's constant and N l the volume-based number density of water molecules in the liquid parent phase. Z is not the same for homogeneous and heterogeneous freezing, because the number of the water molecules in the ice germ, n k,germ , differs. As shown in Pruppacher and Klett (2000) most of the prefactors cancel out in the case of heterogeneous freezing 10 leading to the following expression for the nucleation rate for immersion freezing: Since the energy barrier of immersion freezing is reduced compared to homogeneous freezing, J imm is higher compared to the homogeneous nucleation rate J hom and exhibits a different T dependence (less steep at a certain T ). The shift of the freezing curve is described by the geometric term f and has to be captured by the different CNT formulations. to describe the ice nucleating surface properties of aerosol particles, e.g. Marcolli et al. (2007). Depending on the scheme the IN properties are represented by functions with one or several fit parameters and the complexity for an implementation in a GCM differs accordingly. Note that increasing complexity comes with higher computational costs (Hurrell et al., 2009). Based on the 20 computational complexity we chose three schemes including one or two fit parameters for the following sensitivity analysis (section 3). They are briefly explained here followed by a paragraph about computational costs. A graphical representation of each scheme is shown in Fig. 1. For more details see Marcolli et al. (2007) and Lüönd et al. (2010).
From immersion freezing measurements the frozen fraction FF is obtained, which is the fraction of a droplet population that is frozen at a certain temperature T after a certain time t. To compare CNT based parameterization schemes to measurements, 25 FF is calculated from the nucleation rate J imm . The frozen fraction FF is given by: with A IN being the surface area of the IN. For simplicity particles are assumed to be spherical (A IN = 4πr 2 IN ). Thus, the surface used for the IN of a specific mass represents a lower limit (non-spherical particles would be larger).

Single-α scheme
The single-α scheme assigns one contact angle to the entire surface of all particles. It is based on the assumption that all particles have one uniform surface property responsible for their ice nucleating ability. Consequently all particles have an equal probability to act as IN at given conditions. The scheme requires only one fit parameter (f or α). One variation of the single-α scheme, the single-α + fit ∆g # scheme (see section 3), requires two fit parameters (f or α and ∆g # ).

5
It is the least complex and consequently the cheapest scheme suitable for implementation in GCMs. However, it does not take into account that ice nucleating properties might be variable throughout a particle population. This scheme is used in several models, e.g. Khvorostyanov and Curry (2000, 2004, 2005; Liu et al. (2007);Eidhammer et al. (2009);Hoose et al. (2010); Storelvmo et al. (2011);Ervens and Feingold (2012).
The α-pdf scheme is an extension of the single-α scheme. It assumes a heterogeneity of particles in an aerosol population by using a log-normal probability density function (pdf) for the contact angle α p(α). The log-normal distribution of α within a particle population is expressed by two fit parameters, the logarithmic mean contact angle µ and the variance σ of the distribution: 15 This approach attributes an individual surface property to each particle on the entire particle surface.
The variance σ accounts for the heterogeneity of the particle property within the aerosol population: the larger the variance σ, the larger the heterogeneity among the particles. The frozen fraction FF is derived by integrating the contact angle distribution over all possible contact angles: 20 The approach has been frequently used to interpret freezing data, e.g. Marcolli et al. (2007) issue is closely connected to the time resolution of the GCM, which will be discussed in a future publication.
Another extension of the single-α and frequently used scheme is the active sites scheme, e.g. in Marcolli et al. (2007); 30 Lüönd et al. (2010) The α(T ) scheme is a compromise between the single-α and the α-pdf scheme. It assumes that α is different for different T , which reflects a change of the activated fraction of the α-pdf distribution and with that a change in apparent contact angle.
Efficient IN (with small contact angles) freeze at highest temperatures. The lower the temperature, the higher the chance that particles with larger contact angles can be activated. This leads to a change of the apparent contact angle to larger α with further cooling due to the shift of the activated fraction of the α-pdf distribution with supercooling. The α(T ) scheme is 10 thus representing the shifted mean contact angle of an initial contact angle distribution, but does not take into account how contact angles are distributed among a particle population at a certain temperature. The temperature dependence of α can be approximated to be linear as discussed in Welti et al. (2012). Being capable to describe a variability of the freezing process due to a contact angle distribution without being computationally complex makes the α(T ) scheme attractive for GCMs. However it demands an indirect assumption on how the aerosol population changes with supercooling. It does not circumvent the issue 15 of shifting the contact angle distribution with time.
The frozen fraction FF is estimated analogously to the single-α scheme using a linear function for α(T ): with α(T ) = α 0 + m · (T − 273.15 K) .

20
Note that in general the temperature dependence of α can be interpreted either as a result of the temperature dependence of the interfacial tensions (σ is and σ iw ) or as the apparent contact angle of an ensemble with a diversity of contact angles from particle to particle. In contrast to Zobrist et al. (2007); Alpert et al. (2011); Knopf and Forrester (2011); Rigg et al. (2013) we follow the second interpretation (simplified temperature dependent α-pdf scheme) as explained above. Accounting for a physical dependence of α on T as a result of the temperature dependence of the interfacial tensions leads to a decrease of α, 25 which is contradictory to the assumption we made here.

Computational costs
Since computational costs are strongly linked to the complexity of parameterization schemes (Hurrell et al., 2009), the complexity of a parameterization scheme is an important factor which has to be considered before implementation into a GCM. A general quantification of the computational costs of parameterization schemes depends an many aspects (treatment of aerosol 30 particles and cloud microphysics) of the GCM used. The number and kind of variables needed, e.g. T , size and number of a certain aerosol type, to derive J imm is an indicator for the complexity of the parameterization scheme. As soon as a variable needs to be memorized over several timesteps, e.g. change in contact angle distribution, an extra-tracer is required, which might be computationally expensive.
The cheapest CNT based parameterization schemes are the single-α scheme and the α(T ) scheme. The α-pdf scheme is computationally more expensive because the contact angle distribution changes with time if the contact angles are not replenished from one timetep to the next. An explicit treatment of time evolution of the contact angle distribution requires 5 extratracer to memorize which contact angles from the distribution were already used in the timesteps before and thus make the scheme computationally more expensive. Using extratracer for the contact angle of mineral dust particles would approximately lead to an increase of computational costs of 21% in the GCM ECHAM6-HAM2. If the time evolution of the contact angle distribution is not taken into account, the α-pdf scheme becomes computationally similarly expensive as the single-α and the α(T ) scheme. However, the integral in Eq. 7 can not be solved analytically. Therefore, to minimize computational costs, a 10 look-up table could be used instead of discretized finite sums. Using look-up tables is depending on the size and format of the look-up table more expensive compared to solving an equation with simple constants as in the case of the single-α and α(T ) scheme. The α(T ) scheme is a simplified version of the α-pdf and computationally cheaper, because no integration over a contact angle distribution is necessary. The active sites scheme requires extratracer as well, which would lead to high computational costs (comparable to the α-pdf with explicit change of contact angle distribution with time).

Fitting immersion freezing measurements
In this section, the sensitivity of J imm and FF to different combinations of σ iw and ∆g # (see Ickes et al., 2015 for a discussion of these parameters) in combination with the contact angle schemes discussed in section 2.1 is analyzed by fitting and comparing the different CNT parameterization schemes to experimental data. This helps to understand how fit parameters influence 20 the calculated freezing curves.
The experimental data taken from Welti et al. (2012) consist of optically detected frozen fractions FF of droplets containing single immersed, monodisperse kaolinite (Fluka) particles. The FF was measured as a function of T , the particle radius r IN and the residence time in the measurement setup t. Experiments were performed using IMCA/ZINC [see Welti et al. (2012) for details]. The error bars of the data reflect the uncertainty in the distinction of water droplets and ice crystals in the detection 25 unit. For the sensitivity analysis the dataset measured after 10 s for kaolinite particles with a diameter of 400 nm is used. Note that the size of the particles might be underestimated due to the assumption of sphericity and therefore the calculated nucleation rates J imm from experimental frozen fractions represent a lower limit.
To explore the sensitivity of J imm to thermodynamic and kinetic parameters we use different CNT formulations with thermodynamic and kinetic parameters from different studies. The focus is on σ iw . In the following all CNT formulations, which are 30 used for the analysis, are listed. An overview is given in Table 2. Capital letters in the naming indicate the author from whose publication thermodynamic and kinetic parameters are obtained.
#1: Single-α R&D + Z The first approach is to use a single-α scheme in combination with the thermodynamic and kinetic parameters shown to be in good agreement with homogeneous nucleation rates [see Ickes et al. (2015)]. When using a single-α scheme it might be important that the kinetic and thermodynamic parameters are a combination which reproduces the homogeneous freezing data well as there is only one fit parameter and deviations cannot be compensated by additional parameters. The 5 emerged best fitting combination of σ iw and ∆g # from the analysis of homogeneous data is σ iw from Reinhardt and Doye (2013) and ∆g # from Zobrist et al. (2007). It is tested if these formulations of σ iw and ∆g # are also applicable to reproduce heterogeneous nucleation rates in combination with the single-α scheme.
#2: Single-α R&D + fit ∆g # An second approach which emerged from Chen et al. (2008) is using a constant ∆g # as an additional fit parameter to f /α 10 instead of taking a temperature dependent formulation. This assumption might be wrong in the context of homogeneous freezing especially at very low T (Barahona, 2015). However, it should be applicable for immersion freezing conditions as the change in ∆g # is small in the corresponding temperature range. The approach is used in combination with σ iw from Reinhardt and Doye (2013) and a single-α scheme. To decide if σ iw from Reinhardt and Doye (2013) is the best choice, different expressions for σ iw are tested against a fit of the homogeneous nucleation rate J hom using constant 15 ∆g # (see Fig. 2 analogous to Fig. 17 in Ickes et al.,2015). We find that σ iw from Reinhardt and Doye (2013) remains an appropriate choice even when ∆g # is used as a constant.
#3 and #4: Single-α O + fit ∆g # and single-α E + fit ∆g # To capture the whole possible range, two formulations of σ iw are used. One from Eadie (1971) leading to the lowest homogeneous nucleation rate and a second formulation of σ iw from Ouchi (1954) leading to the highest homogeneous 20 nucleation rate (see Fig. 2). For a summary of the two formulations of σ iw we refer to Ickes et al. (2015).
These two extremes of σ iw are used together with a constant ∆g # (fit parameter) and the single-α scheme to clarify if a fit of α can compensate for a low/high σ iw .
#5: α-pdf R&D + Z This CNT formulation is an α-pdf scheme using the same thermodynamic and kinetic parameter as the CNT formulation 25 for the single-α scheme (CNT #1). The α-pdf increases the complexity and adds an additional fit parameter compared to the single-α R&D + Z formulation (#1). Thus we test the influence of the choice of the contact angle scheme on the fit result. Additionally we examine if the number of free fit parameters has an impact on the result.
The dataset from Welti et al. (2012) is fitted with the previously listed CNT formulations. The fits are done by least-square minimization of FF as a function of T . The results are shown in Fig. 3 and Table 2. Table 2 additionally contains information on ∆g # and σ iw . Overall most CNT formulations with two fit parameters are able to capture the freezing data well with 5 similar root mean square errors (RMSE) of the estimated freezing curve and measured freezing data independently of the thermodynamic and kinetic parameters. In the following the results are discussed in more detail.
The single-α R&D + Z formulation (CNT #1) poorly captures the experimental data and results in a too steep freezing curve. With the single-α scheme it is not possible to reproduce the reduction of the energy barrier in a correct manner and to decrease the temperature dependence of the nucleation rate from the homogeneous to the heterogeneous case. Having only 10 one fit parameter, which in this case is a factor in the exponential term, is not sufficient to shift and flatten the freezing curve compared to homogeneous freezing. Only the T -shift of the freezing curve compared to homogeneous freezing is captured by the fitted single-α scheme. This can be seen in a more general illustration in Fig. 7 (2015)] support the applicability 15 of single-α in certain cases. Possibly a single-α scheme can be used for highly efficient IN triggering ice formation at low supercooling. In the context of the Young equation (Young, 1805) the assumption of a single-α scheme is questionable. If α represents the balance of the three surface tensions it has to change with T since the surface tensions are temperature dependent.
Using ∆g # as an additional fit parameter (CNT #2) can flatten the curve. In this case both fit parameters are factors in the exponential term of the nucleation rate with a similar influence on the fitted FF. The fit parameter which is multiplied with 20 the temperature dependent variable ∆G (f ) mainly shifts the freezing curve but cannot reduce the steepness sufficiently at the same time (see single-α scheme). Using a second fit parameter ∆g # resolves this issue. A simplified view on this is that one fit parameter is responsible for the shift and the other one for the flattening of the immersion freezing curve compared to homogeneous freezing. Using a constant ∆g # might be reasonable based on the results from the homogeneous freezing analysis (Ickes et al., 2015), but fitting ∆g # to immersion freezing data leads to substantially higher ∆g # than those estimated 25 by theoretical calculations (see Ickes et al., 2015). Moreover the fit value of ∆g # is aerosol-specific. This might be an artificial result and it is questionable if the assumption of a temperature independent and aerosol type specific (due to the fitting) ∆g # is a physical valid approach. It contradicts the assumption that the kinetic parameters such as ∆g # are the same for homogeneous and heterogeneous nucleation. The general approach to take the same thermodynamic and kinetic parameters (besides f and the prefactor) for homogeneous and heterogeneous nucleation is based on the assumption that mechanisms (e.g. the diffusion 30 of water molecules across the water-ice boundary) in the supercooled water are not influenced by the immersed aerosol particle.
This hypothesis might not be true. The aerosol particle might influence e.g. the diffusion of water molecules close to the particle which could justify a change in ∆g # depending on aerosol type. A disturbance of the diffusion by an aerosol particle could e.g.
be related to an impact on the hydrogen bond network by IN surface charges or polarizability of the particle surface (Edwards and Evans, 1962).
An alternative to having ∆g # as an additional fit parameter is to use a more sophisticated contact angle scheme, e.g. the α-pdf scheme (CNT #5 and #6) or the α(T ) scheme (CNT #7). Both approaches lead to good fits of the experimental freezing data and can be physically justified because they resemble the natural variability of an IN population. The α(T ) scheme has the disadvantage that it is not inherently known how the apparent α changes with T . A wrong assumption could lead to an unphysical contact angle scheme.

5
The curves resulting from CNT #3-7 (single-α + fit ∆g # schemes, α-pdf schemes and α(T ) scheme) all support the hypothesis that increasing the number of fit parameters from one to two allows to find a reasonable fit, independent of the kinetic and thermodynamic parameters chosen and also independently of the contact angle scheme. This is also supported by the result of the single-α scheme, where one fit parameter alone cannot shift and flatten the freezing curve. However, using a single-α scheme with a different additional fit parameter (e.g. the slope of σ iw instead of a constant ∆g # ) does not lead to a better fit 10 of the freezing curve. This might be due to the formula of the energy barrier preventing a sufficiently large influence of the additional fit parameter on the steepness of the curve. ∆g # as an additional fit parameter is able to reduce the steepness of the freezing curve as it has the opposite temperature dependence compared to the energy barrier ∆G. For a visualization of how the fit parameters influence the freezing curve for each scheme, please see Fig Using one CNT formulation, e.g. the single-α R&D + fit ∆g # formulation (CNT #2) together with fit parameters, e.g. 15 α, emerging from a fit with a different CNT formulation, e.g. the single-α E + fit ∆g # formulation (CNT #4) leads to a wrong freezing curve, which is illustrated in Fig. 3 [solid red line, single-α R&D + fit ∆g # (E) (#2/4)]. It is an example of the implication if fit parameters are used together with a different CNT formulation, that was not the one used to derive the fit parameters. This can unintentionally happen in GCMs if an implemented CNT formulation is later extended, e.g. by another aerosol species, and no care is taken that the fit parameters for the new species from the literature are derived from the same 20 CNT formulation as the one implemented in the model.
Looking at Table 2 one can see that f differs substantially, e.g. when using the single-α + fit ∆g # scheme with different values for σ iw . Comparing σ iw from Reinhardt and Doye (2013) with σ iw from Ouchi (1954) leads to a difference in fitted f of more than 300%, which translates into a difference in contact angle α of approx. 75 • . However, all single-α + fit ∆g # schemes result in a nearly similar freezing curve with the same RMSE. The fit parameters from the contact angle scheme compensate 25 inaccuracies associated with from thermodynamic and kinetic parameters and thus mask potentially wrong assumptions, e.g. Ouchi (1954). This makes it challenging to compare contact angles or fit parameters from different studies if not the same CNT formulation was used. Hence in the next subsection we investigate how fit results vary when thermodynamic and kinetic parameters differ and if there is a possibility to compare fit parameters from different studies using different CNT formulations. 3.2 Uncertainty of fitting α and α-pdf Table 2 shows that, dependent on the choice of σ iw and ∆g # , the estimated fit parameters differ. As an example, a contact angle estimate using CNT with e.g. σ iw from Pruppacher and Klett (2000) compared to using σ iw from Zobrist et al. (2007) leads to different results. In this section the sensitivity of two contact angle schemes to σ iw and ∆g # is investigated.
The two CNT formulations used in this analysis are a single-α R&D + fit ∆g # formulation (CNT #2) and the α-pdf R&D + Z formulation (CNT #5) described in sect. 3.1. We chose these two formulations because CNT #2 is used in GCMs and CNT #5 to interpret our data. Both schemes contain two fit parameters (f and ∆g # in CNT #2, µ and variance σ of the contact angle distribution in CNT #5).

5
We analyze how the two fit parameters depend on a change in thermodynamic and kinetic parameters. For this purpose the thermodynamic and kinetic parameters are varied up to ± 50%. For each variation (e.g. an increase of σ iw by 10%) fitting is done to the same immersion freezing data from Welti et al. (2012) as in the previous section. For the CNT #2 the thermodynamic parameter σ iw is varied, for the CNT #5 the thermodynamic and kinetic parameters σ iw and ∆g # are varied separately. The resulting fits are then compared to the reference fit results of section 3.1 (see Table 2).

10
In both cases (CNT #2 and #5) a similar change in fit parameters can be seen. Changing the thermodynamic parameter σ iw has a stronger impact on the fit parameters than changes in the kinetic parameter ∆g # . This is expected from the nucleation rate formula, where σ iw enters the calculation of the nucleation rate to the power of three and therefore changes the nucleation rate or frozen fraction more drastically than a change in ∆g # .
In case σ iw is increased/overestimated, the fit parameters are decreasing to compensate the change (see dashed arrow in 15 Fig. 4a) and conversely (see dotted arrow in Fig. 4a). The behavior of this compensation is not symmetric but follows the structure of the nucleation rate formula, i.e. 1/x dependence for f or µ and variance σ, respectively. That implies that the change in the fit parameter gets larger the larger the variation of σ iw is. The relative change approaches negative 100% with increasing σ iw . A larger deviation can be seen for the case where σ iw is decreased/underestimated. the variation of σ iw the larger is the deviation in f from the reference fit value, whereas ∆g # remains unchanged. Only one fit parameter (f ) is compensating the change in σ iw . Looking at Eq. 1 and 4 the product of f and σ 3 iw enters the exponential of the thermodynamic part of Eq. 4. As long as f · σ 3 iw = constant the resulting nucleation rate J imm is the same. Therefore ∆g # is not sensitive to a deviation of σ iw . Figure 4b) shows the relative change of the fit parameters as a function of percental change in σ iw and ∆g # for CNT #5. In 25 case ∆g # is changed (Fig. 4b) the compensation is linear for µ, following the structure of the nucleation rate formula. For σ the compensation is linear for small changes of ∆g # (until a change of approximately 30%), but is nonlinear for larger variations of ∆g # .
Summarizing, an over/underestimation of σ iw has a strong effect on the value of the resulting fit parameter, while the effect of over/underestimation of ∆g # is small. If fit parameters were estimated based on fitting different CNT formulations they can 30 not be directly compared. Fig. 4 can be used to estimate how different fit parameters would change due to different assumptions for σ iw or ∆g # . An example how fit parameters change with different CNT formulations is given in Appendix B.

Strategy to evaluate different CNT formulations
The sensitivity analysis in section 3 shows that most CNT formulations with at least two fit parameters are able to reproduce measured freezing data. Contact angle schemes are mostly judged based on the RMSE of the estimated freezing curve and measured freezing data. Only looking at the reproducibility of freezing data, however, is not a conclusive measure of how well immersion freezing is represented, since the fit parameters can mask uncertainties in the thermodynamic and kinetic parameters 5 of CNT. We propose to evaluate different CNT formulations by testing different fit properties against measurements. The evaluation consists of a macroscopic and a microscopic perspective. Three evaluation criteria are suggested.
At a macroscopic level: Is the CNT formulation capable of reproducing the measured T -dependence of FF and how well is dependence on IN size and time captured? According to the sensitivity of immersion freezing on T , r IN and t, the representation of the temperature dependence, followed by the size of the IN and the predicted time dependence of the freezing process should 10 be captured by a suitable CNT formulation if it is used as a function of T , r IN and t in a GCM.
When evaluating fits to a dataset that contains freezing data as a function of T , r IN and t the goodness of the fit implicitly includes all three aspects. However, because T has the strongest effect on the FF , the goodness of fit mostly reflects how well the CNT captures the temperature dependence. Two criteria emerge from the macroscopic level: Crit. 1 How accurately can the overall freezing data be reproduced, i.e. how well is the temperature dependence of the FF 15 captured by the CNT formulation?
Crit. 2 How accurately are the particle size and time dependence of the freezing process captured by the CNT formulation?
Criterion 2 can only be investigated if time and particle size dependent measurements are available.
At a microscopic level: Do the fit parameters match the microphysical assumptions of CNT and are they in general physically reasonable? To evaluate if fit parameters are physically reasonable, the analysis of heterogeneous freezing can be combined 20 with findings from homogeneous freezing. Including homogeneous freezing into the analysis might be useful because of fewer unconstrained parameters in this case. The microscopic criteria to be tested is: Crit. 3 Are the values for the fit parameters reasonable in the context of what we know about the microphysical process of nucleation?
In the following these three criteria are used to decide which CNT formulations are suitable for parameterizing immersion 25 freezing, e.g. in a GCM.

Using experimental data to estimate CNT parameters for different contact angle schemes
In the following a comprehensive dataset of FF (different aerosol species, aerosol particle sizes and residence times in the experiment) is used. Five different mineral dust types were chosen for the analysis: Fluka kaolinite, illite-NX, montmorillonite, microcline (K-feldspar) and ATD (Arizona test dust). Montmorillonite or kaolinite are often used in global models as a surro- 30 gate for mineral dust in terms of ice nucleation, e.g. montmorillonite in ECHAM6-HAM2. Fluka Kaolinite, which was used here, has been widely used to study the mechanism of immersion freezing. Illite-NX was chosen as a the mineral dust reference sample for an instrument intercomparison (Hiranuma et al., 2015). A microcline sample from Namibia and ATD are included to enable sensitivity studies of the freezing parameterization scheme with more efficient IN. The experimental data for kaolinite is taken from Welti et al. (2012). Illite-NX data has been published in Hiranuma et al. (2015). Microcline, montmorillonite and ATD data are new datasets. All measurements were performed using size-selected aerosol particles with diameters of 50, 100, 5 200, 400, 800 and 920 nm and 10 s residence time. The kaolinite dataset contains time dependent measurements for different residence times of 1, 2, 3, 6, 9 and 21 s. Note that the residence times are rounded to full seconds [compared to Welti et al. (2012)] and not all datasets include the smallest and/or largest size (kaolinite: 100 -920 nm, illite 100 -800 nm, montmorillonite 100 -800 nm, microcline 50 -800 nm, ATD 100 -800 nm). The error bars of the data reflect the detection uncertainty and the statistical uncertainty in the measurement by multiple measurements. [single-α scheme with ∆g # as a fit parameter], #5 [α-pdf scheme] and #7 [α(T ) scheme]. For more details see Table 2. CNT #3, #4 and #6 use a σ iw , which was found to not represent homogeneous freezing well and are excluded based on criterion 3.
The wrong assumption of σ iw was chosen on purpose for the sensitivity study in section 3.1 to demonstrate how that influences the fit results and the freezing curves. Note that also the single-α R&D + Z formulation (CNT #1) is not expected to be able 15 to reproduce the experimental freezing data (Crit. 1). It is still included here for comparison of the RMSE value with the other formulations ("bad" reference).
The fit parameters are determined by least square minimization of the calculated versus measured FF . For this purpose the dataset of each dust species, including all measurements as a function of T , r IN and t, is used. To get an impression of the variability of the fit parameters throughout a dataset, the kaolinite dataset is fitted for each size and time separately in Appendix 20 C.
The fit parameters for the different CNT formulations and aerosol types are shown in Table 3 together with the best fit root mean square error (RMSE). The fit curves in comparison to the measured FF are shown in Fig. 5 to Fig. 6 (in the case of kaolinite only a selection of the data is shown).
The geometric term f in Table 3 is smallest for microcline, showing that this is the most efficient IN investigated here. The 25 second lowest value for f is found for ATD. Montmorillonite and kaolinite seem to be quite similar in terms of IN efficiency, whereas illite is the least efficient IN.
Revising the fit results with criterion 1 shows that CNT #1 is too steep and not able to reproduce experimental data, resulting in a high RMSE. One fit parameter is not enough to shift and reduce the steepness of the immersion freezing curve sufficiently compared to homogeneous freezing. Thus the single-α scheme (CNT #1) does not fulfill criterion 1. Reasonable fit results (low 30 RMSE) are obtained with CNT #2 and #5 for all datasets.
In the case of CNT #5 the mean contact angle is very similar to the contact angle of CNT #1, i.e. the medium freezing temperature (FF =0.5) is similar, but the steepness of the freezing curve is reduced by the variance of the contact angle distribution σ. See App. A for additional analysis on the influence of µ and σ on the fit of FF .
When fitting CNT #7 more than one solution for fit parameters is found (no absolute minimum of the fitting function). However, 35 α 0 should be positive and m has to be negative so that α increases with decreasing T . Otherwise criterion 3 is not fulfilled.
Here only the fit parameters that fulfill criterion 3 are reported (local minimum of the fitting function).
For all CNT formulations, the fits with largest RMSE are the ones for ATD which is probably caused by a size dependent mineralogy of ATD. The CNT formulation that best reproduces immersion freezing varies from dust to dust. Therefore, we establish a ranking for each dataset similar to the methodology used in Wheeler et al. (2014). The best CNT formulation gets 5 a ranking of 1, the worst a ranking of 4. From the ranking of the different datasets an average ranking is derived to judge the overall capability to predict FF for each CNT formulation. The ranking (see Table 4) shows that CNT #7 and CNT #2 are the best followed by CNT #5 and CNT #1. Calculating the average RMSE from all fits (as an alternative) leads to the same result, where CNT #7 is the best and CNT #1 the worst (see also Table 4). Note that this ranking does not consider criterion 2 and 3 and is only based on fit statistics. It also does not show directly how good the CNT formulations reproduce time and particle size dependence of the freezing process (Crit. 2).

Testing the time and particle size dependence (Criteria 2)
To test the ability of the single-α R&D + fit ∆g # formulation (CNT #2), the α-pdf R&D + Z formulation (CNT #5) and the α(T ) R&D + Z formulation (CNT #7) to reproduce experimentally observed size and time dependence, the fit parameters for 20 kaolinite (see Table 3) are used to calculate FF s for three different residence times (1, 10 and 21 s) and three different aerosol particle diameters (100, 400, 800 nm). In the case of the size dependent calculation of FF the time is 10 s, in the case of the time dependent calculation the aerosol particle diameter is 400 nm. The calculated FF is compared to measurements of the size and time dependent FF in Fig. 5. The analysis of the RMSE of the fit for each subset of data and CNT formulation revealed marginal differences in the second decimal place and is therefore not reported. 25 Figure 5 shows that CNT #5 is able to capture the time and particle size dependence better than schemes #2 and #7. This leads to an overall smaller RMSE and explains the better ranking for CNT #5 in the case of kaolinite (see Table 3). Looking at Fig. 5, CNT #2 and #7 give very similar results which overpredict both the size and time dependence, while CNT #5 seems to underpredict the particle size dependence but captures the time dependence well. Overpredicting the size dependence translates into an overestimation of FF for particles with an aerosol particle diameter larger than 400 nm and an underestimation of FF 30 for particles with an aerosol particle diameter smaller than 400 nm. Underpredicting the size dependence has the opposite influence on FF . Overpredicting the time dependence means that FF is overestimated when using a timestep larger 10 s as in GCMs. The outcome of the evaluation depends on the dataset used. For different aerosol species the ranking of CNT #5 and #2 differs, e.g. for montmorillonite, microcline and ATD. Due to this limitation, it cannot generally be concluded which contact angle scheme best fulfills criterion 2. All three schemes might thus be chosen for CNT based immersion freezing parameterization schemes in GCMs. The computationally least expensive formulations to use in a GCM would be CNT #2 and #7.
In Fig. 5 the CNT curves are also compared to an empirical immersion freezing parameterization scheme (n s,IN scheme) based on the expression given in Niemand et al. (2012) fitted to the kaolinite measurements (for details please see Appendix D). The 5 n s,IN scheme slightly overestimates the particle size dependence and the scheme does not include any time dependence since it is deterministic. Derived FF curves appear similar to CNT curves. However, due to the general characteristics of empirical relations it is not clear if it can be extrapolated to a wide T -range, which would be mandatory for the use in a GCM.

Testing the physical reasonability of fit parameters (Criteria 3)
Since our knowledge about the microphysical behavior of supercooled water is limited, and measurements on the microscopic 10 level are very difficult, the evaluation of the reasonability of the fit parameters is challenging. Homogeneous freezing measurements or results for molecular dynamics simulations can be used to additionally evaluate certain CNT formulations and contact angle schemes. Molecular simulations recently started to simulate heterogeneous freezing, often using kaolinite surfaces Michaelides, 2007, 2008;Croteau et al., 2010;Šolc et al., 2011). Michaelides (2007, 2008) found, that the ice nucleation ability of kaolinite results from the amphoteric nature of the surface hydroxyl (OH) groups, meaning that these OH 15 groups can act as hydrogen bond acceptors or donators. The surface hydroxyl groups are reacting flexibly to the orientation of the water molecules above the surface. Therefore water molecules bonding more easily with the kaolinite surface compared to other water molecules form an ice-like stable system. Such information can be used to understand the effect a particle surface can have on the kinetics of ice formation. The study of Šolc et al. (2011) is an example for a molecular simulation study that could help to evaluate fit parameters. Šolc et al. used force-field molecular dynamics simulation to investigate water nan-20 odroplets on a kaolinite surface and estimated a microscopic contact angle of approx. 105 • . That is equivalent to a geometric term f of 0.69 (Eq. 3). Unfortunately none of the results for the different CNT formulations in our study is close to this value. In all studies (including the present one) the single-α scheme does not represent the measured freezing data very well (only when ∆g # is used as an additional fit parameter). Even in the recent study of Alpert and Knopf (2016) the single-α scheme is not able to reproduce the freezing curve of a uniform surface (organic monolayer coating). The single-α + fit ∆g # scheme 30 has not been evaluated against other schemes in the above mentioned studies. The ranking of the α(T ) and α-pdf scheme differs from study to study. In Zobrist et al. (2007) Knopf and Forrester (2011); Rigg et al. (2013) the temperature dependence of α is based on the change in interfacial tensions with T (according to Young (1974)).

Comparison to other studies
The α-pdf scheme did not produce good results in the study of Rigg et al. (2013) (aqueous ammonium sulfate droplets). They suggest that in case the α-pdf scheme is used, it should be expressed as a function of water activity. In the study of Wheeler 5 et al. (2014) it leads to good results for ATD, but less so for kaolinite. In both cases the best scheme was the active sites scheme, which was not investigated in the present study. In Welti et al. (2012) as well as in the present study the α-pdf scheme reproduces the data quite well. Note that this study, the study of Lüönd et al. (2010) and Welti et al. (2012) are the only studies where the size and time dependence has been tested separately (Crit. 2) as size and time resolved data was available. From the results presented here, we conclude that the α-pdf scheme is a suitable scheme to represent the time-dependence for kaolinite.

10
This could explain why it achieves a better ranking here than in Wheeler et al. (2014), where no size and time dependence could be tested against measurements. This paper includes the most complete evaluation (taking into account criterion 2 and 3 in addition to criterion 1) of contact angle schemes for immersion freezing of mineral dust types.

Conclusions
In this study the sensitivity of CNT based immersion freezing parameterization schemes to thermodynamic and kinetic param-15 eters is investigated. We discuss their effect on the fit parameters of contact angle schemes when fitting measurement data. For the use in models an assessment of sensitivities is important to estimate uncertainties originating from different parameterization schemes which represent the effect of aerosol particles on the energy barrier of ice nucleation.
Compared to homogeneous freezing, immersion freezing has one more unconstrained parameter, namely the geometric term f (α). Different schemes to represent the contact angle or contact angle distribution with one or two fit parameters based on 20 experimental data are tested. It is found that contact angle schemes containing two fit parameters are able to reproduce experimental FF , while the contact angle scheme with only one fit parameter, the single-α scheme, cannot reproduce the freezing curves.
Analyzing the importance of the choice of σ iw and ∆g # to parameterize immersion freezing revealed that uncertainties in the thermodynamics or kinetics can be compensated by two-parameter contact angle schemes. As a result an under/overestimation 25 of σ iw or ∆g # does not lead to a bad representation of freezing curves as in case of homogeneous freezing. Because the fit parameters compensate inaccuracies or uncertainties of the thermodynamic and kinetic parameters, the absolute value of the found fit parameters is highly dependent on the choice of thermodynamic and kinetic parameters within the formulation of CNT (especially on σ iw ). As a consequence, contact angles for CNT parameterization schemes from different authors should only be applied within the same CNT formulation used to derive the parameters. Implementing a CNT based parameterization scheme 30 into a GCM demands that parameters derived from experiments are calculated using the same CNT formulation. Otherwise an offset of the freezing temperature in clouds might be introduced. The sensitivity of α on the thermodynamic and kinetic parameters used in CNT biases a direct comparison of contact angle values derived in different studies. We emphasize the importance of highlighting which CNT formulation was used for the analysis of experimental data.
Contact angle schemes with two or more fit parameters reproduce freezing curves but can be unphysical or limited to the dataset used for fitting. Contact angle schemes intended to represent immersion freezing properties of a heterogeneous particle population under a variety of environmental conditions and time scales, should reproduce the T -dependence of freezing and 5 the ability to predict the size and time dependence of the freezing process ( Crit. 1 and 2). The schemes should be conform to the microphysical assumptions CNT is based on (Crit. 3). Particle size and nucleation time are implicitly included in the reproducibility of freezing curves of the dataset but should be investigated separately. It should be noted here, that most experimental datasets are not accounting for both, the size and time dependence of freezing. Experimental data without any information on the time or size dependence limit the assessment according to the criteria defined in this study. Here only a 10 full analysis of the CNT parameterization schemes for one mineral dust (kaolinite) was possible. Having limited datasets may lead to unrepresentative conclusions. More size and time dependent measurements of different IN are desirable to compile parameters and find a robust CNT formulation.
In this study we derived fit parameters for five different datasets of mineral dust (kaolinite Fluka, illite-NX, montmorillonite, microcline and ATD) by fitting different CNT formulations to the FF from measurements. Good results in reproducing the 15 freezing curves (Crit. 1: T -dependence and partly 2: size and time dependence) are achieved when using a single-α scheme with fitted constant ∆g # (CNT #2), an α-pdf scheme (CNT #5) or an α(T) scheme (CNT #7).    Simplified semi-singular: All IN have the same contact angle. The contact angle changes with T (α1 < α2 < α3 in the example). α is equivalent to the mean contact angle of the α-pdf scheme (e µ ) or mean activity of the α-pdf distribution, which changes with temperature. Fit q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q qq q q qq q q q qq q q qq q q qqq q q q q q qq q q q qqq q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q qq q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 2. Comparison of the fitted J hom (T ) (solid black line) with calculated nucleation rates using different formulations of σiw and constant values for ∆g # . Grey dots show experimental homogeneous freezing data. σiw from Reinhardt and Doye (2013) Figure 3. Calculated frozen fraction FF as a function of T for different thermodynamic and kinetic parameters in combination with different contact angle schemes for kaolinite with a particle diameter of 400 nm after a residence time of 10 s. More details can be found in Table 2. (e) d=400 nm, t=21 s Figure 5. Calculated FF of kaolinite for three different times and three sizes using the single-α R&D + fit ∆g # formulation (CNT #2), the α-pdf R&D + Z formulation (CNT #5) and the α(T ) R&D + Z formulation (CNT #7) with corresponding fit parameters (see Table 3  for certain sizes using the single-α R&D + fit ∆g # formulation (CNT #2), the α-pdf R&D + Z formulation (CNT #5) and the α(T ) R&D + Z formulation (CNT #7) with corresponding fit parameters (see Table 3). A simplified immersion freezing parameterization scheme based on Niemand et al. (2012) is additionally compared to the dataset.   Decreasing f (reducing the energy barrier) shifts the freezing curve to higher temperatures while the slope slightly increases. Increasing ∆g # shifts the freezing curve to lower temperatures and changes the slope of the freezing curve. A higher activation energy barrier leads to a flattening of the curve. Decreasing m (increasing contact angle with decreasing T because a broader part of the IN populations contact angle distribution causes nucleation) shifts the freezing curve to lower temperatures and changes the slope of the curve (flattening).
Increasing µ (increasing the average contact angle) shifts the freezing curve to lower temperatures and slightly flattens the curve, while an increase in σ (broadening of the contact angle distribution) changes the slope of the freezing curve only (flattening). Note that already small changes in µ lead to a considerable shift of the curve compared to the other schemes.
Appendix B: Estimating the relative uncertainty of fitting α: Example calculations Figure 4 can be used to estimate the deviation of fit parameters from different CNT formulations relative to each other. To show this, we estimate the difference in fit parameter when using σ iw from Eadie (1971) instead of Reinhardt and Doye (2013) in combination with an α-pdf scheme (CNT #6 compared to CNT #5).
Within the 10 K temperature range of the immersion freezing measurements (236-246 K) σ iw is on average 4% higher when 5 using σ iw from Eadie (1971) instead of Reinhardt and Doye (2013). From Fig. 4 b we see that an increase in σ iw by 2.5% (246 K) or 5% (236 K) would lead to a decrease in µ by approximately 7 to 13%. Now we check if that estimated change matches with the real change when fitting the same dataset with the two different σ iw . In Table 2 using σ iw from Eadie (1971) leads to a mean contact angle of 0.44 rad (approx. 25.5 • ) instead of 0.5 rad (approx. 28 • ) when using σ iw from Reinhardt and Doye (2013). This is a difference of 12%, conform with the estimate from Fig. 4 (approx. 7-13%). However, the variance σ of 10 the α-pdf distribution is expected to change less (5 to 9%) but a change by 25% is found for the best fit.
In some cases the predicted change in fit parameters from Fig. 4 deviates from the real change in fit parameters (Table 2) because both parameter µ and σ are changed at the same time. For Fig. 4 one was held constant. Another problem with Fig. 4 is, that the assumption of a constant variation of σ iw over the whole temperature range is invalid in most cases. However, Fig. 4 can be used to illustrate how fit results might change and estimate a rough deviation from the reference when using different 15 thermodynamic and kinetic parameters especially for cases where σ iw changes nearly constant over the fitted temperature range. This can help when comparing fit results to fit results from studies where a different formulation of CNT was used.
Appendix C: Variability of the fit parameters throughout one dataset The variability of the fit parameters throughout the dataset can be seen when fitting the single-α R&D + fit ∆g # formulation (CNT #2), the α-pdf R&D + Z formulation (CNT #5) and the α(T ) R&D + Z formulation (CNT #7) to FF of kaolinite data for different sizes and residence times separately. The resulting fit parameters are compared in Fig. 9 for different sizes in red and for different times in blue and light blue. Each point in Fig. 9 represents the value of the best fit parameter for one subset 5 of the kaolinite dataset. The labels on the x-axis give information which subset of the dataset was fitted. The residence time is 10 s for the data subsets of different sizes and the diameter of the kaolinite particles 400 nm or 800 nm for the data subsets of different times (blue/light blue). The dashed line indicates the mean of the size or time dependent fit parameter. The range is shown as shaded box.
The fit parameters vary depending on the measurement conditions. Omitting the measurement with the shortest residence 10 time (1 s) the variation between the data subsets for different times is small. The variability of the fit parameters is larger for different aerosol particle sizes compared to different residence times, which might be due to the higher sensitivity of the freezing process to particle size compared to time. For CNT #2 the fit parameters seem to be correlated. High values of one fit parameter, e.g. f , correspond to low values of the other fit parameter, e.g. ∆g # . CNT #5 on the other hand does not show a clear correlation for the time dependent subsets, but for the size dependent subsets. The same yields for CNT #7.

15
The different fit results for CNT #5 can be used to study how the shape of the contact angle distribution might change with the size of the particles or the residence time. Whereas the fitted contact angle distribution does not change noticeable with time between 1 and 21 s (Figure not shown here), the mean contact angle (e µ ) and the variance σ changes with particle size (see all mineral dust types except microcline the contact angle distribution broadens with increasing aerosol particle size (neglecting the fit of the 400 nm dataset of kaolinite, which does not fit into the picture). This reflects a larger probability of different α on the aerosol particle population with increasing size. The particle population is more heterogeneous. In the case of microcline the contact angle distribution narrows with increasing aerosol particle size, the ice nucleating properties of the microcline aerosol population seem to get more homogeneous with size. Note that the curves (here shown for kaolinite and microcline) 25 are not considering measurement uncertainties of the fitted data and therefore can only be used to qualitatively interpret the result. In case of idealized measurements the result could be used to derive a relationship between mean contact angle (e µ ) or the width of the contact angle distribution σ and the size of the IN. Using the results of Fig. 9 and developing a size-dependent α-pdf improves the fit results.
C1 Uncertainty of fit parameters due to limited data 30 In many cases no size-and time-dependent measurements are available. Here we investigate the quality of fit parameters if only limited amount of data is available. For that purpose we use the kaolinite dataset (as this is the most extensive dataset available within this study) and use only subsets of the dataset to estimate the fit parameters. The quality of the gained fit parameters is then estimated by using the complete dataset and look how good the freezing curves can be represented (RMSE). We look at four different cases: 1. Reference, the whole dataset is fitted (see also Table 3).
2. Only size dependent measurements are available, time dependence is not known (t=10 s).
3. Only time dependent measurements are available, size dependence is not known (d=400 nm). 5 4. Only one measurement is available (t=10 s, d=400 nm).
The resulting fit parameters and the deviation from the entire kaolinite dataset is shown in Table 5. The fit parameters are not significantly different when the dataset is limited to only size or only time dependent data. Also the deviation from the complete dataset is not significant. This analysis therefore does not allow any conclusion how many dependencies, e.g. size and time, have to be taken into account to successfully fit freezing data. However, if using only one single dataset, the results 10 for the fit parameters are different and the deviation from the measurements is higher. Note that the deviation when fitting only a single dataset could be larger if a dataset is chosen which is not similar to the average values of the dataset as in this case.
This means that there is no guarantee that fits can be extrapolated/used in a universal way across different conditions.  Figure 9. Variability of the fit parameters for subsets of the kaolinite dataset. The variability with the aerosol particle size (d) is shown for a residence time t of 10 s (red). The variability with time (t1 and t2) is shown for a aerosol particle diameter of 400 nm (t1, blue) and 800 nm (t2, light blue). The first row shows the fit results for CNT #2, the second row for CNT #5, the third row for CNT #7. The RMSE value shows the deviation of the fit to the single dataset it was fitted to.   Niemand et al. (2012)]. In this approach it is assumed that ice nucleation occurs on localized sites, called active 5 sites. It is described as a function of T and particle surface area A IN because the amount of active sites is supposed to scale with particle size.
In Connolly et al. (2009) the change of the number of ice active aerosol particles in the size bin j, N i,j , with respect to T is: where N tot,j denotes the total number of aerosol particles in the size bin j and A j the dust particle surface area in the same 10 size bin. The surface site density of ice active sites n s,IN as a function of T can be determined by integrating the factor k(T ) over the whole temperature range: Using Eq. D1 and D2 the frozen fraction FF can be expressed as a function of T : The approximation is valid for A j · n s,IN (T ) 1, which translates into small particles and high temperatures. For low temperatures, e.g. 243.15 K and particles larger than 3 µm the term A j · n s,IN is approximately 1.
The surface site density of ice active sites, n s,IN (T ), is calculated from the total surface area of aerosol particles in the AIDA chamber and the measured ice crystal number concentration during one freezing experiment: with A tot,j the total surface area per unit volume of particles in the size bin j and A tot the total surface area over all size bins. The evaluation of 16 AIDA freezing experiments yields the following fit formula for the ice active surface site density of natural dust: Eq. D6 was used here to fit the dataset for the different mineral dust types. The fits were done in two different ways: By using Eq. D3 and Eq. D6 to fit the measured FF (method 1) or by using Eq. D3 to convert the FF measurements to logarithmic surface site densities and fit ln(n s,IN (T )) directly following Eq. D6 (method 2). The results are shown in Table 6 and in Fig. 11 (method 1). The scheme is labeled "n s,IN " in Fig. 5.  Table 6 shows that the results for the fit parameters are very different depending on whether the FF is fitted directly or the 5 active site density n s,IN is calculated from the FF and fitted afterwards. That is due to the characteristics of the freezing curve.
The very small FF at high temperatures and limited FF (to 1) at low temperatures lead to a flattening of the n s,IN curve.
Calculating n s,IN at low temperatures from FF close to 1 gives the number of active sites, which was needed to freeze all droplets. However it could be that more actives sites were present than needed to freeze all droplets. Therefore the tail (low and high FF ) of the FF dataset is often left out of the fitting. Here we investigate how the fit results and the freezing curves 10 change depending on the share of the dataset accounted for fitting. We use the complete dataset as a first step and then omit FF data higher than 0.9 and lower than 0.1 and 0.8 and 0.2, respectively. Figure 12 shows this exemplary using a dataset of kaolinite particles (d=400 nm, t=10 s). In Fig. 12b it can be seen that the surface site density n s,IN varies depending on the share of the FF dataset. The implication of the different estimations for n s,IN is shown in Fig. 12c on the direct estimated n s,IN captures the data well and falls on the freezing curve of the indirect n s,IN using the same share of FF data. We recommend to cut the tails from the FF data when n s,IN is fitted directly. When estimating n s,IN indirectly by using FF there is no need for cutting the tail. Not cutting the tail increases the amount of data available for the fitting and might therefore be preferable. However, very low/high FF are most susceptible to experimental uncertainties, which could be a legitimation for cutting the tail away from the dataset. Because the results are different depending on the methodology this 5 sensitivity should be taken into account when comparing different fit parameters of n s,IN from literature. We recommend to estimate n s,IN indirectly by fitting FF as a standard procedure.