Modelling mixed-phase clouds with large-eddy model UCLALES-SALSA

The large-eddy model UCLALES-SALSA, with exceptionally detailed aerosol description for both aerosol number and chemical composition, has been extended for ice and mixed-phase clouds. Comparison to a previous mixed-phase cloud model intercomparison study confirmed the accuracy of newly implemented ice microphysics. Further simulation with a heterogeneous ice nucleation scheme, where also ice nucleating particles (INP) are a prognostic variable, captured the typical layered structure of Arctic mid-altitude mixed-phase cloud: a liquid layer near cloud top and ice within and below the liquid 5 layer. In addition, the simulation showed realistic freezing rate of droplets within the vertical cloud structure. The represented detailed sectional ice microphysics with prognostic aerosols is crucially important in reproducing mixed-phase clouds.


Introduction
Clouds are known to have a prominent influence on the hydrological cycle and the atmospheric radiation balance. While significant advances have been made in characterisation of liquid-phase clouds, the microphysical processes, especially heterogeneous ice nucleation, dynamics and radiative effects of mixed-phase and ice clouds remain more poorly constrained. This is mainly because of challenges in obtaining representative observations and a lack of a detailed enough representation of microphysics in climate and numerical weather prediction models. Specific challenges are known to be associated with aerosolcloud interactions (Cox, 1971;Knight and Heymsfield, 1983;Curry, 1995;Solomon et al., 2007;Stevens and Feingold, 2009;Morrison et al., 2011a;Morrison, 2012;Li et al., 2013). 15 What we know about mixed-phase clouds is that by definition supercooled liquid droplets coexist with ice crystals. Such clouds are frequent at temperatures between −10 • C and −25 • C (Filioglou et al., 2019), but can be present from −35 • C to 0 • C and require specific microphysical and dynamical conditions (Andronache, 2017). Ice crystals can form either by homogeneous or heterogeneous freezing (term nucleation used also). In temperatures lower than −38 • C, liquid droplets can freeze homogeneously without the need of ice nucleating particles (INP). In heterogeneous ice nucleation, freezing initiates 20 from the surface of seed particle and can occur at higher temperatures than homogeneous ice nucleation. Droplet freezing processes are not yet fully understood and quantified despite of decades of research (Phillips et al., 2008;Atkinson et al., 2013;DeMott et al., 2011). Kiselev et al. (2017) stated that ice formation on aerosol particles (heterogeneous ice nucleation) is a process of crucial importance to Earth's climate, but it is not understood at the molecular level. However, in Morrison et al. (2011a) it is noted that although many details of droplet freezing are poorly understood, enough knowledge exists to draw 25 first-order (ice water path) conclusions. Furthermore, droplet freezing models and even the representation of cloud structure often require a resolution that is too detailed for large scale models. For instance, the structure of Arctic and mid-altitude clouds is complex with a layered structure with liquid near cloud top and ice within and below the liquid layer (Curry et al., 1997;Hobbs and Rangno, 1998;Pinto, 1998;Rangno and Hobbs, 2001;Zuidema et al., 2005;Shupe et al., 2006;Verlinde et al., 2007;de Boer et al., 2009;McFarquhar et al., 2011;Morrison et al., 2011a). The lack of a proper calculation of ice processes in 30 climate models is seen in the comparisons to observations of mid-and high-latitude mixed-phase clouds. These models tend to underestimate the lifetime of such clouds (Andronache, 2017). Better quantification of droplet freezing processes is expected to narrow the gap between observations and model results.
The importance of including detailed aerosol description is vital in cloud resolving models. Scarcity of INP is an important factor in mixed-phase clouds lifetime and structure, since roughly one in a million particles acts as an ice nucleus and even 35 these particles might have highly different ice-forming activity at different temperatures (Lebo et al., 2008;Morrison et al., 2011a). Therefore, the loss of INP along with precipitating ice crystals limits cloud glaciation and dissipation (Rauber and Tokay, 1991;Harrington et al., 1999;Avramov and Harrington, 2010). Describing this process is not possible without detailed description of aerosols, as is demonstrated in a 1-D cloud model study by Morrison et al. (2005). The significance of aerosols is shown in Filioglou et al. (2019) where high aerosol load was linked with a higher occurrence of mixed-phase clouds. Also, 40 the Norgren et al. (2018) study shows how there is less ice in polluted clouds. Andronache (2017) and Morrison et al. (2011a) are comprehensive review resources to give further details about mixed-phase clouds.
There is a growing number of studies focusing on examining the properties of mixed-phase or ice clouds by combining observations and models, including large eddy simulation (LES) modelling and other cloud-resolving models (CRM) (Jiang et al., 2000;Klein et al., 2009;Morrison et al., 2011b;Ovchinnikov et al., 2014;Andronache, 2017). Large-eddy simulations 45 are particularly attractive for modelling boundary layer clouds since they offer a good compromise between computational cost and accuracy in terms of model resolution (Tonttila et al., 2017;Andronache, 2017). LES models solve the largest eddies in turbulent flows explicitly and use parametrisations for the smallest length scales. In atmospheric applications they are usually coupled with cloud microphysical packages. Recent developments in computational performance of supercomputers have made explicit and detailed description of aerosol-cloud-ice microphysical interactions possible also in LES modelling, allowing the 50 investigation of non-linear cloud phenomena, such as secondary ice production and heterogeneous ice nucleation.
There are several LES models that solve cloud related interactions (Fridlind et al., 2012;Khain et al., 2004;Savre and Ekman, 2015;Fu and Xue, 2017). In comparison to those models, we present a LES model UCLALES-SALSA that brings additional value with a more detailed aerosol description. UCLALES-SALSA explicitly resolves interactions between aerosols, ice crystals and cloud droplets with sectional microphysics for all hydrometeors while keeping track of the aerosol dry size 55 distribution. Sectional description, especially for aerosols, is a clear asset of UCLALES-SALSA and we have now extended this description also for ice crystals. This sectional aerosol description allowed the implementation of a detailed heterogeneous freezing processes. First, the model results are compared with previously published modelling results. Finally, we demonstrate the benefits of this approach to handle heterogeneous freezing over more simplified aerosol-ice-cloud treatments.

Model description 60
The UCLALES-SALSA (Tonttila et al., 2017) model consists of two components: first, the widely used large eddy simulator UCLALES (Stevens et al., 1999(Stevens et al., , 2005, and second the aerosol bin microphysics model SALSA (Sectional Aerosol module for Large-Scale Applications) (Kokkola et al., 2008;Tonttila et al., 2017;Kokkola et al., 2018). UCLALES handles e.g. surface fluxes, transportation of microphysical prognostic variables and atmospheric dynamics including turbulence. The previous version of UCLALES-SALSA incorporated interactions between aerosols, clouds and drizzle (Tonttila et al., 2017). Now we 65 have extended the model with a description for ice crystals. In this study, we focus on how ice crystals and ice nucleating particles (INP) interact with clouds while tracking sectional aerosol size distribution. Figure 1 illustrates the microphysical treatment of different hydrometeor classes and their size distributions in UCLALES-SALSA. All four classes (aerosol, cloud and rain droplets and ice crystals) are tracked with a bin scheme. Bin scheme offers the benefit of greater accuracy in simulating interactions between different sized hydrometeors. Better accuracy is gained 70 by dividing the size distribution into bins. This also enables better flexibility as the shape of the distribution is allowed to evolve. Bulk schemes provide a simpler method and track one or several moments of the size distribution, where the shape of distribution is prescribed. The disadvantage of bin scheme is higher computational cost compared to bulk scheme.
Three of the hydrometeor classes, i.e. aerosol, cloud droplets and ice, are further divided into parallel bins labelled a and b as shown in Fig. 1. This division into a and b bins is done to enable the tracking of externally mixed distributions and to see 75 how different particles affect clouds. For aerosol particles, subrange 1a is an additional feature to describe nucleation mode.
Otherwise, Aitken and accumulation mode size ranges are sufficient to characterise cloud phenomena.
The aerosol, cloud and ice crystal size distributions are discretised into the bins according to the dry aerosol diameter, whereas the rain droplet size distribution is defined by the wet diameter of the droplet. Identical 2a and 2b size bins are used for aerosol, cloud droplets and ice. Such parallel bins are useful for tracking aerosol development through cloud activation, freezing 80 and sublimation. Prognostic variables for each bin include aerosol number and masses of all compounds (water, sulphate, dust, organic carbon, sea salt, nitrate, and ammonium).
In UCLALES-SALSA, recently implemented processes involving ice crystals are droplet freezing, deposition of water vapour, sublimation, melting when T > 0 • C, coagulation between different sized hydrometeors, sedimentation, and interactions with radiation (see also Fig. 1). Most of these processes are included in a similar way as in the previously published 85 version of UCLALES-SALSA (Tonttila et al., 2017).
Regarding the scope of this study, we describe droplet freezing in higher detail. There are five mechanisms for droplet freezing and they are all currently implemented in UCLALES-SALSA.
-Immersion freezing is possible for aqueous droplets that have an insoluble core, which in UCLALES-SALSA is either dust (DU) or black carbon (BC). The rate of heterogeneous germ formation in a supercooled droplet of water or solution 90 is calculated mostly following Khvorostyanov and Curry (2000) and additional parameters are from Jeffery and Austin (1997), Khvorostyanov andK. Sassen (1998), Khvorostyanov andCurry (2004) and Li et al. (2013). See also App. A.
-Homogeneous freezing is possible for any aqueous droplet with or without insoluble particles. This is applied to the model according to Khvorostyanov and K. Sassen (1998). See also App. B.
-Deposition freezing is possible for dry insoluble aerosol at sub-saturated conditions (RH < 100%). This is implemented 95 following Khvorostyanov and Curry (2000) and additional parameters from Hoose et al. (2010). See also App. C.
-Contact freezing is implemented in UCLALES-SALSA following Hoose et al. (2010) so that first the coagulation code is used to calculate collision rates between dry particles and liquid droplets and then immersion freezing code gives the freezing probability.
-Condensation freezing is implemented as a part of immersion freezing, because these droplets can freeze during the 100 modelled condensational growth.
In our simulations (Sect. 3.3), only immersion freezing is active. This is due to high temperatures, when homogeneous freezing is not possible, and mixing state of the INP leading to aqueous droplets, when deposition and contact freezing are not feasible.
Deposition of water, i.e. diffusion limited condensation or evaporation of water vapour, is defined for aerosol when relative 105 humidity (RH) is over 98% and always for other hydrometeors. This is based on the analytical predictor of condensation (APC) scheme by Jacobson (2005)  If ice sublimates, the immersed aerosol nuclei are added back to the aerosol population.
Activation of aerosols to cloud droplets happens when RH is over 100% and aerosol wet diameter exceeds the critical limit 110 corresponding to the resolved supersaturation. At this time, a certain proportion of activated aerosols (i.e. cloud condensation nuclei, CCN) are moved to cloud droplet bins.
Sedimentation is defined as before in Tonttila et al. (2017) and now extended for ice particles. For simulations in this study, a fall rate of ice particles is set as in Ovchinnikov et al. (2014).
Coagulation is implemented with the same way as before and now including also ice particles. Coagulation is affected by 115 diffusion, especially aerosols, and by sedimentation, especially large particles. In a collision, bigger particles absorb smaller particles.
Interaction with radiation is implemented either with the same four-stream radiative transfer solver (Fu and Liou, 1993) as in Tonttila et al. (2017) with extension to include ice particles or parametrised as in Ovchinnikov et al. (2014). We used the latter method in our simulation. In the parametrised radiation, the net upward long-wave radiative flux is computed as a function 120 of liquid water mixing ratio profile. The effect of interaction with radiation can be seen in simulations how radiative cooling weakens after liquid water path decreases below a specific value. Furthermore, UCLALES-SALSA was upgraded with minor bug fixes and improvements. For example, hygroscopicity is now calculated with κ-Köhler (Petters et al., 2006;Petters and Kreidenweis, 2007) instead of previously used ZSR method (Stokes and Robinson, 1966).

Model evaluation
The model performance is evaluated by simulating a well-documented mixed-phase cloud case from the Indirect and Semi-Direct Aerosol Campaign (ISDAC) Arctic study (McFarquhar et al., 2011). This observation case has been used before for comparisons to LES-models (e.g., Savre and Ekman, 2015;Fu and Xue, 2017). Ovchinnikov et al. (2014) presented an in-130 tercomparison of 11 LES models for this same case, where initial profiles were based on aircraft observations in the mixed layer (Flight F31) and idealisation of a sounding on 26th April 2008 at Barrow, AK. Nine of those models had bulk 2-moment microphysics and two of them bin microphysics.
We implemented in UCLALES-SALSA model runs the same semi-idealised simulation setup given in Ovchinnikov et al. (2014) that attempted to minimise intermodel differences by applying identical descriptions for the following processes: surface 135 properties and fluxes (fluxes set to zero), large-scale forcings, radiation, cloud droplet freezing and ice growth processes and sedimentation, and the nudging of horizontal winds, potential temperature and water content above the altitude of 1200 m.
In the simulations ice processes were excluded during the first two hours, i.e. the spinup period, to allow the mixed-layer turbulence and the warm stratus cloud to develop. After the spinup, cloud droplets are allowed to freeze until a specified target ice concentration is reached (Morrison et al., 2011b). Ice shape is described with a mass-diameter parametrisation so that ice 140 can be considered as spherical particles with low effective density (ρ = 84.5 kg m −3 ). Ice fall speed is related to the maximum dimension while capacitance, which is used in the condensation equation, is modified from that of a sphere to C = D/π. Radiation and sedimentation were parametrised similar to Ovchinnikov et al. (2014). For the sake of simplicity, coagulation was switched off as in Ovchinnikov et al. (2014). Warm rain formation was switched off allowing more straightforward model intercomparison. Also, the warm rain mass mixing ratios would have been small due to relatively small cloud droplet size in the 145 simulated case. Aerosol size distribution is given as a sum of lognormal accumulation and coarse modes with concentrations of 159 × 10 6 and 6.5 × 10 6 kg −1 , mode mean diameters of 0.2 and 0.7 µm and geometric standard deviations of 1.5 and 2.45, respectively. Aerosol is composed of ammonium bisulphate. During the simulations this aerosol size distribution provides on average 129 × 10 6 kg −1 cloud droplets. These aerosol distribution parameters provide the best fit to the measured distributions below the liquid cloud layer (Earle et al., 2011;Ovchinnikov et al., 2014). 150 We ran UCLALES-SALSA for the three different simulation setups investigated in the Ovchinnikov et al. (2014) study: no ice (ICE0), average ice (ICE1) and high ice (ICE4) number concentration. The ice number concentration is the only variable that was changed between the evaluation simulations (either 0, 1 or 4 L −1 ). Liquid and ice water paths (marked LWP and IWP from now on), i.e. column integrated mass values averaged over the horizontal domain, in these three cases show how water is distributed between ice and liquid phases depending on the ice crystal concentration.
155 Figure 2 compares the three UCLALES-SALSA simulations to the results presented in the Ovchinnikov et al. (2014) intercomparison paper. In the figure, LES model results from Ovchinnikov et al. (2014) are separated between bulk and bin microphysics to highlight the differences between microphysics schemes. First, Fig. 2a shows LWP for the baseline simulation without any ice (ICE0). It is evident that our model agrees well with the other 11 models. The simulated LWP of UCLALES-SALSA is in the middle of the model spread. Differences are most likely explained by different dynamical cores, which is 160 stated also in Ovchinnikov et al. (2014). A more thorough testing of warm phase cloud microphysics in UCLALES-SALSA was done in the Tonttila et al. (2017), and for the remainder of this work we will concentrate on examining the properties of the ice microphysics implementation.
Second, Figs. 2b and 2c present the LWP and IWP time series, when the target ice number concentration is 1 L −1 , marked with ICE1. Again, LWP in UCLALES-SALSA matches well with the other models being in the lower end of the intermodel 165 spread. As expected, the LWP growth rate is lower than in ICE0 simulation, as some of the water vapour condenses onto ice crystals. Furthermore, IWP matches well especially with other bin models in the Ovchinnikov et al. (2014) study.
Third, Fig. 2d depicts LWP time series with ice number concentration of 4 L −1 , which can be regarded as high ice concentration and is marked with ICE4. Now after spinup, LWP has a decreasing trend since the ice number concentration is so high that it consumes much of the water vapour. Subsequently, IWP in Fig. 2e increases rapidly after the spinup and in UCLALES-170 SALSA reaches its peak value of 15.7 g m −2 just before 4 hours of simulation. It then decreases to a value of 9.4 g m −2 at the end of the simulation. The reduction of IWP is caused by ice crystal precipitation at the surface and evaporation below the cloud.
Compared to the model results in Ovchinnikov et al. (2014), IWP in UCLALES-SALSA declines faster after the peak IWP has been reached in ICE4. One reason for this is that dry particle size is tracked in UCLALES-SALSA and this seem to have 175 an important effect on ice crystal sedimentation. Namely, sedimentation velocities and particle mixing (flux divergency) are here calculated for the dry size bins rather than bins tracking ice particle size. This reduces particle flux especially in the lowest 200 m leaving more particles there to evaporate. Evaporative cooling leads to a surface inversion which prevents the mixing of moist surface air. As such, the higher sensitivity to INP concentration is partly related to the initial conditions of the ISDAC case study. The other reason is related to the model dependent technical details. Our test simulations (not shown) indicate that 180 changing model options, such as flux limiter method, impact IWP and LWP so that the gap between UCLALES-SALSA and the other models decreases. In Ovchinnikov et al. (2014) it was also stated that when the ice number concentration gets higher the differences between models are more caused by discrepancies in microphysics than cloud dynamics. This underlines the sensitivity and significance of microphysics.
To conclude, the spread between models, especially between bin and bulk microphysics models, gets wider as the prescribed 185 ice number concentration gets larger and closer to the limit when the cloud glaciates completely. In UCLALES-SALSA this limit of full glaciation is lower than in other models in Ovchinnikov et al. (2014). This limit is further examined in Sect. 3.2.
Motivated by simulated differences with the 4 L −1 ice concentration, we wanted to further investigate how sensitive cloud glaciation is to changes in ice number concentration. In addition to ICE1 and ICE4 simulations, we performed simulations in 190 which the target ice number concentration was 2, 3, 5 or 6 L −1 (marked with ICE2, ICE3, ICE5 and ICE6 respectively). Figure   3 depicts the LWP and IWP evolution in all six UCLALES-SALSA simulations. The simulation time was extended to 24 hours in those cases where cloud still exists after 8 hours (marked with vertical line in Fig. 3). The simulation time was not extended any further, because we do not see any major changes or trends in the last simulation hours. Figure 3 shows that when the ice number concentration is set to a higher value, LWP decreases faster and cloud glaciates 195 sooner. In simulations ICE4, ICE5 and ICE6, the cloud dissipates totally after glaciation. The cloud glaciation happens because water vapour condenses on the ice crystals at the expense of the cloud droplets. In simulations ICE1, ICE2 and ICE3, IWP stabilises to values of approximately 6.5, 10 and 12 g m −2 respectively towards the end of simulation.
From Fig. 3 we can also see that LWP still increases during the first 8 hours with ICE2 but not anymore with ICE3. With These results show how sensitive the mixed-phase cloud is to ice number concentration either by showing how fast the 205 cloud glaciates or when balance is reached. However, these are highly simplified due to the lack or real aerosol dependent freezing and related feedback processes. These results also show the need for more detailed feedbacks since constant ice number concentration is not a realistic assumption for real clouds.

Prognostic ice simulation
One of the unique features of our model is its ability to keep track of the chemical composition along with sectional aerosol 210 size distribution also in cloud phase. This allows us to model freezing processes related to an ice nucleating compound like dust. Furthermore, parallel bins allow for analysing the relative contribution of e.g. dust particles (INP) on ice formation. We call this prognostic ice because here freezing probability is related to dust aerosol which mass and number concentrations are prognostic variables. We allow interactions between all hydrometeors and ice formation is modelled using the implemented ice nucleation theories, which relate ambient conditions and droplet properties to their freezing rates.

215
To see the difference between fixed and prognostic droplet freezing, we made a prognostic ice simulation that was targeted to have similar IWP during the first 8 hours as in the simulation with ice number concentration of 4 L −1 (ICE4) (see Sections 3.1 and 3.2). This ICE4 simulation was selected for comparison because it is close to the tipping point where cloud either stabilises or glaciates (see Sect. 3.2). number concentration and size distribution of the aerosol remain the same as in the fixed ice number simulations (Sect. 3.1 and 3.2), thus they are the same as in Ovchinnikov et al. (2014) . In the absence of more detailed aerosol observations, INP number concentration and mixing state, and contact angle were considered as adjustable parameters impacting ice nucleation ability.
Here, contact angle represents the angle between the ice embryo and the ice nucleus in an aqueous medium.
First, in order to set the INP number concentration, we incorporated b bins (For more information about bin description, 225 see Sect. 2 and Fig. 1). Proportion x = 0.015 of the total aerosol number concentration was partitioned in b bins as INPs.
Second, the INP mixing state was adjusted so that the particles in the b bins were set to have an insoluble dust core, 50% of the dry mass, and ammonium bisulphate for the other half. Here, dust acts as the INP.

230
Third, the freezing rate was adjusted by setting the cosine of the contact angle of dust to m is = 0.57 (Eq. A3 in Appendix A).
It should be noted that the target IWP could have been reached using different combinations of INP mixing state, x and m is but these simulations showed that the results depend mostly on the resulting ice number concentration rather than the applied parametrisation. These characteristics of aerosol are uniform throughout the whole simulation domain.

235
The simulation time for the prognostic ice run was set to 32 hours. The water paths of the mixed-phase cloud are quite stable after that. The simulation time of ICE4, used to compare with the prognostic run, was not extended any further from 24 hours since the cloud dissipates around 12 hours of simulation.
As in the ICE4 simulation, in the prognostic ice simulation, droplet freezing was set to start after a spinup of 2 hours. Figures   4a and 4b illustrate that the prognostic ice and ICE4 simulations have similar IWP and LWP during the first 8 hours. Hence, the 240 targeting is successful and the initial conditions of the simulations match each other. After that, the prognostic ice simulation diverges from the ICE4 simulation. Figure 4a shows that in the prognostic ice simulation LWP starts to increase after 4.5 hours of simulation. This is caused by a decrease of ice number concentration (Fig. 4c) to such a low level which allows more water vapour for condensation to liquid droplets. The same figure also depicts how the ice number concentration is set to a target value (simulation ICE4) and how 245 the concentration is stable until the cloud dissipates. Figure 4d depicts how droplet number concentration lowers especially right after spinup period when ice number concentration is increasing. However, changes in droplet number concentration is not the driving force behind complete removal of liquid phase. Figure 4e illustrates how the whole cloud with prognostic droplet freezing descends and how the ICE4 is affected by entrainment both below and above the cloud, cloud gets thinner and dissipates. In all simulations, droplet number concentration is specified as prognostic variable.

250
In the beginning of the prognostic ice run, domain mean of dust containing aerosols is approximately 27 L −1 . After 32 hours of simulation the same mean value is about 13 L −1 . Here, the loss of INPs limits the ice number concentration. The mixed-phase cloud persists because the ice number concentration can change. This is so-called self-adjustment of INPs which better reproduces observed evolution of mixed-phase clouds since usually they are more resilient in observations than in models (Andronache, 2017;Morrison et al., 2011a). This is also in line with previous modelling studies, where prognostic INP will 255 reduce the number of ice crystals because of precipitation, thus allowing cloud liquid to sustain (Fridlind et al., 2012;Solomon et al., 2015Solomon et al., , 2018. The decrease of dust (INP) mass concentration in different hydrometeor types is shown in Fig. 5. Dust is an efficient ice nuclei so it will soon end up in ice crystals which are removed from the system by sedimentation (Fig.   5c). Free troposphere is the only source for the boundary layer dust, and the relevant mechanisms are entrainment and large scale subsidence. Subsidence is described with a downward vertical velocity moving mass and energy. Entrainment in this 260 case describes any other kind of mass exchange between cloud top and free troposphere. For instance, subsidence is 0.004 m s −1 at the cloud top and aerosol number concentration in dust containing b bins above the cloud is about 24000 m −3 , so the dust aerosol flux from the free troposphere is approximately 100 m −2 s −1 . Because radiative cooling is strengthening the supersaturation at the cloud top, the most CCN active part of these entrained dust containing particles can be activated immediately as cloud droplets. This can be seen as a higher dust mass concentration within cloud droplets in the upper layer 265 of the cloud (Fig. 5b). If the temperature is low enough, these dust containing droplets will subsequently freeze during the following time steps and therefore takes part in preserving the mixed-phase cloud.
A more detailed examination of droplet activation and ice formation can be done by studying the time evolution of the size distribution. Figure 6 shows how different sized particles are partitioned between different hydrometeor types within the cloud layer. Figures 6c and 6d show how the larger particles freeze first and their number concentration decreases quickly as these 270 particles deposit at the surface within falling ice hydrometeors and are removed from the system. Even though the entrainment from above is providing more particles, this is not fast enough to maintain the original concentration. Removal of the smaller INPs is slower as those are less likely to activate as cloud droplets and the resulting droplets are also less likely to freeze due to the smaller dust core area. However, with time and because of continuous mixing of boundary layer eventually also the smaller INPs are able to form cloud droplet within the strongest updrafts and formed droplets will freeze within the cloud. This will 275 lead to stabilisation of aerosol size distribution. The increase in the total number of particles in bin 1 is a numerical artefact caused by the bin adjustment routine, which can move particles from one bin to another in order to keep the dry size within the predefined bin limits. When a large proportion of particles in bin 2 are activated as cloud droplets, some of the remaining are moved to the smaller bin to avoid numerical problems. However, this numerical artefact does not affect the results.  respectively (McFarquhar et al., 2011;Savre and Ekman, 2015). Ice number concentration is also approximately two orders 285 of magnitude less than the number concentration of efficient INPs above the cloud layer. From that we can estimate that the concentration of INPs entraining from above the cloud should be in the order of 0.1 to 1.0 cm −3 to glaciate the cloud. Figure 7c further illustrates an interesting behaviour of ice particle formation. In the beginning of the simulation ice particles are formed throughout the cloud, but later the most intensive formation takes place at the top of cloud where fresh INPs are entrained into the cloud layer. However, the maximum supersaturation in these entraining downdrafts is so low, that only the 290 largest particles are able to form cloud droplets and consequently freeze. The smaller ones penetrate through the cloud layer as interstitial aerosol particles (i.e. unactivated particle), and are able to form cloud droplets (i.e. activate) and ice particles at the cloud base when they are recirculated back to the cloud with higher supersaturation. This can be seen in Fig. 8. Figure 8a shows how in size bin 2 cloud droplets and ice particles are more frequent in updrafts compared to Fig. 8b which illustrates how aerosols are more favourable in downdrafts. Additionally, ice particles dominate in bigger sizes as aerosols freeze both in down-295 and updrafts (size bin 3 shown Figs. 8c and 8d). Simulated freezing in different vertical velocity conditions in other size bins does not differ from results shown already in Fig. 6. Lower peak at end of simulation in vertical profile of freezing rate in Fig. 7c also shows how recirculated aerosols are frozen in the cloud layer. Such phenomenon has been modelled before e.g. in Solomon et al. (2015), however here the cloud is simulated with explicit calculation of in-cloud supersaturation and representation of aerosol size distribution and chemical composition. If activation is not modelled with this level of details, activation and 300 freezing might happen too early or late and in a wrong part of the cloud. Overall, Figs. 6, 7c, and 8 nicely demonstrate how the relative proportions of particles in different hydrometeors are size dependent and how sectional description for aerosols is required to be able to simulate such processes in LES models.

Conclusions
In this study we have extended our large-eddy model UCLALES-SALSA (Tonttila et al., 2017) for ice and mixed-phase 305 clouds. The model has an exceptionally detailed sectional aerosol description for both aerosol number and chemical composition, which makes this model suitable for examining aerosol-cloud interactions and dynamics. Specifically, this allows the description of an ice nucleation active material such as mineral dust, which can be used in calculating ice formation rates from the nucleation theory.
As the first step, we compared our model predictions with those from a mixed-phase cloud model intercomparison study 310 (Ovchinnikov et al., 2014) to confirm the accuracy of the newly implemented ice microphysics. In this simplified model intercomparison setup, where any cloud droplet will freeze until a specified ice number concentration (from zero up to four particles per liter) is reached, the focus is on cloud dynamics. In agreement with Ovchinnikov et al. (2014) and several other studies (e.g., Klein et al., 2009;Morrison et al., 2011b;Stevens et al., 2018) we conclude that microphysical details such as the fact that dry particle size is tracked in UCLALES-SALSA while most other sectional models track the ice particle size have an 315 impact on predictions. Such details become more important close to the tipping point where the further addition of ice particles leads to the rapid glaciation of the cloud.
In the second part, we constructed a case where ice formation is modelled using a heterogeneous ice nucleation scheme and a prognostic ice nucleating particle population containing mineral dust. This so-called prognostic ice simulation was designed so that it matched with the previous fixed ice number concentration simulation where the cloud was close to the tipping 320 point. When the simulation with fixed ice concentration showed a complete glaciation after about 12 hours, the prognostic ice simulation reached an equilibrium state which lasted up to end of the 32 hour simulation. With this the prognostic simulation showed the importance of the self-adjustment of ice nucleation active particles. This is in good agreement with previous modelling studies (Fridlind et al., 2012;Solomon et al., 2015Solomon et al., , 2018 and a observational study where resilient mixed-phase clouds are seen together with relatively high ice nuclei concentrations (Filioglou et al., 2019).

325
Further examination of the prognostic ice simulation revealed that the efficient INPs entrained from the free troposphere are able to maintain the mixed phase cloud with ice particle number concentration on average 0.1-0.2 % of the INP concentration above the cloud. These entrained particles do not immediately form ice particles in the cloud top. The detailed analysis of the model outputs reveals how particle size and supersaturation-dependent cloud activation eventually control the formation of ice through immersion freezing. Part of entrained INPs penetrate through the cloud as interstitial particles, get mixed within bound-330 ary layer air and contribute ice formation later when recycled back to the cloud. Thus the entrainment process is maintaining INP concentration in the whole boundary layer.
This study emphasises the benefits of the detailed aerosol-cloud-ice module within a LES model. In fact, UCLALES-SALSA is one of the few cloud scale models (Fridlind et al., 2012;Khain et al., 2004;Savre and Ekman, 2015;Fu and Xue, 2017) where details about aerosol and cloud droplet chemical composition can be utilised by using the particle level theoretical 335 understanding about ice nucleation. The model will be a useful tool for the mixed-phase cloud research, which has started to attract more widespread interest.

Appendix A: Immersion freezing
The rate J of heterogeneous germ formation by immersion freezing is a function of T temperature in kelvins, r N radius of insoluble substrate and S w equilibrium saturation ratio at droplet surface based on Kohler theory and is determined as

345
where k and h Boltzmann's and Planck's constants, ∆F act is the activation energy at the solution-ice interface, ∆F cr is the critical energy of germ formation, C het is the normalizing function, r N the radius of an insoluble fraction of an aerosol particle (INP), and c 1s is the concentration of water molecules adsorbed on 1cm −2 of a surface (Eq. 2.1 Khvorostyanov and Curry (2004)). Used parameter values are C = 1.7 × 10.999850 10 N m −2 and c 1s = 1 × 10 19 m −2 .
where T is temperature in kelvin, R is the molar gas constant and N A is the Avogadro constant and parameter values B = 347, T * = 177, D * = 4.17 and D 0 = 349 for p = 1 bar are gained from Table 2 in Khvorostyanov and Curry (2004).
Critical energy is based on Eq. 2.10 in Khvorostyanov and Curry (2000) ∆F is surface tension between ice and solution, where T c is temperature in Celsius degrees (Khvorostyanov and K. Sassen, 1998); ice germ radius:

360
where ρ ice = 900 (kg m −3 is the density of ice, T 0 = 273.15 K; is the effective latent heat of melting (Eq. 6 Khvorostyanov and K. Sassen (1998)); dimensionless parameter where M w is the molar mass of water (Eq. 2.7 in Khvorostyanov and Curry (2000)).

365
The shape factor f is defined as a function of ratio x = r N /r cr and m = m is . It is gained from Eq. 2.9 in Khvorostyanov and Curry (2000), originally from Fletcher (1962), Case dependent parameters , elastic strain produced in ice embryo by the insoluble substrate; α, relative area of active sites; m is , cosine of the contact angle; are defined in our results 3.3 to be The m is was used as targeting parameter since the simulation tests were found to be very sensitive for this parameter. Other case dependent parameters and α were not altered and had their default values.
The number of crystals formed by homogeneous nucleation due to the freezing of supercooled pure water or deliquescent condensation nuclei is described by Eq. 1 in Khvorostyanov and K. Sassen (1998): where k and h Boltzmann's and Planck's constants, ρ w is density of water , ρ ice is density of ice (same as in A) and N c = 380 5.85e16m −2 is the number of water molecules contacting unit area of ice germ (Khvorostyanov and Curry, 2000).
Effective latent heat of melting L ef m is same as A6. Dimensionless parameter G is same as A7. Surface tension between ice and solution σ is is same as A4.

385
Ice germ radius is defined as which is same as A5 with = 0 (Khvorostyanov and Curry, 2000).
Hence we get the critical energy of germ formation (Khvorostyanov and K. Sassen, 1998, Eq. 9b)

Appendix C: Deposition freezing
Deposition freezing is possible for dry insoluble aerosol at sub-saturated conditions (RH < 100%). This is implemented following Khvorostyanov and Curry (2000) and additional parameters from Hoose et al. (2010).
The rate of germ formation J (s −1 ) through deposition freezing is defined as in Khvorostyanov and Curry (2000, Eq 2.13) The pre-exponential factor (kinetic coefficient) is about where r n is the radius of the insoluble portion of the droplet, k is Boltzmann's constant, T is temperature. The pre-exponential factor (kinetic coefficient) is 10 26 (cm −2 )r 2 n (Fletcher, 1962). Here case-dependent activation energy ∆F act is set to zero (Khvorostyanov and Curry, 2000).
Ice germ radius r cr (Khvorostyanov and Curry, 2000, Eq. 2.12) and defined as where S i is water vapor saturation ratio over ice, T is temperature, C is constant 1.7 × 10 10 (N m −2 ).
From previous values we get the critical energy of germ formation : where the shape factor f (m, x) is defined as in A8.