Aerosol-cloud interactions: The representation of heterogeneous ice activation in cloud models

. The homogeneous nucleation of ice in supercooled liquid water clouds is characterized by time-dependent freezing rates. By contrast, water phase transitions induced heterogeneously by ice nucleating particles (INPs) are described by time-independent ice-active fractions depending on ice supersaturation ( s ). Laboratory studies report ice-active particle number fractions (AFs) that are cumulative in s . Cloud models budget INP and ice crystal numbers to conserve total particle number during water phase transitions. Here, we show that ice formation from INPs with time-independent nucleation behavior is 5 overpredicted when models budget particle numbers and at the same time derive ice crystal numbers from s -cumulative AFs. This causes a bias towards heterogeneous ice formation in situations where INPs compete with homogeneous droplet freezing during cloud formation. We resolve this issue by introducing differential AFs, moving us one step closer to more robust simulations of aerosol-cloud interactions.


10
A wide variety of macromolecular or proteinaceous, crystalline, glassy, and solid aerosol particles act as INP in the atmosphere and participate in the formation of cirrus or in the glaciation of supercooled liquid water clouds . Among the various modes of heterogeneous ice formation, immersion freezing caused by INPs present within a volume of supercooled liquid water is considered the most relevant mode in mixed-phase clouds (Vali et al., 2015). Alternative freezing modes include contact freezing where ice forms upon collision of an INP with a cloud droplet, and condensation freezing where ice nucleates 15 while the cloud forms through cloud droplet activation. In conditions below liquid water saturation, deposition nucleation occurring in the absence of liquid water has traditionally been considered the most relevant heterogeneous ice formation mode (Vali et al., 2015). Yet, there is increasing evidence that the loci for ice nucleation on INP surfaces are pores in which water gathers below water saturation through capillary condensation (Marcolli, 2014;Kiselev et al., 2017;Holden et al., 2019).
Pore condensation and freezing (PCF) involves homogeneous ice nucleation within pores in cirrus conditions (air temperature 20 T < 233 K) and may occur heterogeneously through immersion freezing in mixed-phase clouds at higher temperatures (David et al., 2019;Marcolli, 2020).
In laboratory experiments, phase transitions from supercooled liquid water to ice are observed under controlled temperature and relative humidity conditions during set observational times for ice detection . In experiments employing droplet freezing techniques, ice nucleation is detected in arrays of droplets deposited on a substrate. Results are normalized 25 based on total droplet number, surface area, or volume to obtain freezing spectra that are usually reported in terms of cumulative ice-active fractions (Vali, 2019). Laboratory experiments using cloud or continuous flow chambers directly provide number fractions, φ, of ice-activated or frozen particles from a sample of size N 0 as a function of ice supersaturation, s. These fractions vary between 0 at s = 0 (ice saturation) and 1 at sufficiently large s and are cumulative, reflecting measurements in which the ice nucleation ability of a given sample is probed at successively increasing s-values. The total number of ice crystals 30 formed up to a value of s is then determined via N 0 φ(s). In the case of immersion freezing experiments, where an ensemble of water droplets with immersed INPs is cooled, frozen fractions are parameterized as a function of supercooling (temperature, T ) instead of supersaturation such that differential and cumulative AFs are functions of T instead of s.
During immersion freezing, ice nucleates over a wider temperature range compared to homogeneous freezing of pure water droplets (Peckhaus et al., 2016;Tarn et al., 2018). Heterogeneous freezing curves become even broader when a mixture of 35 different INP types is investigated. Yet, when one and the same droplet is repeatedly probed in freezing-thawing cycles during refreeze experiments, freezing occurs in a temperature range that is similarly narrow as the one for homogeneous freezing (Kaufmann et al., 2017).
While the broad range of freezing temperatures observed for an ensemble of droplets with immersed INPs can be ascribed to the deterministic (time-independent) component of freezing given by the characteristic freezing temperature of nucleation sites, 40 the much narrower spread of freezing temperatures observed in refreeze experiments evidences the stochastic (time-dependent) component on specific nucleation sites. Therefore, purely deterministic formulations correctly encompass the broad variability of nucleation sites evidenced in the freezing of particle/droplet ensembles, while neglecting the variability due to stochastic nucleation on specific sites evidenced in refreeze experiments. Applying a deterministic description of immersion freezing in cloud models is therefore justified, as the stochastic component just induces a minor modulation of the characteristic freezing 45 temperatures.
PCF is described by a deterministic parameterization as well, since ice formation in this mode is determined by the relative humidities required either for pore water condensation or ice growth with no stochastic component involved when temperatures are well below the threshold for homogeneous freezing of supercooled solution droplets, which is the case at cirrus temperatures (Marcolli, 2020;Marcolli et al., 2021). At warmer temperatures, PCF is basically immersion freezing in pores: both, the pore 50 filling and immersion freezing process are described deterministically. For the contact and condensation freezing modes, a deterministic description is also appropriate. Since the collision with INPs triggers the glaciation of cloud droplets during contact freezing, the time dependence of ice nucleation can be neglected. Similarly, the time dependence of condensation freezing is determined by the process of cloud droplet activation and ice nucleation can be considered immediate once the INP is immersed in water.

55
For these reasons, a formulation of AF as φ(s) without explicit time dependence is recommended for all modes of ice formation initiated by INPs. Figure 1. In a pool of N0 ice-nucleating particles, ∆N ice crystals form at ice supersaturation s and ∆N+ additional ice crystals form at a higher value, s+. The resulting ice crystals numbers can be directly predicted from time-independent, cumulative ice-active fractions, φ, based on the original INP sample of size N0 ('no budget' approach, black arrows labeled with φ). When already activated INPs are removed at s ('budget' approach, curved arrow), φ can no longer be used at s+ because of the reduced sample size, (N0 − ∆N ). This study derives differential ice-active fractions, ϕ, that can be applied to derive ∆N+ from the smaller sample (blue arrow).

Stating the issue
Treating ice formation as a deterministic process has implications for the use of s-cumulative AFs, φ, from laboratory experiments in cloud models.

60
The following issue arises as illustrated in Fig. 1: After ice has formed on INPs at a given value of s > 0, the latter are budgeted (removed) to ensure that the same INPs are no longer available for nucleation. With ∆N newly formed ice crystals, only (N 0 − ∆N ) INPs are available for further nucleation. The number of crystals formed in a succeeding nucleation event at s + > s must not be diagnosed from (N 0 − ∆N )φ(s + ), since the s-cumulative ice-active fraction is based on a sample of N 0 particles.
We may estimate the differential AF associated with the step process s → s + (Fig. 1). By definition, the number of INPs activating between s and s + is given by (N 0 − ∆N )ϕ(s + ). The number of unactivated INPs at s + is therefore given by the rate equation and the differential AF belonging to the blue arrow in Fig. 1 is expressed solely in terms of the cumulative AF: In the initial step of ice activation, where s increases from a value ≤ 0 to s > 0 for the first time, Eq. (4) simplifies to ϕ(s) = φ(s), because φ(s ≤ 0) = 0.
As we show in section 3, using cumulative instead of differential AFs in the 'budget' approach shifts the outcome of the competition between homogeneous droplet freezing and heterogeneous ice nucleation on INPs artificially towards the latter.

80
3 Solving the issue We derive differential ice-active fractions (section 3.1) and corresponding particle number budget equations (section 3.2) for phase transitions involving INPs with time-independent nucleation behavior and the ice crystals formed from them. The use of differential spectra derived from immersion freezing experiments is discussed by Vali (2019).
3.1 Differential ice-active fractions 85 We define a sequence of ice supersaturation values, {s j } (henceforth s-grid), with grid spacings ∆s j = s j − s j−1 and index j = 1, · · · , j max , starting at s 0 = 0 with φ(s 0 ) = 0. To derive differential AFs, it suffices to assume that s j -values increase.
We view φ j ≡ φ(s j ) as the statistical outcome of many identically prepared laboratory measurements. While φ j describes the fraction of INPs that are ice-active within the interval [0, s j ], the associated differential AF, ϕ j , shall describe only those INPs that are ice-active within [s j−1 , s j ]. Therefore, the probability that INPs remain unactivated at s j , (1 − φ j ), is given by 90 the product of the probabilities for particles not activating in all intervals ∆s prior to s j , (1 − ϕ ): from which ϕ j (j > 1) is obtained by recursion (by definition, ϕ 1 = φ 1 ): generalizing Eq. (4). Note that differential AFs depend on the type of s-grid. Equation (6) tells us that ϕ j equals the fraction of 95 INPs activated within ∆s j in the 'no budget' approach, ∆φ j = φ j − φ j−1 , corrected by a factor accounting for removing INPs that are ice-active below s j−1 .
Depleting INPs from their reservoir after ice activation as done in cloud models is equivalent to using smaller and smaller samples in laboratory experiments. As a result, the correct AFs to be used in such models, ϕ, are smaller than φ, because We model cumulative AFs analytically using: with the 50%-activation point, s (φ(s ) = 0.5), and the slope parameter, δs. Equation (7) 6) and (7), evaluated for a linear s-grid with various constant grid spacings, ∆s: s j = (j − 1)∆s. Consistent with Eq. (6), ϕ approaches φ for s s = 0.35 and ϕ is significantly lower than φ for s s . As ∆s increases and comprises a greater range of s-values, differences between cumulative and differential AFs diminish, but at the cost of resolution. respectively. The dashed curves in the 'budget' approach were obtained by wrongly using the cumulative AF so that the difference to the solid curves indicate the error in simulated particle numbers.

115
In this section, we employ both, a linear and sinusodial s-grid defined by s j = (j − 1)∆s , ∆s > 0 (8) respectively, where ∆s is a constant grid spacing. The linear grid describes a monotonically increasing supersaturation history representing a single ice formation event (s increases linearly due to adiabatic cooling for sufficiently small, constant cooling 120 rates). The wavy grid illustrates an idealized, non-monotonically rising supersaturation history with rising amplitude envelope (set by A j ), such as encountered during gravity wave activity with alternating cooling and heating cycles (controlled by α j ).
Both trajectories are shown in the left panels of Fig. 3.
To simplify the notation, j shall represent a dimensionless time variable. We note that, in general, each grid representation, {s j }, is subject to its own temporal development. For example, s might additionally be affected by latent heat release or ice 125 crystal growth. In cloud models, where grid spacing and temporal evolution cannot be separated, {s j } is determined by the time steps needed to simulate ice formation. The time steps may vary during the simulation depending on accuracy requirements.
The differential AFs from Eq. (6) are then computed based on a variable s-grid.
We denote the number of ice crystals forming from INP that are ice-active at s j as N i,j and the corresponding number of remaining (unactivated) INP as N a,j . We normalize both variables by the initial number of INP (at s = 0), N 0 : η i,j = N i,j /N 0 , 130 η a,j = N a,j /N 0 so that they are bounded by 0 and 1. The equations governing the evolution (j ≥ 1) in the time-independent (deterministic) nucleation framework for both linear and wavy supersaturation histories without budgeting particle numbers are given by: When considering particle number budgets (cf. section 2), we use differential AFs: with η a,0 = 1. The η a,j -values diminish as ice formation progresses while the total particle number, η a,j + η i,j , is conserved 145 (i.e., independent of j). For non-monotonically increasing supersaturation, we modify cumulative AFs byφ j = max{φ j , φ j−1 } to evaluate ϕ j . This ensures thatφ j = φ j−1 stays constant and ϕ j = 0 when s-values descend from (j − 1) to j, as motivated above.
Results for both types of s-grids are presented in Fig. 3, assuming either ∆s = 0.01 in Eq. (8) or variable grid spacing with j max = 101 in Eq. (9). For constant cooling (top panels), s-cumulated, normalized ice crystal numbers η i rise continuously 150 and INP numbers η a stay constant in the 'no budget' approach. When ice crystals are budgeted, using the cumulative AF overpredicts η i , although INPs are removed (η a decreases). When cooling and heating periods alternate (bottom panels), η i again increases at the expense of η a , but using cumulative AF together with budgeting ice crystals again leads to an overprediction of ice crystal numbers.
The impact on cloud properties of wrongly using cumulative AFs in specific simulations cannot be judged based on the 155 results shown in Fig. 3 alone. For example, in cirrus simulations, the change in total nucleated ice crystal numbers is likely small in situations with efficient INPs (with large dφ/ds near s ) and high cooling rates, as most INPs will activate straightaway and the time needed for s to increase above the 50% activation level is short.
A number of models exist to study clouds. Specific cloud processes such as nucleation are simulated in air parcel models on the 160 process level. Cloud-resolving models simulate formation and evolution of clouds with high resolution. Cloud system-resolving models track the life cycles of clouds on regional scales, better accounting for large-scale controls but with poorer resolution and increased need for parameterizations of small-scale processes. Global models with coarse resolution represent much of the atmospheric complexity, but represent clouds only by way of parameterization. In all types of models, ice-forming aerosol particles and cloud ice crystals may be represented by size-integrated properties, such as total particle number, or contain 165 explicit size information via particle size distributions (PSDs). We compare cumulative and differential AFs with size-resolved or size-integrated INP representation using the example of soot particles as INPs.
Soot particles nucleate ice after processing in mixed-phase clouds and aircraft contrails via PCF (Mahrt et al., 2020). Based on laboratory measurements, soot PCF predicts cumulative AFs of soot aggregates as a function of s and mobility diameter, Marcolli et al., 2021). We apply the soot PCF framework to soot particles emitted by aircraft jet engines. We model their with φ(s, D) taken from Marcolli et al. (2021).
175 Figure 4 shows size-resolved and size-integrated cumulative and differential AFs of aircraft soot particles processed in contrails . Size-resolved cumulative AFs decrease strongly with mobility diameter from 400 nm to 100 nm and reduce to zero for D < 40 nm (not shown). Since the soot PSD peaks in the Aitken size range, size-integrated ice activity is low (< 0.01) even at high ice supersaturation (s = 0.5).
A general recommendation on how to include differential AFs in models cannot be given, as this depends on details of the 180 numerical implementation of aerosol-cloud interactions, especially in global models where INP budgets are affected by both microphysics and transport. However, differential AFs are straightforward to implement in cloud models when making use of the budget Eqs. (12) and (13)  AFs that are cumulative in ice supersaturation. In prognostic cloud models, care must be taken to avoid the simulation of multiple ice formation events from the same particles. This is accomplished by the introduction of budget equations for INPs and the ice crystals deriving from them. While straightforward in the case of time-dependent budget equations (suitable for stochastic freezing) containing source and sink terms for the number of aerosol particles (or cloud droplets) and ice crystals, 190 a similar approach is not feasible in the case of INPs with singular ice nucleation behavior. We formulated differential AFs consistent with removal of such INPs after activation and introduced modifications that are necessary when ice supersaturation evolves non-monotonically over time. We discussed the representation of ice activation in cloud models and showed that using differential AFs prevents overestimating INP effects. Finally, we demonstrated the importance of including INP size information in estimations of AF.

195
Our insights help improve cloud simulations and better understand the relative roles of natural and anthropogenic INPs in determining coverage, lifetime, and radiative response of mid-and high-level clouds.

Appendix A: Analytical representation of ice-active fractions
The hyperbolic tangent chosen to represent φ allows to easily compare measured and parameterized s-cumulative ice-active fractions and perform sensitivity studies. It is based on only two parameters, s and δs, with clear physical significance (sec-200 tion 3.1). Here, we apply this function to fit the activation curve from Ullrich et al. (2017) for monodisperse (1 µm) spherical dust particles using s = 0.352 and δs = 0.0175. Figure A1. Comparison of an ice activation parameterization for 1 µm-dust particles (black curve) with an analytical approximation (blue). Figure A1 shows that Eq. (7) approximates the parameterization very well, especially in the crucial part around s , where φ rises steeply from low to significant values. We presume that the hyperbolic tangent provides reasonable fits to activation curves of other INP types as well, which show a similar s-dependence. We note that the parameter values suitable for monodisperse 205 particles change when size-dependent s-cumulative ice-active fractions are integrated over a particle size distribution.
In discussing deterministic ice formation (section 3), we have chosen a larger slope parameter, δs = 0.05. This spreads ice activation over a larger range of s-values (compare black curve from Fig. 2 with the blue curve from Fig. A1) and helps illustrate the ice activation events shown in Fig. 3 more clearly.
Author contributions. B. K. conceived of the study, performed the calculations, and wrote the manuscript draft. C. M. developed the concept of differential ice-active fractions, discussed ice nucleation measurements, and contributed to the final writing.