Analysis of isothermal and cooling rate dependent immersion freezing by a unifying stochastic ice nucleation model

. Immersion freezing is an important ice nucleation pathway involved in the formation of cirrus and mixed-phase clouds. Laboratory immersion freezing experiments are necessary to determine the range in temperature ( T ) and relative humidity (RH) at which ice nucleation occurs and to quantify the associated nucleation kinetics. Typically, isothermal (applying a constant temperature) and cooling rate dependent immersion freezing experiments are conducted. In these experiments it is 5 usually assumed that the droplets containing ice nuclei (IN) all have the same IN surface area (ISA), however the validity of this assumption or the impact it may have on analysis and interpretation of the experimental data is rarely questioned. A stochastic immersion freezing model based on ﬁrst principles of statistics is presented, which accounts for variable ISA per droplet and uses physically observable parameters including the total number of droplets ( N tot ) and the heterogeneous ice nu- 10 cleation rate coefﬁcient, J het ( T ) . This model is applied to address if (i) a time and ISA dependent stochastic immersion freezing process can explain laboratory immersion freezing data for different experimental methods and (ii) the assumption that all droplets contain identical ISA is a valid conjecture with subsequent consequences for analysis and plained by assuming identical ISA in each droplet. When accounting for ISA variability, the cooling rate dependence of ice nucleation kinetics vanishes as expected from classical nucleation theory. The model simulations allow for a quantitative experimental uncertainty analysis for parameters N tot , T , RH, and the ISA variability. In an idealized cloud parcel model applying variability in ISAs for each 25 droplet, the model predicts enhanced immersion freezing temperatures and greater ice crystal production compared to a case when ISAs are uniform in each droplet. The implications of our results for experimental analysis and interpretation of the immersion freezing process are discussed.


Introduction
Ice crystals in tropospheric clouds form at altitudes where temperatures fall below the ice melting point, also known as supercooled temperatures, and for conditions in which water partial pressure exceeds the saturation vapor pressure with respect to ice (Pruppacher and Klett, 1997;Hegg and Baker, 2009).Cirrus or mixed-phase clouds consist entirely of ice crystals or of ice crystals coexisting with supercooled aqueous droplets.These clouds can significantly impact the global radiative budget and the hydrological cycle (Baker, 1997;Rossow and Schiffer, 1999;Chen et al., 2000;Liu et al., 2007;Lohmann and Hoose, 2009;Tao et al., 2012;Rosenfeld et al., 2014); however, their formation is not well understood or constrained in cloud and climate models (Boucher et al., 2013).Ice nucleation precedes the formation of ice crystals.Homogeneous ice nucleation occurs from supercooled aqueous aerosol particles or cloud droplets.Ice formation can also occur at temperatures higher than the homogeneous freezing limit initiated by insoluble particles acting as ice nucleating particles (INPs).Heterogeneous ice nucleation can occur when INPs are immersed in supercooled aqueous droplets (immersion freezing), when INPs make physical contact with supercooled droplets (contact freezing), or when ice nucleates on INPs directly from the supersaturated vapor phase (deposition ice nucleation).It is impossible to observe in situ ice nucleation in the atmosphere and very difficult to infer the ice nucleation pathway (Haag et al., 2003;Hegg and Baker, 2009).Despite the established importance of the impact of heterogeneous ice nucleation on cirrus and mixedphase cloud formation, it is not included in global radiative forcing estimates (Myhre et al., 2013).
Laboratory studies are necessary to investigate at which thermodynamic conditions, i.e., temperature, T , and relative humidity, RH, and by which mode ice nucleation occurs for predictive use in cloud and climate models.This study presents a newly developed model simulation applied for analyses of previously published laboratory immersion freezing data obtained by different experimental methodologies.It allows prediction of atmospheric ice particle production under relevant scales of time and INP surface area (ISA).
Classical nucleation theory (CNT) is currently the only available physical theory to describe ice nucleation.Simply stated, CNT quantifies a maximum Gibbs free energy barrier corresponding to the minimum number of water molecules in a cluster that has to be overcome to initiate ice nucleation (Pruppacher and Klett, 1997).Cluster formation, as well as ice nucleation, occurs stochastically and is dependent on time, t, and in the case of homogeneous ice nucleation, the supercooled liquid volume, V .Koop et al. (2000) parameterized the theoretical homogeneous ice nucleation rate coefficient, J hom , as a function of T and water activity, (a w ; a w = 1.0 for pure water and a w < 1.0 for aqueous solution).This approach yields J hom to be independent of the nature of the solute and avoids the weakness of the capillary approximation in CNT (Pruppacher and Klett, 1997).
Immersion freezing can be described by CNT by reducing the free energy barrier due to the presence of a solid surface.Ice nucleation remains a stochastic process, but is dependent on the available ice nucleating surface area, A, instead of V (Pruppacher and Klett, 1997;Zobrist et al., 2007).The heterogeneous ice nucleation rate coefficient, J het , is a physically and experimentally defined parameter, which gives the rate of nucleation events for given surface area and unit time.Knopf and Alpert (2013) parameterized J het as a function of T and a w following Koop et al. (2000) using direct measurements of J het and J het derived from previous studies (Archuleta et al., 2005;Alpert et al., 2011a, b;Knopf and Forrester, 2011;Murray et al., 2011;Broadley et al., 2012;Iannone et al., 2011;Pinti et al., 2012;Rigg et al., 2013).Known as the a w -based immersion freezing model (ABIFM) (Knopf and Alpert, 2013), J het can be derived for different types of INPs, such as mineral dusts, organic, surfactant and biogenic, applicable for a w ≤ 1.0, and independent of the nature of the solute.The ABIMF is a holistic and computationally efficient description of the immersion freezing process for prediction of ice nucleation for atmospherically relevant conditions and applicable for a variety of experimental methods, including the droplet-on-substrate approach (Zobrist et al., 2007;Knopf and Forrester, 2011;Alpert et al., 2011a, b;Iannone et al., 2011;Murray et al., 2011;Broadley et al., 2012;Rigg et al., 2013), oil-encased droplets (Murray et al., 2011;Broadley et al., 2012;Wright and Petters, 2013), differential scanning calorimetry (Marcolli et al., 2007;Pinti et al., 2012), and continuous-flow diffusion (Rogers et al., 2001;Archuleta et al., 2005;Hartmann et al., 2011;Kulkarni et al., 2012;Wex et al., 2014).These previous studies represent a subset of a much broader selection of experimental methods and designs.
Different parameterizations of J het exist.Zobrist et al. (2007) investigated droplet freezing initiated by organic surfactant monolayers and parameterized experimentally derived J het values using parameterizations of the Gibbs free energy and diffusion activation energy.Reduction in the Gibbs free energy barrier is described by the parameter known as the contact angle, α, defined as the angle of contact between an ice embryo and substrate surrounded by the liquid parent phase and derived by balancing the interfacial surface tensions between the three (Pruppacher and Klett, 1997).It was found that a single value of α could not reproduce the experimental freezing data for organic monolayers, but when one allows α to be a linear function of T the data could be represented (Zobrist et al., 2007), corroborated by others and for different INPs (Knopf and Forrester, 2011;Alpert et al., 2011a, b;Welti et al., 2012;Rigg et al., 2013).It is important to note that a self-assembled organic monolayer is a completely uniform surface down to the molecular level, i.e., the notion of different ice active sites present does not apply (Gavish et al., 1990;Popovitz-Biro et al., 1994;Ma-jewski et al., 1995).For particles that have an uneven morphology, cracks, pits, or ridges such as mineral dust, a single α value has been shown also to not reproduce experimental data (e.g., Marcolli et al., 2007;Lüönd et al., 2010;Niedermeier et al., 2011b;Welti et al., 2012;Rigg et al., 2013;Niedermeier et al., 2014;Wheeler et al., 2015).Unlike Zobrist et al. (2007), some of these studies do not consider that α can vary with T , but instead randomly distribute α values on particles immersed in droplets while using all the other same parameterizations (i.e., the Gibbs free energy and diffusion activation energy) to calculate J het and not experimentally derive J het .This procedure is similar to distributing J het values over different droplets containing INPs as done by Broadley et al. (2012) and Herbert et al. (2014).These are successful approaches to describe the freezing data leading one to think that allowing α to be a linear function of T or distributing α across particles are in principle the same, both resulting in the necessity to change the contact angles to represent the freezing data.However, application of these frameworks are conceptually and mathematically very different and in fact result in very different interpretations of the data and underlying ice nucleation processes.The method of Zobrist et al. (2007) derives a single continuous function of J het (T ) for a specific particle type, while an active site approach randomly distributes multiple J het (T ) functions across particles and their surfaces.In other words, one can either use a single J het (T ) function (i.e., not applying a single α value) or use multiple J het (T ) functions (different only by utilizing different values of α), constituting completely different pictures of ice nucleation.One major advantage of the ABIFM approach chosen here is that it uses a single function of J het (T ) and avoids any use or calculation of α and instead uses water activity as a parameter.It inherently allows α to be both T dependent (Zobrist et al., 2007) or distributed in a particle population (an active site approach).A caveat to our approach is that it cannot make a statement about the active site distribution.Another advantage of our approach is that ABIFM can be used simultaneously for immersion freezing from pure water and aqueous solution droplets.When using the ABIFM, a uniform ice nucleating surface is not assumed; however, a single function of J het (T ) for a single particle type is assumed to describe the experimental data without invoking the presence of different (rare) and non-detectable ice nucleating sites or components present in some but not all droplets.
The major difficulty with a variety of experimental techniques is how accuracy and uncertainty of T , RH, t, and A are assessed and how these uncertainties affect extrapolation of laboratory derived ice nucleation parameterizations to atmospherically relevant conditions.Previous investigations have developed state of the art instrumentation and methods to constrain uncertainties (Connolly et al., 2009;Lüönd et al., 2010;Niedermeier et al., 2010;DeMott et al., 2010;Niedermeier et al., 2011b;Hoose and Möhler, 2012;Niemand et al., 2012;Rigg et al., 2013;Hiranuma et al., 2015;Vali and Snider, 2015).Additionally, progress in understanding ice nucleation behavior has been made by validating empirical parameterizations or models based on the concept of ice active sites, i.e., that surface sites on a particle have variable ice nucleating efficiencies, can be used to reproduce experimental data.However, there is no physical basis for these interpretations (Niedermeier et al., 2010) and they may be inherently constrained to the investigated range of T , RH, t, A and concentration of INPs from which they are derived (Rigg et al., 2013;Knopf and Alpert, 2013).These include the multi-component model (Murray et al., 2011), the timedependent freezing rate parcel model (Vali and Snider, 2015), parameterizations of INPs per liter of air (DeMott et al., 2010), the α-probability density function (PDF) model (Marcolli et al., 2007;Lüönd et al., 2010), the active site model (Marcolli et al., 2007;Lüönd et al., 2010), the singular description (Vali, 1971;Connolly et al., 2009;Alpert et al., 2011a, b;Vali, 2008;Murray et al., 2011;Hiranuma et al., 2015) and the soccer ball model (Niedermeier et al., 2011b).According to the singular hypothesis, the number of active sites, n s (T ), is dependent on T only and neglects ice nucleation kinetics.We suggest that further analytical efforts regarding ice nucleation kinetics can improve our understanding on the governing parameters of immersion freezing.
The immersed ISA per droplet is important for experimental derivation of J het and for deriving empirical quantities such as n s (T ) or fitting functions and their parameters.In previous experimental studies, droplets for ice nucleation experiments were dispensed from a bulk solution containing INPs (Broadley et al., 2012;Rigg et al., 2013;Wright and Petters, 2013;Herbert et al., 2014;Diehl et al., 2014).In other investigations, solid particles were size selected by their electrical mobility and then injected into, or continuously flown through, an ice nucleation chamber where water condensation precedes ice nucleation (Archuleta et al., 2005;Niedermeier et al., 2010;Kulkarni et al., 2012;Welti et al., 2012;Wex et al., 2014).In these studies and those that used polydisperse aerosol (e.g., Niemand et al., 2012), surface area calculations assumed that particles with the same mobility diameter are spherical with identical surface area.Despite this assumption, advancement in accounting for particle size variability considering multiple charged particles in ice nucleation experiments has been made (Lüönd et al., 2010;Augustin-Bauditz et al., 2014).However, extensive theoretical and experimental literature exists on aerosol sizing instrumentation and morphology characterization, which consider particle density, void fraction, shape and electrical charge effects implying their non-sphericity (DeCarlo et al., 2004;Slowik et al., 2004;Zelenyuk et al., 2006;Schmid et al., 2007;Park et al., 2008).In general, neglecting these effects may influence surface area estimates.Surfaces of particles may differ depending on different generation techniques, and previous studies have made progress in understanding why these do or do not result in differences in ice nucleation efficiency for particle types such as hematite (Hiranuma  (Lüönd et al., 2010), KGa-2-and KGa-1b-type kaolinites (Pinti et al., 2012), NX illite (Pinti et al., 2012;Hiranuma et al., 2015), or Arizona test dust (ATD) (Marcolli et al., 2007;Niedermeier et al., 2011a).In the present study, we consider laboratory-generated particles of ATD (Wright and Petters, 2013), NX illite (Broadley et al., 2012;Diehl et al., 2014), KGa-1b kaolinite and K-feldspar (Herbert et al., 2014), Fluka-type kaolinite (Wex et al., 2014) and natural dust (Niemand et al., 2012) each with their own respective generation method and argue that accurate quantification of ISA is crucial to discriminate how surface properties affect ice nucleation efficiency.Frequently, distributions of immersed ISA per droplet are typically assumed to be monodispersed, or in other words, each droplet is assumed to contain identical ISA.Furthermore, the number of droplets applied in an ice nucleation experiment may also affect the significance of the freezing data and thus interpretation of the experiment.It is necessary to question if a potential variability in ISA and/or the assumption of monodisperse ISA and a limited number of observed freezing events become important for interpretation of immersion freezing experiments with subsequent ramifications for the analytical ice nucleation description.
We introduce a newly developed model simulation in which ice nucleation is treated explicitly as a stochastic process applicable for isothermal and cooling-rate experiments.Previous experimental results using different experimental methods are simulated and compared for a wide range of atmospherically relevant conditions.However, this analysis is applied to laboratory-generated particles only and may not be applicable to field or natural samples because of the difficulty to separate INPs from others.Sensitivity studies on frozen fraction data and experimentally derived J het are performed as a function of ISA assumptions, the number of droplets employed in the experiment, T , and RH.The validity of typical assumptions of ISA variability and uncertainty are tested.Then, a detailed analysis of the ability of the model simulation to reproduce experimental results with strict uncertainty estimation is presented for seven independent immersion freezing studies utilizing eight different instrumentation: (i) droplets on a cold stage exposed to air, (ii) droplets on a cold stage covered in oil, (iii) oil droplet emulsions, (iv) droplet acoustic levitation, (v) droplet wind tunnel levitation, (vi) the Leipzig aerosol cloud interaction simulator (LACIS), (vii) a continuous-flow diffusion chamber (CFDC) and (viii) the aerosol interaction and dynamics in the atmosphere (AIDA) cloud chamber.A rigorous uncertainty analysis of the ice nucleation kinetics for typical ranges in experimental conditions is presented and discussed for laboratory application.

Simulation of isothermal freezing experiments
Stochastic immersion freezing simulations (IFSs) are performed to evaluate the effect of variable ISA on droplet immersion freezing experiments conducted in the laboratory.
As discussed above, different droplets in a laboratory experiment will possess different ISA.To account for this fact, ISA in each simulated droplet is sampled from a distribution to mimic this variability.Surface area can be any real positive value and can change by orders of magnitude.For this reason, a lognormal distribution can be assumed with the most probable ISA being A g or a mean distribution parameter µ = ln(A g ).The distribution width parameter is σ = ln(σ g ), where σ g represents the factor by which ISA can vary.A different distribution can also be assumed with knowledge of experimental methods used in particle generation, e.g., assuming a bipolar charge distribution for electrical mobility diameter selected particles.Knowledge of ISA for each droplet can be directly used as an alternative without a need for random sampling.Droplet freezing for isothermal experiments can then be described by where δN ufz represents the change in the number of unfrozen droplets after a certain interval of time, δt, and J het is the heterogeneous ice nucleation rate coefficient.The total available ISA is A tot = A j , where A j is the ISA in the j th droplet.
An assumption typically made is that all droplets contain the same ISA, or A tot = A g N ufz , where A g is the ISA for all droplets (e.g., Marcolli et al., 2007;Lüönd et al., 2010;Niedermeier et al., 2010;Murray et al., 2011;Rigg et al., 2013).Using this assumption and assuming a continuous differential in Eq. ( 1) leads to Integrating Eq. ( 2) further results in the commonly used expression for the fraction of frozen droplets, The form of the expression given in Eq. ( 3) is used in many studies although modified slightly when considering multiple components or contact angle distributions (e.g., Niedermeier et al., 2010;Murray et al., 2011;Broadley et al., 2012;Rigg et al., 2013), and when particle or droplet sizes are discretized or binned (e.g., Marcolli et al., 2007;Lüönd et al., 2010).The major weakness of this exponential form to describe f frz lies entirely in the assumption it is based on; i.e., it is only valid if the ISA is exactly the same for all droplets considered.When taking into account individual droplet ISA for all droplets, this formulation is not valid.Thus, application of this formula to interpret ice nucleation studies, or use in mathematical frameworks, strictly speaking, is also invalid when ISA on a droplet per droplet basis is different.
The ISA in a single droplet is a measurable quantity with a corresponding measurement uncertainty.It is unlikely that every droplet prepared in an immersion freezing experiment has identical ISA.For the same particle type, there will exist a systematic ISA uncertainty with respect to a particular droplet preparation technique.This systematic uncertainty is σ g and can be determined by directly measuring ISA in a population of independently prepared droplets.Since the ISA variability may not be typically resolved in previous experiments, a droplet freezing simulation must be employed to model ice nucleation for interpretation purposes.To accomplish this, freezing of each single droplet is assumed to be stochastic, or in other words, there exists a probability of the j th droplet to freeze, P j,frz , within δt.The probability for a single droplet not to freeze, P j,ufz , is realized as an exponential decay law (Pruppacher and Klett, 1997;Koop et al., 1997) and therefore P j,frz = 1 − P j,ufz = 1 − e −J het A j δt . (4) A time and surface area-dependent immersion freezing process, which follows CNT, is assumed and as a result, all simulations employ J het having units of cm −2 s −1 .However, J het does not explicitly depend on time and ISA, but on T and a w .A droplet can either remain in an unfrozen state or freeze and therefore, is described exactly by a binomial distribution, B(k; n, P j,frz ), with parameters P j,frz given by Eq. ( 4) and n = 1 meaning that only one trial is given for an individual droplet to freeze in δt.A randomly sampled number, k = 0 or 1, is obtained from the distribution B(k; n = 1, P j,frz ) = P k j,frz (1 − P j,frz ) 1−k (5) for each droplet with a normalization prefactor, n!/(k!(n − k)!) = 1.When k = 1, freezing occurs for the j th droplet and if k = 0, the droplet does not freeze and another k is sampled in the next time interval.For a collection of multiple droplets, the number of freezing events that occur in a given time interval is n frz and the cumulative sum as a function of time is N frz (t).For a single IFS starting with N tot liquid droplets, the fraction of unfrozen droplets is A record of n frz and corresponding droplet ISA, i.e., A j , is kept for a single IFS.This record can be thought of as a simulated experimental immersion freezing data set; i.e., it gives a record of droplet freezing time while tracking A j .Due to the stochastic nature of nucleation, repetition of isothermal IFSs will not result in identical values of f ufz over t.Likewise, repetition of a laboratory experiment will not result in exactly the same f ufz (t) curve.Therefore, it is necessary to repeat the simulations in order to reveal a range of f ufz (t) values of which the mean unfrozen fraction, f ufz (t), can be derived from all simulations.We choose an ensemble of 10 5 IFSs to accurately determine f ufz (t).This procedure is a basic form of a Monte Carlo method and yields upper and lower bounds between 5th and 95th percentiles serving as a stochastic uncertainty of the immersion freezing process.We define stochastic uncertainty as the scatter in the data due to the occurrence of random (i.e., stochastic) freezing events upon repeat experiments as a result of a set number of observed freezing events.
An ensemble of IFSs, referred to as a model simulation, requires the selection of parameters, e.g., N tot , A g , σ g , and J het .For demonstration purposes, the parameter choice is arbitrary.However, when reproducing a laboratory derived data set, a parameter selection process is applied.Parameters that can be directly accessed from previous laboratory studies are first selected to mimic experimental conditions.For example, if a study reports that 100 droplets were examined in an immersion freezing experiment, then N tot = 100.Some previous studies report only average ISA per droplet, A avg , and neglect information for estimating σ g .If A avg is reported as 7.1 × 10 −6 cm 2 , then for simplicity we set A g = 7.1 × 10 −6 cm 2 .For all studies in which a parameter or multiple parameters are not available or readily calculated, the model derived f ufz or f frz are fitted simultaneously to experimentally derived f ufz or f frz , and either critically assessed whether or not the parameter best reproduces experimental conditions or the fitted parameter value is compared with independently derived values from other published literature.We define "model derived" to refer to calculated frozen fraction, unfrozen fractions or J het values, which are not model input parameters.Details about the selection of parameters or whether or not parameters are or are not fitted for each simulation will be discussed in Sect.3. In many isothermal immersion freezing laboratory studies, droplet freezing continues over time when all other conditions remain constant, i.e., at constant T (Wright and Petters, 2013;Murray et al., 2011;Broadley et al., 2012;Herbert et al., 2014).Therefore, the J het parameter is selected to be constant for isothermal IFSs.

Experimentally derived J het for model input
When a cooling rate is applied in model simulations, droplet freezing is simulated in discrete temperature intervals and therefore J het at every step is required for deriving P j,frz .In this study, only water droplets are considered and therefore, it is assumed that a w = 1.0 and J het becomes a function of T only.Ideally, experimentally derived J het (T ) should be used for prediction of immersion freezing.However, these data sets are usually limited in T range and are discrete in nature.Knopf and Alpert (2013) where m and c are slope and intercept parameters, respectively, and a w is the independent variable following the formulation of Koop et al. (2000).The a w at which a droplet freezes is calculated by subtracting the a w of the droplet (i.e., 1.0 for pure water) from the water activity point that falls on the ice melting curve, a w, ice (T ), at the same temperature or where and p ice and p • H 2 O are the vapor pressure with respect to planar ice and water, respectively (Murphy and Koop, 2005).
Resulting calculations from Eqs. ( 6) to ( 8) are not computationally demanding and conveniently derive J het (T ) for model input.Details about parameter selection, i.e., m and c, the treatment of ISA variability, or whether or not parameters are or are not fitted will be discussed in Sect.3.

Simulated droplet freezing
Cooling-rate-dependent IFSs are performed to evaluate the effect of stochastic freezing and variable ISA in laboratory immersion freezing experiments.Again, the ISA for a single droplet is sampled; however, Eqs.(1) and (4) are modified to and respectively, where δT is a temperature interval and r = δT /δt is the cooling rate.J het (T ) is calculated from Eq. ( 6) and used in Eq. ( 10).Once the probability for the j th droplet to freeze is calculated for all droplets, freezing is determined by sampling from B(k; n, P j,frz ) (Eq. 5).The number of freezing events that occur in a given δT is n frz , and the cumulative sum as a function of T is N frz (T ) and used to calculate frozen fractions of droplets, f frz (T ) = N frz (T )/N tot .Similar to isothermal freezing, a single r-dependent IFS yields a droplet immersion freezing record analogous to an experimental data set.In this case, the record of droplet freezing and corresponding A j is a function of T .The average frozen fraction for 10 5 simulations, f frz (T ), is calculated along with 5th and 95th percentiles used as a stochastic uncertainty.
It is important to note that application of r-dependent IFSs presented here do not require the ABIFM, as it is only used as a parameterization of previously published immersion freezing data sets to calculate J het (T ).Any other published J het (T ) will work equally as well.The ABIFM parameterization is INP-type dependent and suitable for saturated and subsaturated conditions, i.e., a w ≤ 1, or RH ≤ 100 %, if the droplet is in equilibrium with the water vapor phase.Therefore, the ABIFM is a useful and convenient tool for model input J het (T ).

Isothermal model simulations of individual droplet freezing experiments
Figure 1a shows 5th and 95th percentile bounds of f ufz from four model simulations for different N tot applying either uniformly equal (σ g = 1) or lognormally distributed (σ g = 10) ISA per droplet as given in Table 1.Two of these test cases, Iso1 (isothermal test case 1) and Iso2, have uniform ISA both resulting in f ufz (on a logarithmic scale) linear with t.
However, the spread of the 5th and 95th percentile bounds is much wider for Iso2 having N tot = 30 than for Iso1 having N tot = 1000.It is clear that a larger spread in simulated f ufz is entirely due to applied smaller N tot .This implies that a laboratory experiment using a small N tot , is statistically less significant compared to an experiment with greater N tot .A single experimentally derived f ufz curve under the same conditions as Iso2 will fall anywhere between the upper and lower bounds, and thus may even appear to deviate from a loglinear relationship over time.Therefore, interpretation about the nature of the heterogeneous ice nucleation process from the slope of f ufz over time for an experiment using small N tot should be conducted with care.
Model simulations Iso3 and Iso4 are shown in Fig. 1a where N tot = 1000 and 30, respectively, and the ISA per droplet is sampled from lognormal distribution with σ g = 10.In Iso3, f ufz significantly deviates from a log-linear relationship with t.In Iso4, the same curvature exists; however, the percentile bounds are much wider due to applied smaller N tot .It is important to note that J het is the same and constant for all simulations shown in Fig. 1a.The nucleation rate of each j th droplet can be calculated as, ω het,j = J het A j with units of s −1 .The droplets having a larger or smaller ISA will result in larger or smaller ω het,j , respectively.The fact that f ufz is linear for σ g = 1, the curvature effect in f ufz seen for σ g = 10 must be entirely due to ISA variability.This is because droplets with greater values of ω het,j will tend to nucleate more rapidly than those having smaller ω het,j values.In other words, the curvature of f ufz (t) is entirely due to those droplets having larger and smaller ISA that freeze within shorter and longer timescales, respectively.In addition, the spread in the 5 and 95 percentiles is very similar for Iso1 and Iso3, and for Iso2 and Iso4.This is seen most clearly at the intersection of the blue and green shaded regions (t 1.3 min).In isothermal freezing experiments, variability in ISA will  b A uniform probability density function (U-pdf) was used to define the surface area distribution centered at 2.6 × 10 −7 cm 2 , with distribution endpoints at 9.4 × 10 −8 and 7.5 × 10 −7 cm 2 .See text and Fig. S3 in the Supplement for further details.c A multiple charge distribution (MCD) was used to define the surface area distribution.See text and Fig. S9 for further details.d Isothermal simulations were performed at 0.15 K increments within the stated temperature range.e Values of J het are calculated from the water activity, a w , based immersion freezing model (ABIFM) where m = 53.32 and c = −8.61(Knopf and Alpert, 2013).
not significantly affect stochastic uncertainty estimates, but will cause f ufz (t) to deviate from a log-linear relationship.
In some previous experimental isothermal immersion freezing studies, the number of liquid droplets and an estimate of the average ISA per droplet are provided or can be derived.However, the validity of the assumption that all droplets possess the same ISA is rarely investigated or quantified.Similarly, J het is not often reported.However, laboratory data do provide an opportunity to test our model for robustness while using parameters similar to those reported in the experimental studies.In fact, our model can also provide estimates for parameters typically unreported or unavailable, such as J het and σ g .
Experimental data by Wright and Petters (2013) for isothermal immersion freezing by ATD is very well reproduced by model simulation IsoWR as demonstrated in Fig. 1b.Parameters for IsoWR are given in Table 1 and chosen to mimic experimental conditions in which droplets contained 1 wt % ATD held at 251 K. Bounds at 5th and 95th percentiles of simulated f ufz are shown in Fig. 1b and envelop the laboratory data.A repeat experiment by Wright and Petters (2013) should result in a f ufz curve falling within the percentile bounds 90 % of the time when considering only stochastic uncertainty.
To further evaluate the validity of the simulations, the parameters used are compared with experimental conditions given in Wright and Petters (2013).IsoWR uses N tot = 1000, which agrees with the reported range of 300-1500.Parame-ters A g , σ g and J het used in IsoWR were fitted simultaneously so that the average f ufz from 10 5 simulations (or f ufz ) best reproduced observed f ufz .The first parameter in question is σ g = 9.5, which can be interpreted as a systematic standard error in ISA due to the experimental methods of generating or dispensing droplets containing ATD acting as INPs.We note this is different from an absolute ISA measurement error.Wright and Petters (2013) emulsified a mixture of oil and a bulk solution of water and ATD particles to form droplets with diameters of 50-250 µm.The variability in ISA should scale directly with the variability in droplet volumes (i.e., a factor of 125 or just over 2 orders of magnitude).Additional uncertainty will certainly arise from the variability in ATD particle numbers and the variability in ATD particle size.The surface area distribution and a random sampling is given in Fig. S2 for the IsoWR model simulation.A factor of 125 in range in ISA is a lower estimate of uncertainty, but already accounts for about 75% of the total sampled droplets.While not directly defined by Wright and Petters (2013), we are confident that the overall range in ISA should be well over 2 orders of magnitude and therefore, σ g = 9.5 is a reasonable value for the lognormal distribution width parameter employed in the simulations in Fig. 1b to reproduce the experimental data.The third parameter in question is A g = 6.4 × 10 −3 cm 2 .Unfortunately, an average ISA was not reported by Wright and Petters (2013) Table 1.
tion method (Brunauer et al., 1938).It is important to note that surface area measurements are not unambiguous due to the fact that heterogeneous ice nucleation may involve layers of water molecules interacting with surface molecules (Cox et al., 2013).The BET technique is one of many in which particle surface area is measured, and can be used to represent molecular available surface area.Bedjanian et al. (2013) reported SSA for ATD used in Wright and Petters (2013) as 85 ± 10 m 2 g −1 .The ISA per drop can then be estimated from the drop volume, V drop , and the density of water, ρ w , using the equation Considering only the variability in V drop , average ISA per drop should range between 5.5 × 10 −4 and 7.0 × 10 −2 cm 2 .The A g parameter in model simulation IsoWR falls within this range.
Finally, J het for ATD in water droplets was investigated by Pinti et al. (2012), who reanalyzed ATD immersion freezing data by Marcolli et al. (2007) but did not report J het values.However, estimates can be made following Knopf and Alpert (2013) accounting for f frz = 0.01 and a nucleation time assumed to be 1 s, which yields J het ranging from 5 × 10 6 to 1 × 10 2 cm −2 s −1 between T = 247.4and 252.8 K, in reasonable agreement with J het = 2.6 × 10 3 cm −2 s −1 used in IsoWR at 251 K.
The new model simulation presented here based entirely on CNT can describe freezing experiments by Wright and Petters (2013) accounting for long nucleation timescales and a large number of droplets considering variability in ISA.In addition, all crucial parameters applied are experimentally supported, in particular J het which is in agreement with independent studies (Marcolli et al., 2007;Pinti et al., 2012).Therefore, the isothermal immersion freezing data set of Wright and Petters (2013) may be explained by a time and ISA-dependent stochastic freezing process, in which each droplet contains variable ISA.More experimental investigation and model analysis should be conducted to verify their agreement at different temperatures, timescales and surface areas.Droplet to droplet variability in ice nucleation efficiency is typically parameterized with a variable efficiency of sites to nucleate ice, as done successfully in Wright and Petters (2013) or different contact angles (e.g., Niedermeier et al., 2011b;Broadley et al., 2012).Droplet to droplet variability parameterized in these ways and employing identical ISA can result in a deviation of f ufz from a log-linear relationship, similar to what is seen in Fig. 1.However, using a known ISA variability (Broadley et al., 2012;Wright and Petters, 2013), we reveal that the observed deviation from a log-linear relationship can be accounted for entirely by the ISA distribution.This may imply that the droplet to droplet variability in ice nucleation efficiency is entirely due to variable ISA.
Figure 2 shows results of isothermal freezing experiments by Broadley et al. (2012) for illite compared to model simulation IsoBr and experimental results by Herbert et al. (2014) for the INP types kaolinite and feldspar compared to model simulations IsoHE1 and IsoHE2, respectively (see Table 1).The experimental data and f ufz for all model simulations are in agreement and fall within the percentile bounds.Notice that the scatter in the isothermal immersion freezing data points is much larger than for Wright and Petters (2013) shown in Fig. 1b.As previously discussed, this is entirely due to a smaller number of droplets used in the laboratory experiments by Broadley et al. (2012) (N tot = 63) and Herbert et al. (2014) (N tot = 40) and thus, may be entirely attributed to the stochastic nature of immersion freezing as expected by CNT.The model simulations capture this effect by producing a wide range in f ufz .Only one experiment was performed for each of the laboratory data sets presented in Fig. 2, and if these experiments were repeated, f ufz values would very likely not be the same and may even exhibit either more lin- Orange lines and shading represent f ufz with corresponding 5th and 95th percentile bounds, respectively.Parameter values for model simulations are given in Table 1.
ear or curved behavior with time.Repetition of experiments should provide better estimates of f ufz and σ g , but for any single experiment, f ufz may still fall within the given percentile bounds.In other words, additional experiments would better define the mean of f ufz and the uncertainty in the mean of f ufz , but will not decrease the uncertainty bounds.Only by using more droplets, e.g., Wright and Petters (2013), would a single experiment be more statistically significant.
The selection of parameters and ISA distribution used in IsoBR are discussed.N tot = 63 applied in Broadley et al. (2012) is used here for IsoBR.In the analysis of the experiment "run 20" simulated by IsoBR, the authors sub-selected droplets 10 − 20 µm in diameter from a droplet population.The droplet volume and the ISA variability should scale to the third power of the droplet diameter, i.e., by a factor of about (20/10) 3 = 8.Therefore, the simulated ISA is assumed to follow a uniform probability density function between 9.4 × 10 −8 and 7.5 × 10 −7 cm 2 .This ISA range spans a factor of ∼ 8 with a geometric average of 2.65 × 10 −7 cm 2 as reported in Broadley et al. (2012).We note that a factor of 8 is a lower limit of variability as any additional uncertainty in illite particle size distribution or the numbers of illite particles per droplet is not considered.The parameter J het = 2.82 × 10 3 cm −2 s −1 was not fitted, but instead selected in such a way that this single value resulted in f ufz data falling entirely within the stochastic uncertainty limits.This value is in agreement with the previous ABIFM parameterization, J het = 1.25 × 10 3 cm −2 s −1 , at the same temperature and water activity (Knopf and Alpert, 2013).We note that the decay of simulated f ufz over time appears linear in Fig. 2a although the experimental data appear curved.This curvature led Broadley et al. (2012) to assume an active site model to describe the data.However, we find that when using too small numbers of droplets, experiments may be too uncertain to make any solid claims about the nature of the ice nucleation process.The stochastic uncertainty bounds in Fig. 2a are sufficiently large that the data is still in agree-ment with the model simulations presented here.The IsoBR simulation demonstrates that freezing due to illite can also be described by a stochastic freezing approach with one value of J het .Thus, laboratory derived isothermal immersion freezing of illite may be explained by CNT accounting for the stochastic nature of immersion freezing and variability in ISA.
In the model simulation IsoHE1, parameters A g , σ g and J het are fitted to the experimental data.The parameter A g = 1.2 cm 2 , is in good agreement with experimentally derived A g = 2.4 cm 2 reported in Herbert et al. ( 2014), for kaolinite using SSA = 11.8 m 2 g −1 (Murray et al., 2011), 1.0 wt % concentration and V drop = 1 µL.Herbert et al. ( 2014) did not report sufficient information to estimate an overall variability in ISA; therefore, comparison of σ g to experimental conditions is difficult.As previously discussed, a repeat experiment may result in f ufz exhibiting more linear or non-linear behavior with t within the calculated percentile bounds, i.e., within the stochastic uncertainty.Figure 1a shows that a more linear or non-linear relationship of f ufz with t implies a smaller or larger value of σ g .We note that due to the lack of quantitative information about the variability in ISA, the assumed lognormal surface area distribution may be over or underestimated.Performing more experiments or employing a larger number of droplets will decrease the stochastic uncertainty and better constrain the curvature of f ufz over time.The ABIFM yields J het = 1.75 × 10 −2 cm −2 s −1 at T = 255.15K and a w = 1.0, which is within an order of magnitude of J het used in IsoHE1.The agreement between simulated and experimental parameters implies that CNT may be able to explain observed immersion freezing of kaolinite when variable ISA and stochastic uncertainty is considered.Herbert et al. (2014) and Murray et al. (2011) came to the conclusion that this particular type of kaolinite, KGa-1b, is a "single-component system", which means that a single J het function of T can reproduce the experimental data.Our model simulations lead to the same conclusion, and the derived J het value is in agreement with the independently  1.Those that are fitted are A g , σ g and J het .Average ISA for the data in Fig. 2 is 1.85 × 10 −2 cm 2 , similar to A g = 2.0 × 10 −2 cm 2 used in IsoHE2.Droplets used in Herbert et al. (2014) were dispensed with a digital micropipet with high accuracy; thus, it can be expected that the contribution of droplet volume variability to the σ g parameter is low.However, it is impossible to make any estimate of σ g for comparison with the fitted σ g due to the lack of quantitative experimental information about ISA variability.As in simulation IsoHE1, a more linear or non-linear f ufz curve may imply that our assumed lognormal distribution width is over or underestimated, respectively.However, it will be demonstrated that the same σ g can be used to reproduce the cooling-rate-dependent experiments of feldspar, which gives confidence that this σ g value may be appropriate.To better constrain σ g , more stochastic certainty is required by application of more droplets or conducting multiple experiments.Values of J het for feldspar independent from Herbert et al. ( 2014) in the same temperature range to our knowledge do not exist making comparison difficult.
Model simulations IsoDI1-3 of isothermal immersion freezing experiments by Diehl et al. (2014) for illite in wind tunnel levitation experiments are shown in Fig. 3. Simulation parameters are given in Table 1.Only one droplet was observed in each experiment, and approximately 45 experiments were conducted for each of the three data sets shown in Fig. 3.This is equivalent to one experiment with N tot = 45 droplets, since droplet freezing is independent of the freezing of other droplets.Excellent agreement is observed between simulated and experimental f ufz .The parameter A g for the three simulations match the average surface area per drop reported in Diehl et al. (2014).Parameters σ g , J het (−18 • C) and J het (−21 • C) are simultaneously fitted to experimental data.It can be expected that σ g is the same for all three simulations, due to the fact that Diehl et al. (2014) likely used identical bulk water-illite stock solutions.Therefore, a single fitted value of σ g = 3.2 is used for all simulations.IsoDI2 and IsoDI3 simulate data taken at T = −21 • C and use the fitted parameter J het (−21 • C) = 1.0 × 10 0 cm −2 s −1 .IsoDI1 simulated data at T = −18 • C uses the parameter J het (−18 • C) = 1.8 × 10 −2 cm −2 s −1 .At T = −18 and −21 • C, the ABIFM yields J het = 1.8 × 10 −2 and 2.6 × 10 −1 cm −2 s −1 , respectively, and is in excellent agreement with derived values in IsoDI1-3.An adequate constraint of σ g could not be established due to a lack of information about the ISA variability.However, it is evident that the fitted σ g value may be justified due to the fact that the same value reproduced all three isothermal data sets.We find that a time-dependent and stochastic immersion freezing process may reconcile observations when variable ISA is considered.Depending on ISA variability, trajectories of model derived f ufz over time are significantly altered and thus assuming identical ISA may not be valid.It is well known that immersion freezing depends on surface area; i.e., an increase in ISA translates to an increase in nucleation rate.However, we note that variability in both t and ISA equally affect calculations of droplet freezing probabilities (Eq.4) used in model simulations, and therefore neglecting time dependence will cause erroneous interpretation of immersion freezing data to the same degree as if the surface area dependence is neglected.This simple stochastic immersion freezing model accounting for ISA variability can explain the isothermal ice nucleation data of various experiments without invoking empirical parameterizations, assumptions of particle surface composition, and/or other modifications in parameters and interpretations.

Cooling-rate model simulations of individual droplet freezing experiments
Cooling-rate IFSs were performed to investigate the effects of variable ISA and N tot on experimentally derived J het and f frz as a function of T .For a single cooling-rate IFS, variable ISA per droplet is applied and used to calculate P j,frz from Eq. ( 10) in discrete temperature steps, δT , and then Eq. ( 5) simulates freezing.The IFS stops after some T or when all droplets freeze, and the simulated freezing record is kept detailing which droplets froze or remained liquid at each T and their corresponding ISA.This is analogous to running an immersion freezing experiment in a laboratory setting and recording the observed number of frozen droplets or ice crystals as a function of T .
The simulated freezing record is treated as a freezing data set from which the assumption of identical ISA can be tested.This is accomplished by re-calculating J het from the simulated data.These (re-)calculations use n frz , the length of the time interval, δt = δT /r, and either of two different approaches in determining A tot .For the first approach, A g is assumed to be identical for all droplets, i.e., without the knowledge that immersion freezing was simulated for droplets with variable ISA in the first place.This is equal to assuming a monodisperse INP population in laboratory immersion freezing experiments resulting in an "apparent" heterogeneous ice nucleation rate coefficient, J apparent het (T ), calculated by where n ufz (T ) is the number of unfrozen droplets at T and A tot = n ufz A g .The second approach accounts for the variable ISA present in droplets resulting in the "actual" heterogeneous ice nucleation rate coefficient, J actual het (T ), calculated by and A tot = A j is the total surface area contribution from droplets that remain liquid.Comparing results from Eqs. ( 11) and ( 12) allows for evaluation of the assumption that all droplets have the same ISA, when they actually do not.In this way a null hypothesis is considered, that is if J apparent het (T ) and J actual het (T ) are the same, then the assumption of identical ISA is valid.
Poisson statistics are used to derive upper and lower fiducial limits of J apparent het (T ) and J actual het (T ) at x = 0.999 confidence for n frz following Koop et al. (1997).The upper fiducial limit of the heterogeneous ice nucleation rate coefficient, J up het , accounts for additional freezing events occurring with a probability of x, than observed n frz .Likewise, a lower fiducial limit of the heterogeneous ice nucleation rate coefficient, J low het , accounts for less than the observed n frz occurring with a probability of x.We refer to the upper and lower limits of n frz as n up frz and n low frz , respectively (Koop et al., 1997).The fiducial limits of J apparent het and J actual het for a single simulation can be calculated using Eqs.( 11) and ( 12), but replacing n frz with n up frz or n low frz , respectively.Each simulation results in different J apparent het and J actual het values and different fiducial limits at the same T due to random sampling, therefore, averages are reported.
Figure 4 shows the results of two model simulations, Cr1 and Cr2, having r = 0.5 and 5.0 K min −1 , respectively.For all 10 5 IFSs, J  6) for illite (Knopf and Alpert, 2013).Frozen droplet fractions, f frz , are shown in (c) and (d) where dashed lines and shadings represent f frz and 5 and 95 % bounds, respectively.Parameter values for Cr1 and Cr2 are given in Table 2.
as dashed lines, respectively, along with corresponding f frz curves displayed in Fig. 4c and d.The parameterization of J het (T ) for illite dust (Knopf and Alpert, 2013) with m = 54.5 and c = −10.7 used in Eq. ( 10) for each simulation is shown as the red line in Fig. 4a and b and referred to as the model input J het .Simulation parameters for Cr1 and Cr2 are given in Table 2.
According to CNT, two immersion freezing cooling-rate experiments conducted at different r should result in identical J het values due to the fact that J het is independent of r.CNT is violated if significantly different J het values are derived at different r. Figure 4a shows that values of J apparent het are significantly different between model simulations Cr1 and Cr2.Also J apparent het is overestimated at higher freezing temperatures and underestimated at lower freezing temperatures compared with model input J het (T ).These significant differences do not support the null hypothesis and imply that when experimentally deriving J het , the assumption that ISA per droplet is identical is invalid.Figure 4b shows that accounting for variable ISA, J  shows a cooling-rate dependence; i.e., when the cooling rate varies by an order of magnitude, J apparent het evaluated at one temperature is also about an order of magnitude different.Values of J actual het on the other hand, are not dependent on cooling rate.We reiterate that J apparent het and J actual het are calculated from the same simulated freezing record having lognormally distributed ISA.The only difference between "actual" and "apparent" J het is the surface area assumptions used in their respective calculation (Eqs.10 and 11); i.e., J apparent het intentionally assumes identical ISA as done commonly in experimental analysis and J actual het accounts for the variable ISA.Thus, the apparent cooling-rate dependence in simulations Cr1 and Cr2 is a direct result of assuming identical ISA.This is the case when a broad ISA distribution is simulated; i.e., σ g = 10.When the ISA distribution is very narrow, uniform or σ g is about 1.0, then J actual het will equal J apparent het .This is demonstrated in Fig. S7 by model simulations Cr3 and Cr4 using σ g = 1.If an experimental study succeeds to create a narrow enough distribution, then assuming identical ISA may be applicable.
Towards warmer (T > 248 K) and colder temperatures (T < 238 K), the difference in upper and lower fiducial limits derived in Cr1 and Cr2 are much greater than for the midtemperature range (238 < T < 248 K).In fact the smallest difference occurs at f frz 0.5.This is because calculations are statistically more significant at the median freezing where n frz is largest.Fewer droplets freeze at the beginning and end of a cooling process resulting in a wide fiducial limit range reaching up to 4 orders of magnitude (Fig. 4a and b) in spite of the high number of droplets used (N tot = 1000).The corresponding percentile bounds of f frz shown in Fig. 4c and d do not reflect a considerable uncertainty compared to the upper and lower fiducial limits (Fig. 4a and b).It is important to note that f frz are identical in Fig. 4c and d, because surface area is not used to derive f frz .This analysis suggest that values and uncertainties of f frz are not suited to derive J het and any corresponding error.
Previous immersion freezing experiments by Herbert et al. ( 2014) are modeled in CrHE1 and CrHE2 where r = 0.2 and 2.0 K min −1 , respectively, for the case of feldspar acting as INPs.Herbert et al. (2014) used the same weight fraction of feldspar per droplet in isothermal and cooling-ratedependent experiments.Therefore, it is reasonable to suspect that the parameters should be the same or very close.The parameters σ g and A g are not fitted in model simulations CrHE1 and CrHE2 and instead are taken from the values fitted in IsoHE2.The parameters m and c used to calculate J het as a function of T following the ABIFM are fitted so that experimentally derived f frz is best reproduced by model derived f frz .Since Herbert et al. ( 2014) assumed identical ISA, experimentally derived J het can be directly compared with J apparent het from cooling-rate model simulations.Figure 5 shows experimentally derived f frz and J  6) to reproduce frozen fraction data (Fig. 5a) within 5th and 95th percentile bounds.The laboratory data falls within the percentiles and fiducial limits of f frz and J apparent het , respectively.Previous studies have been successful in interpreting immersion freezing studies to follow an active site approach by considering both cooling-rate and isothermal experiments (Vali, 2014;Herbert et al., 2014).Alternatively, we have found that our isothermal and coolingrate simulations based on a single function of J het (T ) for feldspar are in agreement with experimental results when the same A g and σ g parameters are used.This also gives confidence for the appropriateness of the model parameter val-ues.Although, we reiterate that the ISA distribution may be over or underestimated as the experimental data has a large stochastic uncertainty and thus fitting our model to the data is not well constrained.Figure 5b shows that J apparent het values for different r are not the same; however, they are in agreement with experimental data.Note that simulated and experimentally derived J apparent het for r = 0.2 and 2.0 K min −1 (5b) are different by about 1 order of magnitude at the same T .
Figure 5c shows J actual het with upper and lower fiducial limits derived from CrHE1 and CrHE2.For comparison, J het derived in model simulation IsoHE2 for feldspar is also shown and in agreement with J actual het within our stochastic uncertainty estimates.When accounting for variable ISA, J actual het are in excellent agreement with the ABIFM parameterization derived in this study for feldspar INPs.Furthermore, J actual het calculated for different r are identical as predicted by CNT, a similar finding as in the model simulations Cr1 and Cr2 (Fig. 4b).Therefore, J het (T ) used here can be considered a new J het ( a w ) parameterization for feldspar valid for 0.078 < a w < 0.120.
The differences between J apparent het and J actual het shown in Fig. 5b and c can be attributed to two reasons: (i) a potential misrepresentation of the slope of J het versus T and (ii) a potential misrepresentation of a dependence of J het on r.Regarding the slope of J het vs. T , note that droplets with ISA less than A g likely freeze at colder T compared to droplets with ISA greater than A g , that likely freeze at warmer T .However, assuming identical ISA equal to A g for all droplets either overestimates or underestimates the actual ISA present in droplets that freeze at colder and warmer temperatures, respectively.Due to the inverse relationship between A g and J het , calculations of J apparent het from Eq. ( 11) may then be underestimated and overestimated at colder and warmer temperatures, respectively.When comparing J apparent het against J actual het in model simulations (calculated from the same freezing record generated using lognormally distributed ISA), J apparent het is underestimated and overestimated at colder and warmer temperatures, respectively.Therefore, this analysis clearly shows that assuming identical ISA in each droplet may potentially lead to misrepresentation of the slope J het vs. T .Regarding a dependence of J het on r, we find that J apparent het for simulations CrHE1 (r = 0.2 K min −1 ) and CrHE2 (r = 2.0 K min −1 ) are different by about 1 order of magnitude at the same T .In Herbert et al. (2014), experimentally derived J het applying r = 0.2 and 2.0 K min −1 differ by about 1 order of magnitude at the same T .In separate model simulations not shown here, applying r different by 2 orders of magnitude yields J apparent het values that differ by 2 orders of magnitude.We note that this is the case for a wide distribution width of σ g = 8.5 and that the stochastic uncertainty in J het is large, not considering additional uncertainties, e.g., in temperature or in surface area and its variability.Nevertheless, experimentally derived J het data from Herbert et al. ( 2014) is in good agreement with  6) (Knopf and Alpert, 2013)  indicating a potentially erroneous cooling-rate dependence, which may be caused by the assumption of identical ISA.
Model simulations CrDI1 and CrDI2 of immersion freezing experiments by Diehl et al. (2014) for illite acting as INPs probed in acoustic levitation experiments are shown in Fig. 6.A non-linear r was used in Diehl et al. (2014) and was the same for both experiments and model simulations, but the ISA per droplet was varied.Diehl et al. (2014) reported an ISA per drop of 7.1 × 10 −1 and 7.1 × 10 −3 cm 2 in the two different sets of experiments.These surface areas are not the same as for isothermal experiments, and so we suspect that different stock solutions were prepared for cooling-rate experiments by Diehl et al. (2014).Therefore, parameters from IsoDI1 to IsoDI3 are not used in coolingrate model simulations and instead newly fitted values of simulation parameters σ g and A g for CrDI1 and CrDI2 were derived and are given in Table 2.A continuous function of J het was not fitted, but calculated using the ABIFM for illite (Knopf and Alpert, 2013).When using the A g values reported in Diehl et al. (2014) in conjunction with the other parameters, model simulations cannot reproduce experimental f frz .This is in spite of the excellent performance of IsoBr for reproducing droplet freezing initiated by illite from Broadley et al. (2012).In attempt to reconcile results from Diehl et al. (2014) with previous literature data (Broadley et al., 2012;Knopf and Alpert, 2013), model derived f frz are fit to experimental f frz yielding two different parameter values of A g = 2.94 and 2.91 × 10 −2 cm 2 used in CrDI1 and CrDI2, respectively.We note that fitted A g values differ only by a factor of 4 from values reported by Diehl et al. (2014) and therefore, are in reasonable agreement.However, calculated J apparent het values shown in Fig. 6b still use ISA of 7.1 × 10 −1 and 7.1 × 10 −3 cm 2 as reported by Diehl et al. (2014).
Figure 6a shows that simulated and experimental f frz are in agreement when accounting for ISA variability (σ g = 5.7).Experimental values of J apparent het displayed in Fig. 6b are in agreement with model derived J apparent het .This result is robust since experimental J apparent het data was not used in fitting f frz and the same value of the fitted parameter σ g is used in the two surface area-dependent cooling-rate experiments, which gives confidence that the ABIFM parameterization and the ISA distribution width are appropriate.Accounting for the actual variability in ISA used to simulate freezing, J actual het shown in Fig. 6c is in perfect agreement with the ABIFM parameterization for illite (Knopf and Alpert, 2013).For comparison, J het derived in model simulations IsoDI1-3 for illite are also shown and in agreement with J actual het within our stochastic uncertainty estimates.Again, the data and model supports a stochastic, time-dependent immersion freezing process and may describe laboratory data considering variable ISA.
A major inconsistency between experiment and simulation is shown in Fig. 6b, evident from the agreement of J  (Diehl et al., 2014).The red line in (b) and (c) is calculated from Eq. ( 6) for illite (Knopf and Alpert, 2013).Parameter values for CrDI1 and CrDI2 are given in Table 2. Fitted J het values from model simulations IsoDI1 to IsoDI3 are shown in (c) and their corresponding error derived from Fig. 9 considering N tot = 45, σ g = 3.2 and a temperature error, T = ±0.7 K (see Table 1).
for the case when different ISA but identical cooling rates are applied.According to CNT, J het is independent of surface area.This means that if two experiments are performed with different ISA but use the same r, J het should be the same as a function of T .However, simulated and experimentally derived J apparent het (T ) deviate by more than 1 order of magnitude as the surface area varies by two orders of magnitude.This would indicate that J apparent het values violate CNT, but this is the cause of assuming identical ISA.In fact, this freezing behavior also contradicts all surface-based empirical parameterization of immersion freezing, such as determining n s (T ), or the number of actives sites per particle surface area (Murray et al., 2012;Hoose and Möhler, 2012).This result impacts immersion freezing experiments conducted as a function of ISA that assume identical ISA in each experiment, thereby implicitly imposing a surface area dependence on J apparent het or n s (T ).Accounting for the experimental uncertainty and variability in ISA may reconcile experimental data.

Continuous-flow and cloud chamber immersion freezing experiments
Model simulations IsoCFDC and IsoLACIS (see Table 1) reproduce experimental results of Wex et al. (2014), who used two ice nucleation instrumentation, (i) a continuousflow diffusion chamber (CFDC) (Rogers et al., 2001;De-Mott et al., 2010) and (ii) the Leipzig aerosol cloud interaction simulator (LACIS) (Hartmann et al., 2011), to observe immersion freezing of 300 nm mobility diameter selected kaolinite particles as a function of T and RH > 100 %.
It is important to note that for both instruments, droplet freezing is not observed and instead is optically detected and at the LACIS outlet, a self-built optical particle spectrometer, TOPS-Ice (Clauss et al., 2013), determines if the arriving hydrometeors are liquid droplets or frozen ice crystals, resulting in the determination of a frozen fraction.Thus, f frz is calculated from the ratio between observed ice crystal and aerosol numbers per volume of air.The model simulation parameter N tot is derived from known experimental parameters, including residence time of the CFDC, t r = 5 s, aerosol flow rate, Q = 1.0 L min −1 , and kaolinite particle concentrations, N p = 10 cm −3 (Wex et al., 2014).By defining a single IFS over an interval of time equal to t r , N tot = N p Qt r = 833 particles per IFS.Similarly for LACIS, Q = 0.08 L min −1 , t r = 1.6 s, and N tot = 21 particles per IFS.Note that minimum f frz values for CFDC and LACIS presented in Wex et al. (2014) are approximately equal to 1/N tot .We run 1440 and 6000 isothermal IFSs for IsoCFDC and IsoLACIS, respectively, equivalent to 2 h averages as done in Wex et al. (2014).Simulation parameters for IsoCFDC and IsoLACIS are given in Table 1.
Figure 7 shows that simulated f frz for IsoCFDC and Iso-LACIS agree very well with CFDC and LACIS data by Wex et al. (2014).However, some data points fall outside   as indicated in the legend.Shadings in (b) correspond to upper and lower fiducial limits with x = 0.999 confidence and the red line is calculated from Eq. ( 6) for kaolinite (Knopf and Alpert, 2013).Parameter values for IsoCFDC and IsoLACIS are given in Table 1.
of the 5th and 95th percentiles (Fig. 7a), which may imply that a greater uncertainty exists that cannot be explained by a stochastic freezing process.This may be due, in part, to uncertainty in ice crystal detection, which is not accounted for in model simulations.The surface area for spherical 300 nm particles is A 300nm = 2.8 × 10 −9 cm 2 .However, the assumption that a kaolinite particle with an electrical mobility diameter of 300 nm is equal to a 300 nm diameter sphere is likely not true, due to shape irregularities, variable density, void fractions, multiple charges, and other geometries (De-Carlo et al., 2004;Slowik et al., 2004;Zelenyuk et al., 2006;Schmid et al., 2007;Park et al., 2008) with a tendency for greater surface area than assumed.Additionally, particles of larger diameter, and thus larger surface area, may have the same electrical mobility due to the presence of multiple charges.Therefore, a distribution of particle surface area can be expected.Following Wiedensohler and Fissan (1988), the probability for particles having multiple charges as a function of particle diameter, P (ln D p ), at a constant electrical mobility diameter of 300 nm is shown in Fig. S8.The distribution P (ln D p ) is a probability density function from which particle diameters are sampled in simulations IsoCFDC and IsoLACIS.Individual sampled particle surface area is calculated assuming spherical particles.We note that a lognormal distribution is not used in IsoCFDC and IsoLACIS.Instead, it is assumed that the ISA distribution varies only due to the theoretical multiple charge distribution.Parameters m and c used to calculate J het for Fluka kaolinite are fitted.
Calculations of J apparent het and J actual het assuming constant ISA equal to A 300nm or accounting for variable ISA, respectively, in IsoCFDC and IsoLACIS are shown in Fig. 7b.We find agreement between J apparent het , data by Wex et al. (2014) and J actual het , which accounts for multiple particle charges predicted by Wiedensohler and Fissan (1988).Within the uncertainties presented here, assuming that the electrical mobility diameter corresponds to the physical particle diameter and calculating surface area from a spherical geometry may be a valid assumption.Hence, ice nucleation experiments in which particles are mobility selected may be good examples of cases where ISA variability plays a minor role.Studies which use pre-impactors to remove larger-sized particles, i.e., the selected size is larger than the median size of the total size distribution (Wex et al., 2014;Augustin-Bauditz et al., 2014), may even have a more narrow size distribution than used here.On the other hand, a recent study by Hartmann et al. (2016) derived the numbers of multiple charges on 300 nm mobility diameter size selected particles using simultaneous measurements of cloud condensation nuclei activation curves and total particle counts with a condensation particle counter.The authors found that when utilizing a preimpactor, the multiple charge distribution of mobility diameter selected particles was larger than theoretical predictions (Hartmann et al., 2016).For comparison to the charge distribution used in IsoLACIS and IsoCFDC shown in Fig. S9, we also plot the distribution measured by Hartmann et al. (2016) for the same particle type (Fluka kaolinite) used in (Wex et al., 2014).Hartmann et al. (2016) claimed that when correcting for their measured multiple charge distribution in experiments, values of n s (T ) are shifted by 2 K.We note that J apparent het shown in Fig. 8b is about +1 K shifted from J actual het .This shift is smaller than observed in Hartmann et al. (2016) due to the fact that we applied the narrower theoretical distribution.Despite these issues, the model input J het represents a new parameterization for Fluka kaolinite where m = 53.32 and c = −8.61following the ABIFM applicable for 0.220 < a w < 0.305.and fiducial limits with x = 0.999 confidence, respectively.The red line in (b) is calculated from Eq. ( 6) using new parameters derived for natural dust.Parameter values for CrNI1 and CrNI2 are given in Table 2. Wex et al. (2014) presented a detailed immersion freezing analysis of various kaolinite particle sizes and types of coatings and found that both stochastic and active site approaches can be applied to describe the data.Simulating all these cases using our model simulations is beyond the scope of this paper; however, we are certain that model simulations, which use the same J het (T , a w ), will hold for all systems at all T and RH due to the prediction of immersion freezing kinetics (i.e., using J het ) being independent of experimentally applied ISA, particle size, and particle coating type (assuming the coating dissolves when water is taken up and does not react with the INP surface).These findings demonstrate that our new model simulations and the ABIFM are applicable for ice nucleation studies using a CFDC as previously shown by Knopf and Alpert (2013) and additionally LACIS.
IFSs are used to describe AIDA chamber immersion freezing experiments applying natural dust by Niemand et al. (2012) in model simulations CrNI1 and CrNI2.Among the different types of natural dust investigated, we choose 2 Asian dust experiments at −20.1 < T < −28.1 • C and −14.3 < T < −22.4 • C (see ACI04_19 and ACI04_16 in Tables 2 and 3 in Niemand et al., 2012).A continuous nonlinear cooling rate with time due to adiabatic expansion is fitted to experimental trajectories using a fourth-order polynomial function.In AIDA experiments water saturation is typically reached after cooling begins.To mimic this process, ice particle production in model simulations is allowed after 80 s of cooling (see Fig. 2 in Niemand et al., 2012).Ice crystal concentration in an aerosol sampling flow of 5 L min −1 , from the chamber is observed every 5 s using an optical particle counter (Benz et al., 2005); thus, a volume of 0.42 L of air is simulated.Total particle numbers in the simulated volume are on the order of 10 5 , which agree well with minimum reported f frz of about 10 −5 .Niemand et al. (2012) reported lognormal surface-size distributions with parameters, d S,median and σ g of polydisperse aerosol population.In CrNI1 and CrNI2, A j is derived by sampling particle diameters from the corresponding number-size distributions and assuming spherical particles.We note that a fitted lognormal distribution is not used in CrNI1 and CrNI2, due to the fact that reported size distributions are well defined.Sampling stops when A tot equals total surface area reported by Niemand et al. (2012).Experimentally derived J het is not available and so the ABIFM parameters m and c are fitted to experimentally derived f frz data.Model simulation parameters for CrNI1 and CrNI2 are given in Table 2.
Figure 8 shows simulated f frz and J actual het from CrNI1 and CrNI2 and the time evolution of simulated ice crystal concentration in CrNI1 observed during the experiments.Simulated f frz (Fig. 8a) fall within the experimental uncertainty reported by Niemand et al. (2012) and the scatter in the data for all dust types.Narrow 5th and 95th percentile bounds are attributable to large N tot on the order of 10 5 droplets per cooling simulation.Ice particle concentrations over time in CrNI1 are shown (insert in Fig. 8b) and capture the overall observed trend in observations.This is in spite of the fact that the observed time evolution of ice crystal numbers was not used for fitting parameters m and c. Figure 8b shows J actual het and upper and lower fiducial limits.As frozen fraction decreases the fiducial limits become broader ranging from 0.8 to 2.5 orders of magnitude.We conclude that our model simulations are suitable for describing laboratory immersion freezing in AIDA cloud chamber and further support the necessity of quantification of ISA variability in the derivation of ice nucleation kinetics.
Notice that in Fig. 8a, the vertical scatter in the experimental data increases at warmer T and for low f frz , which implies that uncertainty likely increases as f frz decreases.Since aerosol numbers and surface area in the experiments by Niemand et al. (2012) are relatively the same for the two experiments, decreasing f frz implies fewer detected ice crystals or decreasing numbers of ice nucleation events resulting in an increase in experimental uncertainty.Immersion freezing due to natural dust was parameterized using a deterministic (singular) approach, i.e., using n s (T ), which captured the trend in experimental results (Niemand et al., 2012).However, a deterministic approach for interpretation and analysis of ice crystal production, which inherently ignores stochastic freezing, cannot explain the increase in the data scatter for smaller f frz values at warmer T .These observations can be explained by a stochastic and time-dependent immersion freezing process.We note that other measurement uncertainties may exist which may not be captured either by a deterministic approach or by our model.However, we conclude that stochastic uncertainty is important to consider for future ice nucleation studies.The fiducial limits of J actual het shown in Fig. 8b, in fact, capture this effect of larger scatter as T increases implying the uncertainty in observed ice nucleation kinetics increases.Since the freezing efficiency of Asian dust was shown to be similar for Saharan, Canary Island, and Israeli dust (Niemand et al., 2012), the new ABIFM parameterization of J het (T , a w ) derived here is applicable for natural dust.

Simulation findings and uncertainty analysis
Our results strongly suggest that laboratory immersion freezing studies should provide accurate estimates of ISA variability in droplets.We find that simplified assumptions about ISA can result in misinterpretation and miscalculation of J het values.This includes assuming identical surface area, which implicitly imposes a dependence of J het on both ISA and r.Future laboratory immersion freezing studies should also consider the stochastic nature of ice nucleation following CNT and resulting uncertainties.When only a single ice nucleation experiment is performed or too few droplets are used, stochastic uncertainty can potentially be very large and may limit data interpretation.Once again, stochastic uncertainty refers to large or small expected data scatter from observing small or large numbers of freezing events, respectively.The surface area-based deterministic approach deriving n s (T ) is an alternative to calculating J het , but does not consider stochastic effects or effect of time in analysis of immersion freezing.By design, n s (T ) should therefore, not have any dependence on r.However, this is not supported as n s (T ) has been observed to be dependent on r for feldspar and kaolinite (Herbert et al., 2014).
The model simulation and laboratory data sets investigated here were performed for INPs immersed in pure www.atmos-chem-phys.net/16/2083/2016/water droplets.However, aqueous solution droplets having a w < 1.0 are frequently present in the atmosphere at supercooled temperatures and subsaturated conditions (i.e., RH < 100 %).The ABIFM (Eqs.6-8) inherently and accurately accounts for these conditions and thus provides a complete description of immersion freezing for laboratory experiments, as well as cloud models under atmospherically relevant T and RH.We suggest that future isothermal and cooling-rate-dependent immersion freezing studies investigate aqueous solution droplets in addition to water droplets (e.g., Archuleta et al., 2005;Alpert et al., 2011b;Wex et al., 2014), providing additional data sets to constrain ice nucleation kinetics and to validate and expand ABIFM and other parameterizations.
Uncertainty analysis is crucial for the interpretation of laboratory immersion freezing results.Here we present a quantitative uncertainty analysis of J het , by defining J het as the total uncertainty derived from individual contributions of statistical uncertainty due to N tot , temperature accuracy referred to as T , a w or RH accuracy referred to as RH, ISA variability expressed as σ g , and accuracy of measuring absolute surface area referred to as A g .This uncertainty analysis is applicable to both isothermal and cooling-ratedependent immersion freezing experiments.It is convenient to quantify J het in the form of a × ÷ error instead of a typical ± error due to J het varying exponentially over a linear range in T .If J het = 100 cm −2 s −1 with a factor of ±3 error for example, then J het = × ÷ 3 equivalent to J het = 100 ×3 ÷3 = 100 +200 −67 cm −2 s −1 .In the following analysis, J het is quantified as × ÷ , representing a factor error.The uncertainty due to stochastic freezing is derived by running 10 5 IFSs with different values of N tot and calculating J het where the widths of the fiducial limits are smallest, i.e., at f frz 0.5.Thus, J het derived from N tot yields the smallest error estimate possible or the limit of greatest experimental accuracy.Figure 9a illustrates that smaller N tot results in larger J het .When N tot = 30 for example, J het = ×15 ÷5 , and when N tot = 1000, J het = ×1.3÷1.3 .The uncertainty contribution due to T is calculated using the slope of J het vs. T following a similar procedure as in Riechers et al. (2013).Using the ABIFM at various temperature ranges and for different INP types (Knopf and Alpert, 2013), J het varies by a factor of 7.5 ± 5.5 per degree K.This means that if T = ±1.0K, J het = × ÷ 7.5 on average, but can be × ÷ 2 or × ÷ 13 depending on the INP type and the range in T and RH.For example, T = ±0.5 K translates to J het = × ÷ 3.75 as displayed in Fig. 9a.Considering the uncertainty in RH, Eq. ( 6) is used to derive J het = J het ( a w )/J het ( a w ± RH) = 10 m RH .Values of m in Eq. ( 6) are taken from this study and from Knopf and Alpert (2013) ranging from 15 to 123 and results in 69 on average.The mean and range of J het due to RH are shown in Fig. 9b.For example, if RH = ±3 %, then J het = × ÷ 117 on average.If ISA per droplet varies in an experiment, but is assumed to be uniform, J het is overestimated for f frz < 0.5 and underestimated for f frz > 0.5.This effect is quantified by allowing σ g to vary and calculating the ratio actual het evaluated at f frz = 0.1 and 0.9.The resulting J het is displayed in Fig. 9c as a function of σ g .If σ g = 10, for example, then J het = ×4 ÷20 at f frz = 0.1 and 0.9.Finally, J het is directly proportional to A g shown in Fig. 9c, e.g., if A g = × ÷ 5, then J het = × ÷ 5. Figure 9 demonstrates that each experimental parameter contributes to the uncertainty in J het .The total uncertainty in J het can then be estimated by summing the error contributions due to N tot , T , RH, σ g , and A g , respectively.Figure 9 shows dotted lines serving as example values of experimental uncertainties and corresponding J het .Applying N tot = 30, T = ±0.5 K, RH = ±3 %, σ g = 10, and A g = × ÷ 5, results in J het = ×148 ÷154 .If laboratory immersion freezing studies were to be conducted under these conditions, then the range in experimentally derived J het should be over 4 orders of magnitude.Notice that the uncertainty due to RH alone can potentially dominate the total uncertainty.We hope that Fig. 9 provides guidance in conducting future immersion freezing studies.
We test our analysis to reproduce experimentally derived uncertainty.In Knopf and Alpert (2013), all experimentally derived J het fell within ±2 orders of magnitude as a function of the a w criterion (Eq.7) and as a result, this range was adapted as a conservative uncertainty estimate for the ABIFM model.The root mean square error of over 18 000 droplet freezing events for six different INP types was experimentally derived independent from model simulations, as an alternative uncertainty estimate exhibiting values as high as ±1.3 orders of magnitude.Experimental parameters of studies used in the formulation of the ABIFM for pure water and aqueous solution droplets (Alpert et al., 2011a, b;Knopf and Forrester, 2011;Rigg et al., 2013;Knopf and Alpert, 2013) were about N tot = 300, T = ±0.3K, RH = ±1 %, σ g = 5, and A g = × ÷ 5. Applying the analysis displayed in Fig. 9 results in an uncertainty of J het = ×16 ÷18 (spanning about 2.5 orders of magnitude) for the ABIFM model.This estimate is in excellent agreement with independently derived root mean square errors of J het (Knopf and Alpert, 2013) and demonstrates the accuracy of our uncertainty analysis.
Model simulations reproduced observations of immersion freezing due to illite by Diehl et al. (2014) and Broadley et al. (2012).These experimental data were included in a recent intercomparison study of illite immersion freezing by Hiranuma et al. (2015).Using 17 different instruments, experimentally derived n s (T ) values were observed to increase from 10 −3 to 10 8 cm −2 when T decreased from 263 to 236 K, equivalent to a slope of 0.5 orders of magnitude per 1 K.The instruments used are grouped by common methods and include, (i) cold stage (Broadley et al., 2012;Bingemer et al., 2012;Schill and Tolbert, 2013;Wright and Petters, 2013;O'Sullivan et al., 2014;Budke and Koop, 2015), (ii) liquid aliquots (Hill et al., 2014), (iii) droplet levitation (Szakáll et al., 2009;Diehl et al., 2014;Hoffmann et al., 2013), (iv) cloud chamber (Möhler et al., 2003;Niemand et al., 2012;Tajiri et al., 2013) and (v) continuous flow diffusion chamber (Bundke et al., 2008;Stetzer et al., 2008;Welti et al., 2009;Lüönd et al., 2010;Chou et al., 2011;Friedman et al., 2011;Hartmann et al., 2011;Kanji et al., 2013;Tobo et al., 2013;Wex et al., 2014).The scatter in the n s is roughly 3 orders of magnitude, but depending on T , a n s range of 2 and 4 orders of magnitude can envelop the data.However, the authors provided no quantitative uncertainty analysis to explain this scatter.Since experimental methods and data reproduced by presented model simulations are included in Hiranuma et al. (2015) for illite, we apply the quantitative uncertainty analysis presented in Fig. 9 to provide a potential explanation of the data scatter.We note that the abscissa in Fig. 9 extends to a value of J het equal to a factor of 300, to encompass typical uncertainties of about ±2 orders of magnitude.Although, J het and n s (T ) are different quantities, the contribution to their uncertainties is the same for T , RH, σ g , A g .
Experimental T uncertainty for all methods typically ranged from ±0.2 to ±1.0 K, and hence T = ±0.5 is chosen as a representative value.Considering the slope n s vs. T , T = ±0.5 contributes a factor of ∼ 2 uncertainty to n s (T ), or n s = × ÷ 2. The ISA distribution width parameter of simulated experiments (Tables 1 and 2) is averaged to represent n s (T ) data, yielding a reasonable value of σ g = 7, resulting in n s = ×3 ÷12 .The ISA measurement error is considered to be A g = × ÷ 5, thus n s = × ÷ 5. Calculation of n s (T ) is not stochastic by design, and thus any uncertainty contribution due to N tot on n s (T ) was previously not considered (Hiranuma et al., 2015).Additionally, the intercomparison analysis ignores differences in experimental timescales in n s (T ) derivation.However, this study demonstrates that the stochastic uncertainty may be able to explain immersion freezing data and may contribute to the range of data scatter in n s (T ).Typically, N tot is about 50 which serves as a reasonable representation yielding n s = ×8 ÷4 , although N tot can vary between 10 and more than 1000 depending on the experiment.Previous immersion freezing experiments for illite have shown that when r or residence time differ by 1 order of magnitude, freezing temperatures shift by about 0.75 K on average (Broadley et al., 2012;Welti et al., 2012;Knopf and Alpert, 2013).As discussed in Hiranuma et al. (2015), cooling rates and residence times in the different instruments varied over ±2 orders of magnitude, or t= × ÷ 100, corresponding to T = ±1.5 K, and thus contributing to an error of ±0.75 orders of magnitude or n s = × ÷ 6. Accounting for all uncertainties and making use of ÷(2+12+5+4+6) for a total uncertainty of n s = ×24 ÷29 , or an uncertainty range of 2.8 orders of magnitude.The vast majority of data in Hiranuma et al. (2015) fall within this uncertainty and implies that variability in n s (T ) may be attributed to experimental, time-dependent, and stochastic uncertainties.It is important to note that the uncertainty due to neglecting time, ISA variability and stochastic effect con-tributes more to n s , than T and ISA measurement error.Hiranuma et al. (2015) hypothesized that experimental procedures of droplet or particle preparation, including particle generation, size selection, ice crystal detection, particle loss at instrument sampling inlets, contamination, inhomogeneous temperature, and differences in surface cation concentration between wet dispersed or dry dispersed particles may be the cause in measured scatter in n s (T ) data.These effects are not considered in the uncertainty analysis presented here, but may also contribute.

Summary and conclusions
Immersion freezing simulations based on a droplet resolved stochastic ice nucleation process applicable for various types of INPs and experiments are presented here for both isothermal conditions and applying a cooling rate, r.The parameters in the IFSs are all physically defined and measurable, including the heterogeneous ice nucleation rate coefficient, J het , the number of droplets at the start of an experiment, N tot , and the immersed surface area (ISA) per droplet.When knowledge of ISA per droplet is not known, it may be assumed to be lognormally distributed.For IFSs in which a cooling rate, r, is applied, J het as a function of T and aqueous solution water activity, a w , can be calculated following the water activitybased immersion freezing model (ABIFM) applicable for both pure water (a w = 1.0) and aqueous solution (a w < 1.0) droplets.These IFSs generate frozen and unfrozen droplet fraction data, f ufz and f frz , respectively, and using a Monte Carlo method in which 10 5 IFSs are performed under the same conditions, 5th and 95th percentile bounds are derived as uncertainty estimates.
The sensitivity of f ufz on σ g and N tot was tested using sets of isothermal IFSs, where a single set is referred to as a model simulation.Uniform ISA (i.e., σ g = 1) resulted in f ufz (on a logarithmic scale) being linear with t.When ISA varied lognormally with parameters µ = ln(A g ) and σ = ln(σ g ), where σ g > 1, ln(f ufz ) vs. t exhibit non-linear behavior.When larger or smaller N tot was used, f ufz had a smaller and larger uncertainty, respectively, due to the statistical significance of observing more freezing events.These results demonstrate that in laboratory immersion freezing experiments, variable ISA imposes changes in trajectories of f ufz and f frz over time, and that the number of investigated droplets significantly impacts experimental uncertainty.
Cooling-rate model simulations were used to test the validity of assuming uniform ISA.This was accomplished by recalculating J het after simulation of immersion freezing in two ways, either (i) assuming uniform ISA referred to as the "apparent" ice nucleation rate coefficient, J apparent het , or (ii) accounting for variable ISA referred to as the "actual" ice nucleation rate coefficient, J actual het .When different r were applied in simulations, values of J apparent het were significantly different from each other.When comparing experiments with different ISA but identical r, J apparent het (T ) were again significantly different.For f frz < 0.5 and f frz > 0.5, J apparent het was over and underestimated, respectively, compared to J actual het , yielding an erroneous slope of J apparent het (T ).These results demonstrate that the assumption of identical ISA implicitly imposes a cooling-rate and surface area dependence on experimentally derived J het (T ).However, derivation of J actual het from model simulations accounting for variable ISA were consistent for different r and ISA, supporting a stochastic immersion freezing description as predicted by CNT.
Model simulations in which variable ISA was considered reproduced laboratory experiments using Arizona test dust (ATD) (Wright and Petters, 2013), illite (Broadley et al., 2012;Diehl et al., 2014), kaolinite (Wex et al., 2014;Herbert et al., 2014), feldspar (Herbert et al., 2014), and natural dusts from Asia, Israel, the Sahara desert and Canary Islands (Niemand et al., 2012) acting as INPs.Despite whether isothermal or linear and non-linear cooling rates were applied, modeled and experimental f frz and f ufz were in agreement within the stochastic uncertainty.More importantly, experimentally derived J het (T ) and simulated J apparent het were in agreement for ATD, illite, kaolinite and feldspar possible indicating an imposed bias solely due to the assumption of uniform ISA and not to physical processes governing ice nucleation.On the other hand, variability of ISA in experimental studies using size selected particles from a differential mobility analyzer (Wex et al., 2014), were modeled based on a bipolar charge distribution.There was a 1 K difference between J het assuming uniform ISA and accounting for variable ISA following the bipolar charge distribution, but this was within stochastic uncertainty limits.Thus, the assumption of identical ISA may be valid when size selecting particles.However, the actual ISA distribution for each study should be measured to verify this assumption about ISA in droplets in order to correctly interpret immersion freezing (Hartmann et al., 2016).In general, model simulations can correct for any introduced bias or calculate "actual" values, or J actual het , when not provided.J actual het resulted in consistent agreement between different studies and additionally new a w -based parameterizations of J het ( a w ) for feldspar and natural dusts.
A quantitative uncertainty analysis of J het was presented applicable for experimental studies in which the contribution, due to (i) N tot , (ii) temperature accuracy referred to as T , (iii) a w or RH accuracy referred to as RH, (iv) σ g , and (v) the accuracy of A g referred to as A g , was individually quantified.The following points summarize these error sources and give recommendations for future experimental studies: -Applying too few N tot or performing only a single ice nucleation experiment in laboratory studies results in highly uncertain freezing results.Therefore, repetition of immersion freezing experiments or a statistically significant number of droplets must be applied.We recommend using at least 100 droplets and three independent freezing cycles in order to better quantify data scatter and average J het , f frz , and f ufz values.This contributes to a range of 0.75 orders of magnitude in the uncertainty of experimentally derived J het .
-For different INP types, the slope of J het vs. T is not the same and thus, the uncertainty due to T is INP-type dependent, but can be as high as 1 order of magnitude per 1 K.We recommended that T remain < ±0.5 K to achieve an acceptable uncertainty contribution, i.e., half an order of magnitude.
-The greatest source of error stems from RH, or RH.Immersion freezing experiments for RH < 100 % should aim for RH to be as small as possible.Current and future immersion freezing experiments should be designed to carefully control RH and quantify its uncertainty.
-Droplets in laboratory immersion freezing experiments will not have identical ISA, but will vary from droplet to droplet (σ g ) around some ISA value (A g ).Variability in ISA and corresponding uncertainty should be quantified and accounted for when analyzing ice nucleation experiments.
-Surface area and nucleation timescales clearly affect immersion freezing data.Common assumptions of ISA and neglecting the impact of variable experimental timescales will lead to an incomplete experimental accuracy and uncertainty.Consideration of these effects is recommended to narrow the uncertainty in predicting ice crystal formation.
Considering that INPs have variable ISA may impact atmospheric ice crystal numbers.For a broad surface area distribution of INPs, ice nucleation should occur over a broader range of time and temperature, when compared with a narrow INP surface area distribution.This results in greater ice particle production at warmer temperatures, important for mixed phase cloud formation and their evolution.We suggest that field measurements determine and consider the entire aerosol size distribution as a source of INPs for implementation of a stochastic, time-dependent ice nucleation process characterized by J het , which is easily parameterized following the ABIFM in addition to current methodologies.
Our findings concerning laboratory immersion freezing experiments emphasize the importance of setting constraints on the minimum number of droplets and experimental trials that need to be employed for improved characterization of ISA per droplet.The results presented here resolves commonly used assumptions that contribute to additional uncertainty in predicting immersion freezing data for model implementation.The simulations use ABIFM, shown to be valid for various INP types.We demonstrate that the AB-IFM can reproduce immersion freezing by mineral dust for many vastly different experimental designs and measurement methods.Laboratory derived J het values can aid in testing existing ABIFM parameterizations and formulating new ones.Their application to a very simple stochastic freezing model based on a binomial distribution in accordance with classical nucleation theory, can reconcile immersion freezing data for various INP types and measurement techniques when the applied INP surface areas are treated more realistically.These findings hopefully stimulate further discussion on the analytical procedure and interpretation of immersion freezing and its implementation in atmospheric cloud and climate models.
The Supplement related to this article is available online at doi:10.5194/acp-16-2083-2016-supplement.

Figure 1 .
Figure 1.Sensitivity calculations of the unfrozen droplet fraction, f ufz , as a function of time, t, derived from model simulations for a total number of droplets, N tot , and variability of ice nucleating particle surface area, σ g .(a) Model simulated 5th and 95th percentile bounds of f ufz are shown as dark green (Iso1), light green (Iso2), dark blue (Iso3), and light blue (Iso4) shading.Parameter values are given in the legend.(b) Simulated 5th and 95th percentile bounds of f ufz derived from IsoWR are shown as the orange shading along with experimental data of isothermal immersion freezing by Arizona test dust (Wright and Petters, 2013) shown as black circles.Parameter values for all model simulations in (a) and (b) are given inTable 1.

Figure 2 .
Figure2.Simulated and experimentally(Broadley et al., 2012;Herbert et al., 2014) derived unfrozen droplet fractions, f ufz , as a function time, t.Model simulations and INP types used are (a) IsoBR and illite, (b) IsoHE1 and kaolinite, and (c) IsoHE2 and feldspar, respectively.Orange lines and shading represent f ufz with corresponding 5th and 95th percentile bounds, respectively.Parameter values for model simulations are given in Table 1.

Figure 4 .
Figure 4. Sensitivity calculations of heterogeneous ice nucleation rate coefficients, J het , and frozen droplet fractions, f frz , on cooling rate, r, derived from model simulations Cr1 (orange) and Cr2 (blue) where r = 0.5 and 5.0 K min −1 , respectively.J het as a function of temperature, T , are shown in (a) assuming uniform ice nuclei surface area (ISA) per droplet yielding J apparent het , and (b) accounting for different ISA yielding J actual het .The dashed lines in (a) and (b) are J apparent het and J actual het , respectively.Shadings in (a) and (b) correspond to upper and lower fiducial limits with x = 0.999 confidence and the solid red line is calculated from Eq. (6) for illite(Knopf and Alpert, 2013).Frozen droplet fractions, f frz , are shown in (c) and (d) where dashed lines and shadings represent f frz and 5 and 95 % bounds, respectively.Parameter values for Cr1 and Cr2 are given in Table2.

Figure 5 .
Figure 5. Frozen droplet fractions, f frz , and heterogeneous ice nucleation rate coefficients, J het , from immersion freezing coolingrate, r, dependent model simulations CrHE1 and CrHE2, where r = 0.2 (orange) and 2.0 K min −1 (blue), and experimental data of feldspar acting as immersion INPs (Herbert et al., 2014).Dashed lines and shadings in (a) are f frz and 5th and 95th percentile bounds, respectively.J het as a function of temperature, T , are shown in (b) assuming uniform ice nucleating particle surface area (ISA) per droplet yielding J apparent het and (c) accounting for variable ISA yielding J actual het .The dashed lines in (b) and (c) are J apparent het

Figure 6 .
Figure 6.Frozen droplet fractions, f frz , and heterogeneous ice nucleation rate coefficients, J het , from immersion freezing model simulations CrDI1 (orange) and CrDI2 (blue), and experimental data of illite acting as immersion INPs are shown (Diehl et al., 2014).Dashed lines and shadings in (a) are f frz and 5th and 95th percentile bounds, respectively.J het as a function of temperature, T , are shown in (b) assuming uniform ice nucleating particle surface area (ISA) per droplet yielding J apparent het and (c) accounting for variable ISA yielding J actual het .The dashed lines in (b) and (c) are J apparent het and J actual het , respectively.Shadings in (b) and (c) correspond to upper and lower fiducial limits with x = 0.999 confidence.Experimentally derived f frz and J het are shown as circles in (a) and (b), respectively(Diehl et al., 2014).The red line in (b) and (c) is calculated from Eq. (6) for illite(Knopf and Alpert, 2013).Parameter values for CrDI1 and CrDI2 are given in Table2.Fitted J het values from model simulations IsoDI1 to IsoDI3 are shown in (c) and their corresponding error derived from Fig.9considering N tot = 45, σ g = 3.2 and a temperature error, T = ±0.7 K (see Table1).

Figure 7 .
Figure 7. Frozen droplet fractions, f frz , and heterogeneous ice nucleation rate coefficients, J het , from isothermal model simulations IsoCFDC (orange and black) and IsoLACIS (blue and green), and experimental data of immersion freezing due to kaolinite by Wex et al. (2014) are shown.Dashed lines and shadings in (a) are f frz and 5th and 95th percentile bounds, respectively.J het as a function of temperature, T , are shown in (b) assuming uniform ice nucleating particle surface area (ISA) per droplet yielding J apparent het , and accounting for variable ISA yielding J actual het .The dashed lines in (b) are J apparent het and J actual het

Figure 8 .
Figure 8. Frozen droplet fractions, f frz , and heterogeneous ice nucleation rate coefficients, J het , derived from adiabatic cooling immersion freezing model simulations CrNI1 (blue) and CrNI2 (orange).Simulated and experimentally observed ice crystal concentrations are shown in the insert of panel (a).Dashed lines and shadings in (a) are f frz and 5th and 95th percentile bounds, respectively.Experimentally derived f frz and uncertainties by Niemand et al. (2012) are shown as symbols and error bars.J het as a function of temperature, T , is shown in (b) and accounting for variable ISA yielding J actual het , where dashed lines and shading are J actual het

Figure 9 .
Figure9.Uncertainty analysis derived from immersion freezing model simulations.The relative error in the experimentally derived heterogeneous ice nucleation rate coefficient, J het , is referred to as J het .The y axis indicates J het as a factor error, e.g., J het = 10 indicates an error in J het by a factor of 10 in the positive and negative direction.Stochastic error due to the applied number of droplets, N tot , is shown in (a) where red and blue represent the upper and lower fiducial limits of J het , respectively.The error due to temperature accuracy, T , for a variety of INP types is shown in (a) in orange color where the solid line is average J het as a function of T and the shading is for a range of INP types.The error due to the absolute uncertainty in water activity or equivalently relative humidity, RH, is shown in (b) where the blue line is average J het , and the shading represents the range of values for a variety of INP types.The uncertainty due to variability in INP surface area, σ g , is shown in (c) as black and green lines evaluated at f frz = 0.1 and 0.9, respectively.The uncertainty in measuring absolute surface area, A g , is shown in (c) as the purple line.Further details and example uncertainty values given as dotted lines are described in the text.

Table 1 .
Summary of parameters used in isothermal model simulations.

unifying ice nucleation model formulated
Knopf and Alpert (2013)byKnopf and Alpert (2013).Immersion freezing data of Herbert et al. (2014) for feldspar is reproduced by the model simulation IsoHE2.The parameters for IsoHE2 are given in Table www.atmos-chem-phys.net/16/2083/2016/Atmos.Chem.Phys., 16, 2083-2107, 2016 2092 P. A. Alpert and D. A. Knopf: A Diehl et al. (2014))mentally(Diehl et al., 2014)derived unfrozen droplet fractions, f ufz , as a function time, t, using illite.Model simulated 5th and 95th percentile bounds of f ufz are shown as green, orange and blue shading for IsoDI1, IsoDI2 and IsoDI3, respectively.Temperature and average surface area per droplet reported byDiehl et al. (2014)are given in the legend.Parameter values for model simulations are given in Table1.

Table 2 .
Summary of parameters used in cooling-rate model simulations.
(Niemand et al., 2012)ear cooling rate with time due to adiabatic expansion is fitted to experimental trajectories(Niemand et al., 2012)using a fourth-order polynomial.c Natural dusts from Niemand et al. (2012): Asian, Saharan, Israeli and Canary Island dust.tent.It is important to note that J apparent het using new parameters derived for feldspar.Parameter values for CrHE1 and CrHE2 are given in Table2.The fitted J het value from model simulation IsoHE2 is shown in (c) and its corresponding error derived from Fig.9considering N tot = 40, σ g = 8.5 (see Table1) and a temperature error, T = ±0.4K.