Prognostic simulations of mixed-phase clouds with model AC-1D v1.0: the impact of aerosol types and freezing parameterizations on ice crystal budgets

Mixed-phase clouds at high latitudes contribute to the uncertainty in predicting cloud feedbacks and climate sensitivity, mainly due to the complexity of microphysical processes that influence the partitioning between the supercooled liquid and ice phases, and hence, cloud radiative effects on regional scales. Particularly in Arctic mixed-phase clouds, the activation of ice-nucleating particles (INPs) from various aerosol populations remains a leading source of uncertainty. We developed an aerosol-cloud one-dimensional (AC-1D) model, which provides a novel framework to prognostically treat INP and ice crystal budgets while explicitly accounting for polydisperse and multicomponent aerosol that activate INPs following different freezing parameterizations. The AC-1D model is informed by large-eddy simulations to probe the impact of INP representation on predicted ice crystal number concentrations (N_i) and ice crystal budgets in mixed-phase Arctic stratus. We apply three immersion freezing (IMF) parameterizations, two time-independent (singular) and one time-dependent (classical nucleation theory), to predict the evolution of the INP reservoir and resulting ice crystal budget from polydisperse mineral dust, organic (humic-like substances), and sea spray aerosol particle size distributions. Our analysis focuses on how variations in aerosol number concentration and cloud system parameters such as cloud cooling rate, cloud-top entrainment rate, and ice crystal fall speed influence the INP reservoir and ice crystal budgets. Furthermore, this study investigates the competitive ice nucleation dynamics in mixed aerosol environments and provides a process-level quantification of the INP budget terms, which directly controls ice crystal budgets. For all studied case scenarios, the aerosol types and associated particle size distributions significantly impact INP and N_i, and the choice between a time-dependent and a singular freezing description yields orders-of-magnitude differences in the predicted INP and N_i over the 10 h simulation time, reflecting typical cloud lifetimes. Our results show that the influence of cloud cooling, INP entrainment, and sedimentation varies significantly depending on the chosen freezing parameterization. These findings underscore the critical need for robust IMF parameterizations and precise cloud system observations to enhance the accuracy of models in predicting mixed-phase cloud structure and evolution.

1 Introduction

In the Arctic regions, the rate of warming has been at least twice as large as the global mean since pre-industrial times (Holland and Bitz, 2003; Serreze and Barry, 2011; Gulev et al., 2021), a phenomenon known as “Arctic amplification” (Hahn et al., 2021; Morice et al., 2021) that coincides with a drastic decrease in the cover of sea ice (Stroeve et al., 2012; Richter-Menge et al., 2018). Stratus, a predominant cloud type in polar regions (Shupe et al., 2008; Shupe, 2011; Andronache, 2018; Lubin et al., 2020), plays a crucial role in the surface and top-of-atmosphere radiative budget due to its extensive prevalence (e.g., Dong and Mace, 2003; Zuidema et al., 2005). These clouds control the development of precipitation through their liquid and ice phases (e.g., Field and Heymsfield, 2015; Mülmenstädt et al., 2015; Korolev et al., 2017; Silber et al., 2021). Arctic stratus clouds are predominantly mixed-phase clouds (MPCs) (Curry et al., 2000; Korolev et al., 2003; Shupe, 2011), characterized by the presence of at least one supercooled liquid water layer where ice crystals form and sediment subsequently (Shupe et al., 2006; Morrison, 2012; Silber et al., 2021). The persistence of this mixed-phase state is governed by a delicate interplay between cloud dynamics and microphysics, with the ambient ice concentration determining whether the cloud remains liquid-bearing or glaciates (Khain et al., 2022).

In these Arctic MPCs, ice formation typically occurs from the supercooled liquid phase through heterogeneous ice nucleation. Heterogeneous ice nucleation can commence via several modes including immersion freezing (IMF), where freezing is initiated on the surface of an ice-nucleating particle (INP) that is immersed in a supercooled aqueous solution droplet (Vali et al., 2015). Recent studies have indicated that IMF is the dominant pathway of primary ice production (PIP) in MPCs (e.g., Ansmann et al., 2008; Prenni et al., 2009; Hoose et al., 2010; de Boer et al., 2011; Hoose and Möhler, 2012; Murray et al., 2012; Westbrook and Illingworth, 2013; Kanji et al., 2017; Silber et al., 2021; Burrows et al., 2022).

To accurately predict the number concentration of ice nucleating particles (INPs) from aerosol particles, IMF parameterizations are needed. Over the past decades, IMF parameterizations representing specific particle types have been developed (e.g., Bauer et al., 2008; Murray et al., 2012; Niemand et al., 2012; Knopf and Alpert, 2013; DeMott et al., 2015; Liu et al., 2016; China et al., 2017; Kanji et al., 2017; Knopf et al., 2018; McCluskey et al., 2018; Penner et al., 2018; Alpert et al., 2022; Burrows et al., 2022; Knopf and Alpert, 2023). Typically applied freezing parameterizations include singular schemes based on the singular hypothesis (Levine, 1950; Vali, 1971; Vali, 2014), and stochastic schemes, such as those based on classic nucleation theory (CNT) (Pruppacher and Klett, 2010; Knopf and Alpert, 2023).

Singular (INAS-type) IMF parameterizations represent freezing as a cumulative function of temperature, assuming that a population of nucleation-active sites initiates freezing once sufficiently low temperatures are reached on an immersed particle surface. In this framework, the frozen fraction depends on temperature and aerosol properties (including particle number concentration and available surface area), and freezing is often treated as effectively instantaneous upon reaching the activation threshold. As a result, such schemes do not explicitly represent the observed dependence of freezing on cooling rate or time spent at a given temperature (Bigg, 1953; Pruppacher and Klett, 2010; Alpert and Knopf, 2016; Knopf et al., 2020; Arabas et al., 2025).

CNT-based parameterizations, in contrast, represent immersion freezing as a rate-based process governed by aerosol-specific kinetic parameters (e.g., contact angle distributions or water-activity-dependent coefficients), rather than fixed active sites. In these formulations, the freezing rate scales with the available INP surface area, such that increasing INP surface area increases the rate of ice formation (Pruppacher and Klett, 2010; Knopf et al., 2020; Knopf and Alpert, 2023). Consequently, for constant thermodynamic conditions, the cumulative number of freezing events increases with elapsed time, consistent with laboratory observations (e.g., Biermann et al., 1996; Koop et al., 1997; Alpert and Knopf, 2016; Knopf et al., 2020; Deck et al., 2022).

Although IMF parameterizations are often categorized as deterministic/singular versus stochastic, it is useful to distinguish the microscopic interpretation of nucleation from the numerical implementation used in atmospheric models. At the molecular scale, heterogeneous nucleation is frequently described as a random process; however, most Eulerian models predict the mean evolution of ice formation using deterministic population equations. In this work, the CNT-based ABIFM scheme is implemented as a rate-based (memoryless) freezing formulation, in which the frozen fraction over a model time step depends only on the current thermodynamic conditions and aerosol surface area, not on the prior residence time of individual particles. This distinction helps avoid conflating physical randomness with Monte Carlo sampling approaches used in particle-resolved models (Shima et al., 2020; Arabas et al., 2025), and clarifies how time dependence enters the governing equations across different model hierarchies.

Mineral dust, organic, and sea spray aerosol (SSA) particles are common aerosol types in the Arctic region (Udisti et al., 2020; Schmale et al., 2021). At temperatures below −20 °C, mineral dust has a particularly high ice nucleating efficiency (Kanji et al., 2017) and is assumed to be the best understood source of INPs (Burrows et al., 2022). The variability of atmospheric INP concentrations can, in cases, be attributed to the change in number concentration of long-range transported dust (DeMott et al., 2003; Chou et al., 2011; Boose et al., 2016; Burrows et al., 2022; Shi et al., 2022; Kawai et al., 2023).

Organic aerosol (OA) particles can be directly emitted from fossil fuel combustion and biomass burning, also termed primary organic aerosol (POA). The condensation of oxidized volatile organic compounds (VOCs) can yield secondary organic aerosol (SOA) particles (Hallquist et al., 2009; Shrivastava et al., 2017; Bergman et al., 2022; Srivastava et al., 2022). The organic matter (OM) associated with OA particles that initiate the ice nucleation in the atmosphere is still not well understood (Knopf et al., 2018). OM can exhibit various phase states at the same atmospheric thermodynamic conditions, resulting in different ice nucleation pathways and ice nucleation rates (Knopf et al., 2018). Through recent field measurements and laboratory investigation, it has been pointed out that particulate OM can contribute to atmospheric ice nucleation (e.g., Schnell and Vali, 1975; Knopf et al., 2010, 2014, 2023a; Wang et al., 2012a, b; Hiranuma et al., 2013; Xue et al., 2024).

Several studies have demonstrated that in remote marine regions far from the influence of continental INP sources, SSA particles can act as a significant INP source (e.g., Bigg, 1953; Burrows et al., 2013; Wilson et al., 2015; Vergara-Temprado et al., 2017; Huang et al., 2018; McCluskey et al., 2018; McCluskey et al., 2019; Zhao et al., 2021; Alpert et al., 2022; Raatikainen et al., 2022; Xue et al., 2024).

Soot particles have been widely discussed as a potential source of anthropogenic INPs. Under polluted conditions, relatively high black carbon fraction has been found in both cloud residuals as well as simulated results (Savre and Ekman, 2015). However, recent laboratory studies reveal that soot exhibits minimal efficacy in immersion-mode freezing (Friedman et al., 2011; Schill et al., 2018; Kanji et al., 2020; Schill et al., 2020). Its contribution to atmospheric ice formation at temperatures above −38 °C is likely not significant (Kanji et al., 2020).

In addition to ambient particle types, typical aerosol particle size distributions (PSDs) present in the Arctic are relevant to ice formation rates. The PSDs determine the particle number and surface area concentration, a necessary input for the IMF parameterizations, thus influencing the number concentration of INPs and, ultimately, ice crystals. In this study, we apply the PSDs derived from Indirect and Semi-Direct Aerosol Campaign (ISDAC) (Earle et al., 2011; Hiranuma et al., 2013) and the International Chemistry Experiment in the Arctic LOwer Troposphere (ICEALOT) (Quinn et al., 2017), both carried out in the Arctic region. Detailed information can be found in the methods section.

Our study is motivated by the fact that the inferred strength of ice production within Arctic stratus MPCs has often been challenging to reproduce by models applying a prognostic INP treatment, i.e., a treatment where the budget of INP is tracked (Harrington and Olsson, 2001; Morrison et al., 2005; Fridlind et al., 2007, 2012; Fan et al., 2009; Avramov et al., 2011). Such a treatment contrasts with simpler diagnostic treatments of INP that enable physically unrealistically high ice formation rates to persist indefinitely (e.g., Knopf et al., 2023b), and our previous work indicates that the type of IMF parameterization (singular or CNT-based) plays an important role in the strength of ice formation rates under Arctic stratus conditions (Knopf et al., 2023b; Arabas et al., 2025). In general, when assuming a singular freezing scheme, INP depletion from a turbulently mixed stratus layer occurs within a short period of time, on the order of the mixing time for a well-mixed cloud-topped boundary layer. This rapid depletion results from the INP exposure (via turbulent transport) to conditions that result in their instantaneous activation at the coldest conditions within the cloud layer (within cloud tops). Subsequently, cloud-top entrainment or cooling as drivers of INP activation have often appeared too weak to explain observed ice crystal loading (Fridlind et al., 2012; Westbrook and Illingworth, 2013). Other processes may contribute to sustaining INP under some specific conditions, such as entrainment of INP from below the turbulently mixed layer if the cloud-containing layer is decoupled from the surface (Avramov et al., 2011) or restoration of some of the activated INP via complete ice crystal sublimation within a turbulently mixed layer (i.e., INP recycling) if there is an ice-subsaturated layer (at the base) that is deep enough to enable that (Solomon et al., 2015). Under the simplest conditions without such processes possible, Fridlind et al. (2012) demonstrated how unrealistically high concentrations of INPs from the free troposphere were needed to sustain inferred ice production via cloud-top entrainment when applying a singular IMF parameterization. Unfortunately, direct constraints on INP budgets have not been possible from available data sets owing at least in part to a lack of INP data within cloud layers, and inferring INP budgets from ice crystal number concentrations invokes an additional set of large uncertainties emerging from all of the factors that influence ice crystal PSDs in Arctic MPCs that must also be parameterized (Fridlind et al., 2007, 2012; Avramov et al., 2011; Ovchinnikov et al., 2014; Morrison et al., 2020).

In most previous modelling studies, prognostic ice nucleation has been parameterized with singular IMF schemes (Morrison et al., 2005; Avramov et al., 2011; Fridlind et al., 2012; Solomon et al., 2015; Tully et al., 2023). It has been pointed out that ice nucleation parameterizations derived from laboratory experiments may not be applicable to all atmospherically relevant conditions (temperature and humidity) (e.g., Niemand et al., 2012; Hiranuma et al., 2014; Ullrich et al., 2017), partly due to instrument limitations and the limited amount of data collected (Burrows et al., 2022). Consequently, parameterizations may fail to capture rate-based freezing, which is essential for sustaining continuous ice crystal production (Westbrook and Illingworth, 2013; Yang et al., 2013). Several studies have implemented CNT-based parameterizations (Savre and Ekman, 2015; Raatikainen et al., 2022; Shi et al., 2022) to investigate the effects of different INPs on Arctic MPCs. Savre and Ekman (2015) found that the application of an evolving α-PDF scheme (introducing a contact-angle distribution function that scales the freezing efficacy of the INPs) can support continuous in-cloud ice production controlled mostly by the competition between cloud cooling, cloud-top entrainment, and ice sedimentation in three simulated cases. Raatikainen et al. (2022) showed that the relative significance of marine INP emissions and accounting for INP recycling was crucial for maintaining MPCs in their simulations over water rather than ice surfaces. Knopf et al. (2023b) pointed out another possible explanation that the selection of IMF parameterizations determines the size of the INP reservoir in Arctic stratus, with the CNT-based parameterizations producing a reservoir that is several orders of magnitude larger than singular parameterizations, thereby becoming the dominant factor for sustained strength of ice crystal formation.

While previous studies have implemented different ice nucleation parameterizations in global climate models (e.g., Liu and Penner, 2005; Hoose et al., 2010; Wang et al., 2014), these efforts generally focused on improving agreement with global observations by tuning parameters or refining aerosol composition. Few studies have isolated the structural uncertainty of the parameterization itself (singular vs. CNT-based) within a controlled, prognostic framework for the same cloud scenario, despite recent findings showing that these choices can lead to orders of magnitude differences in ice production (Knopf et al., 2023b; Arabas et al., 2025). Consequently, an “apples-to-apples” comparison of how the fundamental choice of freezing framework dictates the INP reservoir dynamics and subsequent ice crystal budget remains necessary.

Building on these identified uncertainties in primary ice production (PIP), we now focus in detail on the broader characteristics of Arctic aerosol. In this study, we employ the AC-1D model as a prognostic tool designed to isolate the structural uncertainties in PIP. The impact of liquid-ice feedbacks (such as the Wegener-Bergeron-Findeisen process driving ice growth) and secondary ice production (SIP) on the ice crystal budget depends crucially on the accurate description of the PIP. A key feature of this framework is its ability to conduct comprehensive sensitivity analyses by coupling polydisperse and multicomponent aerosol inputs directly to the INP and ice budgets. This setup allows for the simultaneous, prognostic evaluation of fundamentally distinct IMF parameterizations, while permitting user-defined adjustments to thermodynamic profiles and cloud system parameters. By generating detailed process-level data, such as explicit INP and ice crystal budgets, the model serves as a robust testbed to determine how the choice of parameterization dictates the PIP and evolution of the INP reservoir in Arctic mixed-phase clouds (Knopf et al., 2023b; Arabas et al., 2025). A key simplification of the AC-1D model is that the liquid phase is taken as a fixed quantity. This approach decouples ice microphysics from liquid-phase feedbacks, enabling a direct attribution of changes in the ice crystal budget solely to the immersion freezing parameterization under investigation (Knopf et al., 2023b). This idealization is physically justified by observations of long-lived Arctic MPCs, which are often characterized by a quasi-steady liquid water path. This state is maintained by very weak ice precipitation rates that only minimally desiccate the cloud, a feature identified in both specific cases (e.g., Fridlind et al., 2012) and broader statistical analyses (e.g., Silber et al., 2021). The objective of AC-1D is to predict the co-evolution of the size- and composition-distributed aerosol acting as INP and the activated N_i profiles in a simplified framework as a function of specified values of the boundary layer turbulent mixing time scale, cloud-top entrainment rate, cloud cooling rate (CCR), and number-weighted ice crystal fall speed (cf. Fridlind et al., 2012). Knopf et al. (2023b) demonstrate how these specified values exhibit markedly differing controls on ice formation depending on which type of ice nucleation scheme is implemented. Different IMF parameterizations derived from the same freezing experiments will scale differently when predicting INPs for conditions other than those of the experiment (Knopf et al., 2021).

Here we extend the Knopf et al. (2023b) study by investigating the impact of different aerosol types and associated PSDs on INPs, IMF parameterizations, and various microphysical and model parameters on the evolution of INP and N_i. To account for the diversity of IMF parameterizations, including the ice nucleation number (INN)-based parameterization (DeMott et al., 2010, 2015), INAS parameterization (Niemand et al., 2012; China et al., 2017) and water-activity-based IMF model (ABIFM) parameterization derived from CNT (Knopf and Alpert, 2013; Alpert and Knopf, 2016). We further initialize the model with three different aerosol particle types and respective PSDs, including mineral dust, organic (humic-like substances) aerosol particles, and SSA particles, guided by observations in the Arctic regions. Each of these parameterizations requires different sets of parameter inputs, which will be discussed in detail in the methods section. The effect of varying aerosol number concentration and specified cloud system parameters on the INP reservoir and N_i are also assessed. Furthermore, to assess the model's capability in handling more complex atmospheric compositions and to explore the competitive interplay of different INP types, MPC scenarios incorporating co-existing externally-mixed mineral dust, organic, and sea spray aerosol populations are also investigated. The paper is organized as follows. In Sect. 2 we describe our implementation of the AC-1D model and all simulation cases in detail. In Sect. 3 we present the sensitivity simulation results, also the findings from combined aerosol scenarios and a detailed INP budget analysis, followed by a discussion in Sect. 4. Finally, the key findings are summarized and discussed in Sect. 5. Appendix A provides a detailed description of the IMF parameterizations. Appendix B details the INP array calculation in our model for singular approaches. The governing equations are presented in Appendix C. For a comprehensive list of symbols and abbreviations, refer to Appendix D and E, respectively.

2 Methods

The AC-1D model is currently based on cloud system conditions observed in a well-mixed cloud-topped boundary layer over sea ice during the Surface Heat Budget of the Arctic campaign (SHEBA) (Knopf et al., 2023b). In this section, we briefly describe the SHEBA case study. Then we introduce the applied aerosol types and corresponding PSDs and the governing equations of the AC-1D model. Lastly, we summarize the different simulations considered in this study.

2.1 The SHEBA case study

As the initial setup in our AC-1D model, we use an Arctic MPC case study that was well observed during the SHEBA campaign (Curry et al., 2000), and that has been widely used for advancing our understanding of MPC conditions (Morrison et al., 2011; Fridlind et al., 2012). The campaign was carried out at roughly 76° N, 165° W in the Beaufort Sea, Arctic Ocean, approximately 570 km northeast of Prudhoe Bay, Alaska. This case study was characterized by widespread, long-lived, shallow MPCs coupled to an underlying pack ice with cloud tops of approximately 600 to 400 m as indicated by Millimeter Wavelength Cloud Radar (MMCR) echoes (Morrison et al., 2011; Fridlind et al., 2012). Rawinsondes provided profiles of temperature and relative humidity (RH), indicating a well-mixed liquid cloud-topped boundary layer. Two-dimensional cloud optical array probes (2D-C), two-dimensional precipitation optical array probes (2D-P) and the Cloud Particle Imager (CPI) yielded reconstructed cloud particle shapes and sizes (Fridlind et al., 2012). The specific case study from the SHEBA campaign used to define our model's thermodynamic profile is based on observations of a long-lived mixed-phase cloud on 7 May 1998 (Uttal et al., 2002; Fridlind et al., 2012; Fridlind and Ackerman, 2018). This late-spring seasonal context aligns well with the timing of the ISDAC (April 2008) and ICEALOT (March–April 2008) campaigns from which the aerosol data are derived.

Compared to conditions observed during the Indirect and Semi-Direct Aerosol Campaign (McFarquhar et al., 2011) or the Mixed-Phase Arctic Cloud Experiment (Verlinde et al., 2007), this SHEBA case provides a simple starting point owing to ice supersaturated conditions from liquid cloud top to the near-surface as shown in Fig. 1, low liquid water path (approximately 5–10 g m⁻², Fridlind et al., 2012), high droplet number concentration around 200 cm⁻³, an absence of liquid-phase precipitation or active ice aggregation, and relatively sparse concentrations of unrimed and nondendritic ice crystals (Fridlind and Ackerman, 2018). These conditions are not expected to support SIP, which is commonly assumed to proceed by rime splintering (Hallett and Mossop, 1974) or freezing fragmentation (Lauber et al., 2018; Keinert et al., 2020), but is currently debated for its relative role in ice production process in MPCs (Phillips et al., 2017; Luke et al., 2021; Waman et al., 2023; Zhao et al., 2023). Robustly observing SIP in both field and laboratory studies is challenging (Korolev et al., 2020; Silber, 2023), and large uncertainty remains in simulating SIP processes (Phillips et al., 2018; Korolev and Leisner, 2020; Miltenberger and Field, 2021; Burrows et al., 2022; Waman et al., 2023). The apparent occurrence of SIP may be less than 10 % in slightly supercooled (warmer than −10 °C) Arctic clouds (Luke et al., 2021) and is expected to be negligible under conditions for this SHEBA case (Fridlind et al., 2012) in contrast to other cases (Fridlind and Ackerman, 2018). Our setup is therefore chosen to isolate the role of PIP via IMF, which precedes any subsequent SIP in non-seeded Arctic stratus.

Figure 1 Thermodynamic conditions applied in the minimalistic 1D aerosol-cloud model. From left to right: The temperature (T), the relative humidity (RH), supersaturation with respect to the ice (S_ice) and Δa_w. Blue and green lines represent the initial (t = 0 h) thermodynamic conditions and the thermodynamic conditions with cloud cooling rate of 0.3 °C h⁻¹ (CCR = 0.3 °C h⁻¹ sensitivity run) after 10 h, respectively. The blue shaded area denotes the cloud layer.

2.2 Aerosol particle size distributions

To evaluate the impact of different aerosol types on the INP reservoir, N_i and ice crystal formation rate, we examine the effects of mineral dust, organic, and SSA particles. While the thermodynamic profile is derived from the specific SHEBA case study (early May) to represent a typical Arctic mixed-phase boundary layer, the aerosol initialization represents a composite of spring Arctic conditions derived from ISDAC and ICEALOT (April).

Since size-resolved aerosol measurements were not available for the specific SHEBA case study, we approximated the total aerosol inputs based on the observed cloud droplet number concentration (N_d ≈ 200 cm⁻³; Fridlind et al., 2012). Our model assumes that all in-cloud aerosols are activated.

To reconstruct a physically consistent aerosol population, we utilized data from the ISDAC campaign (for mineral dust and organic fractions) and the ICEALOT campaign (for SSA). The ISDAC campaign is an appropriate proxy because it exhibited cloud microphysical properties very similar to SHEBA case, with measured N_d values ranging from 185 to 205 cm⁻³ (Savre and Ekman, 2015).

To achieve quantitative alignment between the prescribed aerosol load and the SHEBA target (N_d ≈ 200 cm⁻³), we selected aerosol inputs corresponding to the “clean case” defined in Earle et al. (2011). This case shows number concentrations of aerosol (N_aer) less than 250 cm⁻³, while N_d is approximately 135 ± 34 cm⁻³ (Earle et al., 2011). Critically, this case exhibits a high activation efficiency, validating our model assumption that all in-cloud aerosols are activated. Within this total number constraint, the relative compositional fractions of mineral dust and organic aerosol were derived from single particle analyses performed during the same field campaign (Hiranuma et al., 2013). Although the PSDs are derived from the subsequent flight (Flight 31, 27 April), Earle et al. (2011) classify both dates as the same meteorological regime, justifying this composite initialization.

While mineral dust is a major component of Arctic INPs during long-range transport events (Creamean et al., 2022; Böö et al., 2023), the background conditions defined for this specific sensitivity study maintain low dust number concentration compared to organic and SSA particles to match the microphysical constraints of the SHEBA case. High-load scenarios, e.g., increase in dust load, are explored separately in sensitivity runs. Similarly, while the SHEBA case occurred over pack ice where local SSA emission is suppressed, we include the ICEALOT SSA distribution (Quinn et al., 2017) to account for potential transport from open water leads and to establish a robust background state for testing competitive nucleation processes.

For each aerosol type, the applied aerosol particle size distributions (PSDs) are polydisperse consisting of two or three lognormal modes (Table 1). In addition to Aitken and accumulation modes, this framework includes a larger accumulation mode for aged aerosols and a source-specific SSA mode.

Table 1 Summary of the applied particle size distribution parameters for different aerosol types.

Figure 2 Aerosol particle size distributions for mineral dust (a), organic particles (b), sea spray aerosol (c) and composite particle size distribution (d) following PSD parameters given in Table 1. Different colors represent different modes for respective aerosol PSDs.

The mineral dust and organic PSDs are based on aerosol samples measured by aircraft for single-layer stratocumulus with below-cloud aerosol concentrations less than 250 cm⁻³ (clean cases) as reported in Earle et al. (2011). We use the same lognormal distribution parameters (geometric mean diameter (D), geometric standard deviation (σ)) for mineral dust and organic particles as provided by Savre and Ekman (2015). The aerosol number concentrations for each mode of mineral dust and organic aerosol particles are derived from micro-spectroscopic single particle analysis of ambient particles that were collected by aircraft (Hiranuma et al., 2013). We average particle types from Flight 30-Substrate 6 (F30-S6) during ISDAC to derive respective particle number concentrations. The applied SSA PSD is based on measurements during the International Chemistry Experiment in the Arctic Lower Troposphere (ICEALOT) (Quinn et al., 2017). SSA particles were sampled 18 m above the sea surface. Lastly, to reflect a more realistic aerosol population we combine the mineral dust, organic, and SSA PSD (composite PSD). Figure 2 displays the lognormal PSDs of the different aerosol particle types and the composite PSD, derived from the modal parameters specified in Table 1.

2.3 AC-1D model setup

LES baseline results for the SHEBA case study (Fridlind et al., 2012) using the Distributed Hydrodynamic Aerosol and Radiative Modeling Application (DHARMA) code (Ackerman et al., 2000; Stevens et al., 2002) serve as the source of the cloud layer conditions for our model setup. As noted above, AC-1D maintains a time-invariant thermodynamic profile, including temperature, relative humidity, and liquid water path. We set the domain height to 390 m with the mixed phase cloud layer located from 262.5 to 390 m determined by thermodynamic conditions obtained from SHEBA LES simulation results (see Fig. 1). The temperature decreases roughly linearly with height from −16.2 to −19.8 °C. In this cloud case, the ice supersaturation ratio is larger than 1 throughout nearly the entire domain (see Fig. 1), so ice sublimation of sedimenting ice crystals is negligible. The governing model baseline parameter values of mixing time scale, entrainment rate, number-weighted ice fall speed, simulation time step, and simulation vertical resolution applied in the model are shown in Table 2. In this setup, PIP and ice crystal sedimentation are the dominant ice microphysical processes tracked. Processes such as the Wegener-Bergeron-Findeisen (WBF) process, riming, and SIP are deliberately excluded to isolate the uncertainty in INP and ice crystal budgets arising directly from the choice of the immersion freezing parameterization. Regarding aerosol physics, the model treats the aerosol population prognostically, explicitly accounting for the depletion of INPs via activation and changes due to transport. However, internal aerosol microphysical processes (e.g., coagulation, condensation growth) are not simulated. IMF parameterizations for respective aerosols are summarized in Table 3 and a detailed description of these three IMF parameterizations is given in Appendix A. To ensure a robust comparison of structural uncertainty, we selected parameterization pairs derived from identical underlying datasets where possible. For mineral dust, the INN (DeMott et al., 2015), the INAS (Niemand et al., 2012) and ABIFM (Alpert and Knopf, 2016) schemes were all derived from the same laboratory experiments conducted at the AIDA cloud chamber. Because these schemes share a common experimental origin, the profound differences we observe in the AC-1D model results isolate the structural uncertainty inherent in how each scheme formulates the freezing process (e.g., time-dependence) rather than differences in the dust samples. Similarly, for organic aerosol, both the INAS and ABIFM parameterizations are derived from the same laboratory dataset of Leonardite (standard humic acid) as described in Knopf and Alpert (2013) and China et al. (2017). For SSA, both schemes are derived from the same mesocosm dataset (Alpert et al., 2022). In contrast, the INN parameterization (DeMott et al., 2010) is included as a benchmark representing averaged ambient measurements. Therefore, comparisons involving INN parameterization are not based on the same datasets, while comparisons between INAS and CNT remain structurally constrained by their common datasets.

Table 2 1D aerosol-cloud model baseline parameters.

^* The sedimentation rate represents a number-weighted average derived from the LES bin-microphysics, corresponding approximately to unrimed plate-like crystals with diameters of around 400–600 µm.

Table 3 Immersion freezing parameterizations used in this study. The table lists the source references, applicable aerosol types, and the valid temperature range derived from the underlying experimental datasets. Input parameters are aerosol number concentration for particles larger than 0.5 µm, N_{aer > 0.5 µm}, aerosol number concentration, N_aer, aerosol surface area, A_aer, temperature, T, and relative humidity, RH. INP array stores the number concentration of activatable INPs using different data structures for different immersion freezing parameterizations. The INP array dimensions include z_i, height at grid cell index i, t, time, d, INP diameter array, T, INP temperature array.

Regarding the temperature validity of these schemes, our simulation conditions (−22.8 to −16.3 °C) fall strictly within the experimentally validated bounds for the Mineral Dust INAS/ABIFM schemes (Niemand et al., 2012; Knopf and Alpert, 2013), all SSA schemes (Alpert et al., 2022), and the ambient INN scheme (DeMott et al., 2010). For the Organic INAS and ABIFM schemes, the underlying laboratory data for Leonardite extends up to ∼ 257 K (−16 °C) (Knopf and Alpert, 2013), effectively covering our simulation range. For the Mineral Dust INN scheme (DeMott et al., 2015), our simulation extends the temperature by approximately 4 K warmer than the original dataset. We consider this minor extrapolation justified as the parameterization is constrained to decay toward zero activity as temperature increases.

Mineral dust is a broad category containing significant mineralogical diversity. As discussed in previous studies (e.g., Murray et al., 2011; Kanji et al., 2017), ambient mineral dust reflects a wide variation of mineral types with different freezing propensities. While this study treats dust as a single species consistent with the bulk dust parameterizations of Niemand et al. (2012) and DeMott et al. (2015), we note that this simplification does not capture the full mineralogical complexity of atmospheric dust.

To define the upper limit of the potential INP reservoir and strictly evaluate the immersion freezing pathway, we assume that all aerosol particles are sufficiently hygroscopic to activate as cloud droplets upon entering the cloud layer. Consequently, above the liquid cloud base, the entire aerosol population is assumed to be immersed within supercooled liquid droplets. To evaluate freezing at subsaturated conditions, i.e., at RH > 90 % in this case (i.e., proceeding from a layer just below the cloud base), assuming the INP is engulfed by an aqueous solution below cloud base, we also include a fourth variation, termed ABIFM*, which permits nucleation to occur in subsaturated conditions (i.e., from diluted aqueous solutions) (Knopf and Alpert, 2013). The calculation of INP arrays for application of the singular freezing parameterizations can be found in Appendix B. These arrays are constructed to discretize the cumulative INP spectrum into differential temperature bins, representing the specific subset of particles that activate within each temperature interval. Prognostic treatments are used for the number concentrations of activatable INPs (the reservoir of INP that can actually be activated within a given cloud layer, refer to N_INP) and N_i and detailed prognostic equations are given in Appendix C. The stability of our simulation calculation is evaluated by ensuring that results are insensitive to a doubling of vertical resolution (5 m) and a much shorter time step (1 s) (Fig. S1 in the Supplement). The model framework allows for the simultaneous tracking of multiple aerosol types, each with distinct PSDs and ice nucleation properties, facilitating the External simulations described in Sect. 2.5. Furthermore, the prognostic equations (Appendix C) enable the diagnosis of individual budget terms for INP number concentration, facilitating the process-level analysis.

2.4 Numerical Implementation of Temperature-Threshold and Rate-Based Freezing

To accurately represent the different physical bases of the parameterizations, we employ distinct numerical implementations. In the following, we adopt the terminology of “singular” (representing temperature-threshold schemes like INN and INAS) and “CNT-based” (representing rate-based schemes like ABIFM). For singular schemes (INN, INAS), the number of activated INPs is determined instantaneously based on the current thermodynamic conditions (temperature, surface area, number concentration).

In contrast, the CNT-based scheme (ABIFM) describes a rate of freezing. To ensure numerical stability while capturing the time-dependence without computationally expensive Lagrangian particle tracking, we calculate the activation source term (S_act) using an implicit form:

$$\begin{array}{}\text{(1)}& \begin{array}{l}{S}_{\mathrm{act}}=\frac{N_{\mathrm{aer}}\left(t\right)}{\mathit{\delta }t+{\mathit{\tau }}_{\mathrm{act}}}\end{array}\end{array}$$

where N_aer(t) is the available aerosol reservoir, δt is the model time step, and the τ_act is the characteristic activation time scale defined by the heterogeneous nucleation rate coefficient (J_het) and particle surface area (πd²):

$$\begin{array}{}\text{(2)}& \begin{array}{l}{\mathit{\tau }}_{\mathrm{act}}=\frac{\mathrm{1}}{J_{\mathrm{het}}\left(T,a_{\mathrm{w}}\right)\mathit{\pi }{d}^{\mathrm{2}}}\end{array}\end{array}$$

The change in ice crystal number concentration (ΔN_i) over one time step is the product of the source rate and the time step. Substituting the definition of τ_act, this can be expressed as:

$$\begin{array}{}\text{(3)}& \begin{array}{l}\mathrm{\Delta }{N}_{\mathrm{i}}=S_{\mathrm{act}}\cdot \mathit{\delta }t=N_{\mathrm{aer}}\left(t\right)\cdot \left(\frac{\mathit{\delta }t}{\mathit{\delta }t+{\mathit{\tau }}_{\mathrm{act}}}\right)\end{array}\end{array}$$

In the conditions relevant to this study, the characteristic freezing time (τ_act) is typically orders of magnitude larger than the model time step (δt= 10 s). For example, mineral dust at −19 °C, τ_act≈ 30 h. Under these conditions (δt≪ τ_act), the denominator is dominated by τ_act, and the equation simplifies to a linear form:

$$\begin{array}{}\text{(4)}& \begin{array}{l}\mathrm{\Delta }{N}_{\mathrm{i}}\approx N_{\mathrm{aer}}\left(t\right)\cdot \left(\frac{\mathit{\delta }t}{{\mathit{\tau }}_{\mathrm{act}}}\right)=N_{\mathrm{aer}}\left(t\right)\cdot \left(J_{\mathrm{het}}\mathit{\pi }{d}^{\mathrm{2}}\right)\cdot \mathit{\delta }t\end{array}\end{array}$$

This allows for the prognostic evolution of the INP reservoir without the need to track the time history of individual aerosol particles.

Finally, to verify that defining the total aerosol population as the INP reservoir does not lead to unrealistic instantaneous depletion of the aerosol population, we calculated the cumulative frozen fraction over the simulation duration. As detailed in Table S1, for representative cloud conditions (−19.0 °C), only ∼ 28 % of the most active reservoir (mineral dust) activates after 10 h, while organic and sea spray activated fractions remain negligible (< 1 %). This confirms that treating the total aerosol population as the potential reservoir is physically robust and does not result in premature exhaustion of INPs.

2.5 Sensitivity case studies

We investigate the sensitivity of N_INP, N_i and ice crystal formation rate (d$N_{\mathrm{i}}/\mathrm{d}t)$ with respect to three immersion parameterizations and several key cloud parameters. A key distinction between the schemes, which governs their fundamentally different behaviors, is the relationship of the INP reservoir (N_INP) to the total aerosol population in its composition class (N_aer). For instance, mineral dust, a canonical INP class, is highly efficient at nucleating ice but often represents a numerically negligible fraction of N_aer compared to more abundant types like organics or sea spray aerosol. In singular schemes, N_INP represents a limited subset of N_aer that are in a composition class based on the temperature at which they may be activated. In contrast, for the CNT-based schemes, all aerosol particles in an INP composition class are treated as potential INPs, making the INP reservoir effectively equivalent to the total aerosol population in that composition class, N_aer. We vary the number concentration of aerosols in each composition class (N_aer), cloud cooling rate (CCR), cloud-top entrainment rate (w_e) and ice crystal fall velocity (v_f) while applying three aerosol particle types including mineral dust, organic aerosol, and SSA PSDs.

We examine various simulation setups as given in Table 4. The thermodynamic conditions and cloud parameters from the LES baseline results of the SHEBA case study serve as the control run (hereafter referred to as CTRL) while applying the three different aerosol PSDs individually. Note that the baseline CTRL simulations are initialized with single aerosol types to isolate their specific freezing behaviors. We note that the CTRL setup differs slightly from the one used in our previous study (Knopf et al., 2023b). Detailed information on these modifications can be found in Appendix B.

Table 4 Parameter choices of the different 1D aerosol-cloud model simulation. Imm_CTRL: The CTRL run with the baseline settings for all IMF parameterizations, no perturbations, used as a reference. h_res_t and h_res_z: Simulation applying higher resolution with doubly refined vertical resolution (5 m) and much smaller time step (1 s).

To examine the sensitivity of cloud properties (N_INP, N_i and d$N_{\mathrm{i}}/\mathrm{d}t)$ towards the parameters N_aer, CCR, w_e and v_f, a series of simulations for each IMF parameterization are carried out by repeating the CTRL case with adjusting the targeted parameters (see Table 4). A detailed summary of all simulations is included in Table 5.

Table 5 Summary of AC-1D model simulations.

We multiply the N_aer by a factor of 10 and 10⁻¹ for CTRL termed experiments N_aer× 10 and N_aer× 01, respectively. A total of 36 cases are obtained, exploring three different aerosol PSDs (mineral dust, organic and SSA particles), four different IMF freezing parameterizations (INN, INAS, ABIFM and ABIFM^*) and three sets of experiments (CTRL, N_aer× 10, N_aer× 01) (see Table 5).

In order to assess the sensitivity of cloud properties to CCR changes, the temperature profiles for CTRL are modified by applying a cooling rate of 0.3 °C h⁻¹ to each layer (see Fig. S2). This means that after 10 hours, the entire layer is 3 °C cooler than at the start of the simulation, termed experiment CCR = 0.3 (see Table 4). This results in 24 cases consisting of three aerosol PSDs, four freezing parameterizations and two sets of cases (CTRL for singular IMF parameterizations and CNT-based IMF parameterizations and CCR = 0.3) (see Table 5). The evolution of the temperature profiles for CTRL, and CCR = 0.3 are presented in Figs. 1 and S2.

Additionally, we investigate the responses of the INPs and N_i evolution to changes in w_e by repeating the CTRL cases while increasing w_e from 0.1 to 1.0 cm s⁻¹, termed experiments w_e= 1.0 which contain 24 cases including three aerosol PSDs (mineral dust, organic, and SSA respectively), four IMF freezing parameterizations and two sets of experiments (CTRL, w_e =  1.0) (see Table 5).

Lastly, the response of INPs and N_i to changes in v_f is evaluated by changing the v_f from 0.3 to 1.0 m s⁻¹, termed experiments v_f= 1.0 containing 24 cases in total, consisting of three aerosol PSDs, four IMF parameterizations, and two sets of experiments (CTRL, v_f= 1.0) (see Table 5).

2.6 Prognostic Evaluation of IMF Parameterizations in the Presence of Different Aerosol PSDs

To move beyond idealized single aerosol type studies and demonstrate the model's capability to simulate aerosol populations more representative of the Arctic atmosphere, we conduct simulations representing an external mixture of mineral dust, organic, and SSA particles each assigned a unique, observationally-based particle size distribution (Table 1) (Riemer et al., 2019). We term this simulation as “External”. This approach ensures that particles from different sources are treated as physically separate entities while coexisting in the same simulation. The primary goal is to prognostically evaluate how competing aerosol types and associated IMF parameterizations impact N_INP and N_i, with freezing initiated by either INN, INAS, or ABIFM as described in Appendix A.

3 Results

3.1 Influence of Immersion Freezing Parameterization on INP and Ice Crystal Evolution

As a crucial step for model evaluation, we first place the chosen IMF parameterizations into an observational context. This provides a baseline for interpreting the prognostic simulations that follow. To this end, we compare the activated INP concentrations predicted by each IMF parameterization (initialized with the specific aerosol PSDs from Table 1) against recent Arctic field observations.

Figure 3 compares these model predictions against a composite of Arctic data. The shaded regions represent the annual variability observed at four ground-based Arctic stations (Alert, Utqiaġvik, Ny-Ålesund, and Villum Research Station) utilizing filter samples for ice nucleation experiments analyzed with a cooling rate of 1 °C min⁻¹ (Wex et al., 2019). Additionally, we show INP measurements recorded during the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition onboard the research vessel Polarstern (Creamean et al., 2022). These particle samples were analyzed using a cooling rate of 0.33 °C min⁻¹. The vertical bars indicate the minimum and maximum INP concentrations observed during the full annual cycle at −17.5, −20, and −22.5 °C, selected to match the temperature range of our simulated cloud layer. For the time-dependent ABIFM (CNT) parameterization, a nucleation time period is required to derive the cumulative activated INP number concentration. We apply t= 1 min (following Alpert et al., 2022) to approximate the nucleation time scales for the experiments.

Figure 3 Predicted activated INP concentrations from immersion freezing parameterizations are shown for mineral dust (black lines), SSA (blue lines), and organic (green lines) aerosols, calculated with their respective PSD from Table 1. Lines indicate the different schemes: INN (solid), INAS (dashed), and ABIFM (dotted). Note that the INN schemes use DeMott et al. (2015) for mineral dust and DeMott et al. (2010) for SSA and organic aerosol particles. Model predictions are compared against a composite of Arctic field observations from Wex et al. (2019) (colored shaded regions reflect Alert, Utqiaġvik, Ny-Ålesund, Villum) and Creamean et al. (2022) (vertical orange bars). The vertical pink background shading indicates the temperature range in the simulation domain.

As shown in Fig. 3, the agreement depends on both the parameterization choice and the assumed aerosol type. The singular INAS scheme, when applied to SSA and organic aerosol PSDs, shows good agreement within the range of field observations. The ABIFM (CNT) predictions for organic aerosol also align well with the annual range observed during the MOSAiC campaign. While the singular INN scheme and dust parameterizations tend to predict concentrations at the upper end of the observational range, they generally fall within the total variability spanning the different stations and seasons. Except for a pure dust case, INAS and CNT parameterizations predict INP concentrations similar to the ones observed in the Arctic regions.

The choice of IMF parameterization fundamentally dictates the simulated temporal evolution of N_INP, N_i, d$N_{\mathrm{i}}/\mathrm{d}t$. Figure 4 illustrates the 10 h domain-averaged time series of N_INP for control (CTRL) simulations under different IMF schemes and aerosol types, with Figs. S3, S4 and S5 providing a detailed view of the initial 0.1 h. The significant differences in predicted INP concentrations when applying different IMF parameterizations to identical particle size distributions are further illustrated in Fig. S6, which demonstrates how parameterization choice can lead to orders-of-magnitude variations in INP predictions. For singular IMF schemes (INN, INAS), where INP activation is treated as instantaneous once temperature criteria are met, a rapid decrease in N_INP is observed, typically by over 90 % from initial values within the first hour (Figs. 4, S3, Table 6).

Table 6 The temporal evolution of the domain-averaged number concentration of activatable ice-nucleating particles (N_INP in L⁻¹) for given aerosol PSDs (dust, organics and SSA described in Table 1) normalized by the initially activatable INP number concentration, that is $N_{\mathrm{INP}}/N_{\mathrm{INP}}^{t=\mathrm{0}}$, using singular ice nucleation number based (INN), the singular ice nucleation active sites (INAS), the classical nucleation theory (CNT) water activity based immersion freezing model (ABIFM) and ABIFM enabling subsaturated freezing (ABIFM^*). N_INP is given for model simulations described in Table 4 for model times 0.5, 5, and 10 h. Note: Values displayed as 1.00 × 10⁰ (particularly for ABIFM cases) indicate that the relative change is less than 1 %, not that the process has ceased.

Figure 4 Time series of simulated domain-averaged activatable INP number concentration (N_INP in L⁻¹) separated by parameterization type: Singular schemes (INN, INAS; left column) and CNT schemes (ABIFM, ABIFM*; right column). The panels display the baseline control (CTRL) simulation (a, b), the sensitivity to cloud cooling rate (c, d), and the sensitivity to cloud-top entrainment rate and sensitivity experiment. Simulations are initialized with different aerosol PSDs (dust, organic, and SSA particles), immersion freezing parameterizations (INN, INAS, ABIFM, ABIFM^*) and cloud parameters (cloud cooling rate, cloud-top entrainment rate (e, f). Brown, green, and blue lines represent the application of aerosol PSDs of mineral dust, organic, and SSA particles, respectively. Immersion freezing parameterizations are distinguished by symbols: INN (no symbols), INAS (cross), ABIFM (triangle), and ABIFM* (circle). In panels (a)–(b), the thin solid lines indicate results with the baseline control cloud parameters (CTRL). The dashed lines (c–d) denote results with the cloud cooling rate (CCR) of 0.3 °C h⁻¹ (CCR = 0.3) and the dash-dotted lines (e–f) show the results with the cloud-top entrainment rate (w_e) of 1 cm s⁻¹ (w_e= 1.0).

This rapid consumption of the initial, limited INP reservoir in singular schemes directly impacts ice formation, as shown by the time series of domain-averaged N_i in Fig. 5. An initial sharp peak in N_i occurs around 0.2 h, followed by a substantial decrease. For the CTRL cases, N_i decreases by 86 %–98 % from this peak value after 10 h.

Figure 5 Temporal evolution of the domain-averaged ice crystal number concentration (N_i in L⁻¹) in response to different cloud system parameters. The nine panels are organized by aerosol type in rows (mineral dust, top; organic, middle; and sea spray aerosol (SSA), bottom) and by sensitivity experiment in columns. The columns from left to right represent simulations with an applied cloud cooling rate (CCR), an enhanced entrainment rate (w_e), and an increased ice crystal fall speed (v_f), respectively. Within each panel, different line styles and colors represent the four immersion freezing (IMF) parameterizations, with legends and styling identical to those used in Fig. 4.

The d$N_{\mathrm{i}}/\mathrm{d}t$ follows a similar trajectory, as depicted in Fig. 6. After an initial burst, the formation rate for singular schemes decreases by 3–4 orders of magnitude over the 10 h period, demonstrating the rapid decline in ice production rates following depletion of the initial INP reservoir.

Figure 6 Temporal evolution of the domain-averaged ice crystal formation rate (d$N_{\mathrm{i}}/\mathrm{d}t$ in L⁻¹ h⁻¹). The six panels are organized by aerosol type in rows (mineral dust, top; organic, middle; and sea spray aerosol (SSA), bottom) and by sensitivity experiment in columns. The columns from left to right show the response to an applied cloud cooling rate (CCR), and enhanced entrainment rate (w_e). Legends and line styles for the four immersion freezing (IMF) parameterizations are identical to those used in Fig. 4.

To provide a more quantitative summary of these trends, the evolution of N_i is normalized and presented in Table 7. Unlike N_INP, which has a defined value at time zero, N_i only begins to form after the first time step when nucleation processes are activated. Therefore, Table 7 normalizes N_i by its value at the first model output time (t = 10 s), providing detailed values that confirm the dramatic decrease observed in the figures.

Table 7 The domain-averaged number concentration of ice crystals (N_i in L⁻¹) normalized by the initial ice crystal number concentration at the first model output time, that is $N_{\mathrm{i}}/N_{\mathrm{i}}^{t=\mathrm{10}\mathrm{s}}$ following Table 6.

Correspondingly, since ice crystal formation rates are only meaningful after nucleation begins, Table 8 normalizes the d$N_{\mathrm{i}}/\mathrm{d}t$ values using the initial rate at t = 10 s as the reference. This allows for a clear, quantitative assessment of the rapid decline in ice production rates over the simulation period.

Table 8 The domain-averaged ice crystal formation rate (d$N_{\mathrm{i}}/\mathrm{d}t$ in L⁻¹ h⁻¹) normalized by the initial ice crystal number concentration at the first model output time, that is d$N_{\mathrm{i}}/\mathrm{d}t/\mathrm{d}{N}_{\mathrm{i}}/\mathrm{d}{t}^{t=\mathrm{10}\mathrm{s}}$ following Table 6.

In contrast, simulations employing CNT-based IMF parameterizations (ABIFM, ABIFM*), where all aerosol particles are potential INPs and activate via a continuous, time dependent freezing process, exhibit a substantially larger and more stable N_INP reservoir. N_INP remain within ∼ 2 % of its initial value at t= 0 throughout the 10 h simulation (Fig. 4b). Table 6 presents these values normalized by their initial value at t= 0. Consequently, N_i in CNT-based simulations reach a quasi-stable plateau after an initial increase and remain orders of magnitude higher compared to the case of singular schemes (Fig. 5 and Table 7). This plateau represents a dynamic equilibrium where the continuous stochastic ice production (sustained by the large INP reservoir) is balanced by the removal of ice crystals via sedimentation. Similarly, d$N_{\mathrm{i}}/\mathrm{d}t$ for CNT schemes is maintained at significantly higher levels (Fig. 6 and Table 8). Simulations with ABIFM^* produce domain-averaged N_INP and N_i trends broadly similar to the ABIFM cases (Figs. 4 and 5).

3.2 Impact of Aerosol Type and Number Size Distribution

While the choice of parameterization establishes the overarching framework for INP and ice evolution, the specific type and size distribution of aerosol particles further modulate ice nucleation efficiency within these frameworks. The impact of aerosol identity (mineral dust, organic, SSA) on N_INP, N_i, and d$N_{\mathrm{i}}/\mathrm{d}t$ varies significantly depending on whether INP activation is parameterized based on INN, INAS, or ABIFM.

For the INAS parameterization, which ties ice nucleation to the ice-active site density (n_s) of particles (Eqs. A4–A7), aerosol composition is paramount. Despite mineral dust having the lowest total aerosol number concentration (N_aer) among the considered types (Table 1), it consistently yields the highest N_INP and subsequently N_i and d$N_{\mathrm{i}}/\mathrm{d}t$ in CTRL simulations (see Figs. 4, 5 and 6, green lines vs. black/blue for INAS). This is directly attributable to the significantly higher n_s values prescribed for mineral dust compared to organic and SSA particles at the relevant cloud temperatures (e.g., Eq. A5 vs. Eqs. A6 and A7). Conversely, for the INN parameterization, which primarily depends on the number concentration of aerosol particles larger than 0.5 µm (N_aer > 0.5 µm, Eqs. A1–A2), the aerosol type itself (beyond its contribution to N_aer > 0.5 µm) plays a lesser role. Consequently, SSA, with the highest N_aer > 0.5 µm in our PSDs (Table 1), results in the largest initial N_INP for INN.

Under the CNT-based ABIFM scheme, where all aerosol particles are potential INPs (N_INP=N_aer), the initial N_INP primarily reflects the total aerosol loading of each type. However, the continuous ice formation rates are governed by the heterogeneous ice nucleation rate coefficient (J_het), which is determined by the aerosol composition (Eqs. A10–A12). Mineral dust, with its generally higher J_het values, tends to produce more ice crystals and at a faster rate compared to organic or SSA particles for a given N_aer and thermodynamic condition (Figs. 5 and 6, ABIFM triangles). Thus, also when applying a time-dependent freezing process, aerosol identity remains a key factor in determining the intensity of ice production. The different PSDs associated with each aerosol type (Fig. 2 and Table 1) inherently influence these outcomes by determining the total number and surface area of particles available for nucleation, regardless of the parameterization.

Beyond these inherent differences between aerosol types, we test the system's sensitivity by scaling the total aerosol loading by factors of 10 and 0.1. The results (see Fig. S7) show a straightforward, linear response: the concentrations of N_INP and N_i, as well as d$N_{\mathrm{i}}/\mathrm{d}t$, scale proportionally with the initial aerosol loading for all parameterizations.

3.3 Sensitivity to Cloud System Parameters

Beyond the intrinsic properties of aerosols and the choice of freezing parameterization, the dynamic and thermodynamic environment of the cloud system, specifically CCR, w_e, and v_f, exerts significant control over the evolution of INPs and ice crystals.

3.3.1 Cloud Cooling Rate (CCR): A Powerful Driver of Ice Production Intensity

Continuous radiative cooling concentrated at cloud top supports turbulent mixing and leads to progressively declining temperatures throughout the cloud domain, which can drive further INP activation by lowering cloud temperatures (Morrison et al., 2011). The impact of an applied CCR of 0.3 °C h⁻¹ is most pronounced for singular parameterizations (Figs. 4c, 5, and 6). For these parameterizations, cooling increases the number concentration of activatable INP by lowering the ambient temperature to meet the fixed activation thresholds of a progressively larger and colder subset of the total INP population (see Eq. B2). For example, with INAS, the N_INP initially declines due to activation as shown in Fig. 4a (lines with cross signs). However, under continuous cooling (Fig. 4c, dashed lines with cross signs), it can recover and even exceed its initial value after several hours, especially for organic aerosols which exhibits high sensitivity of N_INP to temperature changes. This heightened sensitivity arises because the INAS density formulation for organic particles (China et al., 2017) exhibits a significantly steeper slope with respect to temperature compared to the mineral dust parameterization in this temperature range. Table 6 shows that the organic INAS_CCR = 0.3 case reaches over 10 times of initial N_INP after 10 h. This translates to a large increase in N_i (up to a factor of ∼ 300 for organic INAS, see Table 7) and sustains, or even increases, d$N_{\mathrm{i}}/\mathrm{d}t$ compared to CTRL runs where N_i rapidly depletes (Figs. 5 and 6).

For CNT-based approaches, the impact of CCR on the already vast N_INP reservoir is minimal (see Table 6, N_INP changes by < 3 % for ABIFM). However, lower temperatures significantly enhance J_het, leading to substantially increased N_i (e.g., N_i increases by factors of ∼ 3–4 for ABIFM with CCR, Table 7) and d$N_{\mathrm{i}}/\mathrm{d}t$ (Figs. 5 and 6). This highlights that while the potential INP reservoir in CNT-based schemes is less sensitive to cooling-induced expansion, the actual rate of ice formation remains highly sensitive to temperature.

3.3.2 Cloud-Top Entrainment Rate (w_e): A Critical INP Source for Depletion-Prone Singular Schemes

Cloud-top entrainment of free-tropospheric air provides a mechanism for replenishing INPs. In our simulations with instantaneous activation for singular schemes, entrained INPs are immediately converted to ice crystals. Therefore, while the standing N_INP concentration within the cloud layer does not show a sustained increase due to entrainment (Fig. 4a, singular cases, w_e= 1.0 cm s⁻¹ vs. CTRL), this process acts as a continuous source flux (S_ent, Fig. 11) sustaining ice production. This is evident in the significantly higher N_i and d$N_{\mathrm{i}}/\mathrm{d}t$ observed in singular scheme simulations with enhanced entrainment (w_e= 1.0 cm s⁻¹) compared to CTRL runs, particularly after the initial in-cloud INP reservoir is depleted (Figs. 4 and 5). For instance, N_i for INAS with w_e= 1.0 cm s⁻¹ can be 3–7 times higher than CTRL after 10 hours (see Table 7, comparing w_e= 1.0 to CTRL for INAS).

For CNT-based approaches, entrainment directly adds to the total aerosol particle population (N_aer), which constitutes the N_INP reservoir. However, minor effects of entrainment can be found. The impact on N_i and d$N_{\mathrm{i}}/\mathrm{d}t$, while positive, is less significant in relative terms compared to singular schemes (e.g., N_i increases by ∼ 2 %–3 % for ABIFM with w_e= 1.0 cm s⁻¹ after 10 h, Table 7), as the initial reservoir is already substantial and not the primary limiting factor for ice production.

3.3.3 Ice Crystal Fall Speed (v_f): The Primary Sink for Ice Crystals

Ice crystal sedimentation is the sole loss mechanism for ice crystals in our AC-1D model. Increasing the number-weighted ice fall speed (v_f) from 0.3 m s⁻¹ (CTRL) to 1.0 m s⁻¹ leads to a more rapid removal of ice crystals from the cloud layer. This results in substantially lower N_i across all parameterizations and aerosol types (Fig. 6). Specifically, increasing the fall velocity reduces N_i by approximately 70 % (a factor of around 3.5) relative to the control run for both singular and CNT schemes. This consistent relative reduction demonstrates that sedimentation acts as a uniform sink for all parameterizations, while the absolute resulting concentrations are scaled by the strength of the INP source term. The efficiency of sedimentation plays a crucial role in modulating the standing N_i and the lifetime of ice within the mixed-phase cloud.

3.4 Vertical Profiles Reveal Contrasting Sensitivities

While domain-averaged properties provide a valuable overview of the cloud system's response, they can obscure critical height-dependent processes that govern cloud evolution. To deconstruct these mechanisms, we now analyze the vertical profiles of the time-averaged change in INP number concentration (ΔN_INP(z)), ice crystal number concentration (ΔN_i(z)), and ice crystal formation rate ($\mathrm{\Delta }\mathrm{d}{N}_{\mathrm{i}}/\mathrm{d}t\left(z\right)$) in response to changes in CCR, w_e, and v_f. This approach allows us to pinpoint where in the cloud these forcings have their greatest impact and how that impact differs between IMF parameterizations.

3.4.1 INP Concentration Response (ΔN_INP(z))

Figure 7 reveals that the INP reservoir responds in fundamentally different ways to cloud cooling and cloud-top entrainment. For the singular INN and INAS schemes, CCR acts as a powerful volume-wide forcing, increasing N_INP to values up to ∼ 4.6 and ∼ 21.4 times the original, respectively, throughout the cloud layer, with the maximum increase observed in the lower half of the cloud layer (Fig. 7a, b). The CNT-based ABIFM scheme shows a different behavior: its large INP reservoir is less sensitive to both cooling and entrainment, with ΔN_INP remaining less than 1 % throughout the vertical column (Fig. 7c, d).

Figure 7 Vertical profiles of the change in number concentration of activatable INP (ΔN_INP(z) in %) averaged over entire 10 h of simulation time. ΔN_INP differs compared to the respective CTRL runs due to the change of cloud parameters (cloud cooling rate, cloud-top entrainment rate) applying dust, organic and SSA particles, given as brown, green, and blue lines, respectively. Different immersion freezing parameterizations are applied including (a) ice nucleation number based (INN), (b) ice-nucleation active sites (INAS), (c) water-activity based immersion freezing model (ABIFM), and (d) ABIFM enabling subsaturated freezing (d) ABIFM^*. Simulation results for changing cloud cooling rate (solid lines) and cloud-top entrainment rate (dashed lines) are shown. The blue shaded area denotes the cloud layer and the red line in the middle indicates the zero-change reference line (physically representing no deviation from the control simulation).

3.4.2 Ice Crystal Concentration Response (ΔN_i(z))

The response of the ice crystal concentration (ΔN_i(z)), shown in Fig. 8, directly reflects these differing sensitivities. For the INAS scheme, CCR is the dominant factor, driving an enormous ΔN_i(z) of up to 20-fold (Fig. 8b). The ABIFM scheme is similarly dominated by cooling, which yields a ΔN_i(z) of up to tenfold (Fig. 8c). Conversely, the INN scheme is more sensitive to entrainment compared to the ABIFM scheme, which causes a nearly five-fold increase in N_i at the cloud top (Fig. 8a). Increasing the v_f acts as a powerful, universal sink. This leads to a substantial negative change in N_i across all parameterizations. As illustrated in Fig. 8 dotted lines, increasing the fall velocity reduces the ice concentration by over 90 % throughout the domain.

Figure 8 As in Fig. 7 but for the change in number concentration of ice crystals (ΔN_i(z) in %).

3.4.3 Ice Crystal Formation Rate Response ($\mathrm{\Delta }{\mathbf{dN}}_{\mathbf{i}}/\mathbf{dt}\left(\mathbf{z}\right)$)

The change in the ice crystal formation rate ($\mathrm{\Delta }\mathrm{d}{N}_{\mathrm{i}}/\mathrm{d}t\left(z\right)$) in Fig. 9 confirms these findings. Cloud cooling enhances the formation rate throughout the entire cloud layer for all parameterizations, with the most dramatic impact seen in the INAS scheme, where the formation rate shows a peak increase over 25-fold near the cloud top (Fig. 9b). The effect of enhancing cloud-top entrainment is, almost exclusively a cloud-top phenomenon. It provides a significant boost to the formation rate for the singular schemes (increasing it more than fivefold for INN) but has a minor effect on the CNT-based scheme.

Figure 9 As in Fig. 7 but for the change in ice crystal formation rate ($\mathrm{\Delta }\mathrm{d}{N}_{\mathrm{i}}/\mathrm{d}t\left(z\right)$).

3.5 Competing Aerosol in Cloud Ice Formation

The AC-1D model is configured to simulate a more realistic scenario, termed “external”, where mineral dust, organic, and SSA populations coexist, each with their distinct PSDs and INP parameterizations. This approach allows for an examination of the model's ability to handle aerosol complexity and, more importantly, to investigate the relative contributions and potential dominance of different aerosol types (Fig. 10). The results show that the dominant source of ice is governed by a combination of aerosol type and the chosen freezing parameterization. Since the INAS and ABIFM schemes for specific aerosol types are derived from common datasets (Sect. 2.3), these differences can be attributed to the structural formulation of the parameterization rather than discrepancies in the underlying observations.

Figure 10 Initial total aerosol population, activatable INP composition distribution, and resulting ice crystal number contributions (N_i) using different immersion freezing parameterizations for the externally mixed aerosol scenario (refer to Table 4, case “External”). Panel (a) shows the externally mixed aerosol composition. Simulations were run with an aerosol population consisting of dust (brown), organic (green), and SSA (blue) particles. Panels (b)–(d) show the initial (t= 0) fractional composition of the activatable INP number concentration for each parameterization. Panels (e)–(g) show the resulting domain-averaged N_i over 10 h and its contribution from each aerosol type. Each column represents a simulation applying INN_External (b, e), INAS_External (c, f), and ABIFM_External (d, g).

Figure 10a shows the relative contribution of compositionally segregated subpopulations (mineral dust, organic, and sea spray aerosol particles) to the total aerosol number concentration. The aerosol population for these simulations is numerically dominated by organic and SSA particles, which together comprise over 96 % of the total particle number concentration. Mineral dust, in contrast, is a minor component by number.

The simulation INN_External identifies 5.32 L⁻¹ activatable INPs from particles larger than 0.5 µm, with SSA and organic particles comprising > 98 % of the initial activatable INPs (Fig. 10b). In this case, the INP concentration, determined solely by temperature and aerosol number concentration, is controlled by the most abundant particles throughout the simulation. SSA serves as the primary ice source (Fig. 10e), with the combined ice concentration decreasing by approximately two orders of magnitude over 10 h, from ∼ 1 to 10⁻² L⁻¹ .

In contrast, the INAS_External produces an initial activatable INP concentration of only 1.08 L⁻¹, composed almost exclusively (> 96 %) of mineral dust particles (Fig. 10c). Despite dust representing < 2 % of the total aerosol population, its n_s exceeds that of organic aerosols and SSA by several orders of magnitude, thereby dominating ice formation. Accordingly, dust particles control the subsequent ice crystal formation (Fig. 10f). The ice crystal population is composed almost entirely of crystals formed on mineral dust, with negligible contribution from the more abundant but less efficient organic and SSA particles.

The ABIFM_External yields the highest activatable INP concentration (2.31 × 105 L⁻¹), equal to the total aerosol population (Fig. 10a and d). Ice production under ABIFM exhibits markedly different dynamics. Despite starting with the fewest activatable INPs, dust ultimately emerges as the dominant source of ice crystals during the 10 h simulation, followed by SSA and then organic particles (Fig. 10g). Ice crystal number concentration reaching ∼ 2.3 L⁻¹ after 10 h. This is two to three orders of magnitude higher than the final concentrations in the singular frameworks.

3.6 Dominant Processes in the INP Budget

To quantitatively diagnose the controlling factors of the simulated INP budget, a comparative analysis is performed for the singular (INAS) and CNT-based (ABIFM) parameterizations, specifically for the CTRL simulations initialized with the mineral dust aerosol PSD. All budget terms are expressed as domain-averaged rates. The time-averaged budget balances (Fig. 11a, c) provide a concise summary of the relative magnitudes of INP sources (entrainment) and sinks (activation) over the 10 h simulation, while the corresponding temporal evolution plots (Fig. 11b, d) reveal the dynamic mechanisms responsible for these balances.

Figure 11 A comparative budget analysis for the control (CTRL) simulation with the mineral dust PSD, illustrating the fundamentally different controls on INP reservoirs applying either singular parameterization (INAS) (top row) and the CNT-based parameterization (ABIFM) (bottom row). Panels (a) and (c) provide a quantitative summary of the time-averaged domain INP budget balance. Panels (b) and (d) reveal the underlying dynamic evolution of the absolute magnitudes of these budget terms over the 10 h simulation. Note the logarithmic y axis in panels (b) and (d). All budget terms represent domain-averaged rates in units of L⁻¹ h⁻¹.

For the singular (INAS) parameterization, the time-averaged domain-averaged budget (Fig. 11a) shows a net sink, with the activation term (S_act≈ −0.06 L⁻¹ h⁻¹) being an order of magnitude larger than the entrainment source (S_ent ≈ 0.005 L⁻¹ h⁻¹). The temporal evolution of S_ent and S_act (Fig. 11b) reveals the mechanism behind this imbalance: following an initial rapid activation phase that depletes available INPs within the first hour, the system transitions to an entrainment-limited regime. By approximately 2 h, the entrainment source and activation sink converge to nearly equal magnitudes (∼ 0.03 L⁻¹ h⁻¹), establishing a quasi-steady state where ice formation is controlled by the rate at which new INPs are supplied through entrainment.

In stark contrast, the CNT (ABIFM) parameterization exhibits fundamentally different dynamics. The time-averaged budget terms are two orders of magnitude larger, with a massive activation sink (S_act ≈ −7.4 L⁻¹ h⁻¹) continuously overwhelming the substantial entrainment source (S_ent≈ 0.5 L⁻¹ h⁻¹) (Fig. 11c). The temporal evolution (Fig. 11d) demonstrates that this imbalance is established instantaneously and persists throughout the simulation. Unlike the singular case, the activation rate never becomes entrainment-limited but instead maintains relatively high values (∼ 7.7 L⁻¹ h⁻¹), reflecting the continuous availability of potential INPs from the large aerosol reservoir characteristic of the stochastic freezing approach. If, however, this reservoir were to be significantly depleted (e.g., in an aerosol-poor environment or over much longer timescales), we would expect the activation rate to eventually become source-limited by entrainment, converging in principle with the behavior of the singular case.

4 Discussion

4.1 Parameterization Choice Dominates Simulated Ice Phase Evolution

The most striking outcome of this study is the profound impact of the chosen IMF parameterization on the simulated lifecycle of INPs and ice crystals in Arctic mixed-phase stratus over all composition classes and PSDs (Figs. 4, 5, and 6). The orders-of-magnitude difference in the predicted INP reservoir size and N_i between singular and CNT-based approaches highlights a fundamental structural uncertainty in representing ice nucleation. Singular schemes, by defining INPs as a limited, instantaneously activating subset of aerosols (e.g., DeMott et al., 2010; Niemand et al., 2012; DeMott et al., 2015), inherently lead to more rapid INP depletion and a more transient ice phase, rendering them highly dependent on continuous replenishment mechanisms. This aligns with previous modeling studies that have shown the necessity of INP replenishment mechanisms to sustain ice in such frameworks (Fridlind et al., 2012; Fridlind and Ackerman, 2018; Knopf et al., 2023b, and references therein).

Conversely, CNT-based approaches, which treat all aerosol particles in their composition class as potential INPs activating continuously over time (e.g., Koop et al., 2000; Knopf and Alpert, 2013, 2023), maintain a vast INP reservoir governed by time-dependent freezing kinetics. This results in more persistent and significantly higher N_i, suggesting that Arctic MPCs could sustain a more vigorous ice phase if nucleation follows such time-dependent kinetics. The implication is that the perceived need for extremely efficient INP replenishment or additional ice production mechanisms may be reduced if CNT-like processes are dominant. The similarity between ABIFM and ABIFM^* results further suggests that, for the simulated conditions, once saturation is achieved, in-cloud nucleation processes appear to be more critical than slow freezing in subsaturated layers below the cloud base for these simulations and conditions.

Finally, the distinct ice crystal evolution trends observed here highlight the sensitivity of the system to the representation of PIP. By disentangling this process from liquid-phase feedbacks and SIP, we demonstrate that the choice of freezing parameterization constitutes a foundational source of uncertainty in the simulated ice crystal budget.

4.2 Aerosol Characteristics as Modulators of Ice Formation

Within each parameterization framework, aerosol type and PSD further refine ice nucleation predictions. The strong dependence of INAS results on aerosol type (Figs. 4, 5, 6) underscores the importance of accurate n_s parameterizations (Niemand et al., 2012; China et al., 2017; Alpert et al., 2022). The dominance of mineral dust governing N_i, despite lower N_aer, is consistent with laboratory findings of its higher n_s values compared to many organic species or SSA at typical MPC temperatures (Hoose and Möhler, 2012; Kanji et al., 2017). The INN scheme's primary reliance on N_aer > 0.5 µm makes it more sensitive to the aerosol number concentration (Table 1) than to aerosol type differences, a simplification that might not capture the full spectrum of INP behavior (DeMott et al., 2010; Knopf et al., 2021). For CNT schemes, while all particles are potential INPs, the aerosol-specific J_het values (Alpert and Knopf, 2016) ensure that composition and particle size play a role in determining the rate of ice formation, with dust again often showing higher efficacy.

4.3 Interpretation of Cloud System Parameter Sensitivities

4.3.1 Cloud Cooling and Entrainment: Critical Forcing Mechanisms in Singular schemes

The heightened sensitivity of singular schemes to CCR and w_e (Figs. 4, 5, and 6; Tables 6 and 7) reinforces their reliance on external drivers to mitigate the rapid depletion of INPs. Our vertical profile analysis (Figs. 7–9) further illustrates these sensitivities by pinpointing how and where different forcings, cloud cooling, cloud-top entrainment impact the cloud properties. These results highlight a fundamental limitation: singular schemes cannot sustain ice crystal production on a par with CNT schemes without continuous external forcing, resulting in significant differences in ice production within persistent mixed-phase clouds.

In singular schemes, CCR increase the number concentration of activatable INP by lowering temperatures to the activation thresholds of more INPs (Ullrich et al., 2017). This explains the significant increase in ice formation across heights observed in the vertical profiles (Figs. 8 and 9). In contrast, entrainment becomes a direct flux of new INPs that are immediately converted to ice at the cloud top, rather than replenishing a standing INP reservoir. This is why N_INP itself (Fig. 4e) shows no direct change with w_e, while N_i is substantially enhanced, leading to a much higher final N_i than in the CTRL run (Table  7), with the increase in N_i sharply localized at the cloud top (Figs. 8 and 9).

CNT-based schemes exhibit markedly different behavior. Due to the large INP reservoir, their response is governed by nucleation kinetics. Cloud cooling therefore drives a strong, layer-wide increase in ice production by enhancing the nucleation rate, while the minor addition of particles from cloud-top entrainment has a negligible impact on the final ice crystal concentration (Table 7).

The different responses to cloud forcing are a direct consequence of how the INP reservoir is defined. Singular schemes, which assume a scarce INP reservoir derived from a subset of aerosols, are critically dependent on cloud cooling and cloud-top entrainment to mitigate their rapid depletion. In contrast, CNT-based schemes, which treat the entire aerosol population as an abundant potential reservoir, are highly responsive to cooling that accelerates freezing but are largely indifferent to the minor perturbation of entrainment. This foundational difference explains why for the singular approach, even strong, continuous CCR and cloud-top entrainment can only mitigate depletion, not prevent it – a direct consequence of its initial limited INP reservoir.

4.3.2 Sedimentation Controlling Ice Sink

Sedimentation, controlled by v_f, is the primary removal mechanism for ice crystals across all tested IMF parameterizations (Table 7). Our vertical analysis provides confirmation of this unifying role.

Regardless of the IMF parameterization or the specific forcing applied , the vertical profiles of the change in ice concentration (ΔN_i(z)) consistently show their greatest magnitude – whether positive or negative – in the upper portion of the cloud, with the effect diminishing toward the cloud base (Fig. 8). This robust vertical structure is the clear signature of gravitational removal acting on a distributed, in-cloud source of ice. The highest rates of ice production occur in the coldest, uppermost cloud layers. At the same time, sedimentation continuously removes ice crystals from all levels, transporting them downward.

This interplay between in-cloud ice production and removal creates a natural vertical gradient. The impact of any process perturbation is therefore strongest near the region of most active ice formation and is progressively diminished by sedimentation at lower altitudes.

These results demonstrate that ice crystal fall speed emerges as a critical parameter governing cloud ice distribution, regardless of the specific microphysical processes involved. Consequently, accurate representation of ice particle fall speeds becomes essential in models.

4.4 Competitive Ice Nucleation: Insights from Externally-Mixed Aerosol Simulations

Each IMF parameterization uses different physical assumptions that determine how aerosol particles contribute to the activatable INP population, creating distinct relationships between the initial aerosol population, the activatable INP reservoir, and the resulting N_i.

INN parameterization rely on aerosol number concentration. Since SSA and organic aerosols dominate the number concentration in this size range, they predominantly comprise the initial activatable INP pool (Fig. 10b). Consequently, the resulting ice crystal population reflects this abundance pattern, with SSA serving as the primary contributor to ice formation (Fig. 10e).

In contrast, the INAS parameterization is based on aerosol type and surface area. The freezing efficiency for each aerosol type is described by n_s. For the examined external aerosol, our simulations demonstrate that the n_s term governs the ice nucleation response. This is because, the n_s values for mineral dust exceed those of SSA and organic particles by several orders of magnitude. This substantial difference in intrinsic ice nucleation efficiency enables the relatively sparse dust particles to dominate the activatable INP reservoir (Fig. 10c), despite their lower number concentrations. Accordingly, dust particles control the subsequent ice crystal formation (Fig. 10f).

The CNT-based ABIFM parameterization assumes a time-dependent freezing process. Hence, all aerosol particles serve as potential INPs, such that the initial activatable INP reservoir mirrors the total aerosol number distribution, dominated by SSA and organic particles (Fig. 10d). However, the final N_i is dominated by dust (Fig. 10g). Similar to INAS, this is due to the larger nucleation rate coefficient (J_het) of dust compared to SSA and organic particles. Over the 10 h simulation period, the substantially higher nucleation rates of dust particles enable more rapid ice crystal production, ultimately outcompeting other aerosol types.

This study demonstrates our model's capacity to prognostically evaluate competing ice microphysics parameterizations using realistic aerosol PSDs. By simulating competitive nucleation between multiple aerosol populations, we move beyond single-type sensitivity tests toward conditions that better represent atmospheric complexity. These results address a fundamental question in ice microphysics: Is atmospheric ice formation controlled by rare but highly efficient particles, or by the collective contribution of abundant aerosols? Our modeling framework provides the quantitative tools needed to test these competing paradigms against field observations, ultimately constraining which parameterizations most accurately reproduce observed cloud microphysical properties.

4.5 Process-Level Controls on INP Evolution: A Budget Perspective

The budget analysis reveals contrasting behaviors between the INAS and CNT (ABIFM) parameterizations: although their time-averaged activation-to-entrainment ratios are similarly around 10:1 (Fig. 11a, c), the absolute magnitudes differ by two orders of magnitude, and the temporal evolution diverges markedly (Fig. 11b, d). The singular parameterization exemplifies source-limited behavior. Following rapid depletion of the INPs (Fig. 11b), ice production becomes entirely controlled by cloud-top entrainment of INPs. The modest long-term budget balance (Fig. 11a) reflects this supply limitation, where S_ent dictates ongoing ice formation, independent of the freezing parameterization.

In contrast, the CNT parameterization exhibits process-rate-limited dynamics. The sustained, high-magnitude activation sink (Fig. 11d) demonstrates nucleation acting on an abundant aerosol reservoir, with S_act vastly exceeding S_ent (Fig. 11c). Here, ice production is governed by the intrinsic freezing kinetics under ambient thermodynamic conditions, enabling rapid glaciation potential absent in the singular framework.

This process-level comparison demonstrates that choice of freezing parameterization represents a selection between fundamentally different hypotheses for PIP. The stark differences in baseline dynamics (Fig. 11) underscore how this choice directly impacts simulation results.

5 Summary and Conclusions

The objective of this study is to dissect the key processes and uncertainties that govern ice formation in Arctic mixed-phase clouds. First, we explore structural uncertainty in PIP by evaluating how three different immersion freezing ice nucleation parameterizations influence the microphysical evolution of mixed-phase Arctic stratus in the SHEBA MPC case, using three aerosol types with distinct particle size distributions. Second, we assess the sensitivity of INP reservoir dynamics – specifically N_INP and N_i – to variations in aerosol number concentration (N_aer) and key cloud and atmospheric parameters, including cloud cooling rate (CCR), entrainment velocity (w_e), and ice crystal fall speed (v_f). Furthermore, we investigate competitive nucleation among different aerosol types in a more realistic scenario using externally-mixed aerosol populations and conduct a detailed INP budget analysis to diagnose and quantify the dominant processes controlling the INP reservoir. We address these objectives using a version of the AC-1D model extended to include multiple aerosol types and their specific freezing parameterizations, informed by LES of the SHEBA case (Knopf et al., 2023b). The main findings can be briefly summarized as follows:

When treating INPs prognostically as done in this study, the choice of IMF parameterizations can yield contrasting concentrations and trends of simulated N_INP, N_i, and d$N_{\mathrm{i}}/\mathrm{d}t$. Compared to singular approaches, a time-dependent description (CNT) yields an orders of magnitude larger INP reservoir, a finding consistent with previous work (Knopf et al., 2023b). This leads to ice crystal formation at a substantially higher rate without significant INP depletion, resulting in a simulated ice phase that is more persistent and not limited by the same replenishment fluxes as the singular schemes.

Different aerosol types and associated PSDs significantly modify the magnitude of the simulated N_INP, N_i, and d$N_{\mathrm{i}}/\mathrm{d}t$ but do not alter their trends observed in the control run. The relative roles of aerosol types and associated PSDs are primarily dependent on the choice of IMF parameterizations.

Among all microphysical parameters tested, changing the N_aer elicits a directly proportional response in the simulated N_INP, N_i, and d$N_{\mathrm{i}}/\mathrm{d}t$ across all parameterizations. For example, increasing N_aer by a factor of 10 leads to a corresponding ten-fold increase in these variables, and vice-versa.

In contrast to the linear response observed for N_aer, cloud cooling (CCR), w_e, and v_f have non-linear effects, with the dominant process for ice production depending on the chosen parameterization. For singular approaches, where the INP reservoir is limited, ice production is critically dependent on the w_e for replenishment. For CNT-based approaches with an abundant INP reservoir, CCR becomes the dominant control; it can reverse the overall decline in ice crystal concentration and, specifically for organic aerosols, drive N_i to levels that substantially exceed the initial peak.

Initiating IMF under subsaturated conditions (when using ABIFM) leads to a more gradual trend in vertically resolved cloud properties around the cloud base compared to in-cloud only nucleation. These trends do not significantly impact the conclusions of the cases investigated in this study. However, we cannot fully rule out the potential importance of parameterization choices in accounting for IMF under subsaturated conditions more generally.

The AC-1D model simulation results support the following conclusions. Our findings underscore the importance of the IMF parameterization selection in modeling mixed-phase Arctic stratus. Although each examined IMF parameterization was derived from laboratory experiments, the significant differences, particularly between singular and CNT-based approaches, when simulating N_INP, N_i, and d$N_{\mathrm{i}}/\mathrm{d}t$ from aerosol PSDs highlight the profound uncertainty in extrapolating these frameworks to atmospheric conditions and the necessity to constrain them with in-situ observations. Considering that PIP is a ubiquitous phenomenon in Arctic MPCs, characterized by a low liquid water path (Silber et al., 2021), there is an urgent need for more laboratory and field research to better assess the scaling of IMF parameterizations under realistic MPC conditions (Knopf et al., 2020, 2021; Burrows et al., 2022; Knopf and Alpert, 2023).

The starkly different model behaviors emerging from these two different freezing parameterizations have significant implications for climate modeling. The CNT framework predicts sustained ice production without the same dependency on strong, continuous replenishment fluxes seen in the singular framework. As a consequence, interpretation of the evolution of a mixed-phase cloud including cloud lifetime and radiative properties will be different and will crucially depend on which freezing description is more representative of nature. A separate, practical advantage of the CNT approach is its computational simplicity. CNT-based methods do not require an INP array. This simplification removes the necessity of temperature-dependent INP arrays, reducing the complexity of the computational model. By eliminating the need to track additional variables for particle uniqueness, the model becomes more streamlined. Consequently, the time-dependent approach, which provides freezing rates, enhances the feasibility of incorporating prognostic capabilities into global models with size-resolved aerosol modules (Bauer et al., 2008; Liu et al., 2016).

We evaluated three typical aerosol particle types: mineral dust, organic and SSA. Other potentially significant INP types (such as biological particles or biomass-burning aerosols) could also impact ice formation rates in MPCs. In future studies, the impact of other atmospheric INPs on cloud glaciation should be investigated. Continued work is needed to develop accurate ice nucleation parameterizations that cover the temperature range of realistic MPC conditions.

The observed non-linear responses of N_INP, N_i, and d$N_{\mathrm{i}}/\mathrm{d}t$ to changes in CCR, w_e, and v_f highlight the importance of detailed cloud microphysical and dynamical measurements. Additional profiling observations of MPC temperatures, number concentration of free-troposphere aerosol particles, and the characterization of ice crystal fall velocity are essential for robustly evaluating the representation of INP reservoirs in mixed-phase stratus.

Finally, these findings have significant implications for the implementation of ice nucleation in global climate models (GCMs). First, the strong sensitivity of the ice crystal budget to reservoir depletion confirms that GCMs should treat INPs prognostically rather than diagnostically, to avoid unphysical ice production in long-lived clouds.

Second, the choice of parameterization alters the computational burden. Prognostic application of singular schemes generally requires carrying additional tracers to track the subset of particles that have already nucleated at specific temperatures (to prevent “double counting”). In contrast, the CNT approach follows a first-order rate law, eliminating the need for INP-specific history tracers and thus offering a computationally efficient pathway for prognostic implementation.

Regarding the time aspect, the implicit numerical formulation presented in Sect. 2.4 is robust for long GCM time steps (e.g., 30 min), since it integrates the same rate law over the full time interval and yields the same cumulative freezing as many shorter sub-steps. However, this stability relies on strict prognostic subtraction of activated INPs: applying time-dependent freezing rates to diagnostic aerosol fields over long time steps would lead to systematic overestimation of ice formation. Therefore, future GCM developments should prioritize coupling time-dependent nucleation schemes to prognostic aerosol and INP budgets.

Appendix A: Immersion freezing parameterizations

The singular INN parameterization of mineral dust particles acting as INPs is given by (DeMott et al., 2015, D2015)

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\sum _{T=T_{\min }}^{T_{\mathrm{INP}}}{N}_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INN}\right)}\left(T\right)& =\mathrm{cf}{\left(N_{\mathrm{aer}>\mathrm{0.5}\mathrm{\mu }\mathrm{m}}\right)}^{\mathrm{1.25}}\\ & e^{\left(\mathrm{0.46}\left(-T_{\mathrm{INP}}\right)-\mathrm{11.6}\right)}\end{array}\end{array}$$

$N_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INN}\right)}\left(T\right)$ is the number concentration of all activatable INPs which are particles that can serve as INPs and potentially be activated (hereafter referred to as INPs) at temperature T, cf is the calibration factor with a value of 3 suggested when applying this parameterization to atmospheric data, N_{aer > 0.5 µm} is the total number of aerosol particles with diameters larger than 0.5 µm, and T is the temperature in the unit of Celsius.

A similar singular INN parameterization of ambient aerosol particles acting as INPs is given by (DeMott et al., 2010, D2010)

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& \sum _{T=T_{\min }}^{T_{\mathrm{INP}}}{N}_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INN}\right)}\left(T\right)=\mathrm{0.0000594}\\ & {\left(-T\right)}^{\mathrm{3.33}}{\left(N_{\mathrm{aer}>\mathrm{0.5}\mathrm{\mu }\mathrm{m}}\right)}^{\left(\mathrm{0.0264}\left(-T_{\mathrm{INP}}\right)+\mathrm{0.0033}\right)}\end{array}\end{array}$$

We applied this IMF parameterization (Eq. A2) for the case of SSA and organic particles while noting that this does likely not reflect SSA and organic particles but represents ambient atmospheric particles in general.

In the AC-1D model, the INN parameterizations are implemented as follows:

$$\begin{array}{}\text{(A3)}& N_{\mathrm{INP}}^{\mathrm{Imm}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)=N_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INN}\right)}\left(z_i,t,T_{\mathrm{INP}}\right)H\left(T\left(z_i,t\right)\right)\end{array}$$

where $N_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INN}\right)}$ is multiplied by a Heaviside function H(T(z_i,t)), which equals to 1 when $T\left(z_i,t\right)<T_{\mathrm{INP}}$ and 0 otherwise. This function ensures that INP activation occurs only when the temperature is lower than the activation temperature of INPs. In the model configuration used in this study, the array of T_INP designates the activation temperature of INPs per temperature bin using geometric progression bins.

The INAS parameterization of natural dust particles acting as immersion INPs is given as (Niemand et al., 2012, N2012):

$$\begin{array}{}\text{(A4)}& {\displaystyle }\sum _{T=T_{\min }}^{T_{\mathrm{INP}}}{N}_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INAS}\right)}(T,d)=N_{\mathrm{aer}}\left(d\right)(\mathrm{1}-e^{-A_{\mathrm{aer}}\left(d\right)n_{\mathrm{s}}\left(T_{\mathrm{INP}}\right)})\text{(A5)}& {\displaystyle }{n}_{\mathrm{s}}^{\mathrm{dust}}\left(T_{\mathrm{INP}}\right)=e^{\left(-\mathrm{0.517}\left(-T_{\mathrm{INP}}\right)+\mathrm{8.934}\right)}.\end{array}$$

$N_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INAS}\right)}(T,d)$ is the number concentration of all activatable INPs in temperature bin T and diameter bin d, N_aer is the total number concentration of aerosol particles, A_aer is the individual particle surface area, n_s is the INAS density in units of m⁻².

For organic particles (humic-like substances) and SSA, the n_s is given by China et al. (2017) (C2017) and (Alpert, 2022) (A2022), respectively, as

$$\begin{array}{}\text{(A6)}& {\displaystyle }{n}_{\mathrm{s}}^{\mathrm{organic}}\left(T_{\mathrm{INP}}\right)={\mathrm{10}}^{\left(\mathrm{66.90259}\left(\mathrm{\Delta }{a}_{\mathrm{w}}\right)-\mathrm{12.322}\right)}\text{(A7)}& {\displaystyle }{n}_{\mathrm{s}}^{\mathrm{SSA}}\left(T_{\mathrm{INP}}\right)={\mathrm{10}}^{\left(\mathrm{24.02526}\left(\mathrm{\Delta }{a}_{\mathrm{w}}\right)-\mathrm{2.26105}\right)},\end{array}$$

where Δa_w is the water activity criterion, given below.

In this model, the INAS parameterizations are implemented as follows:

$$\begin{array}{}\text{(A8)}& N_{\mathrm{INP}}^{\mathrm{Imm}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)=N_{\mathrm{INP}}^{\mathrm{sing}\left(\mathrm{INAS}\right)}\left(z_i,t,d,T_{\mathrm{INP}}\right)H\left(T\left(z_i,t\right)\right).\end{array}$$

For the INAS case, the T_INP array is set same as in the INN case, while the d array is defined using 50 bins, with a bin to bin mass ratio of 1.5. This implies that each subsequent bin has 1.5 times as much mass as the previous one.

ABIFM parameterizes the heterogeneous ice nucleation rate coefficient J_het, as a function of the water activity criterion, Δa_w (Knopf and Alpert, 2013). Derived from the water activity at the ice melting temperature, $a_{\mathrm{w}}^{\mathrm{i}}$, and water activity at the freezing temperature, a_w, Δa_w is used in ABIFM to describe IMF as a function of ambient T and RH:

$$\begin{array}{}\text{(A9)}& \mathrm{\Delta }{a}_{\mathrm{w}}\left(T\right)=a_{\mathrm{w}}\left(T\right)-a_{\mathrm{w}}^{\mathrm{i}}\left(T\right).\end{array}$$

Notice that, for IMF from water droplets, a_w(T) = 1, while under subsaturated conditions, a_w(T)< 1. The equation for $a_{\mathrm{w}}^{\mathrm{i}}\left(T\right)$ is given in Koop and Zobrist (2009). To evaluate freezing at subsaturated conditions, i.e., at RH > 90 %, we conducted additional model simulations termed ABIFM*.

ABIFM and ABIFM^* derive J_het in units of cm⁻² s⁻¹ for mineral dust, organic particles, and SSA particles in the following way (Alpert and Knopf, 2016; China et al., 2017; Alpert et al., 2022), respectively:

$$\begin{array}{}\text{(A10)}& {\displaystyle }{\log }_{\mathrm{10}}\left(J_{\mathrm{het}}^{\mathrm{dust}}\left(\mathrm{\Delta }{a}_{\mathrm{w}}\left(T\right)\right)\right)=\mathrm{22.62}\mathrm{\Delta }{a}_{\mathrm{w}}\left(T\right)-\mathrm{1.35}\text{(A11)}& {\displaystyle }\begin{array}{rl}{\log }_{\mathrm{10}}\left(J_{\mathrm{het}}^{\mathrm{organic}}\left(\mathrm{\Delta }{a}_{\mathrm{w}}\left(T\right)\right)\right)& =\mathrm{66.90259}\mathrm{\Delta }{a}_{\mathrm{w}}\left(T\right)\\ & -\mathrm{13.40148}\end{array}\text{(A12)}& {\displaystyle }{\log }_{\mathrm{10}}\left(J_{\mathrm{het}}^{\mathrm{SSA}}\left(\mathrm{\Delta }{a}_{\mathrm{w}}\left(T\right)\right)\right)=\mathrm{26.6132}\mathrm{\Delta }{a}_{\mathrm{w}}\left(T\right)-\mathrm{3.9346}.\end{array}$$

In the CNT case, the number concentration of activatable INPs, $N_{\mathrm{INP}}^{\mathrm{CNT}\left(\mathrm{ABIFM}\right)}$, is given by

$$\begin{array}{}\text{(A13)}& N_{\mathrm{INP}}^{\mathrm{CNT}\left(\mathrm{ABIFM}\right)}\left(d\right)=N_{\mathrm{aer}}\left(d\right),\end{array}$$

where N_aer(d) is the number concentration of the applied aerosol particle type per bin.

In the AC-1D model, the CNT-based IMF parameterizations are implemented as follows:

$$\begin{array}{}\text{(A14)}& N_{\mathrm{INP}}^{\mathrm{Imm}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)=N_{\mathrm{INP}}^{\mathrm{CNT}\left(\mathrm{ABIFM}\right)}\left(z_i,t,d\right).\end{array}$$

The d array has the same characteristics as for the INAS case.

Appendix B

B1 INP array calculation for singular approaches

In singular approaches, the model employs a decoupled framework that separates the intrinsic INP activation spectrum from its environmental expression. This design ensures consistent initialization across simulations while allowing dynamic evolution based on atmospheric conditions.

B2 Temperature array construction

The temperature array is constructed to ensure finer resolution at lower temperatures, where ice nucleation is more sensitive. This array is generated using a geometric progression, with the temperature bins increasing exponentially. The freezing array begins at a minimum temperature T_min= −38 °C and the temperature increments are determined by an initial temperature step ΔT₀ (the ΔT between the first and second temperature bin edges) and an exponential factor dT_exp (the ratio of ΔT between consecutive bins).

The progression continues, but the array stops when the next temperature value is greater than or equal to the maximum temperature T_max = −5 °C, meaning the last value in the array may not be reached at exactly T_max. The temperature array is constructed as:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{T}_{\mathrm{array}}& =\left[T_{\min },T_{\min }+\mathrm{\Delta }{T}_{\mathrm{0}},T_{\min }+\mathrm{\Delta }{T}_{\mathrm{0}}\cdot \mathrm{d}{T}_{\exp }^{\left(\mathrm{1}\right)},T_{\min }\right.\\ & \left.+\mathrm{\Delta }{T}_{\mathrm{0}}\cdot \mathrm{d}{T}_{\exp }^{\left(\mathrm{2}\right)},\mathrm{\dots }\right]\mathrm{while}T<T_{\max }\end{array}\end{array}$$

B3 INP reservoir initialization

Using established parameterizations (Appendix A), the model calculates the cumulative INP concentration active at or below each temperature, then discretizes this into a differential activation spectrum N_INP(T) representing INPs that activate specifically within each temperature bin. This creates the baseline INP reservoir array: INP${}_{\mathrm{array}}\left(z,T,t\right)$.

B4 Dynamic evolution and filtering

During simulation, the INP reservoir evolves through entrainment of new INPs from external sources, permanent removal through activation when environmental temperature drops to or below a bin's activation temperature (conversion to ice crystals), and vertical redistribution via turbulent mixing (represented by mixing time scale τ_mix). Note that sedimentation is neglected for the aerosol INP reservoir.

At each timestep, the model applies a dynamic temperature mask based on the minimum column temperature to determine the potentially activatable INP concentration:

$$\begin{array}{}\text{(B2)}& {\mathrm{INP}}_{\mathrm{totarray}}\left(z,t\right)=\sum _{T\ge T_{\mathrm{cloud}\mathrm{top}}}{\mathrm{INP}}_{\mathrm{array}}\left(z,T,t\right)\end{array}$$

This approach ensures that only INPs with activation temperatures warmer than current conditions remain in the activatable reservoir, while consumed INPs are permanently removed. The decoupled framework provides consistent comparisons across simulations while capturing realistic INP dynamics.

Appendix C: Governing equations

Based on the following equations, the simplified 1D aerosol cloud model predicts the budgets for ice-nucleating particles (INPs) and ice particles (Knopf et al., 2023b):

$$\begin{array}{}\text{(C1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{INP}}^{\mathrm{Imm}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)}{\mathrm{d}t}}& =-S_{\mathrm{act}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)\\ & +S_{\mathrm{ent}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)\\ & +S_{\mathrm{mix}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)\end{array}\text{(C2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{i}}\left(z_i,t\right)}{\mathrm{d}t}}& ={\sum }_{k_{*}}^{\mathrm{Imm}}{S}_{\mathrm{act}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)-S_{{\mathrm{i}}_{\mathrm{sed}}}\left(z_i,t\right)\\ & +S_{{\mathrm{i}}_{\mathrm{mix}}}\left(z_i,t\right),\end{array}\end{array}$$

where N_INP presents the number concentration of activatable INPs and N_i denotes the number concentration of ice crystals. The z_i signifies the height at grid cell with index i counted from bottom layer (1, surface) to the top layer (m, PBL top). Notice that the superscript Imm refers to the specific type of IMF parameterization used (INN, INAS, ABIFM or ABIFM*), and $k_{*}^{\mathrm{Imm}}$ represents additional variable dimensions that depend on the chosen IMF parameterizations. For example, for application of INN, $k_{*}^{\mathrm{Imm}}$ is T_INP (reflecting the INP temperature arrays), for INAS, $k_{*}^{\mathrm{Imm}}$ is a combination of T_INP and d (reflecting a 2D array consisting of INP temperature and diameter arrays), while for ABIFM and ABIFM*, $k_{*}^{\mathrm{Imm}}=d$. Moreover, the model utilizes several variables to describe the budget terms. S_act, S_ent and S_mix represent the INP activation, cloud-top entrainment, and turbulent mixing budget terms, respectively. For ice particles, $S_{{\mathrm{i}}_{\mathrm{sed}}}$ and $S_{{\mathrm{i}}_{\mathrm{mix}}}$ denote the ice sedimentation and mixing budget terms. It is essential to note that ice crystals do not retain information about their associated INPs. Therefore, the activation budget is summed over the full T_INP array (in INN), the T_INP and d arrays (in INAS) and the d array (in ABIFM and ABIFM*).

In the model configuration employed in this study, we utilize a time splitting approach to calculate the budgets for ice-nucleating particles (INPs) and ice particles at every time step. The implicit solution for S_act is given as

$$\begin{array}{}\text{(C3)}& S_{\mathrm{act}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)={\displaystyle \frac{N_{\mathrm{INP}}^{\mathrm{Imm}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)}{\mathit{\delta }t+{\mathit{\tau }}_{\mathrm{act}}}},\end{array}$$

where the time step length denoted by δt, is set by default to 10 s. For CNT-based parameterizations (ABIFM, ABIFM*), τ_act is calculated as ${\mathit{\tau }}_{\mathrm{act}}=\frac{\mathrm{1}}{J_{\mathrm{het}}\left(T\left(z_i\right),a_{\mathrm{w}}\left(z_i\right)\right)\mathit{\pi }{d}^{\mathrm{2}}}$. For singular parameterizations (INN, INAS) in this study, INP activation is treated as instantaneous once temperature conditions are met. This is implemented by directly converting the number of activatable INPs (determined by the respective singular parameterization at the given temperature) to ice crystals within the model timestep, effectively bypassing the explicit τ_act-dependent formulation of S_act shown in Eq. (C3) which is used for CNT schemes. The implicit Euler method (Hoffman and Frankel, 2001), is used to ensure numerical stability, where δt≪τ_act  (relevant for CNT schemes), the change in concentration is proportional to $\frac{-\mathit{\delta }t}{\mathit{\delta }t+{\mathit{\tau }}_{\mathrm{act}}}$ , which becomes increasingly similar to the explicit solution, where the change in concentration is proportional to $\frac{-\mathit{\delta }t}{{\mathit{\tau }}_{\mathrm{act}}}$ . However, this implicit Euler method ensures that the loss never exceeds the initial value.

To ensure numerical stability and avoid violating the Courant-Friedrichs-Lewy (CFL) condition, S_ent is computed implicitly:

$$\begin{array}{}\text{(C4)}& \begin{array}{rl}& S_{\mathrm{ent}}\left(z_m,t,k_{*}^{\mathrm{Imm}}\right)=\\ & {\displaystyle \frac{N_{\mathrm{INP},\mathrm{FT}}^{\mathrm{Imm}}\left(k_{*}^{\mathrm{Imm}}\right)-N_{\mathrm{INP}}^{\mathrm{Imm}}\left(z_m,t,k_{*}^{\mathrm{Imm}}\right)}{\mathit{\delta }t+\frac{\mathit{\delta }z}{w_{\mathrm{e}}}}}.\end{array}\end{array}$$

This term represents the source of new aerosol particles and activatable INPs entering the boundary layer from the free troposphere. The number of particles entrained per time step is determined by the concentration gradient between the free troposphere and the cloud top, scaled by the entrainment rate, where w_e is the entrainment rate and $N_{\mathrm{INP},\mathrm{FT}}^{\mathrm{Imm}}$ refers to the free-troposphere immersion freezing INP concentration, which is equal to the initial domain INP size distribution. Notice that z_m represents the height of the domain top. Following Fridlind et al. (2012), w_e is assumed to be 0.1 cm s⁻¹ based on LES results.

The mixing terms for INPs and ice particles are calculated as the following:

$$\begin{array}{}\text{(C5)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{mix}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)& ={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{mix}}}}\left({\stackrel{\mathrm{‾}}{N}}_{\mathrm{INP}}^{\mathrm{Imm}}\left(t,k_{*}^{\mathrm{Imm}}\right)\right.\\ & \left.-N_{\mathrm{INP}}^{\mathrm{Imm}}\left(z_i,t,k_{*}^{\mathrm{Imm}}\right)\right)\end{array}\text{(C6)}& {\displaystyle }{S}_{{\mathrm{i}}_{\mathrm{mix}}}\left(z_i,t\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{mix}}}}\left({\stackrel{\mathrm{‾}}{N}}_i\left(t\right)-N_{\mathrm{i}}\left(z_i,t\right)\right),\end{array}$$

where τ_mix is the PBL mixing time scale. Physically, τ_mix represents the large scale vertical turbulent mixing time scale. This parameter governs the rate at which vertical gradients in scalar concentrations are homogenized toward the boundary layer mean, redistributing entrained particles throughout the cloud layer. ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{INP}}^{\mathrm{Imm}}$ is the PBL mean activatable INP number concentration averaged over the whole domain, and ${\stackrel{\mathrm{‾}}{N}}_i\left(t\right)$ denotes the PBL mean N_i.

Lastly, the sedimentation term for ice crystals is determined as follows:

$$\begin{array}{}\text{(C7)}& S_{{\mathrm{i}}_{\mathrm{sed}}}\left(z_i,t\right)={\displaystyle \frac{v_{\mathrm{f}}}{\mathit{\delta }z}}\left(N_{\mathrm{i}}\left(z_{i+\mathrm{1}},t\right)-N_{\mathrm{i}}\left(z_i,t\right)\right),\end{array}$$

where v_f is the number-weighted ice sedimentation rate (ice crystal fall velocity), maintained as a constant value based on LES estimates of circa 30 cm s⁻¹ (Fridlind et al., 2012). To prevent potential numerical instability, we always set the v_f⋅δt to be smaller than δz, thus leading to CFL < 1. Also, note that in Eq. (C7), the i ranges from 1 to m−1 (index of layer below the PBL top) instead of from 1 to m.

As shown above, we account for the loss and gain of INPs and ice crystals, thus, treating the INPs and ice crystals prognostically.

Appendix D: Nomenclature

aw [–]	Water activity
$a_{\mathrm{w}}^{\mathrm{i}}$ [–]	Ice melting point as a function of water activity
d [µm]	INP diameter array
Δaw [–]	Water activity criterion
δt [s]	Time step length
δz [m]	Model grid cell thickness (vertical resolution)
Jhet [cm−2 s−1]	Heterogeneous ice nucleation rate coefficient
$k_{*}^{\mathrm{Imm}}$ [°C], [°C µm], or [µm]	Additional variable dimensions depending on the applied parameterization
Naer [cm−3]	Aerosol number concentration
Naer>0.5 µm [cm−3]	Aerosol number concentration for particles larger than 0.5 µm
Ni [L−1]	Ice crystal number concentration
${\stackrel{\mathrm{‾}}{N}}_{\mathrm{i}}\left(t\right)$ [L−1]	PBL mean ice crystal number concentration
NINP [L−1]	Activatable INP number concentration
${\stackrel{\mathrm{‾}}{N}}_{\mathrm{INP}}$ [L−1]	PBL mean activatable INP number concentration
$N_{\mathrm{INP}}^{\mathrm{Imm}}$ [L−1]	Activatable INP number concentration for given immersion freezing parameterization
NINP,FT [L−1]	Free-troposphere activatable INP concentration
NINP(T) [L−1]	The number concentration of INPs that activate specifically within temperature bin
ns [cm−2]	Ice nucleation active sites density
RH [%]	Relative humidity expressed in percent
Sact [m−3 s−1]	INP activation budget term
$S_{\mathrm{act}}^{\mathrm{Imm}}$ [m−3]	INP activation budget term calculated based on the selected immersion
	freezing parameterizations
$S_{\mathrm{act}}^{\mathrm{sing}\left(\mathrm{INN}\right)}$ [m−3]	INP activation budget term calculated based on the INN
$S_{\mathrm{act}}^{\mathrm{sing}\left(\mathrm{INAS}\right)}$ [m−3]	INP activation budget term calculated based on the INAS
$S_{\mathrm{act}}^{\mathrm{CNT}\left(\mathrm{ABIFM}\right)}$ [m−3]	INP activation budget term calculated based on the ABIFM
Sent [m−3 s−1]	Cloud-top entrainment budget term
Sice [–]	Supersaturation with respect to ice
$S_{{\mathrm{i}}_{\mathrm{mix}}}$ [m−3 s−1]	Ice mixing budget term
$S_{{\mathrm{i}}_{\mathrm{sed}}}$ [m−3 s−1]	Ice sedimentation budget term
Smix [m−3 s−1]	Turbulent mixing budget term
t [s]	Model time step
τact [s]	Activation time scale
TINP [°C]	INP temperature array
τmix [s]	PBL mixing time scale
vf [m s−1]	ice sedimentation rate
we [cm s−1]	Entrainment rate
zi [m]	Height at grid cell index i
zm [m]	Height at domain top (grid cell index m)
INParray [m−3]	Three-dimensional array that stores calculated INP values.
INPtotarray [m−3]	The sum of all INPs with activation temperatures warmer than the current cloud
	conditions that remain available for ice nucleation.

Appendix E: List of Abbreviations

2D-C	Two-dimensional cloud optical array
	probes
2D-P	Two-dimensional precipitation optical
	array probes
ABIFM	The water-activity-based immersion
	freezing model
ABIFM*	The water-activity-based immersion
	freezing model which enables immersion
	freezing commencing under subsaturated
	conditions
CCR	Cloud cooling rate
CFL	Courant–Friedrichs–Lewy condition
CNT	Classical nucleation theory
CPI	Cloud Particle Imager
DHARMA	Distributed Hydrodynamic Aerosol and
	Radiative Modeling Application
ICEALOT	International Chemistry Experiment in the
	Arctic Lower Troposphere
INAS	Ice nucleation active sites-based
	parameterization
INN	Ice nucleation number-based
	parameterization
INP	Ice nucleating particle
ISDAC	Indirect and Semi-Direct Aerosol Campaign
LES	Large-eddy simulation
MMCR	Millimeter Wavelength Cloud Radar
M-PACE	Mixed-Phase Arctic Cloud Experiment
MPC	Mixed-phase cloud
PBL	Planetary boundary layer
PSD	Particle size distribution
SHEBA	The Surface Heat Budget of the Arctic
	campaign
SSA	Sea spray aerosol