The chance of freezing – a conceptional study to parameterize temperature-dependent freezing by including randomness of ice-nucleating particle concentrations

Ice-nucleating particle concentrations (INPCs) can spread over several orders of magnitude at any given temperature. However, this variability is rarely accounted for in heterogeneous ice-nucleation parameterizations. In this paper, we present an approach to incorporate the random variation in the INPC into the parameterization of immersion freezing and analyze this novel concept with various sensitivity tests. In the new scheme, the INPC is drawn from a relative frequency distribution of cumulative INPCs. At each temperature, this distribution describing the INPCs is expressed as a lognormal frequency distribution. The new parameterization scheme does not require aerosol information from the driving model to represent the heterogeneity of INPCs. The scheme's performance is tested in a large-eddy simulation of a relatively warm Arctic mixed-phase stratocumulus. We find that it leads to reasonable ice masses in the cloud, especially when compared to immersion freezing schemes that yield one fixed INPC per temperature and lead to almost no ice production in the simulated cloud. The scheme is sensitive to the median of the frequency distribution and highly sensitive to the standard deviation of the distribution, as well as to the frequency of drawing a new INPC and the resolution of the model. Generally, a higher probability of drawing large INPCs leads to substantially more ice in the simulated cloud. We expose inherent challenges to introducing such a parameterization and explore possible solutions and potential developments.

1 Introduction

Clouds play an important role in Earth's energy balance by reflecting incoming sunlight and interacting with infrared radiation. The cloud phase influences the extent of a cloud's radiative effect, but more importantly, it determines whether the cloud has a warming or cooling effect. According to Matus and L'Ecuyer (2017), liquid clouds have a global net radiative effect at the top of the atmosphere (TOA) of −11.8 W m⁻², whereas ice clouds exert a warming of 3.5 W m⁻², and mixed-phase clouds cause a net cooling effect of −3.4 W m⁻².

The ice crystals in mixed-phase clouds usually originate from heterogeneous ice nucleation, where ice-nucleating particles (INPs) are necessary to trigger the phase transition. Ice nucleates heterogeneously in the atmosphere at temperatures between 0 and approximately −38 ^∘C. At temperatures below −38 ^∘C, homogeneous ice nucleation occurs spontaneously, and freezing can happen without INPs. Heterogeneous ice nucleation can occur in different so-called modes, namely immersion freezing, contact freezing, deposition nucleation, and condensation freezing (see Vali et al., 2015, for definitions). We focus on immersion freezing, where an INP is immersed in a supercooled cloud droplet and initiates freezing at a specific temperature. Different aerosol types can act as INPs, including mineral dusts and biological and combustion particles. The probability that an aerosol initiates ice nucleation increases with decreasing temperature. Measurements of ambient INP concentrations (INPCs; m⁻³) show that at a given temperature, INPCs can vary by several orders of magnitude over time and space. Examples of spatial INPC variability in different marine locations, from approximately 10⁻¹ to 10³ m⁻³ at −15 ^∘C can be found in Welti et al. (2020). Temporal variability, for example, the annual cycle of INPCs in the Arctic, was reported to be up to 3 orders of magnitude (Wex et al., 2019). However, also within smaller time periods down to single days, INPCs can fluctuate by up to 3 orders of magnitude (Bigg, 1961). It has been found that an INPC at a specific temperature is lognormally distributed. Lognormal frequency distributions of INPC occurrence have been measured at several locations within different environments (e.g., Isaac and Douglas, 1971; Bertrand et al., 1973; Radke et al., 1976; Flyger and Heidam, 1978; Conen et al., 2017; Welti et al., 2018; Hartmann et al., 2019; Schrod et al., 2020; Li et al., 2022). Ott (1990) showed that when an aerosol concentration is observed to have a lognormally distributed occurrence, this likely results because the aerosol was subject to a series of random dilutions subsequent to emission.

Several parameterizations that simulate cloud droplet freezing exist, and they are based on different physical variables. Burrows et al. (2022) distinguish between aerosol-aware and aerosol-unaware parameterizations. One example of the latter is the formulation by Fletcher (1962, hereafter F62), which is a scheme that is based on the observation that the average INP concentration increases exponentially with decreasing temperature for $-\mathrm{10}>T>-\mathrm{30}$ ^∘C. F62 only requires temperature information to calculate INPCs. The scheme by Niemand et al. (2012, hereafter N12) is an example of an aerosol-aware parameterization; it describes immersion freezing based on the active site density of desert dust aerosols observed in laboratory measurements. N12 is valid for $-\mathrm{12}>T>-\mathrm{36}$ ^∘C and requires temperature, the number of dust aerosols, and the average INP surface as input to determine INPCs. Another example of an aerosol-aware scheme is the parameterization by Phillips et al. (2008, hereafter P08), which represents immersion, contact, and deposition freezing on dust or metal aerosol, organic carbon, and biological INPs. P08 uses temperature, water vapor saturation with respect to ice, and aerosol concentrations of the four aerosol species to predict INPCs. P08 is valid for $\mathrm{0}>T>-\mathrm{70}$ ^∘C (or lower). The parameterization of Khvorostyanov and Curry (2000, hereafter K00) is based on the classical nucleation theory (CNT), where each substance is assigned a characteristic contact angle between a nucleating ice cluster and the particle surface. The smaller the contact angle, the more ice-active a substance. K00 utilizes temperature, water vapor saturation with respect to ice, particle radius, and one contact angle per substance to calculate ice-nucleation rates. One method of including the particle-to-particle heterogeneity of INPs in parameterizations is to apply a distribution on the contact angle parameter in CNT to calculate the ice-activity of the INPs (e.g., Marcolli et al., 2007, or Wang et al., 2014). These schemes require temperature and aerosol radius as input to the parameterization and use different values for the mean and standard deviation of the contact angle distribution for different aerosol species.

There are three drawbacks that ice-nucleation parameterizations can have: (i) they can be computationally complex, i.e., require detailed input from models regarding the aerosol type, size, and/or number concentration. (ii) They can be limited to a specific temperature range. (iii) Most of them fail to reproduce INPC variability; i.e., for one set of environmental conditions (temperature, humidity, aerosol type, and concentration) they yield a single fixed INPC value.

The lack of INPC variability in simulations compared to atmospheric observations emerges from non-represented INP types and sources, the large size of model grid boxes, and the use of bulk aerosol concentrations (e.g., dust, soot, and biological elements) as input variables in parameterizations. Even if a model includes information on, for example, the detailed size distribution of dust aerosol, it will not represent all the dust INP types (different minerals) and their variability in the atmosphere. To circumvent these drawbacks, we developed a parameterization of immersion freezing (F23) that simulates observed INPCs and their variability while only using temperature as input variable. F23 is valid for the entire temperature range of immersion freezing ($\mathrm{0}>T>-\mathrm{38}$ ^∘C). The INPCs returned by F23 are drawn randomly from a lognormal distribution of INPCs for each temperature, thereby capturing the natural INP variability, without requiring information about the present ice-active aerosol. The random drawing allows the representation of the INP population by a distribution of INPCs, instead of, for example, one single value (as in F62). F23's lognormal INPC distributions are loosely based upon INPCs observed in the maritime boundary layer (see Sect. 2). The main goal of our study was to test the concept of including INPC variability in a parameterization scheme. We compare the behavior of random drawing from an INPC distribution instead of using fixed values only depending on temperature, for example.

We test the new F23 parameterization in the large-eddy simulation (LES) model MIMICA (MISU MIT Cloud and Aerosol model; Savre et al., 2014) to investigate how it performs at producing ice in a simulated cloud in comparison to two conventional parameterizations (diagnostic ice crystal number concentration and F62). We analyze the sensitivity of the scheme to its characteristics and implementation details. The simulations are initialized based on in situ observations during the ASCOS (Arctic Summer Cloud Ocean Study) Arctic campaign and represent a mixed-phase stratocumulus cloud with in-cloud temperatures between approximately −7 and −10 ^∘C (see Sect. 2.2).

2 Method

2.1 Formulation of the F23 parameterization

The basic idea for the new F23 parameterization is that INPCs have large variability, even on short timescales (Bigg, 1961), which we represent by drawing a new INPC from the given distribution at a certain frequency – if immersion freezing occurs (i.e., T≤0 ^∘C and cloud droplets are present). Drawing a new INPC mimics the evolution of the ice-active aerosol population, for example, by replenishment or time-dependent freezing at each grid point. The latter means that drawing at each time step does not necessarily reflect a change in aerosol population from one time step to another but a change in the present INP population due to the activation of more INPs with time. The main goal of this study is to test the concept of randomly drawing from a distribution of INPCs instead of applying a deterministic value (e.g., the median of observed INPCs), as is done in other freezing parameterization schemes.

We use the conceptional distribution of INPCs shown in Fig. 1, derived from the extensive data sets of a long-term time series of INPC temperature spectra measured in the maritime boundary layer at Cabo Verde (Welti et al., 2018), and widely dispersed ship-based observations (Welti et al., 2020). On average, the INPC field observations show striking consistency in terms of the temperature spectra's shape and variability over time and location. The same consistencies exist when compared with other long-term (e.g., Schrod et al., 2020) or composite INPC data sets (e.g., Petters and Wright, 2015).

Figure 1 Relative frequency distribution (RFD) spectra for INPCs as a function of temperature. Median INPCs are marked by the red line. The loglinear parameterization by Fletcher (1962, hereafter F62) is shown within its validity range (solid black) and extrapolated to higher temperatures (dotted black) for comparison. The inset in the lower-left corner shows the RFD at $T=-\mathrm{16}$ ^∘C (dashed green line).

Based on the aforementioned data sets, we derived a function for the temperature-dependent, lognormally distributed INPC frequency as follows:

$$\begin{array}{}\text{(1)}& D\left(\mathit{\mu },{\mathit{\sigma }}^{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\cdot \mathit{\sigma }}}\exp \left(-{\displaystyle \frac{\left[\ln \right(a\cdot \mathrm{INPC})-\mathit{\mu }(T){]}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

with a=1 m³, INPC in per meters cubed (m⁻³), μ(T) the temperature-dependent mean, and σ² the variance of the lognormal distribution (not the INPC itself). For the marine data sets, we find σ=1.37 and $\mathit{\mu }\left(T\right)=\ln (-(b\cdot T{)}^{\mathrm{9}}\times {\mathrm{10}}^{-\mathrm{9}})$, with $b=\mathrm{1}/$(1 ^∘C) and T given in degrees Celsius (^∘C). By normalizing the distribution, a relative frequency distribution (RFD) as a function of temperature is obtained (Fig. 1). This normalization is necessary, due to the discrete INPC field. Whenever an immersion freezing event occurs in the model (i.e., T≤0 ^∘C and cloud droplets are present), an INPC value is drawn randomly from the RFD. That is, for two grid points with freezing events at the same temperature, the drawn INPC can differ by several orders of magnitude. For a large number of grid points, the relative frequency of the drawn INPC will follow a lognormal function, with the median INPC (red line in Fig. 1) having the highest probability. For example, if all grid points were at a temperature of −16 ^∘C, then the frequency of the drawn INPC would follow the distribution shown as a green curve in the lower-left inset of Fig. 1.

F23 assumes that all INPs are immersed in cloud droplets. To take into account the dynamic evolution of ice formation, the ice crystal number concentration (N_i) at a grid point is subtracted from the drawn INPC. This returns the number of newly frozen cloud droplets in a time step, i.e., hydrometeors moving from the “cloud droplet” class to the “ice crystal” class. Subtracting N_i is one approach to solving the need for time discretization when implementing a deterministic, time-independent scheme. Note that if other frozen hydrometeor species (snow and graupel) are represented in the model (not implemented in our setup), then the sum of their number concentrations should be subtracted from the drawn INPCs. No negative tendencies ($\mathrm{INPC}-N_i<\mathrm{0}$) are allowed in the scheme, since frozen cloud droplets will not melt due to a decrease in the INPCs. The change in the respective mixing ratios is calculated by multiplying the change in the number of frozen droplets by the average cloud droplet mass as follows:

$$\begin{array}{}\text{(2a)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{N}_{\mathrm{i}}=-\mathrm{\Delta }{N}_{\mathrm{c}}=\max \left(\right[\mathrm{INPC}-N_{\mathrm{i}}],\mathrm{0})\text{(2b)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{Q}_{\mathrm{i}}=-\mathrm{\Delta }{Q}_{\mathrm{c}}=\mathrm{\Delta }{N}_{\mathrm{i}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{c}}}}{\stackrel{\mathrm{‾}}{N_{\mathrm{c}}}}}.\end{array}$$

Here, N_i and N_c are the ice crystal and cloud droplet number concentrations per meters cubed (m⁻³). Q_i and Q_c are the ice crystal and cloud droplet mixing ratios (kg m⁻³), respectively, with the mean cloud droplet mixing ratio and number concentration denoted by the bar. The calculations are repeated at each time step, which means that INPs are only partially depleted through the subtraction of N_i. We interpret this as a way to imitate the time-dependent freezing and INP recycling, which has been shown to be crucial for realistic cloud development in LES for Arctic mixed-phase clouds (e.g., Solomon et al., 2015).

The implemented INPC RFD field is discretized into bins of INPCs and temperature. The INPCs differ approximately by a factor of 2 (or Δlog ₁₀(INPC)≈log ₁₀(2)), while the temperature bins have a size of 1 ^∘C. For this reason, we define temperature to the nearest degree when drawing from the INPC RFD, if not stated otherwise. Using the F23 immersion freezing scheme in MIMICA has approximately the same computational expense for the total simulation time as any other interactive ice-nucleation parameterization (for example, F62).

2.1.1 Example of parameterized INP concentrations

To illustrate the INPC distribution from the F23 parameterization, let us assume that we have a uniform −16 ^∘C cloud spreading horizontally over the entire model domain consisting of 1000 grid points. The INPC RFD at all cloud grid points is represented by a lognormal distribution curve (inset in Fig. 1). This results in an INPC of 68.7 m⁻³ being drawn with the highest probability, since this is the median INPC at −16 ^∘C, where $\mathrm{Med}\left[\mathrm{INPC}\right(T=-\mathrm{16}{}^{\circ }\mathrm{C}\left)\right]=\exp \left(\mathit{\mu }\right){\mathrm{m}}^{-\mathrm{3}}=-T^{\mathrm{9}}\times {\mathrm{10}}^{-\mathrm{9}}{\mathrm{m}}^{-\mathrm{3}}\approx \mathrm{68.7}$ m⁻³. In the model, 20.2 % of the grid points will draw the median INPC; that is, 202 of the 1000 cloud grid points will have the median INPC of 54.3 m⁻³ (differing from the theoretical value of 68.7 m⁻³ because of discrete INPC binning). The range of INPC bins with a relative frequency > 0.1 % covers INPCs of 0.8–7543.1 m⁻³, which means that INPCs ≤0.8 or ≥7543.1 m⁻³ will rarely be drawn for the example cloud. If no ice was present previously, then the number of cloud droplets frozen at this time step equals the INPC (Eq. 2a), which is limited by the total number of cloud droplets present.

2.1.2 Representing the RFD

Since the parameterization draws values from a distribution, it needs to be ensured that there are enough random draws in order to represent the distribution well. INPCs in the model should vary according to the distribution, but different model runs should also be reproducible. To investigate how many draws are necessary to represent the INPC distribution, we conducted several drawing tests (drawing, e.g., 50 times from the distribution vs. drawing 100 times) at −16 ^∘C and compared the relative frequencies of the drawn values to the theoretical values by calculating the root mean square error (RMSE) between the drawn and theoretical distributions. Comparing the results of the different drawing tests suggests that 300 random draws lead to a reproducible prediction of the RFD (see Table 1 and Fig. 2); the RMSE converges for ≥ 300, with a constant first derivative (linear slope) of the connecting lines, and the standard deviation decreases only slightly for >300 draws.

Table 1 Average and standard deviation of RMSE for different draw amounts from the INPC RFD at −16 ^∘C compared to theoretical RFD values. Each drawing test was performed 100 times, and the table contains the average and standard deviation of the RMSE over these 100 times.

Figure 2 Average and standard deviation of RMSE for the drawing tests in Table 1 are decreasing for an increasing number of draws at −16 ^∘C.

The domain used in the simulations has 96×96 grid points in the horizontal, leading to 9216 grid points in each layer. The simulated cloud is stratocumulus (see Sect. 2.2), where the temperature field is horizontally uniform. Assuming that at least one layer of grid points in the model has the same temperature due to the uniform structure of the stratus cloud, this translates into a minimum of 9216 grid points with the same temperature. This implies a minimum of 9216 draws from the RFD at one temperature, which is substantially more than the minimum number of draws determined to represent the distribution well (300). Hence, we confirm that a representative INPC distribution is being drawn from the RFD. This has also been verified in MIMICA simulations (not shown).

2.2 Simulation setup

We use the well-established large-eddy simulation (LES) model MIMICA (MISU MIT Cloud and Aerosol model; Savre et al., 2014). For more information on the model, also see Appendix A or, e.g., Savre and Ekman (2015a, b) or Sotiropoulou et al. (2020). The simulated case is based on a mixed-phase Arctic stratocumulus. The case study stratocumulus cloud was observed between 30 and 31 August 2008 during the ship-based ASCOS campaign (Arctic Summer Cloud Ocean Study; Tjernström et al., 2014). At that time, research vessel (R/V) Oden was drifting with an ice floe located at approximately 87^∘ N. The atmospheric conditions were characterized by a high-pressure system, with large-scale subsidence in the free troposphere (for details, see Tjernström et al., 2012). This case has been used to study other microphysical cloud properties like dissipation (Loewe et al., 2017), the influence of cloud condensation nuclei (CCN) hygroscopicity on cloud properties (Christiansen et al., 2020), secondary ice production (Sotiropoulou et al., 2021), and the sustenance (Bulatovic et al., 2021) of an Arctic mixed-phase cloud, as well as in a model intercomparison study (Stevens et al., 2018). We selected this case because it is an established case and because there are large uncertainties about the nature and concentration of INPs in the Arctic, which poses a challenge to modeling mixed-phase clouds in this region. Using other immersion freezing parameterizations in MIMICA, for example, an active site scheme, following Ickes et al. (2017), to simulate this case requires unrealistically active INPs in order to form ice. We focus on immersion freezing in this study because liquid-dependent ice nucleation is dominant in Arctic stratiform clouds (de Boer et al., 2011). Contact freezing can be neglected, since very little interstitial aerosol was present in the case. Furthermore, these aerosols would need to collide with the supercooled cloud droplets in a very stable non-turbulent cloud, making contact freezing unlikely to occur. Deposition nucleation has the same limitation when it comes to interstitial aerosol and can additionally be neglected because of the temperature range which is not favorable for deposition ice nucleation. Secondary ice formation is not explicitly taken into account in this study, since we focus on primary ice formation. Simulations of the case with a diagnostic ice crystal number concentration (DIAG) instead of an interactive representation of ice nucleation assume a minimum ice crystal number concentration of 200 m⁻³ at grid points with a temperature below 0 ^∘C, where there are sufficient cloud droplets. This means that at any given time step in the cloud, if N_i falls below 200 m⁻³ and T<0 ^∘C, then ice crystals are produced to retain N_i=200 m⁻³, no matter the exact temperature below 0 ^∘C. This approach is unspecific to the ice formation mechanisms. It combines immersion and contact freezing (due to the requirement of cloud droplets) and secondary ice processes.

The MIMICA simulations are initialized with the profiles of thermodynamic variables (e.g., potential temperature and pressure) and liquid cloud water measured at approximately 06:00 UTC on 31 August 2008. The initial ice or liquid potential temperature profiles are randomly perturbed in order for convection to develop more quickly. Consequently, any two simulations will yield different results, even if all parameters are held constant. A cloud layer was present between ca. 550 and 900 m a.g.l. (above ground level), capped by a temperature and humidity inversion, and decoupled from the surface (see Sotiropoulou et al., 2021, for profile details). The temperature within the cloud ranged from approximately −7 to −10 ^∘C. The simulation setup follows Sotiropoulou et al. (2021). The domain covers a $\mathrm{96}\times \mathrm{96}\times \mathrm{128}$ grid, with a constant horizontal spacing of $\mathrm{d}x=\mathrm{d}y=\mathrm{62.5}$ m (6 km × 6 km horizontal domain size). The vertical spacing is 7.5 m near the ground and in the cloud layer; between the surface and the cloud, it changes sinusoidally and reaches a maximum dz of 25 m, with a 1.7 km total vertical domain size. The time step is dynamic in order to satisfy the Courant–Friedrichs–Lewy (CFL) condition for the leapfrog time-integration method and ranges from ≈1–3 s to prevent numerical instabilities within the model. Our simulations cover 12 h, with the first 2 h utilized as a spinup period and subsequently omitted from the results. A large-scale steady state is maintained throughout the model runs, and the cloudy layer is present in the initial state of the simulations. However, the initial cloud is liquid, and cloud ice is only formed from the first model time step.

We excluded the hydrometeor categories of snow and graupel from all simulations, since it is known that MIMICA produces rather large volumes of graupel for this case (Stevens et al., 2018), which dominates both the ice water path (IWP) and the number concentration of frozen hydrometeors. Because we are primarily interested in ice formation in our study, having ice crystals as the only frozen hydrometeors simplifies the analysis. For the same reason, snow and graupel were excluded previously from Arctic stratocumulus simulations conducted with MIMICA (Savre and Ekman, 2015b). Excluding snow and graupel means that no cold collection processes are active in our simulations; however, ice crystals can still grow by deposition, be transported, and precipitate. Aerosol is not represented prognostically in the model setup.

Observed liquid water path (LWP) and IWP were derived from measurements by a micrometer radiometer and millimeter cloud radar, respectively (Tjernström et al., 2012). Uncertainties for LWP are 25 g m⁻², while for IWP they are approximately a factor of 2 (Shupe et al., 2013).

3 Results and discussion

We first compare the baseline version of the new F23 parameterization (F23; all parameters are set to the values described in Sect. 2.1) to the standard setup (DIAG; diagnostic ice crystal number concentration), available observations, and F62.

Figure 3 Domain-averaged (a) LWP, (b) IWP, (c) ICNB, and (d) NIRS for 10 simulations of DIAG (red; median with min/max envelope) and with the F23 parameterization (cyan; median with min/max envelope). In both cases, one representative simulation is plotted as a dotted line to be included in Figs. 4, D1, 6, 8, 10, 12, and C1. The interquartile range of observations for LWP and IWP is indicated by the shaded gray areas. Median observations are marked by the dotted white lines. One simulation with F62 is shown in blue. The significance was assessed with a two-sided Kolmogorov–Smirnov test at the 95 % level. Significant differences in the IWPs of F62 to F23 are indicated with an asterisk.

Domain-averaged liquid water path (LWP; g m⁻²), ice water path (IWP; g m⁻²), ice crystal number burden (ICNB; m⁻²), and net infrared radiation at the surface (NIRS; W m⁻²) for DIAG, F23, and F62 are shown, together with the median and the upper and lower quartiles of observational values for LWP and IWP, in Fig. 3. The model was run 10 times using DIAG or F23, respectively. The different “ensemble members” are initialized by randomly perturbed profiles and are represented by their median, maximum, and minimum values in Fig. 3.

Figure 4 Domain-averaged profiles for the simulation time after 2 h spinup time. (a, c) N_i with contours of temperature. (b, d) Ice mixing ratio (Q_i). The upper row (a, b) is the DIAG comparison run (see dotted red lines in Fig. 3), and the lower row (c, d) is the F23 comparison run (see dotted cyan lines in Fig. 3). The dashed lines show the cloud's top and bottom.

3.1 Comparison of F23 to DIAG (diagnostic ice crystal number concentration)

The F23 and DIAG simulations differ only in how ice crystal formation is parameterized. Other related processes like deposition and sublimation are treated in the same manner, and the general setup follows Sect. 2.2. For F23, all parameters are set to the values described in Sect. 2.1, while DIAG ensures a minimum N_i of 200 m⁻³, where the temperature is below 0 ^∘C, and there are sufficient cloud droplets.

LWP is substantially larger for F23 than for DIAG (Fig. 3a). IWP is larger for DIAG than for F23 (Fig. 3b) but does not compensate for the differences in LWP (the total mass of ice and liquid is still larger for F23). The large differences in LWP between F23 (INPC distribution) and DIAG (constant N_i) are due to enhanced freezing in DIAG, which leads to the increased evaporation of liquid droplets and deposition onto ice crystals via the Wegener–Bergeron–Findeisen process (WBF; see Fig. B1).

To analyze the ice crystal number, Fig. 3c shows ICNB, the vertically integrated N_i (Fig. C1a). For DIAG, ICNB increases slightly from 1.88×10⁵ to 1.96×10⁵ m⁻² over 10 h, while it increases dramatically throughout the F23 simulations from 3.2×10⁴ to 1.0×10⁵ m⁻². Several primary processes cause the increase in the ice mass and number, namely (i) the vertical extent of the liquid cloud increases throughout the simulation (see dashed lines in Fig. 4). (ii) The temperature within the cloud decreases throughout the simulation (see white contour lines in Fig. 4a and c). (iii) For F23, new ice crystals form whenever the drawn INPC is larger than the current N_i. Only process (i) is relevant for DIAG, since its requirements for ice formation are T<0 ^∘C and the abundance of cloud droplets. For F23, all three processes are relevant, since it requires cloud droplets (i), the INPC RFD depends on the temperature (ii), and process (iii) is an inherent characteristic.

Figure 3d shows the time series of NIRS (i.e., incoming minus outgoing infrared radiation at the surface). NIRS is larger for F23 compared to DIAG. This can be expected, since LWP is larger for F23. The evolution over time is similar in both cases, with NIRS first increasing and then decreasing from hours 4–5. This pattern emerges from the combination of first increasing, then steady LWP (for F23, the LWP decreases slightly between hours 5 and 8, which causes a larger decrease in NIRS), and continuously decreasing temperature (see white lines in Fig. 4a and c), due to the radiative cooling of the cloud.

The profile of N_i for DIAG in Fig. 4a illustrates the principle of prescribing N_i diagnostically. For DIAG, within the cloud and ca. 100 m below it, N_i has the prescribed minimum value of 200 m⁻³. Down to the surface, the concentration is lower but constant, which is caused by steady sedimentation of ice crystals from above combined with sublimation from ice crystals (see Fig. B1b). Q_i has the maximum values at the cloud bottom (Fig. 4b), which is commonly observed in Arctic mixed-phase clouds (Shupe et al., 2008). Since no cold-phase collection processes are active in these simulations, the mass of single ice crystals can only grow by deposition. The total mass of ice is, however, also affected by transport processes like sedimentation. The vertical distribution of N_i and Q_i in F23 is similar to DIAG (Fig. 4c and d). However, both Q_i and N_i increase by larger rates throughout the simulation time for F23, which can also be seen in the IWP and ICNB values (Fig. 3b and c). Towards the end of the simulation time, F23 ice values are about half of DIAG. Even though ice nucleation depends on the temperature in F23, N_i is quite homogeneously distributed throughout the cloud (Fig. 4c). We explain this by the uniform stratification of the cloud and the subtraction of N_i from INPC (Eq. 2a) in the scheme, which homogenizes N_i vertically over time.

3.2 Comparison to observations

LWP as calculated in MIMICA is substantially higher than the 75th percentile of the observed values, irrespective of the treatment of ice nucleation (Fig. 3a). Stevens et al. (2018) found this already in their model intercomparison study, where MIMICA was amongst the models with the highest LWP. The LWP exceeded the observations for simulations similar to ours, with no prognostic aerosol treatment but a simplified treatment of aerosol activation (note, however, that the MIMICA simulations in Stevens et al., 2018, included snow and graupel). MIMICA simulations, in which aerosol is modeled prognostically, yielded results closer to the observational range in the study by Stevens et al. (2018). The IWP for DIAG lies within the interquartile range of the observations (Fig. 3b), which was found in Stevens et al. (2018) as well. The minimum N_i in DIAG was chosen in order to yield ice masses in the cloud similar to the observations. For the F23 simulations, IWP is lower than the 25th percentile of the observations. F23, however, only represents immersion freezing, leading to a lower ice mass. We can expect that multiplication processes were present in the observed cloud, which would have led to a higher IWP compared to model simulations that do not represent SIP. In fact, Sotiropoulou et al. (2021) simulated the same case, using MIMICA and including a description of ice multiplication from breakup during ice–ice collisions. Their results show an increase in IWP by a factor of 2 to 3 when breakup is activated in the model. Adjusting our simulated IWP accordingly would result in values corresponding to the observed IWP for F23. LWP is decreased by approximately 25–35 g m⁻² in simulations with ice multiplication in Sotiropoulou et al. (2021). This would bring our modeled LWP values closer to the observations, but LWP would still be higher than observed, due to simplified aerosol activation, as explained above.

One aspect that complicates the comparison of our simulations with observations is that we excluded snow and graupel. If snow and graupel were included in the simulations, then IWP might increase as riming converts liquid to frozen mass. On the other hand, this might also lead to faster precipitation and thus depletion of liquid or frozen water.

Note that our goal was to test F23 rather than include all the necessary processes in order to closely model the observations.

3.3 Comparison of F23 to interactive F62 implementation

F62 calculates the number of cloud droplets that freeze to ice crystals from the temperature-dependent INPC (m⁻³) given in F62 (and shown by the black line in Fig. 1), as follows:

$$\begin{array}{}\text{(3)}& \mathrm{INPC}\left(T\right)=\mathrm{0.02}\cdot \exp (-\mathit{\beta }\cdot T),\end{array}$$

with $\mathit{\beta }=\mathrm{0.6}/$(1 ^∘C) and T in degrees Celsius (^∘C). The changes in N_i, N_c, Q_i, and Q_c are calculated according to Eqs. (2a) and (2b). Note that this parameterization is only strictly valid for $-\mathrm{10}>T>-\mathrm{30}$ ^∘C, but we extrapolate it to the temperatures of the simulated cloud (dotted line in Fig. 1). Analogously to F23, freezing happens where T≤0 ^∘C and cloud droplets are present.

Figures 3b and D1 show that very little ice is produced using the F62 ice-nucleation scheme. We explain this by the subtraction of N_i from INPC(T), since this only leads to considerable ice formation at the beginning of the simulation, when no ice is present yet or when the temperature decreases. This comparison between F23 and F62 illustrates the difficulty in simulating reasonable ice masses in warm mixed-phase clouds with conventional schemes that yield one INPC for one set of environmental conditions. Drawing from a distribution of INPCs according to F23 leads to substantially larger ice masses, and these might even be multiplied if, e.g., secondary ice processes would also be considered (see the discussion in Sect. 3.2). In the future, it could be interesting to further analyze if the difference between F23 and F62 could be reduced by tuning F62 and how much it is related to the conceptual differences in the schemes.

3.4 Sensitivity studies of F23

To investigate the sensitivity of the simulated cloud to the characteristics of the F23 scheme, the following parameters (summarized in Table 2) were varied: (i) the median and standard deviation of the INPC distribution (Sect. 3.4.1 and 3.4.2); (ii) the size of the temperature bins (Sect. 3.4.3); (iii) the frequency of drawing (Sect. 3.4.4); and (iv) the resolution of the model domain (Sect. 3.4.5). The sensitivity tests (i) mean that we cover a wider spectrum of INPC RFDs than the one defined in Sect. 2.1 and shown in Fig. 1.

3.4.1 Median of the distribution

The median (exp (μ)) of the INPC RFD is multiplied by a factor (0.5–1.5) shifting the entire distribution in Fig. 1 up or down but without changing the standard deviation. Having more INPs (higher median) or fewer INPs (lower median) impacts the amount of ice formed in the cloud. The effect can be seen from the modeling results of IWP (Table 3; Fig. 5); when increasing or decreasing the median by 25 % or 50 %, IWP increases or decreases accordingly. The vertical profiles exhibit this symmetry even more clearly with linear N_i and Q_i increases or decreases (Fig. 6a and b), and all changes are significant. The significance of differences in IWP between simulations is tested with a two-sided Kolmogorov–Smirnov test at the 95 % level. Changes in the vertical profiles are tested with a two-sided t test at the 95 % level. It is expected that all variables concerning ice crystals increase (decrease) with increased (decreased) median INPCs, since more (fewer) INPs lead to more (fewer) cloud droplets freezing. The change is linear because the median is logarithmized in the formula's exponent (see Eq. 1). No change in the vertical distribution of cloud droplet concentration is apparent (not shown).

Table 2 Overview of sensitivity studies.

^a Simulation with diagnostic N_i and no interactive ice-nucleation parameterization. ^b Time step. ^c In this case, no distribution is used, but the median of the distribution is used for all freezing cases. ^d Simulation in which ice nucleation is started after 2 h, first with a drawing frequency of every 10 s; then after an additional 2 h, the drawing frequency is decreased to 5 min.

Table 3 The three quartiles (25th percentile is P₂₅, etc.) of IWP for all sensitivity studies for the final four simulated hours (8–12 h) and differences relative to F23 (ΔF23) or DIAG (ΔDIAG). Values are given in grams per meter squared (g m⁻²). Simulations with significant differences to F23 (DIAG) are highlighted with an asterisk. Significance was assessed with a two-sided Kolmogorov–Smirnov test at the 95 % level. D-10S-5M is not included because values are similar to 5M or 60M.

3.4.2 Standard deviation of the distribution

The standard deviation (σ) of the RFD determines the variability in the INP concentration. The wider the distribution, the larger the variability in drawn INPCs. For the modeled cloud, a larger variability results in substantially and significantly increased IWP, N_i, and Q_i, while a lower variability results in significantly smaller values (Figs. 7 and 8). The changes are exponential with linear changes of σ, since σ is in the exponent in the INPC RFD (Eq. 1). These results emphasize that it is the large INPCs that dominate ice formation in F23. Increased ice formation triggered by large INPCs could lead to cloud glaciation and subsequent cloud dissipation in colder clouds if ice crystals grow at the expense of liquid droplets due to the WBF process. In other simulations including secondary ice processes, for example, mechanical splintering or break up of ice crystals (e.g., Field et al., 2016), a high N_i can be relevant to even further enhance N_i. For example, Yano et al. (2016) report a critical N_i,crit, which can lead to an explosive enhancement of the ice crystal number. This N_i,crit might only be reached if there is a possibility of drawing large INPCs with F23. If large INPCs become less probable (smaller σ), then IWP, N_i, and Q_i decrease. This is also apparent when analyzing the simulation without an INPC distribution and using the RFD's median value at all time steps (green line in Figs. 7 and 8a and c). All ice variables (IWP, N_i, and Q_i) have the lowest values for this run (see Table 3). Using realistic median INPCs at the rather high temperature of the simulated case leads to almost no ice formation (see also Sect. 3.3, where the INPCs are much larger than the median of the F23 RFD, but almost no ice is produced). Once an INPC distribution is added, considerable ice is formed. Excluding negative ΔN_i (Eq. 2a) causes N_i to progressively increase even without decreasing cloud temperature. Such increases can be expected from time-dependent immersion freezing.

It is remarkable that the large increase in ice for S1.5 even leads to changes in N_c and Q_c (Fig. 8b and d). The cloud bottom is elevated by ca. 100 m in comparison to the F23 case (Fig. 8b). This is probably because more cloud droplets freeze, leading to a depletion of unfrozen cloud droplets. Within the cloud, N_c is very similar for S1.5 and F23. However, Q_c for S1.5 is significantly smaller, indicating an enhanced WBF process resulting in liquid droplet evaporation.

Figure 5 Domain-averaged IWP for 10 F23 simulations (cyan; median with a min/max envelope, and the comparison run used in Fig. 6 is shown as a dotted line). Four sensitivity runs are shown, where the median of the RFD was changed (yellow is increased median INPC; blue is the decreased median INPC). The significance was assessed with a two-sided Kolmogorov–Smirnov test at the 95 % level. Significant differences to F23 are indicated with an asterisk.

3.4.3 Size of the temperature bins

For the F23 simulation, we draw from the INPC RFD defined on 1 ^∘C temperature bins centered on whole degrees (i.e., −0.5 to −1.5 ^∘C is the −1 ^∘C bin). Thus, a change in temperature from, for example, −10.4 to −10.6 ^∘C shifts the INPC from the concentrations at −10 ^∘C to the concentrations at −11 ^∘C of the INPC distribution. Using smaller temperature bins, we expect the parameterization to be less sensitive to small temperature changes and lead to smaller changes in the drawn INPC concentrations. This effect was tested by decreasing the temperature bin range from 1 to 0.5 ^∘C (0.5Deg in Table 3). Figure 9 shows the IWP for a run with 0.5 ^∘C temperature binning. Overall, IWP for the 0.5Deg simulation does not significantly differ from the F23 runs. Until hour 8, the IWP is slightly below the F23 case, and between hours 9 and 11, the IWP is slightly larger. The averaged profiles for the final 4 simulated hours are shown in Fig. 10. None of the variables differs significantly between F23 and 0.5Deg (see Table 3). Judging from these results, there is no compelling benefit in decreasing the size of the temperature bins (from 1 to 0.5 ^∘C).

Figure 6 Profiles averaged over the domain and the simulation period of 8–12 h. (a) N_i. (b) Q_i. The F23 comparison run is plotted in dotted cyan (see the dotted line in Fig. 5), M0.5 and M0.75 are in blue, and M1.25 and M1.5 are in yellow. Simulations with significant differences to F23 are indicated with an asterisk (tested with a two-sided t test at the 95 % level). Horizontal dashed lines indicate the cloud's top and bottom in F23.

Figure 7 Domain-averaged IWP for 10 F23 simulations (cyan; median with a min/max envelope, and the comparison run used in Fig. 8 is shown as a dotted line). Five sensitivity runs are shown where the standard deviation of the RFD was changed (yellow is larger standard deviation; blue is smaller standard deviation; green is no standard deviation (S0)). Significant differences to F23 are indicated with an asterisk.

3.4.4 Frequency of drawing

The frequency of drawing a new value for the INPC sets the length of the time period during which the recently drawn INPC is representative of the INPC at the grid point. The frequency of drawing is coupled to the model's temporal resolution, since the maximum frequency of drawing is restricted by the time step of the model. The sensitivity to the drawing frequency was tested by drawing the INPC once every 5 s (5S; see Table 3), once every 10 s (10S), once every 20 s (20S), once every 5 min (5M), and once every 60 min (60M) instead of at every time step (F23; every 1–3 s). Freezing events still occur at every time step (according to Eqs. 2a and 2b), but within the respective time period (e.g., 5 s or 60 min) the INPC is constant at one grid point.

Figure 8 Profiles averaged over the domain and the simulation period of 8–12 h. (a) N_i. (b) N_c. (c) Q_i. (d) Q_c. The F23 comparison run is plotted in dotted cyan (see the dotted line in Fig. 7), S0 is in green, S0.5 and S0.75 are in blue, and S1.25 and S1.5 are in yellow. Simulations that yielded significant differences to F23 are shown in the respective variable's legend in panels (a), (c), and (d). Horizontal dashed lines indicate the cloud's top and bottom in F23.

Figure 9 Domain-averaged IWP for 10 F23 simulations (cyan; median with a min/max envelope, and the comparison run used in Fig. 10 is shown as a dotted line). Seven sensitivity runs are shown (blue is 0.5Deg, solid yellow is 5S, solid green is 10S, solid magenta is 20S, dashed yellow is 5M, dotted green is 60M, and dash-dotted black is D-10S-5M). Significant differences to F23 are indicated with an asterisk.

IWP, N_i, and Q_i for the runs with lower drawing frequencies exhibit similar relative increases over the simulation time as F23 but have significantly lower absolute values (Fig. 9). The profiles of ice variables also decrease as the frequency of drawing decreases (Fig. 10). This evolution can be explained by the subtraction of N_i at each time step. N_i at a grid point will not change between time steps for which INPC<N_i. Only when a newly drawn INPC>N_i will new ice be formed at the grid point. An exception is if there is a sink for N_i at the grid point (e.g., sedimentation, sublimation, or advection) and N_i becomes smaller than INPC before the next draw of INPC, since INPC−N_i cloud droplets will then freeze. As the overall chance of drawing INPC>N_i is higher with a higher drawing frequency, the frequency of drawing new INPCs changes the amount of new ice formation and ice content in the cloud. Note that since large INPCs and the chance to draw these large INPCs determine the ice in the cloud, this introduces an indirect time-dependency of the scheme. The resulting ice in the cloud depends indirectly on the time until a large INPC value is drawn at a specific grid point. Increasing the drawing frequency leads to increases in the ice variables that do not appear to converge (Figs. 9 and 10; Table 3). Ideally, one expects that the variables should converge for higher drawing frequencies. This could possibly be achieved in future implementations of the scheme by introducing a correction factor that depends on the drawing frequency. As stated above, the maximum drawing frequency is the time step of the model, but the fact that the resulting cloud ice diverges with increased drawing frequency poses a limitation of F23. Note that this does imply that F23 depends on the time step of the model. In order to investigate the behavior of increasing drawing frequencies beyond what has been shown here, the time step of the model would need to be decreased further to be able to decrease the drawing frequency. Even though this could theoretically be done, it would increase the computational costs tremendously, and other aspects of the simulation would potentially change as well, which impedes a comparison. Therefore, we refrained from such a test. Another feature apparent in Figs. 9 and 10, as well as Table 3, is that the ice variables converge for drawing frequencies of 5 min or larger.

Figure 10 Profiles averaged over the domain and the simulation period of 8–12 h. (a) N_i. (b) Q_i. The F23 comparison run is plotted in dotted cyan (see the dotted line in Fig. 9), 0.5Deg is in blue, 5S is in solid yellow, 10S is in solid green, 20S is in solid magenta, 5M is in dashed yellow, and 60M is in dotted green. Simulations with significant differences to F23 are indicated with an asterisk. Horizontal dashed lines indicate the cloud's top and bottom in F23.

The simulation D-10S-5M in Fig. 9 is a run for which no ice nucleation is applied during the first 2 h of the simulation. After that, F23 is called, with a drawing frequency of 10 s for 2 h, and for the final 8 h, the drawing frequency was set to 5 min. The simulation is not included in Fig. 10 and Table 3, since it results in the same values as 5M or 60M. D-10S-5M shows that the current drawing frequency determines the amount of ice in the cloud.

It is important to investigate how the ice content reacts on very high and low drawing frequencies (asymptotic behavior). In this study, however, we are constrained by the model setup. Yet, we can expect that processes other than heterogeneous ice nucleation would limit ice production for high drawing frequencies in more realistic model setups. For example, ice multiplication would lead to a quick increase in N_i, which would result in small or absent ice nucleation due to the subtraction of N_i, according to Eq. (2a). We assume that the converging behavior of the ice mass for low drawing frequencies (no difference between 5M and 60M) is caused by other components of the model than the heterogeneous ice-nucleation scheme. If there were no other processes affecting ice number, then decreasing the drawing frequency can be expected to lead to decreasing ice in the cloud.

3.4.5 Resolution of the domain

The spatial resolution of the model domain affects the F23 scheme in a similar way as the frequency of drawing does, by changing the number of draws of INPC for the same cloud. However, a change in the model resolution impacts all parts of the model, including microphysical processes. The sensitivity of the LES at a lower resolution is tested by doubling the grid spacing in all three dimensions while the domain volume remains constant (DIAG LR and F23 LR). Note that this will also lead to a doubling of the model time step, since it is calculated dynamically to satisfy the CFL criterion. Testing the scheme on a coarser grid is important because one goal for the new parameterization is that it can be adapted to larger-scale models, where the resolution will be much coarser than in large-eddy simulations. Lowering the resolution leads to significantly lower IWP for both the DIAG LR and F23 LR simulations (Fig. 11). IWP is affected more in the simulation with the F23 scheme. The vertical profiles of the cloud droplet variables (Fig. 12b and d) show non-significant but consistent decreases in N_c (Fig. 12b) and Q_c (Fig. 12d) for both runs with a lower resolution. The decrease might be caused by changes in mixing processes, e.g., entrainment at the cloud borders, due to the change in grid size. N_i within the cloud is fixed to 200 m⁻³ for DIAG LR and DIAG (Fig. 12a), as per the definition (see Sect. 2.2). Nevertheless, Q_i is lower for DIAG LR compared to DIAG (Fig. 12c). This can be explained by the decrease in Q_c and thus a less efficient WBF process providing a smaller mass flux from liquid droplets to ice crystals. For F23 LR, all variables shown in Fig. 12 decrease – for N_i and Q_i these changes are significant. The decrease in N_i can be explained by the smaller number of grid points in the domain and a coarser temporal resolution. INPCs will be drawn at only one-eighth of the number of grid points and only about half as often in F23 LR; thus, large INPCs will be drawn less often than for F23. The results shown in Figs. 11 and 12 underline the conclusions from Sect. 3.4.2 and 3.4.4 that the possibility of drawing very large INPCs has the strongest effect on overall N_i.

Figure 11 Domain-averaged IWP for 10 DIAG (red; median with a min/max envelope, and the comparison run used in Fig. 12 is shown as a dotted line) and 10 F23 simulations (cyan; median with min/max envelope, and the comparison run used in Fig. 12 is shown as a dotted line). One run is done with half of the entire resolution, respectively (dashed lines). Significant differences to DIAG or F23 are indicated with an asterisk.

Figure 12 Profiles averaged over the domain and the simulation period of 8–12 h. (a) N_i. (b) N_c. (c) Q_i. (d) Q_c. The DIAG comparison run is plotted in dotted red, the F23 comparison run is plotted in dotted cyan (see the dotted lines in Fig. 11), DIAG LR is in dashed red, and F23 LR is in dashed cyan. Simulations F23 LR in panels (a) and (c) yielded significant differences to F23. Horizontal dashed lines indicate the cloud's top and bottom in DIAG.

3.5 Limitations and challenges

The limitations and challenges of this study can be divided into three aspects, namely the limitations related to the parameterization itself, the data set it is based upon, and the limitations due to the modeling setup used to test the parameterization scheme.

3.5.1 Limitations of the parameterization scheme

Using a scheme with a random component, in our case the random drawing of INPC values, can introduce difficulties when it comes to the technical implementation of such a scheme. Here, we assumed independence between consecutive random draws from the INPC RFD, but a certain degree of autocorrelation both in time and space might need to be added to the parameterization. Another challenge is the high sensitivity of simulated cloud ice to the INPC RFD's standard deviation, the drawing frequency, and the model resolution (i.e., the possibility of drawing large INPCs). This aspect is partially related to the general difficulty in the creation of immersion freezing schemes, where deterministic measurements have to be time-discretized. This leads to the challenge of selecting the best parameters for F23. Further investigations and constraints for the parameters in F23 are needed.

3.5.2 Limitations related to the data set that the parameterization scheme is based upon

We base our parameterization's RFD on an extensive field-based data set. However, this does not necessarily represent a globally valid INPC distribution. Since our study focuses on the conceptual idea of using an INPC distribution for parameterizing immersion freezing, we considered this limitation to be of secondary concern. Additionally, we present sensitivity tests with adjusted distributions that cover a large space of INP concentrations (i.e., the whole range of possible INPC). When using this scheme in the future, region or source-based RFDs can be applied. Another limitation is that the underlying INPC measurements in this study are surface based; we assume the boundary layer to be well-mixed, and thus we assume that concentrations are similar at the cloud level. How surface-based measurements of INPs can represent INPC in the cloud-layer and be used for general freezing parameterizations is still an unanswered question (independent of this study). Finally, INP concentrations might change in a future climate and differ from the INPC RFD assumed here. However, it is easy to use an adjusted distribution in F23.

3.5.3 Limitations of the modeling setup used to test the parameterization scheme

The model used for testing the F23 scheme, MIMICA LES, and the chosen case study, ASCOS, also represent a limitation of the study. The presented results and sensitivities are obtained by this model and based on this one case study, which results in further necessary investigations of F23.

More specifically, the choice of F23's INPC distribution is a limitation. While we simulate an Arctic cloud in this study, F23's INPC distribution is based on measurements at Cabo Verde and other maritime locations. However, Arctic INP measurements (e.g., from Ny-Ålesund in Wex et al., 2019) are covered by our INPC RFD. As mentioned above, another limitation is that F23's INPC distribution is based on surface-based measurements. However, in our case study, the cloud is actually decoupled from the surface; thus, surface INPC may not be representative of INPC in the Arctic cloud layer (e.g., Creamean et al., 2018; Griesche et al., 2021).

Another limitation, when it comes to the modeling setup and especially the comparison of the simulation results to observations, is the general setup of the model (i.e., the representation of microphysical processes). However, we would like to stress that the comparison to observations is not the purpose of this study.

4 Conclusions

A novel parameterization of immersion freezing that takes into account the observed variability in the ice-nucleating particles is presented. The observed INPC variability is reproduced by random drawing from an INPC relative frequency distribution that depends on temperature. This means that the INP population at a specific grid point is represented by a new value at the chosen frequency of drawing. If the newly drawn INPC exceeds the present N_i, then additional ice is formed. The F23 parameterization is valid for the entire temperature range of heterogeneous immersion freezing between 0 and −38 ^∘C. F23 has the additional advantage that it does not require information on the present bulk aerosol from the atmospheric model, which makes it easy to implement and use in many different models. The main goal of this study was to test the general approach of representing immersion freezing by random drawing from an INPC distribution, as opposed to using traditional parameterizations that yield one INPC for the given temperature like, e.g., Fletcher (1962).

We tested the parameterization for the large-eddy simulation of a mixed-phase Arctic stratocumulus cloud case with MIMICA. We used this case as a test bed for F23 because aerosol characteristics in the Arctic are largely unknown, and this lack of information on aerosol bulk properties makes it challenging to use “classical” aerosol-aware freezing schemes. Moreover, it might be important to consider the whole range of possible INPCs in the Arctic summer, where warm mixed-phase clouds are frequently observed. The simulations with F23 lead to less cloud ice than was observed, but our model setup does not include graupel and snow and ice-nucleation modes, other than immersion freezing or secondary ice processes. The latter has been shown to increase IWP by a factor of 2–3, leading to the observed values in MIMICA simulations of the same case by Sotiropoulou et al. (2021). Through our sensitivity tests, we found that the simulated IWP, N_i, and Q_i of the cloud depend linearly on the median of the RFD that describes the INPC variability and exponentially on the distribution's standard deviation. The large dependence on the standard deviation of the distribution is especially interesting, as it implies that the amount of ice in the modeled cloud is particularly sensitive to large INP concentrations. The possibility of drawing large INPC can influence cloud glaciation for colder cases caused by the WBF process and secondary ice processes when N_i,crit is reached (Yano et al., 2016). The relevance of randomly drawing INPC from a distribution is highlighted by the much lower IWP, N_i, and Q_i when applying the F62 parameterization or simply using the median INP concentration (i.e., no INPC distribution). This emphasizes that it is the rare but large INPCs that control the freezing in the cloud, rather than the median INPC values. Additionally, the frequency of drawing a new INPC has a significant impact on the ice variables, since it is when a new, larger INPC is drawn that additional ice is formed (Eq. 2a). The higher the frequency of drawing, the higher the ice content in the cloud. However, drawing every 5 min, or less often (e.g., every 60 min), results in a similar amount of ice in the cloud. The high sensitivity of simulated cloud ice to an increased possibility of drawing large INPC values (as tested by the increased standard deviation of the RFD, high drawing frequency, and higher spatial and temporal model resolution) poses the challenge of choosing the parameters (RFD standard deviation and drawing frequency) for F23, while being in accordance with INPC observations. Further investigations are necessary to specify these parameters in order to apply the parameterization to other cases or in other models (including different model time steps). The unexpected divergence of the cloud ice amount for increased frequencies of drawing needs to be examined as well.

The scheme's independence from aerosol information in the atmospheric model is a strength but can be a limitation. The proposed INPC distribution of F23 may not be representative of distinct scenarios, e.g., a Sahara dust outbreak. The RFD would need to be updated with one that is based on INPC observations specific to an area or event of interest. However, the parameterization is flexible and can be adapted to different INPC observations. To represent several sources or source regions at the same time in a model, different RFDs could be used, depending on the location (remote vs. continental INP as an example). The independence of modeled aerosol might require updating the RFD when simulating future scenarios including changes in aerosol or INP concentrations. Some degree of autocorrelation between subsequent random draws from the INPC RFD in time and space could be added to the scheme in the future. However, it is not clear what degree of autocorrelation would be physically reasonable. Since our study is based on one rather warm Arctic stratocumulus case, the parameterization should be tested and validated for other conditions and case studies in the future, especially for lower in-cloud temperatures. Additionally, F23 should be tested in LES models other than MIMICA and in larger-scale models. The intention of this study was to illuminate what including randomness in INPCs leads to. We here show that representing the INP heterogeneity in an immersion freezing parameterization allows for a realistic simulation of an Arctic stratocumulus cloud, but we are clearly left with remaining challenges.

In order to make it easier to use INPC observations in schemes like the one presented, measurements should be reported either as individual measurement points (e.g., as time series) or, if aggregated, with the average and standard deviation and percentiles of the underlying lognormal distribution.

Appendix A: MIMICA LES model description

The MIMICA LES solves a system of non-hydrostatic anelastic equations and represents cloud microphysics through a two-moment bulk microphysics scheme. The prognostic variables are number concentrations and mass mixing ratios of up to five represented hydrometeors, including cloud droplets, raindrops, cloud ice, snow, and graupel. All hydrometeor mass distributions are regular gamma distributions, and the hydrometeors' terminal fall speeds are calculated from simple power laws dependent on their diameters. Warm microphysical processes are modeled according to Seifert and Beheng (2001, 2006), while collection processes involving frozen hydrometeors and resulting in the hydrometeors sticking together follow Wang and Chang (1993). MIMICA represents the following cold collection processes: (i) riming of ice crystals and graupel by cloud droplets (resulting in graupel); (ii) riming of ice crystals, graupel, and snowflakes by raindrops (resulting in graupel); (iii) autoconversion of ice crystals to snow; (iv) self-collection of snowflakes; (v) collection of cloud droplets by snow (resulting in snow); (vi) growth of a snowflake by aggregating ice crystals; and (vii) collection of snow by graupel (resulting in graupel). Cloud condensation nuclei (CCN) activation is described, following Khvorostyanov and Curry (2006), which is based on a simple power law depending on modeled supersaturation and a prescribed background CCN concentration. It is possible to describe aerosol prognostically in MIMICA, including activation as CCN. The supersaturation is solved pseudo-analytically, following Morrison and Grabowski (2008). In the standard version of MIMICA, primary ice formation is represented by maintaining a constant ice crystal number concentration (N_i) within the cloud. However, it is also possible to choose an interactive ice-nucleation scheme. Radiation is calculated according to Fu and Liou (1993). The initial ice or liquid potential temperature profiles are randomly perturbed in order for convection to develop more quickly. Consequently, any two simulations will yield different results, even if all parameters are held constant.

Appendix B: Averaged Q-tendency profiles for DIAG vs. F23

Figure B1 shows tendencies relevant to the Wegener–Bergeron–Findeisen process (WBF), averaged over the domain and entire simulation time for DIAG and F23, with condensation (positive values) and evaporation (negative values) for liquid water, $Q_{\mathrm{l}}=Q_{\mathrm{c}}+Q_{\mathrm{r}}$, in Fig. B1a and deposition (positive values) and sublimation (negative values) for ice crystals, Q_i, in Fig. B1b. More deposition and evaporation occur in DIAG, illustrating the WBF process with a mass transfer from liquid to frozen water.

Figure B1 Profiles of WBF tendencies averaged over the domain and the entire simulation period, with (a) evaporation and condensation, $\mathrm{d}{Q}_{\mathrm{l}}=\mathrm{d}{Q}_{\mathrm{c}}+\mathrm{d}{Q}_{\mathrm{r}}$, and (b) sublimation and deposition, dQ_i. One DIAG is dotted red, and one F23 run is dotted cyan (see the dotted lines in Fig. 3). Horizontal dashed lines indicate the cloud's top and bottom in F23 during the last 4 h of the simulation.

Appendix C: Averaged N profiles for DIAG vs. F23

Figure C1 shows the time- and domain-averaged profiles of the number concentrations of ice, cloud droplets, and raindrops for DIAG and F23.

Figure C1 Profiles averaged over the domain and the simulation period of 8–12 h. (a) N_i. (b) N_c. (c) N_r. One DIAG is dotted red, and one F23 run is dotted cyan (see the dotted lines in Fig. 3). Horizontal dashed lines indicate the cloud's top and bottom in DIAG.

Appendix D: Averaged N_i and Q_i profiles for F23 vs. F62

Figure D1 shows the time- and domain-averaged profiles of the number concentration and mixing ratio of ice for F62 compared to F23.

Figure D1 Profiles averaged over the domain and the simulation period of 8–12 h. (a) N_i. (b) Q_i. The F23 comparison run is plotted in dotted cyan (see the dotted line in Fig. 3), and one run with F62 is blue. Simulations with significant differences when compared to F23 are indicated with an asterisk (tested with a two-sided t test at the 95 % level). Horizontal dashed lines indicate the cloud's top and bottom in F23.