Studies on the Competition Between Homogeneous and Heterogeneous Ice Nucleation in Cirrus Formation

Cirrus ice crystals are produced heterogeneously on ice‐nucleating particles (INPs) and homogeneously in supercooled liquid solution droplets. They grow by uptake of water molecules from the ice‐supersaturated vapor. The precursor particles, characterized by disparate ice nucleation abilities and number concentrations, compete for available vapor during ice formation events. We investigate cirrus formation events systematically in different temperature and updraft regimes, and for different INP number concentrations and time‐independent nucleation efficiencies. We consider vertical air motion variability due to mesoscale gravity waves and effects of supersaturation‐dependent deposition coefficients for water molecules on ice surfaces. We analyze ice crystal properties to better understand the dynamics of competing nucleation processes. We study the reduction of ice crystal numbers produced by homogeneous freezing due to INPs in both, individual simulations assuming constant updraft speeds and in ensemble simulations based on a stochastic representation of vertical wind speed fluctuations. We simulate and interpret probability distributions of total nucleated ice crystal number concentrations, showing signatures of homogeneous and heterogeneous nucleation. At typically observed, mean updraft speeds (≈15 cm s−1) competing nucleation should occur frequently, even at rather low INP number concentrations (<10 L−1). INPs increase cirrus occurrence and may alter cirrus microphysical properties without entirely suppressing homogeneous freezing events. We suggest to improve ice growth models, especially for low cirrus temperatures (<220 K) and low ice supersaturation (<0.3).

Significant progress has been achieved in quantifying vertical air motion variability due to mesoscale gravity waves based on balloon observations  and in formulating molecular-level physical models for the uptake of H 2 O on ice crystal surfaces based on laboratory measurements (Harrington et al., 2019). Climate models have made great strides in representing cirrus in their cloud schemes (Bardeen et al., 2013;Muench & Lohmann, 2020;Zhu & Penner, 2020). Quantitative evaluation of aerosol-cirrus interactions and man-made cirrus thinning requires full understanding of formation of ice crystals from heterogeneous ice nucleating particles (INPs), as this affects the performance of the cirrus parameterizations used in these models.
A number of atmospheric aerosol types nucleate ice in cirrus conditions (Hoose & Möhler, 2012;Kanji et al., 2017). While ubiquitous supercooled liquid solution (aqueous aerosol) droplets induce homogeneous freezing at high ice supersaturation (defined as s = S − 1, where S is the fractional relative humidity over ice), s > 0.45 below about 235 K (Koop et al., 2000), INPs may form ice heterogeneously at similar or lower s-values. Field observations suggest competition between INPs and solution droplets in the extratropical northern hemisphere (Cziczo et al., 2013;Haag et al., 2003;Jensen et al., 2013) and the tropical tropopause layer (TTL; Ueyama et al., 2015), as indicated by frequency distributions of ICNCs, in-cloud relative humidities, and residual aerosol particle compositions. Heterogeneous ice nucleation may dominate in the polluted Northern Hemisphere (Cziczo et al., 2013) and homogeneous freezing may dominate in the cleaner Southern Hemisphere (Gayet et al., 2006). Homogeneous ice formation dominates in wave clouds due to large updraft speeds, such that any competition for H 2 O by heterogeneous nucleation is overcome as high cooling rates drive high supersaturations (DeMott et al., 1998).
Besides ice-forming aerosol particles, updraft speeds and associated adiabatic cooling rates have long been known to be crucial for cirrus (DeMott et al., 1997;Heymsfield, 1977;Jensen & Toon, 1994). High-flying aircraft and ground-based radar observations pointed to significant enhancements of vertical wind speeds over synoptic values (Hoyle et al., 2005;, but only recently long-duration, quasi-Lagrangian balloon measurements quantified mesoscale temperature fluctuations (MTF) induced by gravity waves. This enabled the development of an approximate probabilistic model describing the Lagrangian evolution of damped MTF at arbitrary latitudes (Kärcher & Podglajen, 2019). This stochastic approach describes background fluctuations away from high-frequency wave sources consistent with observations and is employed here allowing us to perform parametric studies with well-defined mean values of the underlying wind speed fluctuations.
How rapidly H 2 O is incorporated in ice crystals during ice formation events is controlled by the deposition coefficient, α. It is through the impact on supersaturation that INPs may hinder or even suppress homogeneous freezing. The gas-phase diffusion of H 2 O in air toward ice crystals and molecular processes at their surfaces together determine the rate of irreversible vapor uptake in an ice-supersaturated environment. Laboratory measurements have shown that deposition coefficients vary with ambient conditions (i.e., ice supersaturation and temperature), but discrepancies between laboratory data at cirrus temperatures exist and have not been reconciled (Asakawa et al., 2014;Harrington et al., 2019;Nelson, 2001). The development of ice crystal shapes (habits) are caused by variations in deposition coefficients across ice crystal surfaces (Chen & Lamb, 1994); variability in α may also cause microphysical-dynamical feedbacks in simulations of mixed-phase and cirrus clouds relevant for their evolution (Ervens et al., 2011;Kärcher, 2020). Nonetheless, most cloud models resort to constant α in the absence of a complete theory of crystal growth from the vapor encompassing the full cirrus regime, from temperatures, T, below about 233 K at the boundary to mixed-phase clouds to values prevailing in the TTL (<200 K). Here, we follow a hybrid strategy by employing on the one hand variable deposition coefficients that are based on laboratory data for T > 220 K (Harrington et al., 2019) and on the other hand assume ice crystals to be spherical.
With regard to the dependency of deposition coefficient or fresh ice crystal habits on the mode of nucleation, Bailey and Hallett (2002) showed that ice crystals formed on kaolinite and grown in a static diffusion chamber have habits that are quite similar to those observed in the atmosphere, contrary to crystals formed on silver iodide. Thus, this study showed that the manner of nucleation can affect the subsequent crystalline habit. Levitation diffusion chamber studies also indicate that homogeneously frozen droplets have deposition coefficients that are distinctly different from crystals nucleated heterogeneously (Pokrifka et al., 2020). These studies suggest that some of the crystals formed from homogeneously frozen drops have deposition coefficients that evolve from high 10.1029/2021JD035805 3 of 21 to low values during the experiment. Studies with heterogeneously frozen drops do not show a strong decrease of the deposition coefficient as the particle grows, suggesting that the nucleation mechanism played a role in the subsequent growth of the crystals.
Recent balloon-borne in situ sampling of cirrus ice crystals, combined with cold-stage electron microscopy, reveals that polycrystalline ice crystals (polycrystals) dominate the habit of cirrus ice crystals (Magee et al., 2021). Even at small sizes, these crystals are sharply faceted and often with hollowed regions. This indicates that the primary growth mechanism most likely was step nucleation near the facet edges, otherwise symmetric hollows would not be possible. Moreover, many of the hollow regions are circular, indicating that the step edges are strongly roughened. This suggests modeling cirrus ice crystals with step nucleation.
To facilitate process-oriented analyses of cirrus ice formation processes including INPs, we develop a model that only considers the most essential variables and processes operating in ice formation events in Section 2. Our focus is on investigating how INP-derived ice crystals slow the rate of supersaturation increase in cooling air and thereby alter the number of ice crystals produced by homogeneous freezing-a key component in parameterizations of cirrus ice formation. We also delineate the conditions under which INPs quench rising supersaturation and thereby prevent homogeneous freezing. We go one step further than previous approaches and employ both variable deposition coefficients and vertical wind speeds. We introduce basic concepts and quantify conditions in which INPs modify ICNCs in Section 3. Simulation results are presented and discussed in Section 4. Section 5 summarizes our main findings and reiterates important discussion points.

Approach and Approximations
To study INP effects in cirrus clouds, we use an air parcel model framework. Latent heat release due to water phase changes is small at temperatures below 233 K in the absence of cloud droplet freezing and growth (Kärcher, 2017b). Radiative heating of ice crystals is unimportant on the short time scales of cirrus ice formation events.
Neglecting effects of sedimentation is approximately justified for short simulation times. When INP numbers are low causing ice crystals to grow to relatively large sizes and fall out of the nucleation region represented by an air parcel (Murphy, 2014), this may lead to an overestimation of INP effects. We investigate INP effects in two well-separated temperature regimes, 220-225 K ("warm") and 200-205 K ("cold"). INPs are known to operate in various nucleation modes and exhibit stochastic (time-dependent) and deterministic (time-independent) ice nucleation behavior (Vali et al., 2015). Bare (uncoated) INPs, found close to emission sources, nucleate ice deterministically. As such INPs age in the atmosphere, they turn into immersion nuclei by chemical processing and coagulation. A deterministic description also is appropriate for immersion freezing, as the stochastic component does not significantly affect freezing temperatures . Here, we parameterize associated frozen (ice-active) INP fractions, respectively, with a simple analytical function.
In the parcel approach, ice crystal settling is ignored and it is assumed implicitly that INPs cover an entire ice supersaturated layer all the time. A model with vertical resolution would be more realistic, but then the problem would depend on the vertical profiles of moisture and temperature, as well as on the location of the INPs relative to supersaturation maxima. While this approach is suitable for the goal of our study, we may underestimate effects of ice crystal settling and vertical supersaturation variability in real cirrus formation events.
It is not necessary to explicitly simulate homogeneous freezing in a population of size-dispersed solution droplets in order to find the conditions in which INPs suppress homogeneous freezing. Monitoring if and when ice supersaturation reaches a value where homogeneous freezing takes place in the presence of INPs suffices for this purpose. This threshold ice supersaturation depends on the liquid water volume in solution droplets, on air temperature, and on the cooling rate. The dependence on cooling rate arises because the rate of supersaturation increase grows in proportion to the updraft speed. As measured cirrus ICNCs rarely exceed 1 cm −3 , and such numbers are most likely caused by homogeneous freezing, the number of fully soluble aerosol particles is not a limiting factor in homogeneous freezing.
The competition between INPs and solution droplets in ice formation events is controlled by the uptake of H 2 O on INP-derived ice crystals. We represent a given INP type by classes with different ice-active number concentrations belonging to a sequence of ice supersaturation increments. INPs activating at low supersaturation form a 10.1029/2021JD035805 4 of 21 cohort of ice crystals that grows ahead of those forming at higher s-values. In this way, we create a size-resolved ice crystal spectrum accounting for the possibility that deposition growth of early formed ice crystals changes the supersaturation conditions for later nucleation.
We simulate H 2 O mass uptake during the rather short ice formation stage by assuming a simplified ice growth model and approximating the shape of real ice crystals by volume-equivalent spheres (Section 2.4). To date, it is not possible to fundamentally resolve the issue of ice growth in view of the complexity of cirrus ice particle morphology, complexity, and mesoscopic surface roughness on the one hand (Magee et al., 2021) and due to the absence of a consensus on how real-world ice crystal habits can be predicted on the other. However, we drop the common assumption of constant ice growth efficiency and use deposition coefficients on ice that depend on ice supersaturation and other variables. In view of the lack of measurements constraining H 2 O uptake models at temperatures below about 220 K, we approach the extrapolation of such models to lower temperatures with caution.

Model Equations
We first introduce the basic model equations. Equation 1 describes the dry adiabatic cooling of the air parcel by prescribed time series of vertical wind speeds (Section 2.3). Equation 2 tracks the number of H 2 O molecules per ice crystal nucleated in the air parcel determining the mean size of ice crystals in each class. Equation 3 predicts the evolution of ice supersaturation, containing a forcing term that is consistent with Equation 1, and a loss term due to kinetically corrected H 2 O diffusion toward, and deposition onto, ice crystals (Section 2.4): In the above equations, t is the time and w is the vertical wind speed; Γ = g/c p is the absolute value of the dry adiabatic lapse rate with the acceleration of gravity, g, and the isobaric specific heat capacity of dry air, c p ; n sat (T) is the H 2 O number concentration at ice saturation taken from Murphy and Koop (2005), valid for T > 200-210 K (Nachbar et al., 2019); D is the effective H 2 O diffusion coefficient in air (see below) modified to account for kinetic effects of H 2 O uptake on ice crystal surfaces. Air pressure, p, follows directly from the adiabatic relationship, ∝ ( ∕ ) , with the specific gas constant for air, R. The index k denotes an INP class associated with ICNCs, n k , resulting from nucleation of the fraction of INPs that become ice active in a supersaturation interval Δs k (Section 2.5), and N k is the number of H 2 O molecules per ice crystal derived from INPs in each s-class. We infer the corresponding ice crystal radii, r k , from a spherical volume with the INP core at its center: with the volume of one H 2 O molecule in ice, ν. In Equation 4, the radius of the dry aerosol particle core, r c , that is, the initial ice crystal radius, is fixed at a typical INP size, 0.2 μm; results are not sensitive to r c . The thermodynamic parameter a describes the production of supersaturation: where L(T) is the latent heat of sublimation (Murphy & Koop, 2005) and R v is the specific gas constant for H 2 O. Finally, D k is given by: where D v (p, T) is the H 2 O diffusion coefficient in air, ℓ is the H 2 O jump distance (roughly equal to its mean free path), d = 4D v /v is its diffusion length scale with the mean thermal speed of H 2 O, v, and α k is the deposition coefficient belonging to ice crystals generated within Δs k . All p-and T-dependent variables were taken from Lamb and Verlinde (2011).
For α k r k /d ≪ 1, the net uptake rate Equation 2 approaches the Hertz-Knudsen equation, ∕ = 2 sat (free molecular limit). The supersaturation sink term in Equation 3 is proportional to the number concentration of nucleated ice crystals, n k , diagnosed from s based a time-independent (deterministic) nucleation model (Section 2.5).
The system of K + 2 coupled ordinary differential Equations 1-3 is integrated forward in time by means of the initial value solver VODE using an implicit method with functional iteration (Brown et al., 1989). We initialize each simulation in cloud-free, ice-saturated air, prescribing N k,0 = s 0 = 0 along with T 0 = 205 (225) K and p 0 = 100 (250) hPa; the actual homogeneous freezing temperatures are approximately 200 (220) K. ICNCs are initially zero, n k,0 = 0, and evolve over time depending on the nucleation model.

Vertical Air Motions
Updrafts create ice-supersaturated conditions in which aerosol particles nucleate ice. We introduce fundamental concepts and illustrate basic processes based on simulations using constant updrafts. We represent vertical air motions and associated adiabatic fluctuations in T and s caused by mesoscale gravity waves in simulations with more realistic forcings based on a stochastic approach.
To incorporate wave effects, we sample vertical wind speeds, w, from a two-sided exponential distribution , using an autocorrelation time of 2.78 min guided by theoretical arguments to create Lagrangian time series, w(t; Kärcher & Podglajen, 2019). The autocorrelation time of the first order autoregressive time series is based on an average upper tropospheric buoyancy frequency of 0.012 s −1 . Damping of MTF is neglected owing to the short duration of the time series used here. Important for ice nucleation studies, the time series contain contributions from high-frequency waves inducing large cooling/heating rates that are absent in synoptic air parcel trajectories.
The probability distribution of w-values is characterized by a mean updraft speed, w m . To study INP effects for different forcings, we vary w m and simulate a number of w(t)-trajectories, allowing us to compute frequency distributions of nucleated ICNCs.

Depositional Growth
Equation 2 represents the capacitance model for vapor-grown spherical ice crystals (Lamb & Verlinde, 2011). The deposition coefficient, α, accounts for a suite of kinetic processes (formation of surface steps, partially disordered layers, etc.) that control the incorporation of surface H 2 O molecules into the bulk ice crystal lattice. Based on analyses of vapor growth rates encompassing a wide range of droplet properties as well as nucleation and environmental conditions, laboratory measurements show that α ranges from about 0.001 to unity (Harrington et al., 2019). In line with recent balloon-borne measurements that show details of cirrus ice crystal surfaces (Magee et al., 2021), we model depositional growth of newly nucleated ice crystals with step nucleation.
However, it must be kept in mind that step nucleation is considered the dominant growth mode for larger crystals, and especially those that have developed pronounced habit forms (Frank, 1982). Newly nucleated ice crystals initially grow by dislocations and then transition to growth by step nucleation (Nelson & Knight, 1998). Indeed, Harrington and Pokrifka (2021) model the early growth of frozen solution droplets with a combination of facet spreading across the crystal surface and dislocations. Unfortunately, current theories do not provide for a general way to model a transition from one growth mode to another. Another complication arises from the fact that for real (non-spherical) crystals, although one facet may be growing by step nucleation, another facet could be growing by dislocations. It must also be kept in mind that polycrystals, which are often found in cirrus, grow by dislocations that propagate away from the grain boundaries (Pedersen et al., 2011) producing faster growing facets. Stacking faults are prevalent in ice crystals, and can lead to the enhanced growth of facets.
As ice crystals grow and evolve, α changes. However, immediately following nucleation, the development of facets along the crystal surface (lateral growth) alters the growth rate (Harrington & Pokrifka, 2021;Nelson & Swanson, 2019). Lateral growth may be important for the development of complex crystal shapes that occur in cirrus, such as scrolls, stacks of sheaths, and trigonal crystals (Nelson & Swanson, 2019), and does not lead to a significant increase in the maximum dimension of ice crystals.
Since it is not possible to model the above growth modes in detail, we apply a parameteric model of surface processes (Nelson & Baker, 1996): The critical surface supersaturation, s crit (T), in the hyperbolic tangent controls the transition from inefficient to efficient vapor growth, while the parameter μ controls how rapid this transition occurs. The ice supersaturation taken directly above the ice crystal surface is given by Lamb and Verlinde (2011): reminiscent of the gas kinetic correction to D v in Equation 6. Together, these equations mean that α k is determined by diffusion of H 2 O molecules through the air and by surface kinetic processes. We solve Equations 7 + 8 iteratively for α k , that is, separately for each ice crystal class, with and μ = 10 (Harrington et al., 2019), supported by in situ data (Magee et al., 2021).
Prior work demonstrated that the overall H 2 O uptake by nonspherical ice crystals is replicated employing a surface-average s crit -value and the spherical approximation in growth calculations (Zhang & Harrington, 2014). We emphasize that s crit -values for low cirrus temperatures (<220 K) are estimates; we limit T in Equation 9 to −70°C to prevent estimating implausible (too large) α-values. Moreover, as noted above, it is unclear how growth modes are linked at low supersaturation. We refer to further discussion in Section 4.3 and do not recommend using this formulation for low temperatures at this stage.

Heterogeneous Ice Nucleation
Heterogeneous ice nucleation does not occur abruptly at a specific supersaturation, as often assumed in cirrus parameterizations and models, but across a range of s-values, δs, around a characteristic value, s ⋆ . We introduce a function, ϕ(s), bounded by 0 and 1, that represents ice-active INP fractions, and associated ICNCs, cumulated up to s: The total INP-derived ICNC, n ⋆ is defined here as the maximum concentration determined by the maximum ice-active fraction of the INPs that is ever possible. An INP population with a total number concentration n tot may only reach a maximum ice-active fraction ϕ max < 1, in which case n ⋆ = n tot ϕ max . Introducing n ⋆ as a common scaling factor and at the same time using ϕ bounded by unity allows us to study and compare effects for a range of total number concentrations for the various INP types. The latter vary widely in nature. Without loss of generality, we neglect a possible explicit dependence of ϕ on T. Here we use: The location parameter s ⋆ and scale parameter δs have clear physical meanings: s ⋆ describes the 50% activation point (i.e., ϕ(s ⋆ ) = 0.5) and δs is related to the slope of ϕ at s ⋆ . These parameters can be varied parametrically to represent efficiency and abundance of different INP types. We determine their values to model two INP types mimicking the effects of good INP forming ice at low-to-medium s and poor INP forming ice at medium-to-high s (relative to s hom ).
To motivate the use of Equation 11, we compare it to measurement-constrained parameterizations of ice-active fractions for desert dust (Ullrich et al., 2017) and aircraft soot , both taken at 220 K. Figure 1 shows the parameterization results for monodisperse dust (1 μm, orange curve) and soot (0.4 μm, magenta) particles along with our representations of ice-active fractions for good (s ⋆ = 0.3, black) and poor (s ⋆ = 0.4, blue) model INPs. The s ⋆ -value for good INPs represents the lower part of nucleation thresholds for mid-latitude cirrus formation estimated from field data. The s ⋆ -value for poor INPs is chosen such that ϕ approaches unity right before homogeneous freezing would set in. Both cases use δs = 0.03, meaning that ice formation in both INP types overlaps only to a small degree and at the same time approximately bracket the monodisperse dust and soot curves, which are steeper and flatter that the modeled ϕ-curves, respectively. This lets us clearly distinguish between good and poor INPs, but some atmospheric INPs may exhibit an activation behavior that is a combination of both types (i.e., characterized by large s ⋆ and large δs). A similar behavior is mimicked by sea salt/spray particles when investigated at low temperatures (Patnaude et al., 2021;Wagner et al., 2018). The hyperbolic tangent in Equation 11 serves as a reasonable approximation for ice-active fractions once fitted to data. For example, the dust curve in Figure 1 is well reproduced setting s ⋆ = 0.352 and δs = 0.0175 . Changes in s ⋆ and δs occur for both aerosol types when accounting for INP size distributions.
In the numerical simulations, the deterministic INP spectrum, n(s), is divided into various ice supersaturation classes. For this purpose, we define an s-grid with spacing Δs k = s k − s k−1 , allowing us to calculate the associated, partial ice-active number concentrations: once s > s k−1 , INPs turn into ice crystals with ICNCs n k . This allows us to simulate how early nucleating INPs affects the development of supersaturation and thereby the activation and growth of later nucleating INPs. To discretize ϕ, we define a total number of K classes around s ⋆ such that Δs k is constant (Δs). The value of K is chosen such that the activation curve is represented with high resolution around the central value s ⋆ , where ϕ rises most steeply. We recall that the total INP number concentrations, n ⋆ , are prescribed at ice saturation, but nucleation occurs at s > 0. To work with well-defined n ⋆ -values during cirrus ice formation, we neglect the small reduction of n ⋆ due to adiabatic change. Equations 10 and 12 describe deterministic ice nucleation in ice-supersaturated conditions.
It is permissible to determine ICNCs based on s-cumulative ice-active fractions, since we focus on a single ice formation event that does not require removing INP after nucleation . We account for non-monotonically increasing supersaturation histories, s(t), by ensuring that n k -values do not decrease when

Supersaturation for Homogeneous Freezing-Relaxation
We derive a value, s hom , where homogeneous freezing of supercooled solution droplets takes place. This value should ideally quantify the supersaturation where most of the solution droplets freeze and quench the supersaturation, thereby terminating the homogeneous nucleation event (freezing-relaxation; Kärcher & Lohmann, 2002a).
Homogeneous freezing is a stochastic process, where the number of solution droplets diminishes over time in . Here j = VJ is the freezing rate, with the liquid water volume of the droplet population, V, and the rate coefficient, J, taken from a water activity-based formulation (Koop et al., 2000). As the dependence of homogeneous freezing events is insensitive to changes in liquid aerosol particle size distributions (Kärcher & Lohmann, 2002a), we assume a log-normal droplet population with modal radius, r m , and geometric standard deviation, σ. The kth radial moment of a log-normal distribution is given by = exp(0.5 2 ln 2 ) , so that V = 4πM 3 /3. Solution droplets in cirrus conditions are in equilibrium with ambient H 2 O (except for large droplet sizes in vigorous updrafts), so that J = J(S w , T), where S w is the liquid water saturation ratio (fractional relative humidity).
Freezing-relaxation unfolds on a timescale τ defined as: As water vapor uptake is small due to hygroscopic droplet growth during the freezing event, we neglect the volume term in Equation 13 and rewrite the remaining contribution to τ as the first term, when evaluated with the Clausius-Clapeyron equation, is shown to be smaller than the second term, leaving us with with the absolute value of the cooling rate, ̇= | ∕ | . We use the analytical fit ζ [K −1 ] = 304.4 − T(2 − 0.004T) provided by Ren and Mackenzie (2005) to evaluate τ. For given r m , σ, T, and ̇ , the desired values for s hom follow by iterating the relationship: where we neglect the volume of soluble or insoluble material in the solution droplet, justified in so far as the total volume of solution droplets is much larger than that of their dry particle core.
Equation 16 is a thermodynamic relationship, since the water contents of the solution droplets in a given size range are assumed to stay in equilibrium with ambient H 2 O. It ignores that small solution droplets have a smaller liquid water mass fraction, hence water activity and freezing rate, than larger droplets due to the Kelvin effect. Recent measurements showed that onset values for homogeneous freezing may be higher than predicted by Koop et al. (2000) at TTL temperatures (185-205 K; Schneider et al., 2021), but the effect of associated changes in s hom on nucleated ICNCs is small. Figure 2 shows s hom as a function of updraft speed for selected temperatures and monodisperse solution droplets. The threshold s hom varies most notably with T and is less sensitive to variations in cooling rate and aerosol properties owing to the strong T-dependence of J. Here we apply Equation 16 as a function of T and w for r m = 0.25 μm and σ = 1.
Equation 16 does not capture two kinetic effects: (a) large solution droplets need more time to equilibrate with the vapor phase than small droplets due to greater diffusional H 2 O uptake resistance, which introduces a time lag in the freezing rate in the presence of large cooling rates; and (b) vapor depletion by the nucleated ice crystals is associated with a characteristic time scale, which may increase the duration of the freezing event and hence the nucleated ICNCs. Moreover, the probability of freezing for large droplets is higher than for small droplets; depending on their number concentration, early freezing solution droplets grow ahead of the smaller ones and thereby modify the evolution of supersaturation and freezing conditions for the late-freezing droplets. For these reasons, numerical estimates for s hom obtained with aerosol size-resolved, non-equilibrium microphysical models may differ slightly from those obtained by means of Equation 16 to the degree these effects come into play. Again, the impact of these effects on nucleated ICNCs is small due to the self-terminating nature of freezing-relaxation (Kärcher & Lohmann, 2002b).

Suppression of Homogeneous Freezing by INPs
We introduce the quenching velocity, w ↓ , based on Equation 3: as a measure of the H 2 O sink due to growing ice crystals derived from INPs (Kärcher et al., 2006). It encapsulates the dynamical and microphysical aspects of our problem. We may refer to W = w − w ↓ as an effective vertical wind speed, the one relevant for homogeneous freezing in the presence of INPs.
Homogeneous freezing takes place when the condition s = s hom is reached for W > 0. In the other limit, W (initially positive) changes sign at some point, causing s to pass through a maximum and never reaching homogeneous freezing conditions. The two limits are separated by setting s = s hom ; for given w, homogeneous freezing occurs in the presence of INPs when the condition is met. Equation 18 reveals that the role INPs play in the competition for supersaturated vapor in an ice formation event is determined by their intrinsic ice nucleation ability, the mean size of the population of ice crystals deriving from them, and the deposition coefficient via H 2 O volume diffusion and surface kinetic processes. However, even in the case of a constant updraft, Equation 18 cannot be readily evaluated, since for given s(t), D k , and r k depend on the H 2 O uptake history.

Reduction of Homogeneously Produced Ice Crystal Numbers
To explore the impact of INPs on the number concentration of homogeneously nucleated solution droplets, n hom , we employ a cirrus parameterization (Kärcher & Lohmann, 2002a) and apply it to monodisperse droplets (r m = 0.25 μm) in cases where homogeneous freezing is activated (i.e., at the instant where simulated s(t) reaches s hom ). We update the original n hom -parameterization by using a deposition coefficient of 0.7 to represent the homogeneous freezing event (Skrotzki et al., 2013); we use the ice supersaturation thresholds from Section 2.6; and we employ a convenient algebraic fit function proposed by Ren and Mackenzie (2005) (their Equation 24) to speed up the computation of the analytical solution for n hom .
During the simulations with stochastic wave forcing, it may happen that w changes sign from positive to negative when s lies within the small range of values where homogeneous freezing takes place. The reduction of homogeneously produced ice crystals in such relatively rare non-persistent cooling events (Dinh et al., 2016;Jensen, Ueyama, Pfister, Bui, Alexander, et al., 2016;Kärcher & Jensen, 2017) are not represented in this approach.

Regime Definitions
We have already introduced two cirrus regimes categorized by ambient temperature, T, in Section 2.1 and two ice activity schemes depending on the activity location parameter s ⋆ of INPs in Section 2.5. Furthermore, two fundamental ice nucleation regimes, separated by the parameter ω = w /w, describe the competition between heterogeneous ice nucleation and homogeneous freezing in cirrus: the freezing regime (ω < 1) and the quenching regime (ω > 1), wherein INP-derived ice crystals either allow or prevent s from reaching s hom . The ability of INPs to quench s increases as n ⋆ increases (increasing H 2 O uptake rates) and diminishes as T decreases (decreasing H 2 O uptake rates).
Freezing and quenching regimes relate to two distinct factors affecting in situ cirrus formation : (a) Microphysical factor-ice-forming aerosol particles. The ability of INPs to nucleate ice crystals is strong/weak when ice-active fractions increase over a range of s-values (δs) around a low/high characteristic value (s ⋆ ). (b) Dynamical factor-small-scale updraft speeds. The ability of INPs to quench s and prevent homogeneous freezing is for a given updraft speed (w) strong/weak for high/low total INP number concentrations (n ⋆ ) and for high/low temperatures (T).
The presence of INPs may not be consequential when their impact on supersaturation is weak, for example, due to low number concentrations or slow depositional growth. In that case, even good INPs are inefficient when they are not capable of preventing homogeneous freezing. However, even poor or inefficient INPs may exert a notable effect on cirrus microphysical properties, since homogeneous freezing takes place at a lower effective updraft speed than in the unperturbed cloud, producing fewer but on average larger ice crystals. Table 1 summarizes the regimes addressed in our study along with defining parameter values. INPs are regarded as efficient when they are abundant enough to quench ice supersaturation before homogeneous freezing occurs, depending on the updraft speed for a given set of values {T 0 , n ⋆ , s ⋆ }. Otherwise they are viewed as inefficient.
Here, we set out and explore the evolution of the vapor-ice system in different regimes.

Results and Discussion
We present and discuss numerical solutions to the K + 2 model Equations 1-3 for representative values of w m and n ⋆ . Most w m -values found in field observations lie in the range 10-20 cm s −1 , while individual values may cover a broader range (Kärcher & Podglajen, 2019). Total INP number concentrations n ⋆ can vary widely, mostly around a few 10 L −1 up to 100 L −1 and peak values of several 100 L −1 in anvil cirrus (Prenni et al., 2007). In this section, one set of simulations represents the freezing regime (n ⋆ = 1 L −1 , inefficient INPs), another set with higher number concentration (50 L −1 , efficient INPs) illustrates the quenching regime.

Analysis of Ice Microphysics
We illustrate the microphysical processes in action with constant updraft speeds. Figure 3 shows the temporal evolution of s, r, and α for the "good" INP type (s ⋆ = 0.3) in both warm (220-225 K) and cold (200-205 K) cirrus regimes in a constant updraft (w = 15 cm s −1 ). We show two sets of curves for r and α, representing Lagrangian tracking of ice crystals that formed at early (class k = 10, s 10 = 0.25) and late (s 30 = 0.35) stages of the evolution; the K = 41 individual activation classes are separated by Δs = 0.005). The underlying ice activation curve is shown in Figure 1 (black curve). We first discuss the warm cirrus case (top panel in Figure 3). In the freezing regime due to low ICNCs up to 1 L −1 (solid curves), ice supersaturation increases approximately linearly up to the time where homogeneous freezing relaxation sets in and the simulation is terminated (at t ≃ 45 min, s hom ≃ 0.52). Ice crystals forming in classes 10 (30) at approximately 23 (31) min, experience different supersaturation conditions during deposition growth and acquire different radii, 38 (33) μm, at the point where homogeneous freezing sets in. In both cases, the deposition coefficient peaks close to unity for the very short time where the crystals remain small (μm-sized), because s-values exceed the critical supersaturation, s crit . During the growth phase with s still increasing, α decreases to values slightly above 0.03 prior to homogeneous freezing.
In the quenching regime with ICNCs up to 50 L −1 (dashed curves), ice crystals derived from the INPs are abundant enough to prevent homogeneous freezing. Ice crystal formation from INPs in the high s-class is slightly delayed due to the growth of INPs already active at lower s, compare to solid and dashed green curves. The ice supersaturation peaks at about 0.35 and then diminishes due to ongoing ice crystal growth. Deposition coefficients of ice crystals from the two s-classes evolve similarly initially owing to comparable supersaturation histories, but they reduce to values near 0.01 at 80 min. Deposition growth slows down accordingly; a quasi-steady state of the ice-vapor system is not reached as long as deposition coefficients change.
The fact that sizes of INP-derived ice crystals differ depending on the ice supersaturation at which the associated INPs activate underscores the importance of defining multiple supersaturation classes in order to accurately simulate the supersaturation history and, by inference, homogeneous freezing and the competition between the various ice nucleation modes. The same holds true for pure homogeneous freezing events, where early freezing large solution droplets influence later freezing smaller droplets, ultimately self-terminating the freezing event (Kärcher & Lohmann, 2002a). This represents a great challenge for numerical cloud models that are based on bulk microphysics, where ice nucleation across a particle size or supersaturation spectrum is parameterized.
We turn to the cold cirrus case (bottom panel). In the freezing regime, the supersaturation history is qualitatively similar to the warm case, with three exceptions. The rate of increase is slightly faster and homogeneous freezing commences slightly later at a higher supersaturation (t ≃ 40 min, s hom ≃ 0.58). Ice crystals remain considerably smaller with radii ≈10 μm when homogeneous freezing sets in due to the lower available H 2 O content in colder air. As a result, s peaks at a larger value (≈0.51) in this regime, indicating a smaller supersaturation sink (smaller quenching velocity, w ↓ ) compared to warm cirrus. The different initial behavior of the deposition coefficient in the early and late activating cases is due to the ice supersaturation. When ice initially forms in the early activation case, s-values are low and close to s crit causing low α-values initially; as s rises, α also increases. On the other hand, α declines with time in the late activating case, because s is higher when ice forms leading to deposition coefficients near unity. The ice supersaturation remains high after ice formation in this case, and therefore does not greatly affect α. However, ice crystal size increases with time causing α to decrease. As in the warm case, α-values converge for early and late activating INPs due to the combination of decreasing s and increasing r, approaching a higher value ≈ 0.1 owing to lower T. In the quenching regime, the subsequent evolution of s, r, and α is analogous to the warm cirrus regime. Figure 4 shows results of simulations in the two temperature regimes carried out with good INPs (s ⋆ = 0.3) to illustrate the evolution of the supersaturation sink represented by w ↓ . We choose INPs with a total number concentration n ⋆ = 1 L −1 and prescribe updraft speeds in a range around observed mean values, so that in all cases the freezing regime is realized (ω < 1). For both temperatures at the lowest w, w ↓ is very close to w; for slightly smaller w we would enter the quenching regime. Furthermore, it takes longer for w ↓ to become significant as w decreases; peak w ↓ -values increase, since s hom is reached later and ice crystals have more time to grow to larger sizes. Figure 5 shows results of simulations in the warm cirrus regime, illustrating the evolution of the ice water content (IWC) of the INP-derived ice crystals: where ϱ is the mass density of bulk ice. Two sets of curves illustrate the freezing regime realized for the two total INP number concentrations and for good and poor INPs. In the case of high n ⋆ , ice crystals are efficient vapor sinks and IWC increases rapidly to values greater or equal to 28 mg m −3 . In the case of low n ⋆ , deposition growth is kinetically limited (IWC stays below 8 mg m −3 ) and the difference between good and poor INPs is more discernible. A notable difference between good and poor INP is only visible for the low updraft speed.
While ice crystal sedimentation may modify the above results at low updraft speeds, they suggest that n ⋆ is the primary control for any heterogeneous INP type, poor INPs in sufficient numbers being as important as good INPs. This is relevant, for example, for situations like biomass burning plumes, where the efficiency of INPs could be low but enhancements of INP number concentrations are large.
From a fundamental perspective, the frequent use of constant α in models of ice growth from the vapor phase is an approximation that is only valid for a small range of conditions. From a practical perspective, this tends to overestimate the depositional growth rate, especially for small ice crystals. Relating to this study, doing so may either under-or overstate the role of INPs in cirrus formation. As to the consequences of the impact of cirrus on the moisture budget or for simulated cirrus properties, a general answer cannot be given, because the consequences of the constant-α approximation have to be judged against other cloud-controlling processes. For example, setting α to a constant value is unlikely to cause large discrepancies in simulated IWC when updraft speeds are high, because the time available for quenching of supersaturation is minimized.

Competing Nucleation
We explore in more detail the impact of INPs on homogeneously nucleated ICNCs. Figure 6 shows total ICNC (INP-derived plus homogeneously nucleated) as a function of updraft speed in the warm and cold cirrus regimes and for good and poor INPs with low and high number concentrations. The results are taken from a series of numerical simulations with prescribed updraft speeds. The ICNCs reflect the ice activation curves of the good and poor INP types shown in Figure 1. However, for updraft speeds <0.4-0.75 cm s −1 , simulated total ICNCs fall off more rapidly than the hyperbolic tangent ice activation curve would indicate (not shown in Figure 6). This is an artifact, because the time it takes to activate INPs then exceeds the overall simulation time. Regardless, ICNCs increase with updraft speed until a plateau is reached that corresponds to the prescribed total INP number concentrations, n ⋆ . The width of these regions (in terms of w) increases with n ⋆ and is somewhat smaller for poor INPs (larger s ⋆ ). Up to the end of the plateau regions, any further increase in w merely enhances the growth of already nucleated ice crystals preventing s to rise to values where homogeneous freezing sets in. For even higher w, homogeneous freezing sets in and ICNCs increase correspondingly to values above n ⋆ ; this happens earlier for the poor INPs. For the highest updraft speeds, the ICNC curves converge to the pure homogeneous freezing results, which are attained when the deposition sinks get weak enough (w ↓ ≪ w). This occurs at much higher updraft speeds (50-100 cm s −1 ) in the high INP number case compared to 5-10 cm s −1 in the low INP number case.
These results demonstrate that competition between INPs and solution droplets can take place across a wide range of updraft speeds, depending on n ⋆ , and even for poor INPs. Importantly, this w-range covers values of some 10 cm s −1 , making competing nucleation atmospherically relevant. While these results are qualitatively comparable to those obtained with the help of a simplified parametrization based on sharp ice nucleation thresholds ( Kärcher et al., 2006), they are based on a more realistic treatment of INP properties, where ice nucleation occurs over a range of supersaturation and are thus more accurate.

Simulations With Wave-Driven Vertical Wind Speeds
We discuss results from the full simulations with stochastic (mesoscale) temperature forcing. We begin by illustrating the temporal evolution of the activation process for one specific stochastic realization of vertical wind speeds and resulting fluctuations in T and s based on a mean updraft speed w m = 15 cm s −1 . We recall from Section 2.3 that the vertical wind speeds are sampled from an exponential distribution and alternate every 2.78 min. it takes almost 3 hr for s to decrease from about 0.15 to 0.05. This effect becomes even more pronounced in the cold cirrus regime due to slower growth rates. However, in situ data indicate that in most cirrus, most relative humidities measured within cirrus scatter around ice saturation (Krämer et al., 2009). The scatter might be explained by the presence of wave-driven variability for cirrus with low ICNCs (Kärcher et al., 2014). In our simulations, the long time it takes for s to finally reach zero is tied to the modeling of α by means of Equations 7 + 9. These equations imply that growth to the facets is inefficient for s < s crit ≈ 0.19 (0.36) for 220 (200) K according to Equation 9. Given the predominance of cirrus observed near equilibrium and that both the value for s crit and the physical nature of the underlying surface growth mechanism (embodied in the growth parameter μ) are poorly known or unconstrained in the case of cold cirrus, this issue warrants further analysis.
Modeling ice crystals as solid spheres may also be causing this issue, at least in part. For real ice crystals, deposition coefficients can vary substantially over different facets, potentially accelerating overall growth rates. In addition, complex polycrystals may be growing by dislocations for a substantial period of time (meaning μ = 1) or by stacking defaults (μ = 3). When keeping the parametric growth model (Section 2.4), reducing the values of the parameters μ and s crit would allow more substantial H 2 O uptake. More refined studies consider habit-predicting growth models accounting for habit mixes across a given ice crystal population constrained by observations.
We present results from ensemble simulations of competitive homogeneous freezing and heterogeneous ice nucleation for the low and high total INP number concentration cases in Figure 8. The frequency distributions of total nucleated ICNCs were obtained for a range of typical wave-driven updraft speeds added to a slow background uplift (1 mm s −1 ) to speed up the simulations. We run the warm cirrus case only because the rather high s crit -values predicted by Equation 9 for T < 220 K are only estimates (Section 2.4). All simulations are based on 1,000 realizations of stochastic vertical wind speed time series in each case. We do not show the distributions below 0.1 L −1 , as they are created only in the weakest updrafts and associated ICNCs are not captured due to runtime limitations. A prominent feature in the low INP number concentration case (green curves) is the very broad distribution of nucleated ICNCs. While the distributions show a peak at the corresponding ice-activated INP concentration (near 1 L −1 ), they extend to much larger values for all mean updraft speeds. The peak is absent in frequency distributions of ICNCs formed by homogeneous freezing . The broad distribution tails toward much higher ICNCs are caused by additional homogeneous freezing events in stronger updrafts. As the mean updraft speed increases, high values are more likely to occur in the exponential distribution of vertical wind speeds, leading to further broadening and a shift of the homogeneous nucleation mode toward greater ICNCs, while at the same time the frequency of occurrence of peak ICNCs due to INPs diminishes. Moreover, a shadowing effect occurs right above the total INP number that is closely related to the plateau region discussed above by means of Figure 6, best visible as a local distribution minimum in the simulations with 25 cm s −1 where updraft speeds are insufficient to overcome the deposition sink of the fully activated INP population.
Another salient feature is the cut-off near the prescribed high INP number concentration of 50 L −1 (orange curves). Homogeneous freezing does not occur in the 10 and 15 cm s −1 cases due to efficient supersaturation quenching, as indicated in Figure 7. However, some homogeneous freezing events are simulated for 25 cm s −1 owing to the sampling of updraft speeds that are large enough to overcome the INP deposition sink. The highest ICNCs are related to the highest updraft speeds, a phenomenon coined as preferential freezing .
After normalization of the frequency distributions, we analyze the mean ICNC and associated standard deviation for the most representative mean updraft speed (15 cm s −1 ): 16 ± 30 L −1 for n ⋆ = 1 L −1 and 4 ± 10 L −1 for n ⋆ = 50 L −1 (rounded values), which seems counter-intuitive at first glance. In the low INP case, mean ICNC is significantly higher than n ⋆ due to the frequently activated homogeneous freezing mode generating significant distribution variance. In the high INP case, mean ICNC is significantly lower than n ⋆ , and here the large variance is caused by variability in updraft speeds only along with the characteristic ice activation behavior shown in Figure 1. For comparison, we estimate the mean value for pure homogeneous freezing in the presence of waves: 313 L −1 , substantially exceeding the mean values including the INPs. Similar comparisons can be made for the other w m -cases. The solid magenta curves in Figure 8 represent the simulated ICNC statistics in the case of pure homogeneous freezing (n ⋆ = 0). The features discussed above are readily confirmed by comparing the curves with n ⋆ > 0 to these results and confirms that the highest n-values are correlated with the highest updraft speeds; we recall that homogeneously nucleated ice numbers increase strongly with increasing w. The dashed curves are analytical approximations of the frequency distributions for pure homogeneous freezing discussed in ; while the regions around the mean ice concentrations compare very well with the simulation results, the latter show poorer agreement toward the low concentration tails and at the highest ice concentrations. As the agreement worsens the higher the mean updraft speed, this is most likely a statistical abberration caused by undersampling in the simulations.
The strong reduction of homogeneously nucleated ICNCs for INP concentrations as low as 1 L −1 together with the fact that ICNC frequency distributions have been found to extend above values of 100 L −1 even in the TTL Jensen, Ueyama, Pfister, Bui, Lawson, et al., 2016), is strong evidence that competing nucleation is the rule rather than the exception in the formation of cirrus clouds. However, our analysis also hints at the difficulty of inferring nucleation mechanisms from airborne measurements.
We expect observed frequency distributions to be broader than simulated here toward their low tails and overall shifted toward lower ICNCs due mainly to ice crystal sedimentation, in particular for low mean updraft speeds. Sedimentation may also produce a local maximum at low concentrations unrelated to n ⋆ . With full lifecycles of cirrus (including sedimentation) included, the contrast between sets of simulated frequency distributions with a wide range of INP assumptions are rather subtle (Jensen et al., 2018). Taken together, we caution the reader against comparing the frequency distributions of nucleated concentrations, which predict the maximum ICNCs for a given cirrus cloud, to measurements capturing the full cirrus life cycles.

Summary and Outlook
We develop a parcel model enabling detailed analyses and enhancing understanding of the impact of INPs on the formation of pure ice clouds, including mesoscale gravity wave-driven vertical wind speeds and ice supersaturation-dependent deposition coefficients. In particular, we isolate and study factors controlling competitive homogeneous freezing and heterogeneous ice nucleation in cirrus. Heterogeneous ice nucleation is treated with time-independent ice-active fractions which are varied parametrically to cover a wide range of total number concentrations and nucleation behavior seen in the laboratory and field. Homogeneous freezing is parameterized based on an analytic treatment of the time-dependent ice formation process from supercooled aqueous aerosol droplets. Our methodology and findings pave the way for further improvements of cirrus parameterizations used in global models.
An important focus of our study is to obtain a lucid conceptual understanding of how ice crystal activation on INPs competes with homogeneous freezing of solution droplets in cirrus cloud formation. For this purpose, we define three different regimes. The first characterizes the temperatures in which cirrus form, affecting macroscopic cloud properties such as IWC and cloud optical depth. The second divides INPs into good and poor categories according to the central value and range of ice supersaturations around this value where most of them activate. Most importantly, we define the nucleation regime, wherein INP-derived ice crystals either allow or prevent homogeneous freezing, respectively. We introduce the quenching velocity as a key parameter separating the two cases.
Reinforcing earlier findings, cirrus formation is to be interpreted on the basis of the underlying dynamical forcing. Knowledge of updrafts speeds and their variability during ice formation is crucial to characterize the nucleation regime. The mean updraft speed determines the probability of occurrence of actual vertical wind speeds along Lagrangian air parcel trajectories that eventually leads to a frequency distribution of nucleated ICNCs from which respective mean values and variances can be inferred.
For the representative range of wave-driven mean updraft speeds and INP number concentrations considered here, homogeneous freezing events dominate distributions of nucleated ice crystal numbers toward large number concentrations, as seen in field observations. Among the most important finding of this study is the strong potential of INPs to lower homogeneously produced ice crystals numbers, and at times suppress homogeneous freezing events altogether, even for moderately high mean updraft speeds and rather poor ice activity. Besides the mean updraft speed, the next most important variable we need to know in simulations of cirrus formation is the total number concentration of potentially ice-active INPs, which we expect to be highly variable in space and time at cirrus levels.
Our study lets us dismiss the notion of either pure homogeneous freezing or heterogeneous ice nucleation during cirrus formation in many situations. The latter scenarios appear to be rather idealized and only expected to be realized in situations with high INP number concentrations, low mean updraft speeds, and low ambient temperatures (or combinations thereof). In situ measurements of frequency distributions of ICNCs, if representative for cirrus formation events, contain important signatures due to ice nucleation and growth. However, as measurements capture various stages of cirrus cloud development where a number of processes, including sedimentation, affect ICNCs, it is not straightforward to make inferences related to the impact of INPs.
To our knowledge, our study is the first to discuss effects of variable deposition coefficients in cirrus formation with multiple nucleation pathways. It is desirable to include such effects in models, as the use of constant deposition coefficients lacks basic physics and prevents habit prediction. In a first step, we consider facet growth by step nucleation based on spherical ice crystal geometry. While providing meaningful results for cirrus properties due to competing nucleation at moderate to high supersaturation (s = 0.3-0.6), our simulations suggest that supersaturation quenching times become very long once ice growth becomes inefficient at low supersaturation. This implies that cirrus may not stay close to ice equilibrium as frequently observed, preventing us from recommending the default use of variable α in cirrus simulations, especially for TTL cirrus.
It is currently unclear if the current parameterization of ice crystal growth from the vapor can be modified or if the growth theory itself needs to be revised. We have offered pointers to problems that require further study, including: the growth mode for complex polycrystals; the transition between various growth mechanisms while ice supersaturation increases or declines; and the fact that many of the observed, complex ice crystal features are influenced by lateral growth. More research is needed to advance predictive ice growth models for pure ice clouds in all temperature regimes.
We hope that our effort to bring together important developments in cirrus research that emerged in recent years-stochastic dynamical forcing of ice supersaturation, deterministic heterogeneous ice nucleation, and ice crystal growth from the vapor including surface kinetic effects-helps advance detailed numerical simulations and parameterizations of cirrus formation and development.

Data Availability Statement
The data that support the findings of this study are available on Zenodo: https://doi.org/10.5281/zenodo.5806109. The solver VODE is a collection of subroutines for the numerical solution of the initial-value problem for systems of first-order ordinary differential equations. The source files are publicly available under https://computing.llnl. gov/projects/odepack/software.