Immersion Freezing in Particle‐Based Aerosol‐Cloud Microphysics: A Probabilistic Perspective on Singular and Time‐Dependent Models

Abstract

Cloud droplets containing immersed ice‐nucleating particles (INPs) may freeze at temperatures above the homogeneous freezing threshold temperature in a process referred to as immersion freezing. In modeling studies, immersion freezing is often described using either so‐called “singular” or “time‐dependent” parameterizations. Here, we compare both approaches and discuss them in the context of probabilistic particle‐based (super‐droplet) cloud microphysics modeling. First, using a box model, we contrast how both parameterizations respond to idealized ambient cooling rate profiles and quantify the impact of the polydispersity of the immersed surface spectrum on the frozen fraction evolution. Presented simulations highlight that the singular approach, constituting a time‐integrated form of a more general time‐dependent approach, is only accurate under a limited range of ambient cooling rates. The time‐dependent approach is free from this limitation. Second, using a prescribed‐flow two‐dimensional cloud model, we illustrate the macroscopic differences in the evolution in time of ice particle concentrations in simulations with flow regimes relevant to ambient cloud conditions. The flow‐coupled aerosol‐budget‐resolving simulations highlight the benefits and challenges of modeling cloud condensation nuclei activation and immersion freezing on insoluble ice nuclei with super‐particle methods. The challenges stem, on the one hand, from heterogeneous ice nucleation being contingent on the presence of relatively sparse immersed INPs, and on the other hand, from the need to represent a vast population of particles with relatively few so‐called super particles (each representing a multiplicity of real particles). We discuss the critical role of the sampling strategy for particle attributes, including the INP size, the freezing temperature (for singular scheme) and the multiplicity.

Introduction

Mixed-phase clouds consist of supercooled droplets, ice particles and the moist air surrounding them. Nucleation of both water drops and ice crystals predominantly relies on the presence of aerosol particles serving as nuclei facilitating the onset of phase transition. Ice nucleating particles (INPs) include, for example, dust and sea spray (see Kanji et al., 2017; Knopf & Alpert, 2023, for reviews). Cloud heights overlap with the topmost layer of the biosphere, and inorganic, biogenic, and organic matter, while less abundant, can also serve as INPs (see, e.g., Després et al., 2012; Fröhlich-Nowoisky et al., 2016; Knopf et al., 2018, and references therein).

Clouds affect the Earth's energy budget and climate (e.g., Pincus & Chepfer, 2020). Particles of both natural and anthropogenic origin contribute to cloud condensation nuclei (CCN) and INP populations. Clouds thus play a key role in the indirect effects of anthropogenic particle emissions on the radiative properties of the atmosphere (see Bellouin et al., 2020, for a review). There are numerous open questions linked with the mixed-phase aspects of cloud evolution. Examples range from the impact of aerosol on the persistence of shallow Arctic clouds (Fridlind & Ackerman, 2018; Morrison et al., 2012) to the role of aerosol in determining the dynamics of deep convection (Marinescu et al., 2021).

Tackling research questions (and climate-change policy questions alike) related to mixed-phase clouds continues to pose significant challenges originating from knowledge gaps in both observational and modeling domains (see, e.g., Burrows et al., 2022). Among the modeling challenges that are particularly relevant to aerosol-cloud interactions in mixed-phase clouds are: (a) the representation of interconnected aerosol, droplet and ice particle concentration budgets and the CCN-activation and ice-nucleation processes, for example, immersion freezing, dependent on CCN and INP concentrations (e.g., Stevens et al., 2018); (b) the diversity of both size and composition of atmospheric aerosol serving as CCN and INPs (e.g., Andreae & Rosenfeld, 2008; Knopf et al., 2021); (c) the two-way nature of aerosol-cloud interactions in which ambient CCN and INPs shape cloud droplet size spectrum and the likelihood of glaciation, while cloud and precipitation processes shape ambient aerosol through scavenging (Wood et al., 2012) as well as resuspension (Solomon et al., 2015) of chemically or physically processed nuclei (Kilchhofer et al., 2021).

Improving our ability to accurately model mixed-phase cloud processes in high-resolution models, which are the focus of this work, can enhance our understanding of these clouds and their feedback within climate systems. This knowledge is crucial for improving global circulation models, which currently have uncertainties related to the representation of mixed-phase clouds (see, e.g., Ceppi et al., 2017, for a review). A modeling technique which, by-design, can consider all of the above-listed challenges is referred to as super-particle or particle-based cloud microphysics, a probabilistic approach gaining significant momentum in high-resolution atmospheric cloud modeling (Morrison et al., 2020). In this work, we focus on representation in particle-based models of the immersion freezing phenomenon upon which the presence of INPs within supercooled droplets allows freezing to occur at temperatures higher than the homogeneous freezing threshold. Our contribution builds upon earlier developments detailed in Alpert and Knopf (2016) and Shima et al. (2020) and compares how the so-called “singular” and the alternative “time-dependent” models of immersion freezing can be used in Monte-Carlo particle-resolved simulations of mixed-phase clouds.

Particle-based cloud microphysics modeling is based on the concept of splitting the simulation of a particle-laden flow into two bi-directionally coupled simulation components: an Eulerian continuous-phase component and a Lagrangian dispersed-phase component (for an overview, see Grabowski et al., 2019; Morrison et al., 2020). The Eulerian component is a computational fluid dynamics solver handling the motion of moist air, as well as mass and heat budget corresponding to phase changes of water. The Eulerian component employs gridded representation of scalar fields such as water vapor density. The Lagrangian component solves for the location in both physical and attribute space of computational particles, each representing a large multiplicity of real aerosol or cloud particles (hence the terms simulation-/computational-/notional-/super-particles and super-droplets). The Lagrangian component employs discrete particle representation with each particle being assigned a number of attributed (e.g., mass of water, phase of water, surface of immersed aerosol). The key advantages of particle-resolved microphysics models, as opposed to continuous-field formulations also termed bulk or bin models, include (see, e.g., Grabowski et al., 2019): (a) maintaining the identity of particles throughout the simulation, (b) favorable scaling with the attribute space dimensionality and (c) by-design, strict preservation of the positivity of derived density fields, while their numerical diffusion can be either zero (trajectory integration) or adapted to match the Eulerian component numerical diffusivity (Curtis et al., 2024). Noteworthily, both particle-based and bin-resolved models of the evolution of aerosol and cloud droplet size spectrum are bound by spectral resolution limitations stemming from the limited number of super-particles or bins, both of which represent large numbers of real particles.

The identity preservation is particularly useful for research on aerosol-cloud interactions. Unlike bin or bulk cloud microphysics models, particle-based techniques do not differentiate, at the level of simulation, if a given super-particle represents an aerosol particle, a cloud droplet or a rain drop—such categorization is done only at the level of simulation output analysis based on values of particle attributes (i.e., the wet size in the case of differentiation between aerosol or droplet). For instance, all liquid particles, regardless of their size, are subject to the same set of basic processes such as sedimentation, transport by the flow, condensational growth or evaporation and collisions. Phenomena such as CCN activation, aerosol washout or resuspension are effectively simulated through the combination of the above basic processes. The aerosol reservoir dynamics are thus inherently resolved, which enables studies on aerosol processing by clouds (understood as a combination of both chemical processing and physical processing caused by particle collisions, for a recent discussion, see e.g., Hoffmann & Feingold, 2023). This capability has been exemplified in particle-resolved studies of in-cloud chemistry followed by aerosol re-suspension on cloud droplet evaporation (Jaruga & Pawlowska, 2018; Yao et al., 2021). Furthermore, what is also relevant in the context of triggering of freezing, particle-based models offer flexibility in terms of handling particle response to ambient temperature and supersaturation fluctuations (Abade & Albuquerque, 2024; Abade et al., 2018; Hoffmann & Feingold, 2019) or even tracking the heat content of each particle and resolving the heat accommodation inertia (Richter et al., 2021). In general however, particle-based microphysics models capturing aerosol-cloud interactions impose demanding constraints on the temporal and spatial resolution of the host flow solver, and typically are used in Large Eddy Simulation (LES) frameworks. From the point of view of resolving aerosol composition diversity, particle-resolved models can track the aerosol mixing state (Riemer et al., 2019, Section 6.3.2) which is also relevant for determining the ice formation potential of atmospheric aerosol (Burrows et al., 2022; Knopf, Wang, et al., 2023; Knopf et al., 2018, 2022; Lata et al., 2021). An objective of this work is to inform the future development of climate model aerosol physics schemes such as MATRIX (Bauer et al., 2008). The MATRIX mixing state scheme has been evaluated against single-particle mass spectrometry measurements of composition (Bauer et al., 2013), but the results of such an evaluation for ice nucleation potential would depend on the selection between parameterization that may differ as investigated in this work.

Several authors have used particle-based concepts in studies on ice-phase and mixed-phase cloud microphysics modeling (e.g., Brdar & Seifert, 2018; Jensen & Pfister, 2004; Paoli et al., 2004; Seifert et al., 2019; Shima et al., 2020; Shirgaonkar & Lele, 2012; Sölch & Kärcher, 2012; Welß et al., 2022). This work builds upon the work of Shima et al. (2020) where an aerosol-budget coupled representation of immersion freezing was proposed. The model formulated in Shima et al. (2020) features a probabilistic formulation of the singular representation of immersion freezing which embodies the assumption that ice nucleation only depends on INP-characteristic freezing temperatures (which are sampled taking into account the distribution of immersed insoluble surface area) and thus is time-independent. The singular (time-independent) approximation of the immersion freezing process has long been challenged by evidence supporting the stochastic (and thus time-dependent) nature of the process (see Knopf & Alpert, 2023; Knopf et al., 2020). Moreover, from an implementation point of view, the singular models are found to be grossly inadequate for diagnostically representing the INP reservoir dynamics owing to the inherent requirement that activation temperature-dependence within each INP class be separately tracked in order to accurately account for losses (Burrows et al., 2022; Fridlind et al., 2012; Kärcher & Marcolli, 2021; Knopf, Silber, et al., 2023). By contrast, the simplest stochastic schemes treat all aerosol surfaces uniformly, thus requiring no additional prognostic variables (e.g., Knopf, Silber, et al., 2023). Such findings motivate this work to inform large-scale models that currently require a balance of physical accuracy and numerical efficiency.

In this work, both the time-dependent as well as the singular freezing models are cast in a probabilistic, nuclei-reservoir-resolving and super-particle-number-conserving form. We highlight how a simple Poissonian model of the rate of heterogeneous freezing in time constitutes a common base for both approaches, and how integrating out time, to derive the singular model, entails limitations in robustness of singular schemes to varying cooling rates and flow regimes. We also comment on the importance of representing the diversity of immersed surface areas (as opposed to monodispersity) of which both approaches are capable of representing.

Models of Immersion Freezing

Origins of Models and Nomenclature

Among the seminal works for the discussion presented herein, there is Bernard Vonnegut's (1948) paper reporting on a series of quasi-isothermal experiments exploring freezing of droplets suspended in oil, carried out at different temperatures. Vonnegut's work depicted (Figure 1 therein) both strong temperature dependence of the overall nucleation rate, as well as persistent time-dependent nature of the process. Quoting Vonnegut (1949), “time required for these drops to freeze could be best explained on the basis of the chance formation of stable nuclei on the foreign surfaces.” In Levine (1950), a seminal statistical theory of heterogeneous freezing was presented, and this work is considered the origin of the “singular” hypothesis. It veils the role of time and puts forward the ansatz of each singularity (mote) causing freezing to be associated with a characteristic freezing temperature. For further references and a detailed recount on the earlier works in both the meteorological as well as other domains, see for example, Vali (1971) and Vali (2014, Appendix A therein). In what follows, we give a brief revisit to the key concepts in both the time-dependent and singular approaches, with emphasis on the links with the shared Poissonian base model.

In an attempt to revisit Levine's theoretical considerations with inclusion of dependence on time, in Bigg (1953a, 1953b), a unifying probabilistic description of the process was presented which admits both the time-dependent and the singular characteristics. Albeit Bigg's theory was developed “without appealing to the action of foreign ice-forming nuclei” (and in general without involving a description of any microscopic mechanism), its reinterpretation and applicability to heterogeneous freezing was highlighted already in Mossop (1955), Langham and Mason (1958) and henceforth had been referred to as Mason-Biggs theory (e.g., Michel, 1967). An alternative, simpler derivation of the results obtained by Bigg was subsequently presented by Carte (1959). The “stochastic” label to describe these developments was first used by Stansbury (1961) and Marshall (1961). A review of these developments in the meteorological context is presented in Pruppacher and Klett (2010, Section 9.2.5).

The early works in the development of the theory and nomenclature of heterogeneous freezing outlined above do not feature the later-defined terminology (e.g., Isaac & Douglas, 1972) which in particular differentiates: (a) condensation-, (b) immersion- and (c) contact modes of freezing as three different pathways of triggering heterogeneous nucleation of supercooled droplets (see, e.g., Laaksonen & Malila, 2022, for a recent overview). Herein, we focus on the immersion freezing mode, and contrast the two different approaches of time-dependent and singular conceptualization of the process, neglecting hybrid models that bridge both descriptions (e.g., Vali & Snider, 2015; Wright & Petters, 2013).

Poissonian Counting and Active Sites

The starting point for description of heterogeneous freezing in the treatments akin to Bigg (1953a) and Carte (1959) is the Poisson counting process implying the following form of cumulative probability of $$k$$ freezing-triggering events occurring in time $$t$$ in a population of INPs (see, e.g., Section 7.2.3 Pruppacher & Klett, 2010, presented in the context of homogeneous freezing): 1 $${P}^{\ast }(k\hspace*{.5em}\text{events}\hspace*{.5em}\text{in}\hspace*{.5em}\text{time}\hspace*{.5em}t)=\frac{{(rt)}^{k}\mathrm{exp}(-rt)}{k!}$$ where $$r$$ is a process rate (in arriving at the Poisson limit, a large number of independent events is considered with the probability of each event being small and likelihood of simultaneous events negligible). By evaluating the complement of the probability of zero events occurring, the cumulative probability of freezing can be defined as: 2 $$P(\text{one}\text{or}\text{more}\text{events}\text{in}\text{time}\hspace*{.5em}t)=1-{P}^{\ast }(k=0,t)$$ or: 3 $$\ln(1-P)=-rt$$ In Bigg (1953b), the following factorization of the right-hand-side $${-}rt$$ term was postulated, incorporating a particle-size-related variable labeled $$X$$ (notation following Vali, 2019, see discussion of Equation 1 therein) and the ambient temperature-related variable $$T$$ (cf. Equation 2 in Bigg (1953b) or Equation 1 in Carte (1959)): 4 $$\ln\left(1-P(X,t)\right)=-X\underset{{n}_{X}(T)}{\underbrace{\underset{0{}}{\overset{t}{\int }}{J}_{X}\left(T\left({t}^{\prime }\right)\right)d{t}^{\prime }}}$$ where $$T(t)$$ denotes the ambient temperature evolution, and $${J}_{X}(T)$$ is the nucleation rate expressed per size $$X$$ (for generalization to both surface and volume dependent rate, see discussion of Equations 1–3 in Leonard & Im, 1999). Commenting on the assumption inherent in the Poisson counting process, a prerequisite for application of the model is that the macroscopic size $$X$$ (total over all particles in a population) is large with respect to the microscopic length scale on which freezing is initiated (on a single particle), and that the freezing events are uncorrelated in space or time, essentially instantaneous, hence non-overlapping. Hereinafter, the model is applied to atmospheric aerosol particles populating macroscopic spatial domains and described collectively by a probability density defined over $$X$$.

In Bigg's 1953 works, $$X=V$$ was conceptualized as the drop volume and $${J}_{X}$$ as the homogeneous nucleation rate. Starting with Mossop (1955) and Langham and Mason (1958), an alternative interpretation was provided leading to the presently commonly used model where $$X=S$$ is the immersed insoluble surface area and $${J}_{\mathrm{S}}$$ is the heterogeneous (immersion) freezing rate (which provides a connection to the classical nucleation theory formalism, see Knopf & Alpert, 2023, and references therein). Nevertheless, Bigg's model proved seminal, and it allows to capture the link between time-dependent and singular views of the process through the relation embodied in Equation 4. Bigg's analysis includes a derivation of a temperature-dependent probability-of-freezing (theoretical curves in Figures 2 and 3 in Bigg (1953a)) employing an exponentially linear form of $${n}_{X}(T)$$: 5 $${n}_{X}(T)={c}_{0}+{c}_{1}\mathrm{exp}\left({c}_{2}\cdot T+{c}_{3}\right),$$ where $${c}_{k}$$ $$(k\in 0,1,2,3)$$ are coefficients constrained by vanishing of the probability of freezing at $$T=0{}^{\circ}$$C.

Bigg's formula given here as Equation 5 is used widely in large-scale models for representation of droplet freezing, embracing the $$X=V$$ interpretation, which implies that the freezing probability is modeled as dependent on volume of water rather than on the surface of immersed particles. Examples include: Wisner et al. (1972, see discussion of Equation 26 where parameter uncertainty is attributed with a 10 K freezing temperature uncertainty), Reisner et al. (1998, Equation A.22), Khain et al. (2000, Section 4.2.2), Seifert and Beheng (2006, discussion of Equation 44), Sullivan et al. (2018, passage between Equations A2 and A3) and Barrett and Hoose (2023). Note that, in a number of instances, the formula is given with a typo (e.g., Seifert & Beheng, 2006) with a “$${-}1$$” term included under the exponent instead of next to it to match Bigg's argument of vanishing probability at zero temperature, that is, “$$\mathrm{exp}(\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\cdot {\Delta }T)-1=0$$.”

Even though, Equation 5 lacks explicit dependence on time, the $${c}_{k}$$ coefficients do in fact depend on the cooling rate (Herbert et al., 2014; Murray et al., 2011; Niedermeier et al., 2015; Vali, 1994). This thus, in principle, renders $${n}_{X}$$ applicable only for the cooling rate for which the coefficients $${c}_{k}$$ are fitted (see Section 2.4). An intuitive parameter interpretation links the median $$\left(P=\frac{1}{2}\right)$$ freezing temperature with the size parameter, which was embraced in the measurement data analysis in Bigg (1953a).

Building on the heterogeneous freezing interpretation $$(X=S)$$, the time-integrated (hence time-independent) density $${n}_{\mathrm{S}}(T)$$ as a function of temperature mathematically corresponds to what is presently referred to as the active sites density, a term used as early as in Fletcher (1969, Section 4 and references to earlier works therein) and more recently popularized in the modeling community with the “IASSD” (ice-active surface site density, Connolly et al., 2009) or “INAS” (ice nucleation active surface site, e.g., Hoose & Möhler, 2012) acronyms. The INAS parameterization was used also in Shima et al. (2020) with the following form: 6 $${n}_{\mathrm{S}}(T)={n}_{0}\cdot \mathrm{exp}\left({a}_{\text{INAS}}\cdot \left(T-{T}_{0{}^{\circ}\mathrm{C}}\right)+{b}_{\text{INAS}}\right)$$ where $${n}_{0}=1\hspace*{.5em}{\mathrm{m}}^{-2}$$, $${a}_{\text{INAS}}=-0.517\hspace*{.5em}{\mathrm{K}}^{-1}$$ and $${b}_{\text{INAS}}=8.934$$ (measurement fit from Niemand et al., 2012). In Section 2.4 herein, an implicit dependence on cooling rate embedded in the above coefficients will be further discussed.

Monte-Carlo Simulations of Immersion Freezing

Numerical simulations of the immersion freezing process can be realized using Monte-Carlo techniques which imply random sampling of the particle attribute spectrum, and consequently a probabilistic look on the process (with each model run yielding a unique realization). In the next two subsections, we outline two previously developed aerosol-microphysics-coupled Monte-Carlo frameworks for representing immersion freezing using time-dependent (Alpert & Knopf, 2016) and singular (Shima et al., 2020) parameterizations. Noteworthily, a Monte-Carlo technique has recently been introduced in Frostenberg et al. (2023) for representing immersion freezing but without an explicit link to aerosol properties.

We describe both aerosol-coupled models using particle-based cloud microphysics nomenclature revolving around super-particles and their attributes. One of the key super-particle attributes is its multiplicity (also termed weighting factor) which describes the number of real particles represented by a given super-particle. Whenever a freezing event occurs, all particles represented by a given super-particle crystallize, and regardless of the employed immersion freezing scheme the INP represented by the newly frozen super-particle are removed from the INP reservoir (and reinstated upon melting). There is no splitting of super-particles, and the number of computational particles is constant to ensure practical realizability of such simulations.

The concept of multiplicity is depicted in Figure 1. Without particle coalescence or fragmentation, the multiplicities for each super-particle are maintained throughout the simulation. A degree of freedom exists in the choice of initial multiplicities for super-particles. Figure 1 depicts a case where only a subset of particles contains insoluble immersed surfaces, and the abundances of super-particles with and without insoluble immersed surfaces are equal. In contrast, in the population of particles represented by super-particles, those without immersed surfaces are more abundant. This is taken into account by assigning larger multiplicities $$\xi $$ to the more abundant super-particles.

[IMAGE OMITTED. SEE PDF]

Time-Dependent Scheme Using ABIFM Parameterization

The time-dependent Monte-Carlo model explored herein uses the stochastic water-activity based immersion freezing model (ABIFM, Knopf & Alpert, 2013) applied following Alpert and Knopf (2016) where it was used in analysis of laboratory experiments. Particle freezing is triggered by comparing a uniform random number in the interval (0,1) (shortened as “u01” in the schematic) with the probability of freezing evaluated at instantaneous ambient conditions in each timestep. The probability is evaluated by assuming constant temperature within a timestep, thus Equation 4 leads to (cf. Equation 1 in Bigg (1953b)): 7 $${P}_{i}=1-\mathrm{exp}\left(-{J}_{\text{S}}(T)\cdot {S}_{i}\cdot {\Delta }t\right),$$ where $$i$$ denotes super-particle index, and $${J}_{\mathrm{S}}$$ is expressed using the ABIFM parameterization (Knopf & Alpert, 2013): 8 $${\log }_{10}\left({J}_{\mathrm{S}}(T)\right)={c}_{\text{ABIFM}}+{m}_{\text{ABIFM}}\cdot {\Delta }{a}_{\text{w}},$$ with $${\Delta }{a}_{\text{w}}$$ denoting the difference between droplet's water activity and its value along the ice melting curve here taken as (assuming pure water, ambient conditions at liquid saturation and a single equilibrium value of temperature for all phases): 9 $${\left.{\Delta }{a}_{\text{w}}\right\vert }_{RH\approx 100\%}=1-\frac{{p}_{\text{s}}^{\text{ice}}(T)}{{p}_{\text{s}}^{\text{liq}}(T)},$$ where RH denotes the relative humidity and $${p}_{\text{s}}^{\text{ice}}{\fracslash}{p}_{\text{s}}^{\text{liq}}$$ is the ratio of saturation vapor pressures over solid and liquid phases of water (Koop, 2002). The pure-water assumption allows for comparison of ABIFM-based simulations with the INAS-based runs for the latter does not feature dependence on water activity. The neglected effects of solutes on the water activity are readily representable using particle-based aerosol-cloud microphysics and using ABIFM; if represented they allow for simulation of freezing of unactivated solution droplets (see Figure 1).

Consequently, for a particle-based cloud microphysics model, each super-particle is required to include the immersed surface area $${S}_{i}$$ of the modeled particle as an attribute (here, $${S}_{i}$$ denotes the value for a single particle, not multiplied by the multiplicity). Since, unlike in the case of coagulation, evaluation of the probability does not involve particle concentration, its value can be applied for super-particles of any multiplicity (see also the formulation of the probabilistic transport model in Curtis et al., 2016). However, the longer the timestep and the larger the multiplicities (stemming from coarser size-spectral resolution), the larger will be the spread among different Monte-Carlo realizations of the process.

Singular Scheme Using INAS Parameterization

The singular formulation follows Shima et al. (2020). As in the case of the above time-dependent formulation, it also uses a probabilistic approach, however, here, it is limited to random attribute sampling at initialization, which reflects the time independence of $${n}_{\mathrm{S}}(T)$$. The randomly sampled attribute is the freezing temperature. It is sampled from a probability density function $$p$$ based on the cumulative $$P(S,T)$$ given in Equation 4 by employing a parameterized INAS density function Equation 6. This yields (Equation 1 in Shima et al., 2020): 10 $${\left.p(T)\right\vert }_{S}=-S\frac{d{n}_{\mathrm{S}}(T)}{dT}\mathrm{exp}\left(-S\cdot {n}_{\mathrm{S}}(T)\right),$$ which can also be found in an approximated form in Bigg (1953a) expressed using present notation through $$p={\left.{\partial }_{T}P\right\vert }_{X}\approx (1-P)\ln(1-P)$$. Effectively, such singular formulation may be regarded as Poissonian-in-space (where the nucleus surface $$S$$ stands for space, see Equation 11) as opposed to the space-time Poissonian model (Equation 4) it is derived from.

In the singular formulation, the freezing of a super-particle is triggered by comparing the ambient temperature with the freezing temperature associated with a given particle at initialization. Phase-change-triggering is thus deterministic. Thus, the immersed surface area $$S$$ does not need to be retained as a particle attribute and its values can be discarded after particle attribute initialization (or retained as an attribute which is not used by the singular freezing scheme).

In a way, this bears analogy to the simplified aerosol activation representation in particle-based microphysics models using the so-called Twomey CCN activation model (Grabowski et al., 2018). In such a model, the super-particles are assigned a critical supersaturation by sampling from a distribution at initialization, and the aerosol characteristics (such as mass or hygroscopicity) do not need to be retained as particle attributes.

The singular approach offers a significant performance gain, for there is no need to re-evaluate the probability of freezing or draw random numbers at every timestep. This holds even if the immersed surface area is retained as particle attribute, whereas in so-called “bin” (non-super-particle) modeling approaches (e.g., Knopf, Silber, et al., 2023), resolving both the immersed surface area and the freezing temperature increases the dimensionality of the problem making the singular scheme significantly more expensive.

INAS-Embedded Cooling-Rate Signature Issue

The performance gain and appealing simplicity of the singular scheme comes with numerous tradeoffs, though. As noted hereinbefore discussing Equation 5, and highlighted by several other authors (Herbert et al., 2014; Murray et al., 2011; Niedermeier et al., 2015; Vali, 1994), the coefficients defining $${n}_{\mathrm{S}}(T)$$ measurement fits embed a signature of the cooling rate characteristic of the laboratory freezing experiments the fits are based upon. Noteworthily, coefficients derived for $${J}_{\text{het}}$$ and $${n}_{\mathrm{S}}$$ parameterizations share several other sources of uncertainties (e.g., Murray et al., 2011; Vali & Stansbury, 1966) stemming from: challenges in obtaining statistically significant sample of freezing events, characterization of ambient thermodynamic state, and characterization of the immersed surface area spectrum (Alpert & Knopf, 2016)—all of which apply to both methods of parameterizing the INP activity.

Using the simple Poissonian model introduced above and substituting a constant cooling/heating rate $$c=dT/dt$$ into Equation 4 yields (cf. Equation 3 in Bigg (1953a)): 11 $$\ln(1-P(S,t))=-\frac{S}{c}\underset{{T}_{0}{}}{\overset{{T}_{0}+ct}{\int }}{J}_{\text{het}}\left({T}^{\prime }\right)d{T}^{\prime }=-S\cdot {n}_{\mathrm{S}}(T).$$ Taking temperature derivatives of both sides yields (Equations 1 and 4 in Hoose & Möhler, 2012): 12 $${-}\frac{d{n}_{\mathrm{S}}(T)}{dT}=-{a}_{\text{INAS}}\cdot {n}_{\mathrm{S}}(T)=-\frac{1}{c}{J}_{\text{het}}(T),$$ which may be used to mathematically link the INAS and ABIFM parameterizations, namely the $${n}_{\mathrm{S}}(T)$$ and the $${J}_{\text{het}}(T)$$ functions. This is done without a link to the physics attributed to the central notions of both models (i.e., surface site density and nucleation rate) and implicitly assuming identical material type of particles in the modeled population. We note also that this is only defined for a cooling process, whereas natural turbulent clouds subject parcels to chaotic cycles of cooling and warming. Under such alternating cycles, this equation has limited practical modeling utility if the INP reservoir budget is to be closed (cf. Knopf, Silber, et al., 2023), but it nonetheless serves as a useful benchmark.

For such a monotonic cooling process, Figure 2 presents the relevant mathematically comparable quantities of ABIFM and INAS immersion freezing parameterizations (for the same material type). The terms $${-}d{n}_{\mathrm{S}}(T)/dT$$ and $${-}{J}_{\text{het}}(T)/c$$ are plotted as ordinate in logarithmic scale with temperature on the abscissa. Given the form of Equation 12 in which the cooling rate is featured only on one side of the equation, expectedly only one value of the cooling rate allows for a match between INAS and ABIFM. An approximate match, within the 245–240 K temperature range, is depicted for $$c=-0.75$$ K min⁻¹. The above analysis highlights that a discrepancy between singular and time-dependent simulations is expected. It stems from the INAS coefficients, which only apply to a single cooling rate matching the one characteristic of the lab experiment data to which the coefficients were fitted. For this reason, INAS parameterization is considered applicable only for specific cooling-rate range corresponding to updrafts of $${\sim} 1.5-3$$ m s⁻¹ (Kanji et al., 2017, p. 1.22). Because mixing in Arctic stratiform clouds is typically driven by cloud-top cooling, downdrafts exceed updrafts that are typically weaker than 1 m s⁻¹ as inferred from observation-constrained mean Doppler velocities in Avramov et al. (2011), Fridlind et al. (2012, note opposite Doppler sign convention) and Silber (2023, Figure 5, top panel in column 4). This implies cooling rates outside of INAS applicability range.

[IMAGE OMITTED. SEE PDF]

Early discussion of the relation of singular and time-dependent formulations can be found in Fletcher (1958) and Vali and Stansbury (1966). Recent works covering it include: Ervens and Feingold (2013, Equation 6 therein), Kubota (2019, Equation 25 therein) and Cornwell et al. (2021, Equations 1 and 2 in Supplement therein), albeit providing a less general relationship than here or in Hoose and Möhler (2012). Depiction of the cooling-rate dependence of $${n}_{s}$$ based on experimental data can be found in Herbert et al. (2014, Figures 5 and 6) and in Niedermeier et al. (2015, Figure 3), while the robustness of ABIFM $${J}_{\text{het}}$$ to the cooling rate is discussed in Knopf and Alpert (2013) and Alpert and Knopf (2016).

This line of research led to introduction of modified singular models that account for cooling rate dependence in INAS-like parameterizations (Broadley et al., 2012; Murray et al., 2011; Vali, 1994). While such parameterizations do address the issue highlighted in simulations presented herein, that is, the experimental cooling rate “signature” in INAS fits, the lack of time-dependence in models employing an INAS fit renders any such parameterization still less general compared with a nucleation-rate based approach. From the herein presented particle-based perspective, this can be noted in the assignment of $${T}_{\text{fz}}$$ to each super particle at initialization. Such singular ansatz still entails prescription of a cooling rate, which limits application to flow regimes with a wide spectrum of cooling rates. Moreover, any kind of singular approach will yield no nucleation in the limit of quiescent flow or within downdrafts where temperature likely increases along particle trajectory. Another advantage of the time-dependent freezing event triggering is that the immersed surface area is not fixed for a particle and may be coupled with dynamically varying physical properties of the particles (e.g., coalescence-triggered increase in total immersed surface area or water-content dependence in case of computing the area covered by surfactants) as well as chemical characteristics of the solution (i.e., changes in the water activity and hence $${\Delta }{a}_{w}$$ depending on the presence of solutes).

Handling of Freezing-Related Attributes Upon Particle Collisions

While in this study, particle collisions causing coalescence or fragmentation are not simulated, the choice of super-particle attributes stemming from the formulation of the immersion freezing model has implications on the collisional growth/breakup handling. In the case of above-presented singular probabilistic freezing model, the randomly sampled freezing temperatures are super-particle attributes. In Shima et al. (2020), it was proposed to select the higher of two freezing temperatures when assigning a freezing temperature for a freshly coalesced super-particle. The time-dependent model inherently allows a more general treatment because the immersed surface is an extensive attribute and scales additively as any other extensive attribute upon coalescence. In the case of collisions of particles with similarly sized immersions, likely having similar freezing temperatures in the singular scheme, the freshly coalesced particle has close-to-unchanged likelihood of freezing in the singular scheme, while in the time-dependent scheme the likelihood doubles by summing the surface areas of immersed material. Analogous limitations apply to the Twomey CCN activation scheme (Grabowski et al., 2018) where activation supersaturation is used as an attribute leading to potentially unphysical attribute evolution upon coalescence or breakup.

It is an open question how to handle freezing-attribute value assignment upon super-particle breakup. One of the challenges is linked with the insolubility of the immersed surfaces which conflicts with the equipartition assumption implied, for example, in representation of hygroscopicity-related aerosol attributes upon breakup. Intuitively, only a subset of fragmented super-particles should retain an immersed surface, however representing it would require an additional super-particle to be introduced into the system which is essentially unfeasible in large-scale simulations for which fragmentation algorithms are required to maintain constant super-particle count (de Jong, Mackay, et al., 2023; Kotalczyk et al., 2017; Lee & Matsoukas, 2000).

Summary: Probabilistic Perspective on Singular and Time-Dependent Models

From descriptions of particle-based immersion-freezing models given in Sections 2.3.1 and 2.3.2, it follows that not only the time-dependent scheme but also the singular scheme generate different realizations if the simulation is repeated. This is obvious for the time-dependent scheme because it simulates the time evolution with probabilistic sampling, but it may come as a surprise for the singular scheme since it is also referred to as “deterministic” in literature.

To understand why the singular scheme also generates different realizations, we need to consider that the particle-based modeling approach requires assignment of the initial particle attributes to a finite set of computational particles. That is, in the case of the time-dependent scheme, we need to assign a surface area to each computational particle, while in case of the singular scheme, we need to assign both a surface area and a freezing temperature. How can this be accomplished?

This per-particle information could, in principle, come from single-particle measurements, but these do not exist in practice. We therefore adopt the standard approach of sampling to assign the initial particle attributes, which is explained in detail in Sections 3.2, 4 and Appendix C. This random initialization is the reason why the singular scheme also exhibits probabilistic features even though the freezing process itself is entirely deterministic. We will provide more detail on the sampling process in Section 3. Figure 4 summarizes the two approaches schematically.

The initial sampling of freezing temperatures for the singular scheme corresponds to assigning active sites to the computational particles. Equation 10 tells us how to do this—the INAS density function $${n}_{\mathrm{S}}(T)$$ is key here. An implication is that even if two particles have the exact same surface area, the number of active sites and their associated freezing temperatures may be different, hence they will freeze at different temperatures.

In summary, it is helpful to think of the initialization process and the time evolution process of the simulations as two separate issues. For the singular scheme, this results in a probabilistic assignment of freezing temperatures to each computational particle, which in turn will result in a different outcome for repeated simulations even though the time evolution for this scheme is strictly deterministic. For the time-dependent scheme, particles are initialized with randomly sampled surface areas but are otherwise identical, but then evolve with a stochastic freezing process, necessarily resulting in different outcomes for each repeated simulation.

Ensemble means (and ensemble spread) from multiple realizations of the singular and the time-dependent models are expected to match if: (a) the ambient cooling rate matches the conditions corresponding to the data set from which the INAS coefficients were derived, (b) the distribution of sampled immersed surface areas is consistent, (c) the freezing temperatures are sampled from a distribution matching condition Equation 12 and (d) auxiliary conditions for freezing triggering are consistent (RH$${ >} $$100% in this case).

Box-Model Simulations

Simulation Framework and Initial Condition Attribute Sampling

The box model simulations lack any spatial context, and can thus be referred to as zero-dimensional or pertaining to an unspecified volume of air. The only ambient dynamics considered is an imposed temperature evolution.

The aerosol parameters are the median surface area (here always set to the surface area of a sphere of diameter $$0.74\ {\upmu }\mathrm{m}$$ corresponding to the large mode of the “Desert Dust” spectrum used in Knopf, Silber, et al., 2023) and the geometric standard deviation of a lognormal spectrum of the immersed insoluble surface (unless indicated otherwise, set to an arbitrary value of $$\sqrt[4]{e}$$). Two parameters defining the immersed material type are needed for both singular and time-dependent models, and these are related with the coefficients in the INAS (Equation 6) and ABIFM (Equation 8) formulae, respectively. For INAS, we use the experimental fits from Niemand et al. (2012) which correspond to mineral dust (see discussion of Equation 6 for values). For ABIFM, we use coefficients for mineral dust given in Alpert and Knopf (2016, Table 2): $${m}_{\text{ABIFM}}=22.62$$ and $${c}_{\text{ABIFM}}=-1.35$$ (note that $${J}_{\mathrm{S}}$$ is expressed in cm⁻² s⁻¹ there). Both fits were obtained using the same AIDA chamber experiment. Note however, that procedures for deriving the coefficients for INAS and ABIFM fundamentally differ (Rigg et al., 2013). INAS coefficients are fitted to aggregated data in frozen fraction versus temperature space. ABIFM parameters are calculated from $${J}_{\text{het}}$$ versus water activity space. $${J}_{\text{het}}$$ is directly calculated from experimentally measured temperature and time of freezing, and immersed surface area as detailed in Zobrist et al. (2007) and Alpert et al. (2012). $${J}_{\text{het}}$$ is derived for a given temperature interval from the number of observed nucleation events and associated observation time and total INP surface area (Alpert et al., 2012). Frozen fraction is thus not an input for fitting, but rather can be calculated from the derived $${J}_{\text{het}}$$ values, and used for validation of the ABIFM parameters (see discussion in Alpert et al., 2011; Rigg et al., 2013).

For each temperature profile (discussed in Section 3.2) or each aerosol spectrum shape (Section 3.3), a set of multiple simulations is performed—an ensemble of runs for the singular model and a second ensemble for the time-dependent model, with each member of the ensemble using a different random number generator seed. The sampling procedure for the initial conditions is as follows. Every simulation uses super-particles of equal multiplicity each (for discussion of other attribute sampling approaches, see Appendix C). The value of multiplicity is determined by assigning each super-particle a consecutive quantile of the dry aerosol mass distribution. To match the equal-multiplicity attribute-space sampling with freezing-related attributes, in the case of the time-dependent simulations, the values of the immersed surface attribute are initialized by inverting the cumulative lognormal distribution, and randomly sampling the resultant quantile function. In the case of the singular simulations, the immersed surface areas are sampled in the same way, however, they are not used as particle attributes; instead the immersed surface area values are used to evaluate the cumulative freezing probability as a function of temperature, and this cumulative distribution is inverted into a quantile function $$Q$$ used for random sampling (see discussion of Equation 1 in Shima et al., 2020): 13 $${\left.{Q}_{\text{INAS}}(p)\right\vert }_{S}={T}_{0}+\frac{1}{{a}_{\text{INAS}}}\left[\ln\left(\frac{\ln(1-p)+\mathrm{exp}\left(-S{n}_{0}\,\mathrm{exp}\left(-{a}_{\text{INAS}}{T}_{0}+{b}_{\text{INAS}}\right)\right)}{-S{n}_{0}}\right)-{b}_{\text{INAS}}\right]$$ where $$p\in [0,1]$$. Effectively, the above procedure corresponds to sampling a two-dimensional probability density function as in Figure 3 (depicted for three different values of the geometric standard deviation) and defined through the chain rule by: 14 $${p}_{(T,S)}(T,S)={p}_{T\vert S}(T)\cdot {p}_{S}(S),$$ where $${p}_{S}(S)$$ is a probability density function of the immersed insoluble surface (e.g., a lognormal distribution) and $${p}_{T\vert S}$$ is the conditional probability labeled in Equation 10 as $${\left.p(T)\right\vert }_{S}$$. Noteworthily, despite the singular scheme being associated with determinism (e.g., Vali et al., 2015), the hereby embraced particle-based formulation involves random sampling at initialization and thus leads to different realization in each simulation (see Figure 4).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

From a statistics perspective, Figure 3 presents dependency among two variables which can in turn be described by two marginal distributions (plotted in Figure 3 with teal and black curves) and a so-called copula, the latter containing solely the information on the dependency structure (see, e.g., Schölzel & Friederichs, 2008). From a cloud physics perspective, the figure depicts a two-dimensional probability density of which marginal distributions are commonly measurable quantities, namely: the size spectrum (per surface) and the freezing spectrum (per temperature). In our example, both the INAS parameterization parameters as well as the parameter of the lognormal size spectrum are based on laboratory measurements, while the dependency structure is defined by the $${n}_{\mathrm{X}}$$ function.

Response of Simulated Frozen Fraction to Ambient Cooling Rates

In order to highlight the differences between the singular INAS-based and the time-dependent ABIFM-based schemes, here we present a set of simulations in which the two models are driven by a set of diverse idealized temperature profiles, with all other settings unchanged (approach akin to the analysis depicted in Figure 3 in Kärcher & Marcolli, 2021). The particle attribute sampling, particle dynamics and analysis logic are schematically illustrated in Figure 4. The leftmost section outlines sampling of immersed surface areas for the time-dependent scheme, and of the freezing temperatures for the singular scheme. The middle section depicts how freezing is triggered: using a path-independent probability of transition evaluated in each timestep for the time-dependent scheme (as in a discrete-time Markov chain), or using a deterministic transition for the singular scheme (as in a finite state machine). The right section illustrates that the same procedure is used for both schemes to derive the fraction of frozen particles.

The results are depicted in Figure 5 and are based on simulations performed with 32 super-particles, each simulation covering a time period of 1 hr divided into 120 timesteps. The immersed insoluble surface area is sampled from a lognormal spectrum (see Section 3.1 for parameters). The ensemble size is set to three realizations per model. To emphasize differences in the temporal evolution of the systems depending on the model choice, we present the simulated frozen fraction as a function of time rather than as a function of temperature. Note that immersion freezing is the only process modeled here; while in a natural environment, the frozen fraction evolution would be governed by an interplay of immersion freezing with other processes (particle sedimentation, diffusional growth).

[IMAGE OMITTED. SEE PDF]

Panel (a) in Figure 5 features a linear temperature gradient of $${-}\frac{2}{3}\mathrm{K}\hspace*{.5em}{\min }^{-1}$$ which is comparable to the cooling rate in the AIDA chamber experiment on which both INAS and ABIFM parameterization fits were based (Niemand et al., 2012, Figure 2). The temperature profile is given by the red line with axis markings on the right-hand side. The teal and black line-connected points correspond to the time-dependent and singular simulations, respectively, and are plotted as frozen fraction versus time. The results from the two models match well, with the discrepancies being comparable in magnitude to the differences between realizations for a particular scheme.

Panel (b) in Figure 5 complements the analysis using a lower magnitude of cooling rate with $$c=-\frac{1}{6}\mathrm{K}\hspace*{.5em}{\min }^{-1}$$. Here, the discrepancy between singular and time-dependent simulation is pronounced and clearly larger than the inter-realization spread. The singular scheme triggers fewer freezing events compared to the time-dependent one. Compared to the results in panel (a), as expected, the singular scheme reaches the same frozen fraction at $$240\hspace*{.5em}\mathrm{K}$$ in both panels—regardless of the ambient cooling rate, while the time-dependent scheme is sensitive to the temperature time series.

Panel (c) in Figure 5 depicts how the two schemes respond to a steep temperature drop with a temperature gradient of $${-}\frac{8}{3}\mathrm{K}\hspace*{.5em}{\min }^{-1}$$ followed by a constant-temperature leg. First, complementing the results plotted in panel (b), during the initial steep temperature drop, the singular scheme triggers more freezing than the time-dependent one. Second, during the constant-temperature leg, the two models differ qualitatively as the singular scheme yields no freezing at all and the time-dependent scheme continues to trigger freezing up until the frozen fraction reaches unity.

The scenario depicted in panel (d) in Figure 5 is constructed with a twofold aim. First, it confirms that during constant-temperature legs at the temperature of $$270\hspace*{.5em}\mathrm{K}$$ where freezing is unlikely, neither of the schemes triggers any increase in the frozen fraction. Second, it can be observed that the impulse-like drop in temperature around $$t=15\hspace*{.5em}\min $$ leads to almost four-fold higher frozen fraction in the singular simulation than in the time-dependent one. The behavior depicted in the latter corresponds to a kinetic limitation, connected with the nucleation rate $${J}_{\text{het}}$$ which limits the amount of ice produced, even though a low temperature is reached for a short instance (with the length of this instance multiplying $${J}_{\text{het}}$$, see Equations 4 and 7).

Panel (e) in Figure 5 presents results from a simulation in which a sharp initial temperature drop as in panel (c) is followed by a steady slow increase in temperature. For the INAS-based scheme, this scenario does not differ from the one in panel (c) as freezing can only be triggered in singular schemes while the temperature gradient is negative. For the ABIFM-based time-dependent scheme, there is no such limitation, and while the droplets remain supercooled and the temperature is low enough for the freezing probability to be non-negligible, a freezing-while-warming behavior is observed (see also top panel in Figure 6). Note this is contrary to typical paradigms. For example, CCN activation in a downdraft is not expected to occur. This reflects one of the fundamental differences which makes prediction of ice formation difficult.

[IMAGE OMITTED. SEE PDF]

Panel (f) in Figure 5, complements the analysis with a scenario composed of repeated freezing cycles (in a way in the spirit of refreezing studies of INP as in, e.g., Fornea et al., 2009; Kaufmann et al., 2017; Wright et al., 2013). As in panel (d), each sharp drop in temperature leads the singular scheme to triggering much more freezing events than the time-dependent one. As in panel (e), freezing-while-warming can be observed only for the time-dependent scheme. After the temperature reaches back to $$0{}^{\circ}\mathrm{C}$$, instantaneous melting is happening (for simplicity, as in Shima et al. (2020), the phase change back into liquid state is represented without any inertia). Subsequent freezing cycles follow exactly the same pattern for the singular scheme, while the time-dependent simulation features different realizations in each cycle. This depicts where the stochastic nature of immersion freezing proceeds: for the singular scheme, only the initial sampling of insoluble ice nucleus sizes is probabilistic, while the time-dependent scheme additionally performs Monte-Carlo triggering in each time step.

The simulation presented in panel (f) of Figure 5 is presented in more detail in Figure 6. The top graph is plotted in the same manner as panel (f), while the middle and bottom graphs depict the evolution of the model state vector (particle attributes). All three plots share the abscissa scale. The ordinate scales correspond to freezing attribute values: immersed surface for the time-dependent scheme (middle graph), and the freezing temperature for the singular one (bottom graph). Whenever a line turns from thin to thick, freezing occurs; melting is indicated by thinning of a given line. Data plotted in teal correspond to the state vector of the time-dependent simulation for which different super-particles freeze in each cycle (due to the random sampling procedure). In contrast, the state vector of the singular simulation depicted with black lines features a repeating pattern in which the same super-particles freeze in the same order in each cycle. The probability density functions used for attribute sampling are shown with dashed lines in the right sub-panels of both the middle and bottom graphs, and correspond to the marginal distributions plotted in Figure 3. While the state vector evolution is shown for single realization per model only, the initial sampling for subsequent realizations is presented with unconnected dots.

Response of Simulated Frozen Fraction to Immersed Surface Spectrum

After exploring the impact of cooling rate on the frozen fraction evolution, we now extend the discussion to cover also the impact of the surface area of the immersed particles (Figure 7) and its spectrum width (Figure 8). In Figure 7, two sets of simulation are presented for three different cooling rates: $$c=-0.15$$ K min⁻¹, $$c=-0.75$$ K min⁻¹ and $$c=-3.75$$ K min⁻¹. The lognormal immersed surface area distribution has a width defined by $${\sigma }_{\mathrm{g}}=\mathrm{exp}(1.5)$$. The simulations are carried out with 64 super-particles, five ensemble members are plotted. The panels depict the fraction of frozen particles as a function of the temperature. Monte-Carlo simulation results are plotted with connected dots: teal for the time-dependent simulations and black for simulations using the singular scheme. The spread across different realizations originates from the probabilistic treatment. In addition, analytic INAS-derived cumulative distribution functions (i.e., inverse of Equation 13) corresponding to monodisperse particle populations are plotted with dashed lines for three values of the median surface diameter: $$0.074\ {\upmu }{\mathrm{m}}^{2}$$ (yellow, $$0.1\times $$ label), $$0.74\ {\upmu }{\mathrm{m}}^{2}$$ (burgundy) and $$7.4\ {\upmu }{\mathrm{m}}^{2}$$ (orange, $$10\times $$ label).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

When plotted as a function of temperature, the frozen-fraction profiles are insensitive to the cooling rate for the singular scheme (black connected dots) and match the INAS-derived cumulative probability distribution for the median surface (burgundy dashed line). The time-dependent simulations (teal connected dots) match the theoretical INAS cumulative curve only for the case of $$c=-0.75$$ K min⁻¹ (roughly corresponding to the AIDA chamber conditions used for the derivation of this parameterization). This is in line with the remarks above pointing out the role of the cooling rate for the coefficients defining time-integrated quantities based on the Poissonian model of freezing events. It is also in line with the behavior depicted in panels (a)–(c) in Figure 5 herein and in previous results (Alpert & Knopf, 2016, Figure 5 therein). The singular scheme yields lower ice concentrations than the time-dependent scheme for the slower cooling case, and higher ice concentrations for faster cooling. The behavior depicted in Figure 7 can be reconciled with the expected behavior of the two models in the limiting cases. For an instantaneous temperature drop, the time-dependent scheme will yield negligible ice concentrations (infinitely lower ice concentrations than the singular scheme). In the opposite limit, for a constant temperature in the supercooled regime though above the freezing temperatures of all droplets, the singular scheme would not yield any ice and hence would have infinitely lower concentrations than the time-dependent simulation.

Comparing the results obtained for the three different cooling rates plotted in Figure 7, it is evident that for the presented case, a discrepancy associated with a five-fold increase/decrease in the cooling rate is roughly comparable to that originating from a tenfold change in the median immersed surface area. The 25-fold decrease in the magnitude of cooling rate corresponding to the change from $$c=-3.75$$ K min⁻¹ (panel c) to $$c=-0.15$$ K min⁻¹ (panel a) results in over 5 K shift of the frozen fraction curve, that is, over 2 K shift for a tenfold decrease (cf. Bigg, 1953a, Section 2 iii).

Figure 8 presents two sets of simulations (five singular and five time-dependent) for the following values of the immersed surface spectrum width expressed through the geometric standard deviation $${\sigma }_{\mathrm{g}}$$: $$\mathrm{exp}(1.25)$$ and $$\mathrm{exp}(2.5)$$ (i.e., increasingly larger than in simulations summarized in Figure 7 where $${\sigma }_{\mathrm{g}}=\mathrm{exp}(0.25)$$ was used). In all simulations, the cooling rate is set to $$c=-0.75$$ K min⁻¹, which is the value corresponding to the best fit across INAS and ABIFM parameterizations depicted in Section 2.4. Hence, regardless of the geometric standard deviation value, the singular and the time-dependent simulations agree with each other. However, compared to the monodisperse analytic frozen-fraction profiles, increasing $${\sigma }_{\mathrm{g}}$$ causes a flattening of the frozen fraction profiles.

The significance of the polydispersity of immersion freezing nuclei has been previously highlighted, for example, in Alpert and Knopf (2016) where analogous Monte-Carlo simulations were employed to quantify the bias in laboratory data analyses stemming from assuming monodispersity where actual spectra were polydisperse. Noteworthily, the simulations presented in Shima et al. (2020) and Abade and Albuquerque (2024), using the singular particle-based scheme, were only performed with monodisperse spectra. Our analysis confirms that both the INAS- and ABIFM-based formulations explored herein are capable of capturing the impacts of polydispersity of the immersed surface spectrum.

2D Prescribed-Flow Super-Particle Simulations

To extend the discussion beyond simple box-model considerations, we present an analysis of a set of idealized flow-coupled particle-based simulations. The simulations are performed in two spatial dimensions (2D) and are driven by a prescribed flow field. The goal is to qualitatively explore the macroscopic impact of the choice of either the singular or the time-dependent immersion freezing model in a system characterized by a range of cooling rates (which vary in time from the perspective of a droplet); to depict challenges stemming from simultaneous simulation of aerosol-constrained CCN activation processes; and to relate the preceding discussion with model resolution parameters characteristic of flow-coupled applications (such as LES). We use a minimal framework consisting of (a) an Eulerian fluid-flow component solving for conservation of water vapor and heat in the domain and (b) a Lagrangian super-particle component solving for transport, diffusional growth/evaporation and immersion freezing of the particles. Despite the prescribed-flow simplification, the two components are bidirectionally coupled: (a) particle transport is driven by the fluid flow, diffusional growth, evaporation and immersion freezing is driven by ambient thermodynamic conditions; (b) vapor and heat budget solved by the Eulerian component feature sink/source terms representing uptake/release of heat and moisture by the particles. Collisions, ice sedimentation and ice diffusional growth are not represented in the simulations. Unlike in preceding box-model simulations and unlike in analytically solvable setups, each super-particle is exposed to different cooling-rate history which stems from random sampling of particle locations at initialization (the flow is laminar, particle trajectories are not perturbed).

The concept of using such kinematic simulation framework for cloud microphysics schemes development can be traced back to the work of Kessler (1969, Section 3C) and was subsequently embraced in multiple studies (Arabas et al., 2015; Gedzelman & Arnold, 1993, 1994; Grabowski, 1998, 1999; Jaruga & Pawlowska, 2018; Lebo & Morrison, 2013; Morrison & Grabowski, 2007; Muhlbauer et al., 2013; Rasinski et al., 2011; Schmeller & Geresdi, 2019; Slawinska et al., 2009; Sulia et al., 2013; Szabó-Takács, 2011; Szumowski et al., 1998; Yang et al., 2015). The setup employed here, which mimics a stratiform cloud deck and features periodic horizontal boundary condition and vanishing flow at vertical boundaries, was introduced in Morrison and Grabowski (2007) and later adopted for particle-based simulations in Arabas et al. (2015). Herein, we modify the setup parameters to roughly resemble thermodynamic and aerosol conditions of an Arctic mixed-phase cloud. The dry-air potential temperature and water-vapor mixing ratio fields are initialized with constant values throughout the domain: $$\theta =(289-33.3)\hspace*{.5em}\mathrm{K}=255.7\hspace*{.5em}\mathrm{K}$$ and $${q}_{\mathrm{v}}=(7.5-6.66)\hspace*{.5em}\mathrm{g}\hspace*{.5em}{\mathrm{kg}}^{-1}=0.84\hspace*{.5em}\mathrm{g}\hspace*{.5em}{\mathrm{kg}}^{-1}$$ (where the minuends are the values from the original setup of Morrison and Grabowski (2007) used for warm-rain simulations, and the subtrahends are arbitrarily chosen for the relative humidity profile to roughly match). The dry-air density profile is initialized by integrating the hydrostatic equilibrium equation for dry air. These conditions result in water supersaturation in the upper part of the domain (demonstrated in Figure 9) and supercooled conditions throughout the domain (hence melting never occurs in the simulation). The setup parameters are further outlined below, for details on the simulation framework, see Appendices B and C for Eulerian and Lagrangian components, respectively.

[IMAGE OMITTED. SEE PDF]

We note that in the case of dominant control of dynamics by the liquid-phase properties and radiative processes in Arctic stratiform clouds (e.g., Silber et al., 2020), the main effect of ice is a relatively slow glaciation of coupled or decoupled layers as a whole via net sedimentation of vapor-grown crystals (e.g., Fridlind et al., 2012; Morrison et al., 2011). Thus, the approach introduced in Morrison and Grabowski (2007) serves here to represent such typical quasi-steady states in Arctic stratiform clouds with weak vertical motions $${< } $$1 m s⁻¹. Under such conditions, the strongest parcel-wise cooling rates will follow moist adiabats at cloud base.

The domain extents are set to $$X=1500$$ m and $$Z=500$$ m discretized on a $$60\times 20$$ grid $$({\Delta }x={\Delta }z=25\hspace*{.5em}\mathrm{m})$$. The domain is populated with $$32\times 60\times $$ 20 super-particles with random spatial locations (i.e., 32 per grid cell on average, 38,400 of super-particles in total in the domain).

The aerosol composition, sizes and concentrations are specified as in the mixed-phase particle-based simulations reported in Shima et al. (2020, Section 6.1.3). The super-particles are split into two subpopulations: soluble particles with insoluble immersed surfaces (1,000 per liter) and immersed-surface-free soluble particles (315 cm⁻³). Both subpopulations have the same soluble mass spectrum defined by lognormal dry-radius distribution with two modes having concentrations of $$270\hspace*{.5em}{\mathrm{cm}}^{-3}$$ and $$45\hspace*{.5em}{\mathrm{cm}}^{-3}$$, geometric mean radii of 0.03 $${\upmu }$$m and 0.14 $${\upmu }$$m and geometric standard deviations of 1.28 and 1.75, respectively. The hygroscopicity parameter for the soluble masses is set to $$\kappa =0.61$$ corresponding to ammonium bisulfate (Petters & Kreidenweis, 2007, Table 1). While the 1,000 per liter concentration of immersed-surface-rich particles matches the setup in Shima et al. (2020), instead of a monodisperse population, here a lognormal spectrum of immersed surfaces is used with the same geometric mean surface and geometric standard deviation as used for the box-model simulations presented above (i.e., corresponding to the surface of a sphere with diameter of 0.74 $${\upmu }$$m and with geometric standard deviation of 2.55 following the large mode of the “Desert Dust” spectrum used in Knopf, Silber, et al., 2023).

All given concentrations are interpreted as corresponding to the standard atmosphere conditions at zero height (i.e., $$T=15\hspace*{.5em}{}^{\circ}$$C, $$p=1013.25$$ hPa and $$\mathrm{R}\mathrm{H}=0\%$$) labeled here as STP. The actual volume concentrations have a vertical gradient due to the stratification of the dry-air density.

The INAS fit parameters used to compute freezing temperatures for the singular scheme are taken from Niemand et al. (2012) and thus correspond to mineral dust. Consistently, for the runs with the time-dependent scheme, the ABIFM parameters for mineral dust are used (see Section 3).

In the present study, both subpopulations are set to have equal super-particle counts. However, those subpopulations effectively represent contrasting concentrations of immersed-surface-rich (out of which only a subset may take part in freezing) and immersed-surface-free particles of 1 and 315 particles per cubic centimeter, respectively. This results in multiplicities ranging from $$1.2\times 1{0}^{6}$$ to $$1.7\times 1{0}^{10}$$. For comparison, in the Arctic-stratiform-cloud modeling study of Fridlind et al. (2012), the concentration discrepancy between INPs and non-INP aerosol was set at 1.7 per liter (all taking part in freezing) versus 352 cm⁻³. The employed super-particle sampling strategy is thus motivated by the challenge of resolving rare particles (and rare freezing events) with relatively low number of super-particles. The concept is depicted in the schematic in Figure 1 and its implementation is outlined in Appendix C. Note that without employing the subpopulation split, taking an indicative LES grid cell volume of (100 m)³, 100 super-particles per grid cell, and employing common constant-multiplicity sampling of the aerosol distribution with concentration of 1,000 cm⁻³ results in each super-particle having multiplicity of 10¹³. Such setup would result in most grid cells having no INP representation at all.

The simulations are performed with either singular or time-dependent representations of immersion freezing described in Section 2.3. For both schemes, freezing is set to be contingent on ambient vapor supersaturation, and applies to all particles, regardless of their wet size (i.e., both aerosol-sized and droplet-sized particles may freeze). While the ABIFM scheme is capable of resolving the influence of dissolved soluble material on immersion freezing activity, and the particle-based representation of particle attributes can be set to track the relevant aerosol-composition and water-content parameters, for simplicity and for allowing comparison with INAS results, it is not taken into account.

With the main goal of the simulations being to explore and compare characteristics of the immersion freezing schemes, no other ice-phase processes are enabled. Upon freezing, the super-particles start to act as tracers, vapor transfer to or from the ice surface is not represented, hence neither is the Wegener–Bergeron–Findeisen (WBF) process. Only condensation/evaporation-related latent heat exchange with the environment is simulated.

The wet radii of the particles are initialized by solving for equilibrium size at ambient humidity or at RH = 100%—whichever is smaller (i.e., solving for wet radius matching zero growth rate). Due to the initial supersaturation in the domain, the simulations are carried out with an initial spinup period of 10 min during which freezing is disabled. During the spinup period, the initially activated droplets formed under unrealistically high supersaturation are first deactivated in the downdraft region of the domain, and then reactivated in realistic supersaturation within the updraft region. As a result of the spinup, the cloud deck becomes abundant with interstitial aerosol particles while the size spectrum of the activated particles ceases to be related to the initial unphysically high supersaturation. The spinup procedure is analogous to that described in Arabas et al. (2015, Section 2.2) where sedimentation and coalescence were switched off during the spin-up.

Figure 10 depicts snapshots of the simulation state at 600 s, 1,800 s (i.e., 1,200 s after the spinup period) and 6,000 s. The figure is meant to depict the simulation resolution–both the grid of the Eulerian component and the size of the population of super-particles of the Lagrangian component. Each super-particle in the domain is rendered as a macroscopic sphere with its color corresponding to: (a) the value of the wet radius attribute of the super particle, (b) the immersed surface area (for the time-dependent scheme) and (c) the freezing temperature (for the singular scheme). The range of freezing temperatures stems from the choice of the INP size spectrum (same spectrum and sampling are used for both time-dependent and singular runs). Concentration of liquid and ice particles is depicted with greyscale shading of grid cells. In panel (a), the cloud base is visible at an altitude of around 350 m above which most of the super-particles have sizes indicating activation into cloud droplets; interstitial aerosol is visible in the cloud layer. In panels (b) and (c), the yellow-colored particles correspond to INP-void super droplets (either sampled from INP-void subpopulation, or having negligible immersed surface area or very low freezing temperature); distribution of ice concentration throughout the domain is due to lack of representation in the simulations of ice growth/shrinkage mechanisms; immersion freezing is the only ice process represented, there are supercooled conditions throughout the domain, hence, once frozen, particles remain as frozen tracers.

[IMAGE OMITTED. SEE PDF]

In order to assess the impact of differing ambient cooling rates on the modeled efficiency of immersion freezing, a set of $$2\times 6$$ simulations is analyzed (six settings, each repeated twice for different random seed). Figure 11 summarizes the simulations. The simulations are run for three different settings of the eddy frequency parameter $$A$$ (different line colors in Figure 11) and with either singular (solid lines) or time-dependent (dashed lines) representation of immersion freezing. The different eddy frequencies correspond to maximal air velocities of ca. 2, 1 and 0.5 m s⁻¹ (see legend in panel b in Figure 11). Each line in the plots depicts a single realization of the system evolution, there are two realizations for each setting shown. In principle, these are the ensemble-mean quantities that can be a basis for conclusions on the system behavior (and the inter-realization spread can be a basis for uncertainty estimation).

[IMAGE OMITTED. SEE PDF]

Panel (a) in Figure 11 presents the time evolution of the domain-mean time-accumulated ice concentration. The whole-domain averaging implies counting cells with no cloud or no ice (see gray shading in Figure 10) (local in-cloud droplet concentrations are depicted in panel (c) discussed below). The somewhat unphysical choice of domain-wide statistics does not affect the comparison between the singular and time-dependent schemes. Given the simplistic treatment of ice particles as tracers after freezing (no ice growth or ice sedimentation physics is featured) spatial subsetting of the particles would not result in more realistic concentration values, for this setup.

The 10-min period of spinup when freezing is disabled is indicated with a gray bar at the top of the plot. After the spinup concludes, it is evident that singular simulations (solid lines) consistently predict an order of magnitude lower concentrations of ice than the time-dependent simulations. Such differences are notable because precipitation and desiccation rates can be expected to scale with ice number concentration (e.g., Fridlind et al., 2012). Furthermore, the singular curves level off within 30 min after the end of spinup, while the time-dependent curves exhibit continuous upward trend. This is in line with the model formulation as the singular scheme limits attainable ice concentrations to those corresponding to the number of particles with freezing temperatures above the minimal temperature in the domain, while in the time-dependent scheme there is no such limitation and the concentration of ice monotonically increases with time. It is also evident, in agreement with the box-model results, that the different ambient cooling rates (see discussion of panel b below) stemming from different air velocities trigger differences in singular scheme behavior visible in the first minutes after spinup, while the time-dependent freezing rates are more robust to the flow regime.

Plots in panels (b) and (c) are constructed by analyzing occurrence counts of values of selected parameters in the domain. To ensure that transient spin-up effects at the beginning of the simulation are excluded from the statistics (e.g., concentration of interstitial aerosol stabilizes only after a full eddy revolution), counting is carried out over a period of 40 min from 60 till 100 min of the simulation as indicated with a gray bar at the top of panel (a). The occurrence counts are recorded in bins and plotted using 3-bin moving average for plot clarity.

Panel (b) in Figure 11 depicts frequency of occurrence of cooling rate values in the domain. Since the super-particle model assumes uniform thermodynamic state within each grid cell (see discussion in Arabas et al., 2015, Section 5.1.3), the particles encounter changes in ambient temperature only when crossing grid cell boundaries; however, the derived per-grid-cell cooling rate reported in the plots is a multiplicity-weighted average of cooling rates for all super particles within each grid cell (while the cooling rates for each particle are computed taking into account ambient state properties of both the current and the previous cell for each particle). This approach is consistent with the stochastic super-particle transport representation proposed in Curtis et al. (2016) in which particle positions within a cell are not part of the model state at all. It is evident from the plotted data that even for a simple flow regime governed by the single-eddy stream function (Equation B6), there is a range of cooling rates encountered in the domain. The magnitudes of the cooling rates grow with the eddy spinning rate. The simulated magnitudes of up to 0.025 K min⁻¹ are smaller (by more than an order of magnitude) than those for which agreement among ABIFM and INAS parameterizations is expected given the cooling rates observed in experiments. Only in a few instances (note the logarithmic scale for occurrence counts), the cooling rate magnitudes exceed the $$\vert c\vert > 0.01$$ K min⁻¹ threshold above which the validity of the singular approximation was explored in Wright et al. (2013, see the last row of panels in Figures 1 and 2 therein), and in general the values are well beyond the $$1-2$$ K min⁻¹ INAS applicability range reported in Kanji et al. (2017).

Figure 11c depicts analogous occurrence statistics for values of the concentration of particles larger than 1 $${\upmu }$$m in diameter. Similarly as in the case of the cooling rates, the choice of immersion freezing parameterization has a negligible effect on the presented data, as expected. The vast difference in ice concentrations (few per liter) and droplet concentrations (few hundreds per cubic centimeter) makes the different outcomes of freezing for singular and time-dependent schemes not discernible from the data plotted in panel (c) (note that no ice growth mechanisms are represented here, in particular there is no representation of the WBF process). The gradual shift toward larger droplet concentrations with increasing air velocities is consistent with the expectation of CCN activation in more vigorous updrafts resulting in higher concentration of activated droplets. The spread in concentration values is consistent with the range of vertical velocities encountered at cloud base (analogous to the range of cooling rates depicted in panel b). Quantitatively, in the case of the slow-spinning eddy (maximal velocity of ca. 0.5 m s⁻¹), the occurrence counts peak roughly at 200–250 cm⁻³ (with domain-mean aerosol concentration of 316 cm⁻³), while in the case of the fast-spinning eddy, almost all particles activate into cloud droplets. This confirms that the setup captures the influence of the flow regime on aerosol activation efficiency.

Summary, Limitations, and Future Directions

Here, we contrasted two probabilistic representations of immersion freezing for aerosol-cloud interaction models: a singular INAS-based scheme proposed in Shima et al. (2020) and an alternative ABIFM-based time-dependent scheme akin to the Monte-Carlo approach used in Knopf and Alpert (2013). We depicted that both schemes can be used in particle-based constant-state-vector-size microphysics frameworks. Both approaches involve Monte-Carlo-type random sampling of the attribute space and both admit arbitrary initial attribute distributions (including polydisperse immersed surface spectra). The singular scheme employs deterministic triggering of freezing (finite state machine), while the time-dependent scheme employs Monte-Carlo triggering (discrete time Markov chain). Both schemes keep track of the INP reservoir budget. We explored these two contrasting approaches, using both box-model and flow-coupled simulations, portraying the key trade-off inherent to the singular description in particle-based cloud microphysics modeling, namely: lower computational cost traded off for reduced robustness to diverse cooling-rate regimes.

The key takeaways from the herein discussion, confirming and extending findings from earlier works, are:

• Both the singular and time-dependent schemes can be viewed as constituting numerical representations of the same Poissonian model of the heterogeneous nucleation. Such perspective, employed herein, helps to discern why integrating out the time for a given temperature profile, needed to formulate the singular parameterizations, effectively embeds the cooling rate characteristics within the parameterization coefficients. As a result, as reported in literature, the INAS-based and ABIFM-based simulations match only for one specific cooling rate which is characteristic for the laboratory measurements employed for fitting the INAS coefficients. Only parameterizations featuring time dependence can be robust to ambient flow regimes which might substantially differ from laboratory settings or feature flow patterns precluding nucleation with the singular scheme (downdrafts or quiescent flow). Whereas time-dependent and singular parameterizations derived from the same laboratory data provide a match to that data in a laboratory scenario, they can result in order-of-magnitude and greater differences in ice formation on timescales that are more typical of natural mixed-phase clouds (see also Knopf, Silber, et al., 2023; Knopf et al., 2021).

• Despite the time-dependent approach being computationally costlier than the singular approach when cast in the particle-based framework, it provides the added value of robustness to flow regimes; moreover since the time-dependent approach features immersed insoluble surface area as the freezing-relevant super-particle attribute (as opposed to freezing temperature for the singular approach), it opens up possibilities (not explored herein) of online coupling of the immersion freezing scheme with aerosol physio-chemical dynamics that lead to evolution of the immersed surface. Being an extensive attribute, the immersed surface dynamics are also more readily representable upon particle collisions.

• For both singular and time-dependent schemes, one of the key challenges from the viewpoint of particle-based mixed-phase cloud microphysics modeling, is the sampling of attribute space in which rarely occurring INPs (outnumbered by CCN) need to be represented by admitting much lower multiplicities for immersed-surface-rich super-particles. Furthermore, neglecting polydispersity of the immersed surface area spectrum was shown to have an effect on the frozen fraction predictions comparable to multi-fold changes in cooling rate or median surface values.

Appendix

Appendix A - Convergence Properties of the ABIFM/Monte-Carlo Scheme

The fidelity of the time-dependent Monte-Carlo scheme depends on three resolution-related aspects, namely: (a) temporal resolution (i.e., the timestep length), (b) size-spectral resolution (i.e., the number of super-particles) and (c) ensemble size (i.e., the number of realizations simulated). To depict the above and quantify the convergence, we are analyzing the departures from expected analytic solution of exponential decay for a simple case of constant event rate $$r=1\times 1{0}^{-7}$$s⁻¹ (i.e., corresponding to constant $$T$$ and constant $$X$$ in Equation 4). In such case, since evaluation of the probability integral in Equation 4 is exact regardless of the timestep, the temporal resolution is not expected to play a role (while, in a realistic simulation, the rate changes in time, hence the smaller the timestep the more accurate is the representation). For the particle population sampling, from the law of large numbers, an inverse square root scaling of the error is expected with the number of independent samples (super-particles). Analogous scaling is expected with the size of the simulation realizations ensemble (see, e.g., a recent discussion in the context of particle-based modeling in Nordam et al., 2023, Sections 2.6 and 3.2.3 therein).

To quantify the departure from the expected analytical solution, we use a simple L2 root-mean-square norm defined as $${E}_{L2}=\sqrt{\overline{{\left(A[i]-E[i]\right)}^{2}}}$$ where $$E[i]$$ are the “expected” values (analytic solution: $$1-\mathrm{exp}(-rt)$$), $$A[i]$$ are the “actual” numerical solution results, and the mean is computed over all timesteps indexed by $$i$$.

Figure A1 depicts the data set using a small subset of the simulations the analysis is based on: eight sample simulations carried out with two different random number generator seeds, two different timesteps and two different numbers of super-particles. Time is given on the abscissa and the ordinate represent the frozen fraction. The thick-dashed-red line in the plot depicts the analytic solution (see caption for details). Each freezing event results in the entirety of particles represented by a single super-droplet to become frozen, hence the step-wise evolution of the frozen fraction. The more super droplets (thin lines), the smaller is their multiplicity, and the smaller is the magnitude of the steps. Since constant-multiplicity sampling is employed, each step is of the same “height.”

[IMAGE OMITTED. SEE PDF]

The whole data set consists of $$32\times 6\times 4\times 2=1536$$ simulations performed for: 32 different seeds, six different numbers of super-particles (powers of two, from $${2}^{8}$$ through $${2}^{18}$$), four different timestep settings, and single or double floating point precision. Figure A2 depicts the computed error measures plotted in base-2 logarithmic scale as a function of base-2 logarithm of the number of super-particles. In addition, for the simulations with longest timesteps (5,000 s), overlaid error bars depict values of the standard error computed from 4, 8, 16, and 32 realizations (different random seeds). The plot depicting computations in double precision shows that, as expected, the errors are essentially insensitive to the timestep length (due to assumed constant decay rate). The slope of the dependence of the error on super-particle count matches the theoretical inverse square root lines. The realization spread expectedly diminishes with the number of runs considered.

[IMAGE OMITTED. SEE PDF]

The single-precision calculations in general match, with the intriguing exception of results for shortest timesteps and largest numbers of super particles. This is attributed to the accumulation of arithmetic errors in computation of the frozen fraction (the more super particles, the more numbers being summed up; the smaller the multiplicities and shorter the timesteps, the more contrasting the individual vs. total numbers are). The captured effect of error accumulation is reported here as a cautionary remark in the context of ongoing discussion of introduction of reduced-precision (i.e., below single precision) arithmetics in atmospheric modeling (e.g., Banderier et al., 2024), including in super-particle simulations (Matsushima et al., 2023).

Appendix

Appendix B - Kinematic Framework: Eulerian Component

The conservation of moisture and heat solved by the Eulerian component of the simulation is expressed through two heterogeneous advection equations in a stratified incompressible flow (of a compressible fluid—dry air): B1 $${\partial }_{t}{\theta }_{\mathrm{d}}+\frac{1}{{\rho }_{\mathrm{d}}}\nabla \cdot \left(\vec{u}{\rho }_{\mathrm{d}}{\theta }_{\mathrm{d}}\right)=-\frac{{\theta }_{\mathrm{d}}{l}_{\mathrm{v}}(T)}{{c}_{\text{pd}}T}{\overset{\cdot }{q}}_{\text{v}}$$ B2 $${\partial }_{t}{q}_{\mathrm{v}}+\frac{1}{{\rho }_{\mathrm{d}}}\nabla \cdot \left(\vec{u}{\rho }_{\mathrm{d}}{q}_{\mathrm{v}}\right)={\overset{\cdot }{q}}_{\text{v}}$$ B3 $$\nabla \cdot \left(\vec{u}{\rho }_{\mathrm{d}}\right)=0,$$ where $${c}_{\text{pd}}$$ is the dry-air specific heat at constant pressure, $${\theta }_{\text{d}}=T{\left(1000\mathrm{hPa}/{p}_{\mathrm{d}}\right)}^{{R}_{\mathrm{d}}/{c}_{\text{pd}}}$$ is the dry-air potential temperature and $${q}_{\text{v}}={\rho }_{\mathrm{v}}/{\rho }_{\mathrm{d}}$$ is the water vapor mixing ratio (with $${p}_{\mathrm{d}}$$, $${\rho }_{\mathrm{v}}$$ and $${R}_{\mathrm{d}}$$ denoting partial pressure of dry air, water vapor density and gas constant for dry air, respectively). The density profile with altitude $${\rho }_{\text{d}}(z)$$ is assumed constant in time and constant in the horizontal dimension. The velocity field $$\vec{u}$$ is also constant in time and is defined by: B4 $${u}_{x}=-\frac{1}{{\rho }_{\mathrm{d}}}{\partial }_{z}\psi $$ B5 $${u}_{z}=\frac{1}{{\rho }_{\mathrm{d}}}{\partial }_{x}\psi ,$$ where $$\psi $$ is a stream function constructed to mimic a single eddy spanning the periodic domain (see example in Figure 9) and specified in the following form (for discussion, see also Maxey & Corrsin, 1986, Equation 2.3): B6 $$\psi (x,z)=-A\frac{X}{\pi }\sin \left(\pi \frac{z}{Z}\right)\cos \left(2\pi \frac{x}{X}\right),$$ where $$A$$ is a constant, with the dimension of density times velocity, controlling the eddy frequency (and hence air velocities); $$X$$ and $$Z$$ are domain extents in horizontal and vertical directions, respectively. The choice of $${\theta }_{\mathrm{d}}$$, $${q}_{\mathrm{v}}$$, $${\rho }_{\mathrm{d}}$$ state variable triplet to describe moist air is discussed in Arabas et al. (2015, Appendix A therein). The advective terms in the transport equations B1 and B2 are solved numerically using the MPDATA solver (see, e.g., Smolarkiewicz, 2006, for an overview) using the PyMPDATA package (Bartman, Banaśkiewicz, et al., 2022). The timestep for the Eulerian component (MPDATA) is set to 2.5 s.

Appendix

Appendix C - Kinematic Framework: Lagrangian Component and Attribute Sampling

The source terms in Equations B1 and B2 are evaluated by summing over the changes in super-particle wet volumes within a given timestep as in Equations 25 and 26 in Arabas et al. (2015). The numerics of particle attribute dynamics are solved as follows: particle displacement is integrated with an implicit-in-space scheme (see Section 5.1.2 in Arabas et al., 2015); condensational growth/evaporation (including CCN activation and deactivation) is integrated with an explicit-in-time/implicit-in-supersaturation scheme (see Section 5.1.3 in Arabas et al., 2015). The Lagrangian component uses adaptive substepping for condensational growth/evaporation and for particle transport; freezing is calculated using 2.5 s timestep and double-precision arithmetics (see Appendix A for discussion).

There are multiple possible sampling strategies for aerosol attributes to be applied in context of particle-based cloud microphysics modeling. First, the multiplicity can either vary across the super-particle population or can be set constant (or constant across super-particle subpopulations). Second, there is degree of freedom in the specification of the multidimensional probability density across particle attribute space.

In Unterstrasser et al. (2017), an analysis was presented highlighting how the choice of constant-multiplicity versus variable-multiplicity sampling has a profound effect on the performance of coagulation algorithms. A sampling strategy particularly relevant to the present work, where concentration of particles that drive the population evolution is small relative to the background density, was presented in DeVille et al. (2019). In Dziekan and Pawlowska (2017) and Li et al. (2022), the different multiplicity choices were analyzed in the particular context of sampling for the so-called “lucky droplet” system with contrasting sizes of particles. In Shima et al. (2020), a uniform sampling per subpopulation (immersed-surface-rich vs. immersed-surface-free particles) was used, while in Section 5.3 therein alternative strategies are discussed. Shima et al. (2020) pointed out that the sampling strategy analyses of Unterstrasser et al. (2017) and Dziekan and Pawlowska (2017) were limited to zero-dimensional coagulation-only setups, and generally it is expected that in higher-dimensional attribute space a uniform sampling strategy potentially reinforced with a quasi-random numbers is optimal. Noteworthily, simulations involving processes such as condensation and freezing do increase the dimensionality of the attribute space by incorporating such particle properties as hygroscopicity, soluble mass and insoluble surface, and particle habit information. Clearly the challenge in representation of INP is very relevant as these particles are often rare in terms of concentration compared to other aerosol, while the fidelity of the representation of their collisions also influences the potential for representing the contact-freezing mechanism.

Here, the super-droplets are initialized in pairs, with both super-particles within a pair sharing their location in space and their soluble-substance dry radius. Locations in space are shuffled from a uniform distribution separately in horizontal and vertical dimensions. Soluble-substance dry radii are sampled from the lognormal distribution by inverting its cumulative distribution function (to compute quantiles) and assigning each super-particle pair with an equal multiplicity. Within each super-particle pair, one super-droplet belongs to the immersed-surface-rich subpopulation, while the other to the immersed-surface-free subpopulation. The multiplicities are split among super-particle within each pair according to the immersed-surface-free to immersed-surface-rich concentration ratio. Super-particles belonging to the immersed-surface-rich subpopulation are initialized by sampling from the lognormal distribution of immersed surface areas by inverting its cumulative distribution (consistently with constant multiplicity sampling of the soluble mass distribution). In the case of simulations with the time-dependent scheme, the insoluble immersed surface is used as the particle attribute. In the case of singular simulations, the insoluble immersed surface is used to evaluate the freezing temperature which is used as the particle attribute (i.e., insoluble immersed surface area is not retained as super-particle attribute). Super-particles belonging to the immersed-surface-free subpopulation have their freezing attribute (either freezing temperature or immersed insoluble surface) set to zero (Kelvins or meters squared) precluding triggering of freezing.

Acknowledgments

This study was supported by the Atmospheric System Research Program sponsored by the US Department of Energy (DOE), Office of Science, Office of Biological and Environmental Research (OBER), Climate and Environmental Sciences Division (CESD), US DOE Grants DE-SC0021034 (SA, IS, DAK, and NR) and DE-SC0022130 (JHC and MW). AMF was supported by the NASA Modeling, Analysis, and Prediction Program. SA acknowledges support from the Polish National Science Centre (Grant 2020/39/D/ST10/01220) and AGH Excellence Initiative—Research University (Grant IDUB 9056). Prior to joining AGH University of Kraków, SA conducted the research in part at the University of Illinois at Urbana-Champaign, IL, USA, and at Jagiellonian University in Kraków, Poland. SA is grateful to Shin-ichiro Shima and Thomas Nagler for helpful exchanges on the project. We thank three anonymous reviewers for their feedback on the initial version of the manuscript.

Data Availability Statement

All simulations performed for the study and discussed in the text were carried out using free and open-source Python packages PySDM (Bartman, Bulenok, et al., 2022; de Jong, Singer, et al., 2023) and PyMPDATA (Bartman, Banaśkiewicz, et al., 2022). All presented figures can be reproduced using Jupyter notebooks shipped in the PySDM-examples package (v2.80 and newer) available on () as well as persistently archived at Zenodo (PySDM developers, 2022-2025). The notebooks can be executed in the cloud using platforms such as Google Colab, ARM JupyterHub or anonymously through (links provided within each notebook). In addition, PySDM-examples includes notebooks reproducing using PySDM several simulation plots from Alpert and Knopf (2016).

References

Abade, G., & Albuquerque, D. (2024). Persistent mixed‐phase states in adiabatic cloud parcels under idealised conditions. The Quarterly Journal of the Royal Meteorological Society, 150(763), 3450–3474. https://doi.org/10.1002/qj.4775

Abade, G., Grabowski, W., & Pawlowska, H. (2018). Broadening of cloud droplet spectra through eddy hopping: Turbulent entraining parcel simulations. Journal of the Atmospheric Sciences, 75(10), 3365–3379. https://doi.org/10.1175/JAS‐D‐18‐0078.1

Alpert, P., Aller, J., & Knopf, D. (2011). Ice nucleation from aqueous NaCl droplets with and without marine diatoms. Atmospheric Chemistry and Physics, 11(12), 5539–5555. https://doi.org/10.5194/acp‐11‐5539‐2011

Alpert, P., Aller, J., & Knopf, D. (2012). Initiation of the ice phase by marine biogenic surfaces in supersaturated gas and supercooled aqueous phases. Physical Chemistry Chemical Physics, 13(44), 19882. https://doi.org/10.1039/C1CP21844A

Alpert, P., & Knopf, D. (2016). Analysis of isothermal and cooling‐rate‐dependent immersion freezing by a unifying stochastic ice nucleation model. Atmospheric Chemistry and Physics, 16(4), 2083–2107. https://doi.org/10.5194/acp‐16‐2083‐2016

Andreae, M., & Rosenfeld, D. (2008). Aerosol–cloud–precipitation interactions. Part 1. The nature and sources of cloud‐active aerosols. Earth‐Science Reviews, 89(1–2), 13–41. https://doi.org/10.1016/j.earscirev.2008.03.001

Arabas, S., Jaruga, A., Pawlowska, H., & Grabowski, W. (2015). libcloudph++ 1.0: A single‐moment bulk, double‐moment bulk, and particle‐based warm‐rain microphysics library in C++. Geoscientific Model Development, 8(6), 1677–1707. https://doi.org/10.5194/gmd‐8‐1677‐2015

Avramov, A., Ackerman, A., Fridlind, A., van Diedenhoven, B., Botta, G., Aydin, K., et al. (2011). Towards ice formation closure in arctic mixed‐phase boundary layer clouds during ISDAC. Journal of Geophysical Research, 116, D00T08. https://doi.org/10.1029/2011JD015910

Banderier, H., Zeman, C., Leutwyler, D., Rüdisühli, S., & Schär, C. (2024). Reduced floating‐point precision in regional climate simulations: An ensemble‐based statistical verification. Geoscientific Model Development,17(14), 5573–5586. https://doi.org/10.5194/gmd‐17‐5573‐2024

Barrett, A., & Hoose, C. (2023). Microphysical pathways active within thunderstorms and their sensitivity to CCN concentration and wind shear. Journal of Geophysical Research: Atmospheres, 128(5). https://doi.org/10.1029/2022JD036965

Bartman, P., Banaśkiewicz, J., Drenda, S., Manna, M., Olesik, M., Rozwoda, P., et al. (2022). PyMPDATA v1: Numba‐accelerated implementation of MPDATA with examples in Python, Julia and Matlab. Journal of Open Source Software, 7(77), 3896. https://doi.org/10.21105/joss.03896

Bartman, P., Bulenok, O., Górski, K., Jaruga, A., Łazarski, G., Olesik, M., et al. (2022). PySDM v1: Particle‐based cloud modelling package for warm‐rain microphysics and aqueous chemistry. Journal of Open Source Software, 7(72), 3219. https://doi.org/10.21105/joss.03219

Bauer, S., Ault, A., & Prather, K. (2013). Evaluation of aerosol mixing state classes in the GISS‐modelE‐MATRIX climate model using single particle mass spectrometry measurements. Journal of Geophysical Research: Atmospheres, 118(17), 9834–9844. https://doi.org/10.1002/jgrd.50700

Bauer, S., Wright, D., Koch, D., Lewis, E., McGraw, R., Chang, L.‐S., et al. (2008). MATRIX (Multiconfiguration Aerosol TRacker of mIXing state): An aerosol microphysical module for global atmospheric models. Atmospheric Chemistry and Physics, 8(20), 6003–6035. https://doi.org/10.5194/acp‐8‐6003‐2008

Bellouin, N., Quaas, J., Gryspeerdt, E., Kinne, S., Stier, P., Watson‐Parris, D., et al. (2020). Bounding global aerosol radiative forcing of climate change. Reviews of Geophysics, 58(1). https://doi.org/10.1029/2019RG000660

Bigg, E. (1953a). The supercooling of water. Proceedings of the Physical Society Section B, 66(8), 688–694. https://doi.org/10.1088/0370‐1301/66/8/309

Bigg, E. (1953b). The formation of atmospheric ice crystals by the freezing of droplets. The Quarterly Journal of the Royal Meteorological Society, 79(342), 510–519. https://doi.org/10.1002/qj.49707934207

Brdar, S., & Seifert, A. (2018). McSnow: A Monte‐Carlo particle model for riming and aggregation of ice particles in a multidimensional microphysical phase space. Journal of Advances in Modeling Earth Systems, 10(1), 187–206. https://doi.org/10.1002/2017MS001167

Broadley, S., Murray, B., Herbert, R., Atkinson, J., Dobbie, S., Malkin, T., et al. (2012). Immersion mode heterogeneous ice nucleation by an illite rich powder representative of atmospheric mineral dust. Atmospheric Chemistry and Physics, 12(1), 287–307. https://doi.org/10.5194/acp‐12‐287‐2012

Burrows, S., McCluskey, C., Cornwell, G., Steinke, I., Zhang, K., Zhao, B., et al. (2022). Ice‐nucleating particles that impact clouds and climate: Observational and modeling research needs. Reviews of Geophysics, 60(2). https://doi.org/10.1029/2021RG000745

Carte, A. (1959). Probability of freezing. Proceedings of the Physical Society, 73(2), 324. https://doi.org/10.1088/0370‐1328/73/2/126

Ceppi, P., Brient, F., Zelinka, M., & Hartmann, D. (2017). Cloud feedback mechanisms and their representation in global climate models. WIREs Climate Change, 8(4), 8. https://doi.org/10.1002/wcc.465

Connolly, P., Möhler, O., Field, P., Saathoff, H., Burgess, R., Choularton, T., & Gallagher, M. (2009). Studies of heterogeneous freezing by three different desert dust samples. Atmospheric Chemistry and Physics, 9(8), 2805–2824. https://doi.org/10.5194/acp‐9‐2805‐2009

Cornwell, G., McCluskey, C., DeMott, P., Prather, K., & Burrows, S. (2021). Development of heterogeneous ice nucleation rate coefficient parameterizations from ambient measurements. Geophysical Research Letters, 48(23). https://doi.org/10.1029/2021GL095359

Curtis, J., Michelotti, M., Riemer, N., Heath, M., & West, M. (2016). Accelerated simulation of stochastic particle removal processes in particle‐resolved aerosol models. Journal of Comparative Physiology, 322, 21–32. https://doi.org/10.1016/j.jcp.2016.06.029

Curtis, J., Riemer, N., & West, M. (2024). Explicit stochastic advection algorithms for the regional scale particle‐resolved atmospheric aerosol model WRF‐PartMC (v1.0). EGUsphere. https://doi.org/10.5194/egusphere‐2024‐825

de Jong, E., Mackay, J., Bulenok, O., Jaruga, A., & Arabas, S. (2023). Breakups are complicated: An efficient representation of collisional breakup in the superdroplet method. Geoscientific Model Development, 16(14), 4193–4211. https://doi.org/10.5194/gmd‐16‐4193‐2023

de Jong, E., Singer, C., Azimi, S., Bartman, P., Bulenok, O., Derlatka, K., et al. (2023). New developments in PySDM and PySDM‐examples v2: Collisional breakup, immersion freezing, dry aerosol initialization, and adaptive time‐stepping. Journal of Open Source Software, 8(84), 4968. https://doi.org/10.21105/joss.04968

Després, V., Huffman, J., Burrows, S., Hoose, C., Safatov, A., Buryak, G., et al. (2012). Primary biological aerosol particles in the atmosphere: A review. Tellus B: Chemical and Physical Meteorology, 64(1), 15598. https://doi.org/10.3402/tellusb.v64i0.15598

DeVille, L., Riemer, N., & West, M. (2019). Convergence of a generalized weighted flow algorithm for stochastic particle coagulation. Journal of Computational Dynamics, 2019(1), 69–94. https://doi.org/10.3934/jcd.2019003

Dziekan, P., & Pawlowska, H. (2017). Stochastic coalescence in Lagrangian cloud microphysics. Atmospheric Chemistry and Physics, 17(22), 13509–13520. https://doi.org/10.5194/acp‐17‐13509‐2017

Ervens, B., & Feingold, G. (2013). Sensitivities of immersion freezing: Reconciling classical nucleation theory and deterministic expressions. Geophysical Research Letters, 40(12), 3320–3324. https://doi.org/10.1002/grl.50580

Fletcher, N. (1958). Time lag in ice crystal nucleation in the atmosphere. Part II theoretical (Vol. 1). BULLETIN DE L'OBSERVATOIRE DU PUY DE DOME. Retrieved from https://web.archive.org/web/20231213115021/https://www.phys.unsw.edu.au/music/people/publications/Fletcher1958b.pdf

Fletcher, N. (1969). Active sites and ice crystal nucleation. Journal of the Atmospheric Sciences, 26(6), 1266–1271. https://doi.org/10.1175/1520‐0469(1969)026<1266:ASAICN>2.0.CO;2

Fornea, A., Brooks, S., Dooley, J., & Saha, A. (2009). Heterogeneous freezing of ice on atmospheric aerosols containing ash, soot, and soil. Journal of Geophysical Research, 114(D13). https://doi.org/10.1029/2009JD011958

Fridlind, A., & Ackerman, A. (2018). Simulations of arctic mixed‐phase boundary layer clouds: Advances in understanding and outstanding questions. In C. Andronache (Ed.), Mixed‐phase clouds: Observations and modeling (pp. 153–183). Elsevier. https://doi.org/10.1016/B978‐0‐12‐810549‐8.00007‐6

Fridlind, A., Van Diedenhoven, B., Ackerman, A., Avramov, A., Mrowiec, A., Morrison, H., et al. (2012). A FIRE‐ACE/SHEBA case study of mixed‐phase arctic boundary layer clouds: Entrainment rate limitations on rapid primary ice nucleation processes. Journal of the Atmospheric Sciences, 69(1), 365–389. https://doi.org/10.1175/JAS‐D‐11‐052.1

Fröhlich‐Nowoisky, J., Kampf, C., Weber, B., Huffman, J., Pöhlker, C., Andreae, M., et al. (2016). Bioaerosols in the earth system: Climate, health, and ecosystem interactions. Atmospheric Research, 182, 346–376. https://doi.org/10.1016/j.atmosres.2016.07.018

Frostenberg, H., Welti, A., Luhr, M., Savre, J., Thomson, E., & Ickes, L. (2023). The chance of freezing—A conceptional study to parameterize temperature‐dependent freezing by including randomness of ice‐nucleating particle concentrations. Atmospheric Chemistry and Physics, 23(19), 10883–10900. https://doi.org/10.5194/acp‐23‐10883‐2023

Gedzelman, S., & Arnold, R. (1993). The form of cyclonic precipitation and its thermal impact. Monthly Weather Review, 121(7), 1957–1978. https://doi.org/10.1175/1520‐0493(1993)121<1957:TFOCPA>2.0.CO;2

Gedzelman, S., & Arnold, R. (1994). Modeling the isotopic composition of precipitation. Journal of Geophysical Research, 99(D5). https://doi.org/10.1029/93JD03518

Grabowski, W. (1998). Toward cloud resolving modeling of large‐scale tropical circulations: A simple cloud microphysics parameterization. Journal of the Atmospheric Sciences, 55(21), 3283–3298. https://doi.org/10.1175/1520‐0469(1998)055<3283:TCRMOL>2.0.CO;2

Grabowski, W. (1999). A parameterization of cloud microphysics for long‐term cloud‐resolving modeling of tropical convection. Atmospheric Research, 52(1–2), 17–41. https://doi.org/10.1016/S0169‐8095(99)00029‐0

Grabowski, W., Dziekan, P., & Pawlowska, H. (2018). Lagrangian condensation microphysics with Twomey CCN activation. Geoscientific Model Development, 11(1), 103–120. https://doi.org/10.5194/gmd‐11‐103‐2018

Grabowski, W., Morrison, H., Shima, S.‐I., Abade, G., Dziekan, P., & Pawlowska, H. (2019). Modeling of cloud microphysics: Can we do better? Bulletin of the American Meteorological Society, 100(4), 655–672. https://doi.org/10.1175/BAMS‐D‐18‐0005.1

Herbert, R., Murray, B., Whale, T., Dobbie, S., & Atkinson, J. (2014). Representing time‐dependent freezing behaviour in immersion mode ice nucleation. Atmospheric Chemistry and Physics, 14(16), 8501–8520. https://doi.org/10.5194/acp‐14‐8501‐2014

Hoffmann, F., & Feingold, G. (2019). Entrainment and mixing in stratocumulus: Effects of a new explicit subgrid‐scale scheme for large‐eddy simulations with particle‐based microphysics. Journal of the Atmospheric Sciences, 76(7), 1955–1973. https://doi.org/10.1175/JAS‐D‐18‐0318.1

Hoffmann, F., & Feingold, G. (2023). A note on aerosol processing by droplet collision‐coalescence. Geophysical Research Letters, 50(11). https://doi.org/10.1029/2023GL103716

Hoose, C., & Möhler, O. (2012). Heterogeneous ice nucleation on atmospheric aerosols: A review of results from laboratory experiments. Atmospheric Chemistry and Physics, 12(20), 9817–9854. https://doi.org/10.5194/acp‐12‐9817‐2012

Isaac, G., & Douglas, R. (1972). Another “time lag” in the activation of atmospheric ice nuclei. Journal of Applied Meteorology, 11(3), 490–493. https://doi.org/10.1175/1520‐0450(1972)011<0490:ALITAO>2.0.CO;2

Jaruga, A., & Pawlowska, H. (2018). libcloudph++ 2.0: Aqueous‐phase chemistry extension of the particle‐based cloud microphysics scheme. Geoscientific Model Development, 11(9), 3623–3645. https://doi.org/10.5194/gmd‐11‐3623‐2018

Jensen, E., & Pfister, L. (2004). Transport and freeze‐drying in the tropical tropopause layer. Journal of Geophysical Research, 109(D2). https://doi.org/10.1029/2003JD004022

Kanji, Z., Ladino, L., Wex, H., Boose, Y., Burkert‐Kohn, M., Cziczo, D., & Krämer, M. (2017). Overview of ice nucleating particles. Meteorological Monographs, 58, 1.1–1.33. https://doi.org/10.1175/AMSMONOGRAPHS‐D‐16‐0006.1

Kärcher, B., & Marcolli, C. (2021). Aerosol–cloud interactions: The representation of heterogeneous ice activation in cloud models. Atmospheric Chemistry and Physics, 21(19), 15213–15220. https://doi.org/10.5194/acp‐21‐15213‐2021

Kaufmann, L., Marcolli, C., Luo, B., & Peter, T. (2017). Refreeze experiments with water droplets containing different types of ice nuclei interpreted by classical nucleation theory. Atmospheric Chemistry and Physics, 17(5), 3525–3552. https://doi.org/10.5194/acp‐17‐3525‐2017

Kessler, E. (1969). On the distribution and continuity of water substance in atmospheric circulations. Meteorological Monographs, 10, 1–84. https://doi.org/10.1007/978‐1‐935704‐36‐2_1

Khain, A., Ovtchinnikov, M., Pinsky, M., Pokrovsky, A., & Krugliak, H. (2000). Notes on the state‐of‐the‐art numerical modeling of cloud microphysics. Atmospheric Research, 55(3–4), 159–224. https://doi.org/10.1016/S0169‐8095(00)00064‐8

Kilchhofer, K., Mahrt, F., & Kanji, Z. (2021). The role of cloud processing for the ice nucleating ability of organic aerosol and coal fly ash particles. Journal of Geophysical Research: Atmospheres, 126(10). https://doi.org/10.1029/2020JD033338

Knopf, D., & Alpert, P. (2013). A water activity based model of heterogeneous ice nucleation kinetics for freezing of water and aqueous solution droplets. Faraday Discussions, 165, 513. https://doi.org/10.1039/C3FD00035D

Knopf, D., & Alpert, P. (2023). Atmospheric ice nucleation. Nature Reviews Physics, 5(4), 203–217. https://doi.org/10.1038/s42254‐023‐00570‐7

Knopf, D., Alpert, P., & Wang, B. (2018). The role of organic aerosol in atmospheric ice nucleation: A review. ACS Earth and Space Chemistry, 2(3), 168–202. https://doi.org/10.1021/acsearthspacechem.7b00120

Knopf, D., Alpert, P., Zipori, A., Reicher, N., & Rudich, Y. (2020). Stochastic nucleation processes and substrate abundance explain time‐dependent freezing in supercooled droplets. npj Climate and Atmospheric Science, 3(1), 2. https://doi.org/10.1038/s41612‐020‐0106‐4

Knopf, D., Barry, K., Brubaker, T., Jahl, L., Jankowski, K., Li, J., et al. (2021). Aerosol‐ice formation closure: A southern great plains field campaign. Bulletin of the American Meteorological Society, 102(10), E1952–E1971. https://doi.org/10.1175/BAMS‐D‐20‐0151.1

Knopf, D., Charnawskas, J., Wang, P., Wong, B., Tomlin, J., Jankowski, K., et al. (2022). Micro‐spectroscopic and freezing characterization of ice‐nucleating particles collected in the marine boundary layer in the eastern North Atlantic. Atmospheric Chemistry and Physics, 22(8), 5377–5398. https://doi.org/10.5194/acp‐22‐5377‐2022

Knopf, D., Silber, I., Riemer, N., Fridlind, A., & Ackerman, A. (2023). A 1D model for nucleation of ice from aerosol particles: An application to a mixed‐phase arctic stratus cloud layer. Journal of Advances in Modeling Earth Systems, 15(10). https://doi.org/10.1029/2023MS003663

Knopf, D., Wang, P., Wong, B., Tomlin, J., Veghte, D., Lata, N., et al. (2023). Physicochemical characterization of free troposphere and marine boundary layer ice‐nucleating particles collected by aircraft in the eastern North Atlantic. Atmospheric Chemistry and Physics, 23(15), 8659–8681. https://doi.org/10.5194/acp‐23‐8659‐2023

Koop, T. (2002). The water activity of aqueous solutions in equilibrium with ice. Bulletin of the Chemical Society of Japan, 75(12), 2587–2588. https://doi.org/10.1246/bcsj.75.2587

Kotalczyk, G., Devi, J., & Kruis, F. (2017). A time‐driven constant‐number Monte Carlo method for the GPU‐simulation of particle breakage based on weighted simulation particles. Powder Technology, 317, 417–429. https://doi.org/10.1016/j.powtec.2017.05.002

Kubota, N. (2019). Random distribution active site model for ice nucleation in water droplets. CrystEngComm, 21(25), 3810–3821. https://doi.org/10.1039/C9CE00246D

Laaksonen, A., & Malila, J. (2022). Ice nucleation. In Nucleation of water (pp. 209–248). Elsevier. https://doi.org/10.1016/B978‐0‐12‐814321‐6.00018‐X

Langham, E., & Mason, B. (1958). The heterogeneous and homogeneous nucleation of supercooled water. Proceedings of the Royal Society of London. Series A. https://doi.org/10.1098/rspa.1958.0207

Lata, N., Zhang, B., Schum, S., Mazzoleni, L., Brimberry, R., Marcus, M., et al. (2021). Aerosol composition, mixing state, and phase state of free tropospheric particles and their role in ice cloud formation. ACS Earth and Space Chemistry, 5(12), 3499–3510. https://doi.org/10.1021/acsearthspacechem.1c00315

Lebo, Z. J., & Morrison, H. (2013). A novel scheme for parameterizing aerosol processing in warm clouds. Journal of the Atmospheric Sciences, 70(11), 3576–3598. https://doi.org/10.1175/JAS‐D‐13‐045.1

Lee, K., & Matsoukas, T. (2000). Simultaneous coagulation and break‐up using constant‐N Monte Carlo. Powder Technology, 110(1–2), 82–89. https://doi.org/10.1016/S0032‐5910(99)00270‐3

Leonard, J., & Im, J. (1999). Modelling solid nucleation and growth in supercooled liquid. Materials Research Society Symposia Proceedings, 580, 233. https://doi.org/10.1557/proc‐580‐233

Levine, J. (1950). Statistical explanation of spontaneous freezing of water droplets. Washington (NACATech.Note2234). Retrieved from https://core.ac.uk/download/pdf/42803258.pdf

Li, X.‐Y., Mehlig, B., Svensson, G., Brandenburg, A., & Haugen, N. E. L. (2022). Collision fluctuations of lucky droplets with superdroplets. Journal of the Atmospheric Sciences, 79(7), 1821–1835. https://doi.org/10.1175/JAS‐D‐20‐0371.1

Marinescu, P. J., vanden Heever, S. C., Heikenfeld, M., Barrett, A. I., Barthlott, C., Hoose, C., et al. (2021). Impacts of varying concentrations of cloud condensation nuclei on deep convective cloud updrafts—A multimodel assessment. Journal of the Atmospheric Sciences, 78(4), 1147–1172. https://doi.org/10.1175/JAS‐D‐20‐0200.1

Marshall, J. (1961). Heterogeneous nucleations is a stochastic process (Vol. 4). Nubila: rivista di fisica delle nubi. Retrieved from https://web.archive.org/web/20220404141323/https://cma.entecra.it/Astro2_sito/doc/Nubila_1_1961.pdf

Matsushima, T., Nishizawa, S., & Shima, S. (2023). Overcoming computational challenges to realize meter‐ to submeter‐scale resolution in cloud simulations using the super‐droplet method. Geoscientific Model Development, 16(21), 6211–6245. https://doi.org/10.5194/gmd‐16‐6211‐2023

Maxey, M., & Corrsin, S. (1986). Gravitational settling of aerosol particles in randomly oriented cellular flow fields. Journal of the Atmospheric Sciences, 43(11), 1112–1134. https://doi.org/10.1175/1520‐0469(1986)043<1112:GSOAPI>2.0.CO;2

Michel, B. (1967). From the nucleation of ice crystals in clouds to the formation of frazil ice in rivers. In International Conference on Low Temperature Science, 1966, Sapporo, Japan. Retrieved from https://web.archive.org/web/20230810010945/https://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/20291/1/1_p129‐136.pdf

Morrison, H., de Boer, G., Feingold, G., Harrington, J., Shupe, M., & Sulia, K. (2012). Resilience of persistent arctic mixed‐phase clouds. Nature Geoscience, 5(1), 11–17. https://doi.org/10.1038/ngeo1332

Morrison, H., & Grabowski, W. (2007). Comparison of bulk and bin warm‐rain microphysics models using a kinematic framework. Journal of the Atmospheric Sciences, 64(8), 2839–2861. https://doi.org/10.1175/JAS3980

Morrison, H., van Lier‐Walqui, M., Fridlind, A., Grabowski, W., Harrington, J., Hoose, C., et al. (2020). Confronting the challenge of modeling cloud and precipitation microphysics. Journal of Advances in Modeling Earth Systems, 12(8). https://doi.org/10.1029/2019MS001689

Morrison, H., Zuidema, P., Ackerman, A., Avramov, A., de Boer, G., Fan, J., et al. (2011). Intercomparison of cloud model simulations of Arctic mixed‐phase boundary layer clouds observed during SHEBA/FIRE‐ACE. Journal of Advances in Modeling Earth Systems, 3(2). https://doi.org/10.1029/2011MS000066

Mossop, S. (1955). The freezing of supercooled water. Proceedings of the Physical Society Section B, 68(4), 193–208. https://doi.org/10.1088/0370‐1301/68/4/301

Muhlbauer, A., Grabowski, W., Malinowski, S. P., Ackerman, T. P., Bryan, G. H., Lebo, Z. J., et al. (2013). Reexamination of the state of the art of cloud modeling shows real improvements. Bulletin America Meteorology Social, 94(5), ES45–ES48. https://doi.org/10.1175/BAMS‐D‐12‐00188.1

Murray, B., Broadley, S., Wilson, T., Atkinson, J., & Wills, R. (2011). Heterogeneous freezing of water droplets containing kaolinite particles. Atmospheric Chemistry and Physics, 11(9), 4191–4207. https://doi.org/10.5194/acp‐11‐4191‐2011

Niedermeier, D., Augustin‐Bauditz, S., Hartmann, S., Wex, H., Ignatius, K., & Stratmann, F. (2015). Can we define an asymptotic value for the ice active surface site density for heterogeneous ice nucleation? Journal of Geophysical Research: Atmospheres, 120(10), 5036–5046. https://doi.org/10.1002/2014JD022814

Niemand, M., Möhler, O., Vogel, B., Vogel, H., Hoose, C., Connolly, P., et al. (2012). A particle‐surface‐area‐based parameterization of immersion freezing on desert dust particles. Journal of the Atmospheric Sciences, 69(10), 3077–3092. https://doi.org/10.1175/JAS‐D‐11‐0249.1

Nordam, T., Kristiansen, R., Nepstad, R., van Sebille, E., & Booth, A. (2023). A comparison of Eulerian and Lagrangian methods for vertical particle transport in the water column. Geoscientific Model Development, 16(18), 5339–5363. https://doi.org/10.5194/gmd‐16‐5339‐2023

Paoli, R., Hélie, J., & Poinstot, T. (2004). Contrail formation in aircraft wakes. Journal of Fluid Mechanics, 502, 361–373. https://doi.org/10.1017/S0022112003007808

Petters, M., & Kreidenweis, S. (2007). A single parameter representation of hygroscopic growth and cloud condensation nucleus activity. Atmospheric Chemistry and Physics, 7(8), 1961–1971. https://doi.org/10.5194/acp‐7‐1961‐2007

Pincus, R., & Chepfer, H. (2020). Clouds as light. In A. Siebesma, S. Bony, C. Jakob, & B. Stevens (Eds.), Clouds and climate: Climate science’s greatest challenge (pp. 99–122). https://doi.org/10.1017/9781107447738.005

Pruppacher, H., & Klett, J. (2010). Heterogeneous nucleation. In Microphysics of clouds and precipitation. Springer. https://doi.org/10.1007/978‐0‐306‐48100‐0_9

PySDM developers. (2022‐2025). PySDM Zenodo archive (DOI represents all permanently archived versions, resolves to the latest one, previous releases selectable from the same site). https://doi.org/10.5281/zenodo.6321262

Rasinski, P., Pawlowska, H., & Grabowski, W. (2011). Observations and kinematic modeling of drizzling marine stratocumulus. Atmospheric Research, 102(1–2), 120–135. https://doi.org/10.1016/j.atmosres.2011.06.020

Reisner, J., Rasmussen, R., & Bruintjes, R. (1998). Explicit forecasting of supercooled liquid water in winter storms using the MM5 mesoscale model. The Quarterly Journal of the Royal Meteorological Society, 124(548), 1071–1107. https://doi.org/10.1002/qj.49712454804

Richter, D., MacMillan, T., & Wainwright, C. (2021). A Lagrangian cloud model for the study of marine fog. Boundary‐Layer Meteorology, 181(2–3), 523–542. https://doi.org/10.1007/s10546‐020‐00595‐w

Riemer, N., Ault, A., West, M., Craig, R., & Curtis, J. (2019). Aerosol mixing state: Measurements, modeling, and impacts. Reviews of Geophysics, 57(2), 187–249. https://doi.org/10.1029/2018RG000615

Rigg, Y., Alpert, P., & Knopf, D. (2013). Immersion freezing of water and aqueous ammonium sulfate droplets initiated by humic‐like substances as a function of water activity. Atmospheric Chemistry and Physics, 13, 6603–6622. https://doi.org/10.5194/acp‐13‐6603‐2013

Schmeller, G., & Geresdi, I. (2019). Study of interaction between cloud microphysics and chemistry using coupled bin microphysics and bin aqueous chemistry scheme. Atmospheric Environment, 198, 366–380. https://doi.org/10.1016/j.atmosenv.2018.10.064

Schölzel, C., & Friederichs, P. (2008). Multivariate non‐normally distributed random variables in climate research—Introduction to the copula approach. Nonlinear Processes in Geophysics, 15(5), 761–772. https://doi.org/10.5194/npg‐15‐761‐2008

Seifert, A., & Beheng, K. (2006). A two‐moment cloud microphysics parameterization for mixed‐phase clouds. Part 1: Model description. Meteorology and Atmospheric Physics, 92(1–2), 45–66. https://doi.org/10.1007/s00703‐005‐0112‐4

Seifert, A., Leinonen, J., Siewert, C., & Kneifel, S. (2019). The geometry of rimed aggregate snowflakes: A modeling study. Journal of Advances in Modeling Earth Systems, 11(3), 712–731. https://doi.org/10.1029/2018MS001519

Shima, S., Sato, Y., Hashimoto, A., & Misumi, R. (2020). Predicting the morphology of ice particles in deep convection using the super‐droplet method: Development and evaluation of SCALE‐SDM 0.2.5‐2.2.0, ‐2.2.1, and ‐2.2.2. Geoscientific Model Development, 13(9), 4107–4157. https://doi.org/10.5194/gmd‐13‐4107‐2020

Shirgaonkar, A., & Lele, S. (2012). Large eddy simulation of early stage contrails: Effect of atmospheric properties. In 44th AIAA Aerospace Sciences Meeting and Exhibit. https://doi.org/10.2514/6.2006‐1414

Silber, I. (2023). Arctic cloud‐base ice precipitation properties retrieved using bayesian inference. Journal of Geophysical Research: Atmospheres, 128(16). https://doi.org/10.1029/2022JD038202

Silber, I., Fridlind, A., Verlinde, J., Russell, L., & Ackerman, A. (2020). Nonturbulent liquid‐bearing polar clouds: Observed frequency of occurrence and simulated sensitivity to gravity waves. Geophysical Research Letters, 47(10). https://doi.org/10.1029/2020GL087099

Slawinska, J., Grabowski, W., & Morrison, H. (2009). The impact of atmospheric aerosols on precipitation from deep organized convection: A prescribed‐flow model study using double‐moment bulk microphysics. The Quarterly Journal of the Royal Meteorological Society, 135(644), 1906–1913. https://doi.org/10.1002/qj.450

Smolarkiewicz, P. (2006). Multidimensional positive definite advection transport algorithm: An overview. International Journal for Numerical Methods in Fluids, 50(10), 1123–1144. https://doi.org/10.1002/fld.1071

Sölch, I., & Kärcher, B. (2012). A large‐eddy model for cirrus clouds with explicit aerosol and ice microphysics and Lagrangian ice particle tracking. The Quarterly Journal of the Royal Meteorological Society, 136(653), 2074–2093. https://doi.org/10.1002/qj.689

Solomon, A., Feingold, G., & Shupe, M. (2015). The role of ice nuclei recycling in the maintenance of cloud ice in Arctic mixed‐phase stratocumulus. Atmospheric Chemistry and Physics, 15(18), 10631–10643. https://doi.org/10.5194/acp‐15‐10631‐2015

Stansbury, E. (1961). Stochastic freezing (Tech. Rep. No. MW‐35). McGill University.

Stevens, R., Loewe, K., Dearden, C., Dimitrelos, A., Possner, A., Eirund, G., et al. (2018). A model intercomparison of CCN‐limited tenuous clouds in the high Arctic. Atmospheric Chemistry and Physics, 18(15), 11041–11071. https://doi.org/10.5194/acp‐18‐11041‐2018

Sulia, K. J., Harrington, J. Y., & Morrison, H. (2013). A method for adaptive habit prediction in bulk microphysical models. Part III: Applications and studies within a two‐dimensional kinematic model. Journal of the Atmospheric Sciences, 70(10), 3302–3320. https://doi.org/10.1175/JAS‐D‐12‐0316.1

Sullivan, S., Barthlott, C., Crosier, J., Zhukov, I., Nenes, A., & Hoose, C. (2018). The effect of secondary ice production parameterization on the simulation of a cold frontal rainband. Atmospheric Chemistry and Physics, 18(22), 16461–16480. https://doi.org/10.5194/acp‐18‐16461‐2018

Szabó‐Takács, B. (2011). Numerical simulation of the cycle of aerosol particles in stratocumulus clouds with a two‐dimensional kinematic model. Quarterly Journal of the Hungarian Meteorological Service, 115(3), 147–165. Retrieved from https://web.archive.org/web/20220617054242/https://www.met.hu/downloads.php?fn=/metadmin/newspaper/2012/04/115‐3‐2‐szabo.pdf

Szumowski, M., Grabowski, W., & Ochs, H. T., III. (1998). Simple two‐dimensional kinematic framework designed to test warm rain microphysical models. Atmospheric Research, 45(4), 299–326. https://doi.org/10.1016/S0169‐8095(97)00082‐3

Unterstrasser, S., Hoffmann, F., & Lerch, M. (2017). Collection/aggregation algorithms in Lagrangian cloud microphysical models: Rigorous evaluation in box model simulations. Geoscientific Model Development, 10(4), 1521–1548. https://doi.org/10.5194/gmd‐10‐1521‐2017

Vali, G. (1971). Quantitative evaluation of experimental results an the heterogeneous freezing nucleation of supercooled liquids. Journal of the Atmospheric Sciences, 28(3), 402–409. https://doi.org/10.1175/1520‐0469(1971)028<0402:QEOERA>2.0.CO;2

Vali, G. (1994). Freezing rate due to heterogeneous nucleation. Journal of the Atmospheric Sciences, 51(13), 1843–1856. https://doi.org/10.1175/1520‐0469(1994)051<1843:FRDTHN>2.0.CO;2

Vali, G. (2014). Interpretation of freezing nucleation experiments: Singular and stochastic; sites and surfaces. Atmospheric Chemistry and Physics, 14(11), 5271–5294. https://doi.org/10.5194/acp‐14‐5271‐2014

Vali, G. (2019). Revisiting the differential freezing nucleus spectra derived from drop‐freezing experiments: Methods of calculation, applications, and confidence limits. Atmospheric Measurement Techniques, 12(2), 1219–1231. https://doi.org/10.5194/amt‐12‐1219‐2019

Vali, G., DeMott, P., Möhler, O., & Whale, T. (2015). Technical note: A proposal for ice nucleation terminology. Atmospheric Chemistry and Physics, 15(18), 10263–10270. https://doi.org/10.5194/acp‐15‐10263‐2015

Vali, G., & Snider, J. (2015). Time‐dependent freezing rate parcel model. Atmospheric Chemistry and Physics, 15(4), 2071–2079. https://doi.org/10.5194/acp‐15‐2071‐2015

Vali, G., & Stansbury, E. (1966). Time‐dependent characteristics of the heterogeneous nucleation of ice. Canadian Journal of Physics, 44(3), 477–502. https://doi.org/10.1139/p66‐044

Vonnegut, B. (1948). Variation with temperature of the nucleation rate of supercooled liquid tin and water drops. Journal of Colloid Science, 3(6), 563–569. https://doi.org/10.1016/S0095‐8522(48)90049‐X

Vonnegut, B. (1949). Nucleation of supercooled water clouds by silver iodide smokes. Chemical Reviews, 44(2), 277–289. https://doi.org/10.1021/cr60138a003

Welß, J.‐N., von Terzi, L., Kneifel, S., Seifert, A., & Siewert, C. (2022). Exploring the origin of increasing ice particle number in the dendritic growth zone combining polarimetric radar observations and novel Lagrangian particle modeling. In EGU22, the 24th EGU General Assembly, Held 23–27 May, 2022 in Vienna, Austria and Online. https://doi.org/10.5194/egusphere‐egu22‐5159

Wisner, C., Orville, H., & Myers, C. (1972). A numerical model of a hail‐bearing cloud. Journal of the Atmospheric Sciences, 29(6), 1160–1181. https://doi.org/10.1175/1520‐0469(1972)029<1160:ANMOAH>2.0.CO;2

Wood, R., Leon, D., Lebsock, M., Snider, J., & Clarke, A. (2012). Precipitation driving of droplet concentration variability in marine low clouds. Journal of Geophysical Research, 117(D19). https://doi.org/10.1029/2012JD018305

Wright, T., & Petters, M. (2013). The role of time in heterogeneous freezing nucleation. Journal of Geophysical Research: Atmospheres, 118(9), 3731–3743. https://doi.org/10.1002/jgrd.50365

Wright, T., Petters, M., Hader, J., Morton, T., & Holder, A. (2013). Minimal cooling rate dependence of ice nuclei activity in the immersion mode. Journal of Geophysical Research: Atmospheres, 118(18). https://doi.org/10.1002/jgrd.50810

Yang, F., Ovchinnikov, M., & Shaw, R. (2015). Long‐lifetime ice particles in mixed‐phase stratiform clouds: Quasi‐steady and recycled growth. Journal of Geophysical Research: Atmospheres, 120(22). https://doi.org/10.1002/2015JD023679

Yao, Y., Dawson, M., Dabdub, D., & Riemer, N. (2021). Evaluating the impacts of cloud processing on resuspended aerosol particles after cloud evaporation using a particle‐resolved model. Journal of Geophysical Research: Atmospheres, 126(24). https://doi.org/10.1029/2021JD034992

Zobrist, B., Koop, T., Luo, B., Marcolli, C., & Peter, T. (2007). Heterogeneous ice nucleation rate coefficient of water droplets coated by a nonadecanol monolayer. Journal of Physical Chemistry C, 111(5), 2149–2155. https://doi.org/10.1021/jp066080w