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 1Relative 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).

Download

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

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:

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.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 1Average 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.

Download Print Version | Download XLSX

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

Download

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.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:

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 2Overview 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.

Download Print Version | Download XLSX

Table 3The 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.

Download Print Version | Download XLSX

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 5Domain-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.

Download

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 6Profiles 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.

Download

Figure 7Domain-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.

Download

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 8Profiles 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.

Download

Figure 9Domain-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.

Download

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 10Profiles 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.

Download

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 11Domain-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.

Download

Figure 12Profiles 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.

Download

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.