Simulating mixed-phase Arctic stratus clouds: sensitivity to ice initiation mechanisms

Abstract. The importance of Arctic mixed-phase clouds on radiation and the Arctic climate is well known. However, the development of mixed-phase cloud parameterization for use in large scale models is limited by lack of both related observations and numerical studies using multidimensional models with advanced microphysics that provide the basis for understanding the relative importance of different microphysical processes that take place in mixed-phase clouds. To improve the representation of mixed-phase cloud processes in the GISS GCM we use the GISS single-column model coupled to a bin resolved microphysics (BRM) scheme that was specially designed to simulate mixed-phase clouds and aerosol-cloud interactions. Using this model with the microphysical measurements obtained from the DOE ARM Mixed-Phase Arctic Cloud Experiment (MPACE) campaign in October 2004 at the North Slope of Alaska, we investigate the effect of ice initiation processes and Bergeron-Findeisen process (BFP) on glaciation time and longevity of single-layer stratiform mixed-phase clouds. We focus on observations taken during 9–10 October, which indicated the presence of a single-layer mixed-phase clouds. We performed several sets of 12-h simulations to examine model sensitivity to different ice initiation mechanisms and evaluate model output (hydrometeors' concentrations, contents, effective radii, precipitation fluxes, and radar reflectivity) against measurements from the MPACE Intensive Observing Period. Overall, the model qualitatively simulates ice crystal concentration and hydrometeors content, but it fails to predict quantitatively the effective radii of ice particles and their vertical profiles. In particular, the ice effective radii are overestimated by at least 50%. However, using the same definition as used for observations, the effective radii simulated and that observed were more comparable. We find that for the single-layer stratiform mixed-phase clouds simulated, process of ice phase initiation due to freezing of supercooled water in both saturated and subsaturated (w.r.t. water) environments is as important as primary ice crystal origination from water vapor. We also find that the BFP is a process mainly responsible for the rates of glaciation of simulated clouds. These glaciation rates cannot be adequately represented by a water-ice saturation adjustment scheme that only depends on temperature and liquid and solid hydrometeors' contents as is widely used in bulk microphysics schemes and are better represented by processes that also account for supersaturation changes as the hydrometeors grow.

These glaciation rates cannot be adequately represented by a water-ice saturation adjustment scheme that only depends on temperature and liquid and solid hydrometeors' contents as is widely used in bulk microphysics schemes and are better represented by processes that also account for supersaturation changes as the hydrometeors grow.

Introduction
The surface energy budget over the Arctic ice pack is determined to a large extent by radiative fluxes that in turn are strongly dependent on the presence of clouds. Lowlevel Arctic clouds contribute about the half of the total cloud fraction throughout the year due to their persistence and horizontal extent (e.g., Curry and Ebert, 1992). The vertical structure and radiative properties of these persistent low-level clouds depend on their microphysics, and thus estimation of the relative significance of the microphysical processes that occur in these clouds is important.
An efficient method to investigate the role of different microphysical processes that determine the microphysical and radiative cloud properties is the utilization of multidimensional cloud models with bin resolved microphysics (BRM). There are many BRM cloud models that are "warm" rain models (e.g., Cotton, 1972a;Ogura and Takahashi, 1973;Clark, 1973;Tzvion et al., 1989;Bott et al., 1990;Kogan, 1991;Kogan et al., 1995;Stevens et al., 1996). To account for the ice phase several BRM cloud models use one size distribution function to describe cloud ice (e.g., Hall, 1980;Sednev and Khain, 1994;Khvorostyanov and Sassen, 1998;Ovtchinnikov and Kogan, 2000a). In these models only one type of solid hydrometeors that is ice crystals is considered or the bins with the smallest ice sizes are assigned to ice crystals while the rest of bins are interpreted as graupel. The ability of these models to simulate realistically microphysical processes in mixed-phase clouds is limited to situations when the processes of precipitation formation do not play a key role. It is difficult to expect that solid cloud hydrometeors, whose bulk densities and terminal velocities vary widely (Macklin, 1962), can be described adequately by one distribution function. Only a few BRM models use designated distribution functions for different types of ice hydrometeors and calculate growth rates of microphysical processes due to several transformations of liquid and solid hydrometeors in mixed-phase clouds (Cotton, 1972b;Young, 1974;Scott and Hobbs, 1977;Chen and Lamb, 1994;Khain and Sednev, 1996;Reisin et al., 1996a;Takahashi and Shimura, 2004). As a rule these models also include a budget equation for the size distribution function for aerosol particles (AP) that can be of different chemical composition. The values of supersaturation calculated in the course of model integration are used to determine the size of APs to be activated, and the corresponding spectrum of cloud droplets just nucleated are directly calculated. It is thought that despite high computational costs these mixed-phase cloud microphysical models provide more accurate simulations of cloud-aerosol interactions and processes of precipitation formation in mixed-phase clouds than models with simpler ice schemes (Lynn et al., 2005).
Although Arctic mixed-phase low level clouds due to their existence throughout much of the year have important climatic impacts, the number of studies, in which BRM models are used for the investigation of microphysical characteristics of these clouds, is quite limited. Using a mixed-phase BRM scheme (Reisin et al., 1996a), which utilizes three distribution functions for the ice phase (crystals, snow, and graupel), coupled to a 2D cloud resolving model, Harrington et al. (1999) studied model performance for idealized situations that mimic environmental conditions typical for the transition (spring and fall) Arctic season. A set of sensitivity runs was performed to reveal the impact of different microphysical processes on glaciation time, longevity, stability, and radiative properties of simulated mixed-phase clouds. It was shown that environmental conditions have a strong impact on modeled cloud properties.
Recently an advanced 3D LES model was used to quantify the role of different ice formation mechanisms in mixed-phase stratocumulus observed during MPACE . The BRM scheme utilized in this model includes sophisticated equations for ice nuclei (IN) that can be activated in the heterogeneous modes (contact, deposition, condensation, and immersion nucleation), formed due to water drop evaporation and scavenged by water droplets. Additional ice origination mechanisms such as rime splintering, drop freezing during evaporation, shattering of drops during freezing, and fragmentation during ice-ice collisions are also considered. The model is able to reproduce persistent mixed-phase stratocumulus cloud decks as well as cloud microphysical properties (liquid and ice water content, droplet, and ice nuclei concentration profiles) within the observed ranges for particular combinations of ice formation mechanisms mentioned above. They found that glaciation time and longevity of mixed-phase MPACE clouds are determined by formation of ice nuclei due to water drop evaporation and drop freezing during evaporation, whereas processes of ice multiplication were less important.
There is a broad consensus that the ice initiation process is of crucial importance for the longevity of mixed-phase clouds. To study the ice initiation processes (IIP) through nucleation from water vapor and transformation of super-cooled liquid water, as well as the transformation of water vapor due to condensation/deposition, evaporation/sublimation, and the Bergeron-Findeisen process (BFP) in Arctic mixed-phase clouds, we use the BRM scheme Sednev, 1995, 1996) coupled to the GISS SCM (Menon et al., 2003) called the GISS-LBL SCM. In our simulations of singlelayer stratiform mixed-phase clouds observed during the DOE ARM Mixed-Phase Arctic Cloud Experiment (MPACE) Intensive Observing Period (IOP) in October 2004 at the North Slope of Alaska Verlinde et al., 2007) with the GISS-LBL SCM, we consider two mechanisms of ice initiation. If liquid phase is not involved in IIP, we parameterize nucleation of ice crystals from water vapor as a function of supersaturation w.r.t. ice. Otherwise, ice crystal origination is considered to proceed via drop freezing, and its rate is a function of the shape of droplet distribution, water droplet mass and temperature. Once nucleated small crystals grow rapidly due to deposition/BFP. To treat the BFP bulk microphysics schemes use various modifications of the "saturation adjustment" assumption that mainly only temperature dependent and are not account for hydrometeors' shapes and size distributions. As opposed to bulk microphysics schemes that use this oversimplified approach and, in fact, are not able to treat the BFP process adequately, the BRM scheme uses analytical solutions to solve equations for supersaturation w.r.t. water (SSW) and ice (SSI) taking into account the hydrometeors' size distributions, densities, and shapes. Moreover, the BRM scheme takes into account supersaturation (SS) changes for the liquid/solid particle growth equations during the microphysical time step, thus providing a better representation of the BFP.
We describe the GISS-LBL SCM bin-resolved microphysics in Section 2. Simulation setup is given in Section 3. The results of several sensitivity experiments and comparison to observations are presented in Section 4, and a summary and discussion are provided in Section 5. Finally, in Appendix A we outline some details of the BRM scheme that are relevant to this study.

Model Description
In this study we use the GISS-LBL SCM that is a modified version of the GISS SCM (Menon et al., 2003) adapted from the GISS GCM. The cloud physics parameterizations in the GISS GCM is based on an assumption that only liquid or ice phase is permitted to exist at temperatures below freezing (Del Genio et al., 1996. This cloud scheme also has limited ability to represent cloud-aerosol interactions, relying on diagnostic calculations of cloud droplet concentration as a function of aerosol mass (Menon et al., 2002). To better account for effects of environmental conditions and microphysical processes on cloud development and persistence several modifications to the GISS SCM have been done. These modifications include: 1) redesign of numerical algorithms used in the turbulence scheme and reformulation of the numerical surface boundary conditions; 2) implementation of a bin resolved microphysical (BRM) scheme that was specially designed to represent mixed-phase clouds. The BRM scheme takes into account numerous microphysical processes, some of which (IIP and BFP) are of special interest in this study.

Equations for size distribution functions
The modified BRM scheme originally developed by Sednev (1995, 1996) directly solves integro-differential equations for or mass (size) distribution functions describing seven types of hydrometeors: f 1 is the mass distribution function for water droplets and rain drops; f 2 , f 3 , and f 4 are the size (mass) distribution functions for columnar crystals, plate crystals, and dendrites, respectively; f 5 , f 6 , and f 7 are the mass distribution functions for snowflakes, graupel particles, and frozen drops/hail and does not use assumptions regarding the shape of distribution functions. The distribution functions are defined on mass grids that can contain different number of bins. The changes in distribution functions for liquid and solid phases are governed by the following equations: where k = 1. . .7 denotes the type of hydrometeor, V k is its terminal velocity, w is vertical velocity, and ρ a is air density. [∂f k /∂t] act/nucl are the rates of changes of f k due to activation/nucleation processes; [∂f k /∂t] cond/ev are the rate of changes of f due to 7 the condensational growth or evaporation of droplets (for k = 1 ) or due to deposition/sublimation of ice particles (for k > 1). [∂f k /∂t] coal are the rates of change of f k due to coalescence between hydrometeors of any type including type k; [∂f k /∂t] f rz/mlt and [∂f k /∂t] brk are the rates of change of f k due to the freezing of droplets and melting of ice particles and breakup processes; [∂f k /∂t] mult describes ice multiplication process, and operator ∆() denotes the contribution of small scale turbulence. The BRM scheme provides calculation of precipitation amount, concentration, mass contents and precipitation fluxes of different hydrometeors, radar reflectivity from water and ice, the mean and effective radii of droplets and ice particles as well as provides information for calculation of cloud optical properties such as single scatter albedo, optical depth and asymmetry parameter. For example, concentrations n k , water/ice contents q k , and precipitation fluxes P k of hydrometeors are determined by means of distribution functions as follows: In the above, f k are given in number of particles per kg of water in kg of air, and n k , q k , and P k are in number of particles in kg of air, kg of condensate in kg of air, and kg of condensate per m 2 per second, respectively.

Initiation of liquid phase
Liquid phase initiation parameterization in BRM scheme is based on solving the supersaturation equations, which predict SSW and SSI, Köhler theory, and a prognostic spectrum of aerosol particles that can be of different chemical composition. The equation for distribution functions for AP f ccn that is defined on separate mass grid m c is as follows: where [∂f ccn /∂t] act is the rate of change of f ccn due to activation. Sedimentation velocities of AP V ccn (m c ) are set to be zero, and wet removal of AP is not considered in this study. Aerosol particles of a certain size can be activated when the supersaturation calculated exceeds the critical value determined by the Köhler equation. Köhler theory is also used to calculate so called critical and equilibrium aerosol radii. In the case the AP distribution contains aerosols with dry radius greater than critical radius at a given point, these APs can be activated and transformed into droplets. The size of new nucleated droplets is equal to equilibrium radius if dry AP radii is less than 0.03 microns; otherwise, the radii of cloud droplets are five times as much as dry AP radii since large CCN does not reach their equilibrium sizes. This approach prevents origination of unrealistically large droplets and too fast warm rain formation (Kogan, 1991;Khain and Sednev, 1996). This droplet nucleation scheme starts with activation of the APs, whose critical supersaturation is the smallest, calculates corresponding droplet sizes and liquid content increase, and assures that total vapor and cloud water content is conserved and there is still enough water vapor exceed to activate APs in the next bin.

Initiation of solid phase
Solid phase initiation parameterization in BRM scheme accounts for two general mechanisms distinguishable according to the involvement of liquid phase in the ice generation process.
If liquid phase is not involved in ice initiation process, we parameterize nucleation of ice crystals from water vapor as a function of SSI (Meyers et al., 1992: where A ms and B ms are set to -0.639 and 12.96, respectively, N ms is ice nuclei (IN) concentration in 1/L (Table 2), and N mc determines the upper limit of concentration, up to which ice crystals can be nucleated from water vapor at a particular point. Nucleation is not permitted if concentration of ice hydrometeors is already greater than that determined by (6). Only the number of crystals needed to reach the concentration given by (6) is nucleated. All ice crystals just nucleated have the minimum size permitted by mass grid, and basic crystal habits, which depend on temperature (Takahashi et al., 1991), are plates (-4 If liquid phase is involved in ice initiation process, water drop freezing that contribute to ice initiation is considered separately, and is treated as a main mechanism of primary graupel formation. The freezing probability is given by: A similar dependence was used in Takahashi (1976), Alheit et al. (1990), and Khain and Sednev (1995). In Eqs. (7) m is the drop mass, T sup =T f -T is the degree of supercooling, T f = 273.16 K is temperature threshold value, B f = 0.66 K −1 , and values for arbitrary constant A f are in Table 2. It is arbitrarily assumed that frozen droplets with radii less than 100 µm are transformed into plate crystals with the density of 0.9 g cm −3 , and drops with greater radii become graupel particles. It should be noted that opposite to the parameterization of the activation of aerosol particles to form cloud droplets that depends on aerosol chemical composition (among others), parameterization of solid phase initiation does not depend on chemical composition of IN. To account for this effect two-dimensional size distributions functions are needed (Bott, 2000). In this approach the particles are classified according to their water and total aerosol mass on a two-dimensional grid. For example, Diehl and Wurzler (2004) and Diehl et al. (2006) studied heterogeneous drop freezing employing the two-dimensional treatment of cloud physics that allows the coexistence of similarly sized drops with different contents of soluble and insoluble particles. Because the two-dimensional approach is very computationally expensive, in our BRM scheme the one-dimensional approach is used, and freezing probability of equally sized droplets remains the same and does not depend on the origination mechanism.
The first IIP is only active in supersaturated w.r.t. water environment, whereas the second IIP operates in both saturated and subsaturated (w.r..t. water) environment. Even if liquid particles evaporate and ice particles sublimate, this transformation takes place. If SSW > 0 transformation of liquid phase into ice phase is accompanied with droplet condensational and ice particle depositional growth. The same transformation occurs when droplets are evaporated supplying additional water vapor for ice particle depositional growth (Bergeron-Findeisen process).

Treatment of Bergeron -Findeisen process
The parameterization of the Bergeron -Findeisen process is a challenging problem because of the necessity to partition water vapor excess between solid and liquid phases. For this we have to answer a question on how we can calculate the amount of water vapor, water, and ice if their total amount, temperature, and pressure are known. Usually this partitioning is assumed to be only a function of temperature or "saturation adjustment" schemes of different degrees of sophistication are used in atmospheric models of different scales (Lord et al., 1984;Tao et al., 1989;Ferrier, 1994;Fowler et al., 1996;Del Genio et al., 1996). Advanced schemes (Rotstayn et al., 2000;Morrison et al., 2005) solve growth equation for ice particles using some additional assumptions. Essential features of the Rotstayn et al. (2000) scheme are (a) saturation w.r.t. water; 11 and (b) supersaturation w.r.t ice is constant during a microphysical time step. The assumption (a) works well for situations when only liquid phase exists. To some extent assumption (b) works for only icy clouds. If both assumption are used simultaneously it leads to the following consequences: 1) liquid phase becomes completely decoupled from ice phase (the amount of water vapor deposited on ice particles will be the same if LWC is equal to 0.001 g kg −1 , 1 g kg −1 or 10 g kg −1 ; 2) BFP parameterization might be valid for environmental conditions without liquid phase. Such kind of parameterizations are able to treat only icy clouds. Moreover, it contradicts the definition of BFP process ("depositional growth of ice particles on expense of evaporated cloud droplets). It's well known that in mixed-phase clouds supersaturation w.r.t water and ice might coexist at the same time. In this situation we have condensational growth of droplets and depositional growth of ice particles at the same time but not BFP. In Morrison et al. (2005) equation for SSI that account for coexistence of liquid and solid phases is solved, but equation for ice phase growth rate due to sublimation/deposition relies on assumption that supersaturation w.r.t ice remains unchanged during the microphysical time step. It might lead to excessive ice phase growth rates and the atmosphere drying at different altitudes (especially those where the temperature varies in a the range -15 • C < T < -10 • C, in which the differences between saturated water pressure w.r.t. water and ice have a maximum).
To improve BFP representation in BRM scheme we utilize a treatment for condensational processes in mixed-phase clouds considering different growth regimes: (1) condensation for liquid phase (SSW > 0) and deposition for ice phase (SSI > 0); (2) evaporation for liquid phase (SSW < 0) and sublimation for ice phase (SSI < 0); and (3) Bergeron-Findeisen process or evaporation for liquid phase (SSW < 0) and deposition for ice phase (SSI > 0). This treatment accounts for effect of hydrometeors size distribution as well as supersaturation changes during the time step on liquid phase growth rate due to condensation (evaporation) and ice phase growth rates due to deposition (sublimation). Because of the significant impact of the BFP on mixed-phase clouds glaciation time, we outline the numerical procedure that is used to calculate the condensation / deposition and evaporation / sublimation rates of liquid and solid hydrometeors in saturated/subsaturated (w.r.t. water/ice) environments in Appendix A.
It should be noted that this approach is applicable in any cloud model that resolves supersaturation. It can be used in bin-resolved and bulk microphysical schemes. In the latter case it can be significantly simplified using prescribed a priori liquid/solid phase size distributions. The approach outlined in the Appendix can be used for the numerical modeling of cloud processes using cloud resolving models and developing parameterization of processes of vapor/liquid/solid phase transformations for use in large-scale models.

Simulation setup
The initial vertical profiles used to drive the SCM (36 levels with 25 mb resolution near the surface) are given by the idealized profiles from observations during the MPACE IOP. We use the large-scale forcing, horizontal velocity components, subsidence velocity, surface pressure, temperature, and fluxes as defined by Klein et al. (2007). We focus on Period B (17Z October 9 to 5Z October 10) when single-layer stratocumulus mixed phase clouds with temperatures varying between -5 • C and -20 • C were observed. These clouds are characterized by persistent liquid phase with liquid water content (LWC) that increases with altitude reaching its maximum at the top of a well mixed boundary layer. Ice phase exists both in clouds and under the liquid cloud base.
The typical values of liquid water path (LWP) and ice water path (IWP) are 200 g m −2 and 20 g m −2 , respectively. A detailed description of the environmental conditions and microphysical characteristics of observed clouds can be found in Klein et al. (2007).
In all our 12-hour simulations (from 17Z October 9 to 5Z October 10) the pressure P s near the surface is 1010 hPa with constant sea surface temperature T s = 0.85 • C. Surface sensible and latent heat fluxes are 138 W g m −2 and 108 W g m −2 , respectively, and vertical profiles of horizontal velocities are also prescribed (Klein et al., 2007). There is no ice phase at all altitudes, and no liquid phase exists above the inversion 13 (P inv = 850 hPa) at the initial time. Idealized vertical profiles of total water mixing ratio q t = q + q w and liquid water potential temperature Θ l are defined as (Klein et al., 2007): where A q , B q , C q , P q and A t , B t , C t , P t are set to be equal to 1.950 g kg −1 , 0.291 g kg −1 , 0.00204 g kg −1 hPa −1 , 590 hPa and 269.20 K, 275.33 K, 0.07910 K hPa −1 , 815 hPa, respectively. Adiabatic LWC q w , vapor content q, and potential temperature Θ derived from Eqs. (8)-(9) are shown in Fig. 1. The initial bimodal distribution of dry APs are assumed to be composed of ammonium sulfate and independent of altitude. We do not simulate cloud origination and development, but use idealized vertical profile of LWC and droplet effective radii from flight measurements (Klein et al., 2007) to initialize the BRM. These characteristics as well as the derived vertical profile of droplet concentration and droplet size distributions at different altitudes at the initial time are shown in Fig. 2. In all our runs we use the BRM scheme with a 10 s time step. With this small time step, prescribed dynamics, and idealized forcing there is no significant spin up time, and we account for all data obtained during the course of our 12-hour simulations.
We perform four sets of simulations for the "warm" and "ice" microphysics cases. We present the microphysical cloud properties obtained from the sensitivity experiments described in Table 1 and Table 2 for simulations with only "warm" and "ice" microphysics, respectively, as moments (concentration, content, and effective radius) of size distribution functions for liquid/solid hydrometeors. In simulations with "warm" microphysics, listed as W1, W2, W3, and W4, "ice" microphysics is inactive, and we switch on/off processes of water-water interactions (coagulation) and processes responsible for changes of initial CCN distribution as described in Table 1. Differences between 14 W1 and W2 indicate the effects of coagulation and those between W1 and W4 indicate the effects of changes in the CCN distribution when coagulation is switched off. W3 includes both coagulation and an updated CCN distribution. For simulations with "ice" microphysics switched on, listed as I1, I2, I3, and I4, we use different ice crystal origination rates due to the IIP under consideration as described in Table 2. For I1 and I2, the first IIP is active, and the IN concentrations differ by a factor of 10 as shown in Table 2. Only the second IIP is active in I3, and both IIP are active in I4. Both I2 and I4 have the same IN concentrations. For similar IN concentrations, differences between I2 and I4 indicate the relative effect of the second IIP.

Results
As described in the previous section we perform a series of simulations to evaluate the impact of idealized forcing on modeled SS, the importance of the CCN spectrum shape for droplet activation and ice initiation processes. These runs are described in Table 1 and Table 2. We compare cloud microphysical properties from these simulations with observed values as obtained from MPACE IOP shown in Table 3  . These simulations are discussed below.

Sensitivity runs with warm microphysics
We perform "warm" sensitivity runs to ensure that simulated cloud characteristics for liquid phase (droplet concentration, content, and effective radius) are in range with observations. In addition, these runs permit us to evaluate influence of initial setup (temperature, water vapor mixing ratio, and microphysical characteristics vertical profiles and bimodal distribution of dry aerosols) and applied forcing (temperature and water vapor mixing ratio tendencies, prescribed subsidence and surface sensible and latent fluxes) proposed for MPACE intercomparison study (Klein et al., 2007) on simulated clouds. Here we describe our results for simulations without ice microphysics. Table 4 shows the average values of liquid-phase microphysical properties during the course of model integration. Figure 3 shows the SSW for these experiments. At altitudes where SSW is negative, no activation of new cloud droplets is permitted, and cloud droplets instantly evaporate and sediment due to their own terminal velocities and applied large-scale subsidence at all levels. At altitudes where SSW is positive, activation of new cloud droplets can occur. The BRM droplet activation scheme is sensitive to modeled SS that determines critical CCN radius, which is the cut off radius for the CCN spectrum, and the number of droplets just nucleated.
In W1 coagulation is switched off, and droplet activation, condensation, evaporation, and sedimentation are the only active microphysical processes. Droplet activation at a particular level mainly occurs when SSW exceeds its value at the previous time steps because if an activation event takes place, the corresponding bins in the CCN spectrum are likely to be empty. In all our experiments we do not model processes of new AP formation as well as their growth due to condensation/coagulation. There is no AP supply due to the large-scale horizontal processes in W1 and the only physical mechanisms that supply AP at a particular altitude are large-scale subsidence and vertical turbulent diffusion. The implied large-scale tendencies of temperature and water vapor mixing ratio together with the prescribed subsidence velocity result in mainly negative tendencies of SSW in cloudy regions. The balance between tendencies, turbulence, and radiation are such that SSW rarely increases, and the critical SS and CCN critical radius remain practically unchanged. This means that the amount of water droplets just activated is negligibly small. Figure 4 shows the droplet concentration N w , and Figure 5 shows the LWC for all the "warm" microphysics simulations. Both droplet concentration and LWC diminish with time due to sedimentation and evaporation at all levels during the first six hours in W1. After this time, in sub-cloud layers SSW becomes positive due to the instantaneous vapor supply from the surface and droplet evaporation just below the initial cloud base. Starting with the lowest layer and propagating upward, SSW remains positive determining the existence of non-dissipated warm clouds near the surface. In these clouds droplet effective radii R ew shown in Fig. 6 and droplet pre-cipitation flux P w (not shown) reach about 30 µm and 2.2 mm d −1 , respectively, and maximum droplet concentration and LWC never exceed their initial values (Fig. 2).
As in W1, in simulation W2, in which coagulation is active, N w and LWC have maximum value at the initial time (Fig. 2) and diminish with time as can be seen in Fig. 4 and Fig. 5. However, the process of rebuilding of SSW starts early, and SSW reaches very high values (about 3.5 %, Fig. 3) because coagulation effectively reduces droplet concentration. The "warm" rain formation process determines the increase in R ew (Fig. 6) and P w (not shown), whose average values are about 28.5 µm and 0.6 mm d −1 , respectively ( Table 4). The R ew values are significantly greater than those in W1. To prevent unrealistically high values of supersaturation and very short glaciation time in experiments with ice microphysics, we update the CCN spectrum after each time step with its initial values assuming that air masses with similar aerosol properties travel through the domain considered.
Supersaturation and microphysical characteristics (N w , LWC, and R ew ) for runs without "ice" microphysics using the CCN spectrum update assumption, W3 (with coagulation) and W4 (without coagulation) are shown in Fig. 3 and Fig's. 4 -6, respectively. In both W3 and W4 there are no areas of largely positive SSW as in W2 Fig. 3). At the same time average N w in W3 and W4 reaches 80 cm −3 and 66 cm −3 , respectively, as compared to 4 cm −3 and 2 cm −3 in W1 and W2, respectively. With coagulation turned on, as in W3, precipitation flux is reduced, and the LWC is higher compared to W4 in which coagulation is turned off.
Although the measurements of cloud droplets by a one-dimensional cloud probe (1DC, 20 − 640 µm maximum particle dimension) show drizzle development at the top of some of the MPACE single-layer clouds , MPACE observations indicate that the spectrum of water droplets remains relatively narrow, and there is no remarkable precipitation during October 10-12. Cloud microphysical values for W3 and W4, shown in Table 4, are in better agreement with observations (Table 3) compared to those obtained for W1 and W2. For example, W3 and W4 have R ew average values that are within the R ew observed range of 9 to 10.9 µm, shown in Table 3. How-ever, compared to observations, W1 and W2 overestimate average R ew values by a factor of 2 and 3, repsectively. Average values of N w and LWC for W1 and W2 are severely underestimated compared to observed ranges of 23 to 72 cm −3 and 154 to 193 mg m −3 for N w and LWC, respectively. On the other hand, values of N w for both W3 and W4 are within uncertainties in observations for N w and LWC for W4 is closer to the observed range in LWC than are values simulated for W3. These results indicate that regardless of the warm rain formation process, the CCN spectrum update assumption is crucial to maintain a persistent liquid phase with values of LWC that are comparable with observations. Based on numerous supplemental runs (not shown) and analysis of observation we conclude that water-water interactions (coagulation) is relatively minor if we believe that applied forcing and AP distribution are typical for MPACE period B conditions. Based on differences between the four sets of simulations shown in Table 4 and observations shown in Table 3 we suggest that the CCN spectrum shape for droplet activation is more important than is the process of water-water interactions (coagulation). Moreover, because of small concentration of ice particles and significantly reduced collision efficiencies for water-ice and ice-ice interactions at low temperatures, we conclude that coagulation is relative unimportant as compared to the BFP.
All these validate to some extend the assumptions used in experiments with "ice" microphysics that processes of water-water, ice-water, and ice-ice interactions may be relatively minor for the MPACE single-layer mixed-phase clouds (coagulation is switched off).

Sensitivity runs with ice microphysics
To evaluate the impact of the rates of the different IIP's on single-layer cloud evolution, we perform a set of runs I1, I2, I3, and I4 with ice microphysics Table 2. As shown in the previous section, simulations with the updated CCN spectrum, (W3 and W4) show more realistic cloud properties than do those without the updated CCN spectrum (W1 and W2). In the runs with ice microphysics we restore the CCN spectrum to its initial 18 values after each droplet activation event to prevent the cloud glaciating in unrealistically short time-scales. Tables 5 -6 show the average values of cloud microphysical properties for droplets and individual ice crystals, and Fig's. 7 -8 show SSW and SSI evolution for runs with ice microphysics. In these runs we consider two mechanisms of ice initiation. The fundamental difference between the two ice origination processes is the involvement of the liquid phase in the IIP. If the liquid phase is not involved in the IIP, we parameterize the origination of ice crystals from water vapor as a function of the SSI as shown in Eq. (6). It is assumed that this function provides the maximum concentration, up to which ice crystals can be nucleated at a particular point. We assume that all ice crystals just nucleated, whose shape (plates, columns, or dendrites) depends on temperature, have the minimal size permitted by the mass grid (of about the average values associated with a cloud droplet of 2 µm). This process operates for temperatures T< -2 • C. When the liquid phase is involved in the IIP, ice origination is considered to proceed via drop freezing, and its rate is a function of the shape of the droplet distribution, droplet mass and temperature, as shown in Eq. (7). Nucleated ice crystals of different sizes are assumed to be plate-like crystals. This process is active at negative temperatures in both saturated and subsaturated (w.r.t. water) conditions.
In all experiments with ice microphysics I1, I2, I3, and I4 the implied forcing assures the existence of high (∼ 20%) SSI (Fig. 8), and crystals thus formed grow rapidly reaching sizes of hundreds microns due to deposition and the BFP in mainly subsaturated (w.r.t. water) environments (Fig. 7).
In I1, only the first IIP is active. Figure 9 shows cloud droplet concentration N w , and Figure 10 shows LWC for the liquid phase for all simulations I1, I2, I3, and I4. Figures  11 -13 show the microphysical properties for the ice phase (concentration N i , IWC, and effective radius R ei ) for the same simulations. For I1, N w (Fig. 9) and LWC (Fig. 10) have maximum values of 73 cm −3 and 468 mg m −3 at the initial time (Fig. 2) and are continuously diminished due to evaporation and the BFP. Cloud glaciation time is ∼ three hours in I1. Simulated fields of SS show that initially intensively glaciated clouds con-tinue their development as icy clouds in sub-saturated (w.r.t. water) conditions (Fig. 7). Ice-phase concentration and content have maximum values in this experiment as compared to other experiments with ice microphysics as shown in Fig's. Fig. 11 and Fig. 12. For example, the maximum value of the concentration/content is 9.2 L −1 / 11 mg m −3 and 7.4 L −1 / 152 mg m −3 for plate and dendrite crystals, respectively. We note that the total ice-phase concentration, content and effective radii in I1 are significantly higher than those observed during MPACE.
As in I1 only the first IIP is active in I2, but the maximum concentration of ice crystals which can be nucleated at a particular point for the same SSI is reduced by an order of magnitude due to the assumption made for N ms , which is an order of magnitude smaller than that used for I1. As a result, the liquid phase in I2 exists for the course of model integration (12 hrs) supplying water vapor due to droplet evaporation for ice crystal depositional growth. Activation of new droplets also takes place because the maximum value of N w is greater that its value at the initial time. N w (Fig. 9) and LWC ( Fig. 10) have average values of 29 cm −3 and 123 mg m −3 , respectively (Table 5). LWC is reduced by an order of magnitude in ∼ nine hours (Fig. 10). The maximum ice concentration is 7 times less in I2 than in I1, and R ei in I2 is larger compared to that in I1 (Fig. 11).
I3, in which only the second IIP is active, is characterized by persistent liquid phase with maximum values of droplet concentration (Fig. 9) and LWC ( Fig. 10) near cloud top, significantly higher crystal concentration, and minimum values of ice precipitation flux (not shown) as compared with I2. Crystals effective radii R ei in I3 also have minimum values (Table 5).
Both IIP are active in I4, that combines some microphysical features of I2 and I3. Its main features are reduced droplet concentration and LWC as compared to I3 (Table 5), increased ice concentration and reduced effective radii as compared to I2, with about the same precipitation fluxes (not shown) for both runs (Table 6). I4 also agrees qualitatively with M-PACE data  that show the typical vertical structure of single layer clouds: existence of mainly liquid and ice phases at cloud top and near cloud base, respectively, with mixed phase in the middle of cloudy region. We expect that the relative importance of the second IIP will increase for long-lasting Arctic stratocumulus clouds within the temperature range -5 • C and -20 • C in less supersaturated (w.r.t. ice) environments than used in our runs.
It should be noted that using MPACE ice nuclei measurements Prenni et al. (2007) (P2007) proposed non-temperature dependent formula that has the same functional form as formula (2.4) in M1992, which we use to parameterize nucleation of ice crystals from water vapor (our formula (6)). Values for N ms , A ms , and B ms used in P2007 are 1.0 1/L, -1.488, and 0.0187, respectively. In our sensitivity experiments we use different values for N ms (Table 2). For example, these value are 1.0 1/L and 0.1 1/L for runs I1 and I2, respectively. N mc values calculated using our formula (6) in I2 and P2007 formula coincide (0.285 1/L and 0.288 1/L, respectively) when supersaturation w.r.t. ice is about 13%. For SSI between 2% − 13% and 13% − 25% these two formulae provide mean values equal to 0.154 1/L and 0.260 1/L and 0.663 1/L and 0.320 1/L, respectively. Mean N mc ratios for our formula (1) as used in I2 and P2007 formulation are 0.591, 2.071, and 1.506 for SSI between 2% − 13%, 13% − 25%, and 2% − 25%, respectively.
Both formulae give N mc that has the same order of magnitude for wide SSI range if N ms is set equal to 0.1 1/L in our formula (6). It can be demonstrated on P2007 Fig.  3 by parallel shifting of blue line that depicted M1992 formulation. N ms is set equal to 0.1 1/L in our sensitivity runs I2 and I4 (Table 2). If formula (6) were replaced by P2007 formulation in these runs, we would expect increasing significance of the second IIP. We would not expect to get dramatical changes if P2007 formulation were used. The relative significance of expected differences is determined by how often simulated SSI is less than 13% (P2007 formulation overpredicts observed MPACE values) and greater than 13% (our formulation overpredicts observed MPACE values). In this paper we use in formula (6) exactly the same values for A ms and B ms as in M1992. As opposed to P2007 who derived new A ms and B ms values for M1992 type formula, we changed N ms value. Moreover, N ms that is constant in this study can be a function of altitude and geographical location (among others) if used in global models. In this case it would be better to treat A ms and B ms as constants and make changes to N ms in a manner that it is done to account for maritime and continental CCN concentration differences in some GCMs.

Comparison with observations
To facilitate a comparison between observations shown in Table 3 and simulations, we show the same averaged characteristics for experiments with ice microphysics -I1, I2, I3, and I4 -in Table 7. Comparison of these tables indicate that observed and simulated microphysical characteristics (concentration of liquid and solid particles, LWC, and IWC) are quite similar. The R ew calculated using observed and simulated data are also comparable.
At the same time ice crystal effective radii R ei calculated from observations and simulations differ significantly. The R ei calculated from observations are about 25 microns for October 9-10 flights, whereas the R ei calculated from simulations are systematically greater. For example, values of R ei for I2 and I3 are 8 and 5 times greater than that from observations. Possible reasons for these differences are from numerical diffusion and different techniques used for R ei calculations. Numerical diffusion is an unavoidable feature of any numerical scheme used to solve equations for distribution functions for condensation/evaporation, deposition/sublimation, and BFP. Because favorable conditions for the BFP exist in modeled clouds during glaciation, depositional growth of ice crystals at the expense of evaporated cloud droplets is a reason that might determine the artificial spectra broadening in numerical simulations (see the Appendix A for details). A second reason for possible differences between observed and simulated R ei is different techniques used to calculate its values.
To calculate R ei from the observations  the following definition based on the ice water content (IWC) and cross-sectional area of the particle distributions (A c ) is used (Fu, 1996): where ρ i = 0.9 g cm −3 is the bulk density (mass divided by volume) of the ice crystals. The R ei calculated from the observations are highly dependent on the mass-diameter (m-D) relation that is assumed to characterize the observed size distributions (for details see McFarquhar and Heymsfield (1998)). The R ei calculated from the simulations correspond to a "composite" crystal distribution because more than one type of crystals with different shapes and densities are used. These R ei are provided in Table 3 (R ei for individual ice crystals are listed in Table 6). The "composite" ice phase effective radius shown in Table 7 is calculated as where r k (m) are bulk radius for columns (k = 2), plates (k = 3) and dendrites (k = 4), respectively. Definition (11) is useful for analysis of radar data providing information about ice particles sizes. As Table 7 and Table 3 show, R ei calculated using (11) reflect the contribution of large crystals to size distribution and are significantly greater than those calculated using (10) chosen in such a way that if a lot of large ice crystals exist the R ei are actually small. To compare R ei calculated from the observations and simulations using definition (10) a "composite" m-D relation is needed. It is not easy to determine what m-D relation might apply to the "composite" crystal distribution from the simulations. It becomes evident that techniques used to calculate different microphysical characteristics from observations and essential BRM scheme characteristics (mass grids, m-D relations, hydrometeor densities, capacitances, and terminal velocities among others) should be interrelated. Otherwise, direct comparison of data derived from observations and simulations is not logically based. 23 To determine if the differences between observed and simulated R ei arise due to different definitions, we use the formula that mimics (10) for individual ice crystals: where R ek , ρ ik , and A ck are ice crystal effective radius, density, and projected area, respectively, for columns (k = 2), plates (k = 3), and dendrites (k = 4). Corresponding "composite" ice phase effective radius R ei then calculated as Another definition of the "composite" ice phase effective radius shown in Table 7 and Table 8 reads as where r m (m) is melted radius (radius of sphere that has the same mass as ice particle and whose density is equal to water density). Both (11) and (14) definitions are useful for analysis of radar data because they provide information about crystal sizes. Table 8 shows the composite ice effective radii calculated from the simulations using (11), (13), (14), respectively, and effective radius for individual ice crystals (plates R ep and dendrites R ed ) calculated according to (12). As can be seen from Table 8, R ei calculated using (13) from simulations are within observed ranges (Table 3) indicating comparability of observed and simulated ice crystal distributions. But these radii show relatively small variability from experiment to experiment as compared to R ei calculated using definition based on melted radius. Definition based on melted radius requires only distributions of ice particles on mass grids, and additional knowledge of ice crystals m-D relations, projected areas, bulk radii and bulk densities is not necessary. Thus, ice crystal effective radius definition base on melted radius should be recommended for evaluation of relative importance of different microphysical processes such as different ice initiation mechanisms in intercomparison studies.
The differences between effective radius calculated using different definitions for individual ice crystal (Table 6 and Table 8) as well as "composite" ice phase effective radius R ei (Table 7 and Table 8) highlight the necessity to standardize calculation of ice effective radii since these are ultimately provided as input for radiation calculations.

Discussion
To improve the representation of mixed-phase cloud processes in the GISS GCM and facilitate the improvement of bulk microphysics parameterizations that do not use known a priory shape of hydrometeors' distribution functions, we couple a mixed-phase BRM scheme to the GISS SCM. We perform sensitivity simulations with and without ice microphysics to evaluate the impact of the CCN spectrum shape, process of warm rain formation, different ice initiation mechanisms, and the Bergeron-Fendeisen process on glaciation time and longevity of mixed-phase clouds observed during the ARM MPACE IOP.
Based on differences between our sensitivity simulations that do not include ice microphysics, we find that the process of water-water interaction may be relatively minor compared to that of the CCN spectrum shape for droplet activation for the MPACE single-layer mixed-phase clouds.
For the ice phase initiation we consider two main mechanisms. The first mechanism is active in cold supersaturated (w.r.t. ice) environments and determines the number of small ice crystals originating from water vapor, whose shapes depend on temperature. The second mechanism of ice initiation is active at negative temperatures in both saturated and under-saturated (w.r.t. water) environments due to the transformation of super-cooled droplets, whose spectrum and masses as well as degree of supercooling determine the rate of origination of bigger (up to 100 µm) plate-like crystals. Because the freezing rate depends on the droplet mass, the bigger droplets are likely to freeze faster. These two ice initiation mechanisms act quite differently. The first IIP is responsible for the supply of small ice crystals with different shapes. These crystals grow fast at different rates in a highly supersaturated (w.r.t. ice) environment at the expense of evaporated cloud droplets. The second IIP is responsible for the supply of bigger (assuming the droplet spectrum is broad enough) ice crystals that continue to grow mainly due to riming, reducing droplet concentration and water vapor supply for the ice phase due to droplet evaporation. The second mechanism indicates the importance of the AP spectrum for the ice initiation process. It crucially depends on the shape of the AP distribution and not only on the concentration of cloud droplets but also on the broadness of the spectrum of cloud droplets just activated. We speculate that in maritime stratiform clouds with broader droplet spectra the second IIP might be of greater importance. In our simulations with prescribed large-scale forcing that assures the existence of high supersaturation (w.r.t. ice) (up to 20 %) and coalescence processes switched off, the net supply of new ice particles due to the two ice initiation mechanisms has the same order of magnitude.
The differences between ice effective radii calculated using ice crystal cumulative cross-sectional area and melted radius definitions indicate importance of the first definition for radiation calculations and the second definition for analysis of precipitation formation process in mixed-phase clouds. Because of the relatively small variability of ice effective radius calculated using cross-sectional area definition, ice effective radius definition based on melted radius should be used as additional microphysical characteristic for evaluation of relative importance of different microphysical processes such as different ice initiation modes in intercomparison studies.
Recently, a 2-D CRM was used to obtain differences in cloud properties in simulations with one and two-moment bulk microphysics (BLK) for MPACE conditions (Luo et al., 2007). MPACE mixed-phase clouds were also simulated with a 3-D Arctic ver-sion of MM5 with a two-moment bulk microphysics scheme to evaluate sensitivity of clouds properties to cloud condensation and ice nuclei concentration . Although BLK schemes are usually able to represent adequately the variations of droplet concentration for maritime and continental clouds, their ability to represent the process of droplet activation for maritime and continental clouds with respect to broadness of spectrum of cloud droplets just activated is limited. Accounting only for the variations of the droplet concentration under different aerosol conditions is necessary, but not sufficient, for the appropriate representation of ice initiation processes in mixed-phase clouds. This fact has to be taken into account if bulk microphysics schemes are used to investigate relative importance of different ice initiation modes in mixed-phase clouds.
Analytical considerations highlight the effect of the BFP on the longevity of mixedphase clouds (Korolev, 2007;Korolev and Field, 2008). In our sensitivity runs, originated ice crystals continue to grow in simulated clouds mainly due to the BFP that is identified as a process responsible for the rate of glaciation of single layer mixed-phase MPACE clouds. An adequate treatment of Bergeron-Findeisen process is important for models that use BRM or BLK schemes to investigate these types of Arctic clouds. Despite the high computational cost, our calculations of hydrometeors' growth rates due to the BFP are based on analytical solution to equations for supersaturation (w.r.t. water and ice), and the changes of supersaturation during the microphysical time step in liquid/solid particle growth equation are also taken into account. It is difficult to expect that the utilization of different modifications of "saturation adjustment" that is widely used in BLK schemes can represent the simultaneous growth rate of cloud particles due to the BFP. Since the droplet nucleation process (w.r.t. broadness of spectrum of cloud droplets just nucleated) and the BFP (w.r.t. calculation of simultaneous evaporation rates for droplets and deposition rates for ice particles) are difficult to be reliably represented in bulk schemes, the interpretation of the results with these schemes in the case of mixed-phase clouds as observed during MPACE has to be done very carefully.
One of the possible ways to improve the creditability of mixed-phase bulk micro-physics schemes is the creation of a unified modeling framework that includes a computationally expensive BRM-type scheme and a computationally efficient but less sophisticated microphysics scheme. Development of such a scheme should be based on observations and numerical simulations obtained using the BRM scheme that is considered as a benchmark. This work is underway. Our future study will focus on the investigation of the impact of different environmental conditions and processes of water-ice and ice-ice interaction on the longevity and glaciation time of mixed-phase MPACE clouds using the BRM scheme and a two-moment BLK scheme (Morrison et al., 2005) coupled to the GISS-LBL SCM. Numerical implementation of Bergeron-Findeisen process outlined in Appendix A is applicable in any cloud model that resolves supersaturation and utilizes bin-resolved or bulk microphysical schemes. In the latter case it can be significantly simplified because of prescribed a priori hydrometeors' size distributions. This approach can also be used for developing parameterization of processes of vapor/liquid/solid phase transformations for use in large-scale models.

Appendix A
The rate of changes of distribution function f 1 for liquid phase due to condensation (dm 1 /dt > 0) or evaporation (dm 1 /dt < 0) is written as Equation (A1) provides two useful computational constrains for condensation or evaporation processes 1) Integrating equation (A1) with respect to mass m 1 from 0 to ∞ where n we is total number of evaporated droplets.
The first equation (A2) [ has the simple physical meaning that in the condensation process concentration of droplets is constant. The second one expresses the fact that in the evaporation process the total number of existing and evaporated cloud droplets remains unchanged 2) Multiplying equation (A1) by mass m 1 and integrating resulting equation with respect to m 1 and using definition (3) for k = 1, we get The last equation determines the increase in liquid water content (LWC) q w due to condensed water vapor supply or decrease in LWC due to evaporation.
The rate of change of the water vapor mixing ratio q due to condensation/evaporation in ice free environment can be written as Adding (A5) and (A6), it follows that Since the mass conservation law has to be satisfied, we obtain The rate of change of the temperature T can be written as where L w is the specific latent heat of evaporation and c p is specific heat of air at constant pressure. Combining (A6) and (A10), we get energy conservation law The rates of changes of distribution functions f k for solid hydrometeors (k=2. . . 7) due to deposition (dm k /dt > 0) or sublimation (dm k /dt < 0) are given as where k is the type of hydrometeor (k = 2 . . . 4, ice crystals; 5, aggregates; 6, graupel; and 7, frozen drops/hail).
Equations (A3), (A5), (A6), (A8)-(A11) for ice phase can be written as where L i is the specific latent heat of sublimation, n i = 7 k=2 ∞ 0 f k (m k )dm k and q i = 7 k=2 ∞ 0 m k f k (m k )dm k are ice concentration and ice water content (IWC), respectively.
In mixed-phase cloud the rates of changes of water vapor mixing ratio and temperature due to diffusional processes are governed by where ε w and ε i are rates of changes of LWC and IWC, which are defined by (A9) and (A17), respectively. Both ε w and ε i depend among other characteristics on supersat-31 uration w.r.t. water S w and ice S i that change during one microphysical time step. To account in this fact and calculate ε w and ε i , we define size distribution function for each type of hydrometeors on the mass grids. The mass grid for each type of hydrometeor is represented by different numbers of mass bins N k : where j is the mass bin number, m k0 is the minimal mass for hydrometeor of type k, J k0 and a k > 1 are parameters that characterize the mass grid. For example, N k = 33, J k0 = 1 and a k = 2 were used in Khain and Sednev (1996) (KS96). Diffusional growth (evaporation) of water droplets of mass m 1j in (A9) is expressed (Pruppacher and Klett, 1978) The changes of ice particles mass m kj (k > 1) due to deposition (sublimation) in (A17) is written as (PK78): In the above, D v , k a , and R v are the water and air diffusivity coefficients and the moist air gas constant, respectively; expressions for the "electrostatic capacitance" of particles of different shape C kj are taken from PK78 (see also KS96). The method used for the calculation of supersaturation (SS) is similar to that used by Tzvion et al. (1989) and KS96 with some additional modifications. The calculation of SS w.r.t. water S w = (e/e sw − 1) and ice S i = (e/e sw − 1) (where e, e sw , and e si are water vapor pressure and its saturated values w.r.t. water and ice, respectively), are performed in two steps. First, the equations for the advection of potential temperature Θ and water vapor mixing ratio q are integrated during a dynamical time step 32 ∆t dyn without microphysical terms. As a result, the values of supersaturations S * w and S * i , as well as the non-microphysical tendencies of (δS w /δt) dyn = (S * w − S 0 w )/∆t dyn and (δS i /δt) dyn = (S * i − S 0 i )/∆t dyn are calculated at each grid point. The dynamical time step is divided into several microphysical time steps, ∆t dif . The change of supersaturation at each microphysical time step is calculated as the sum of the nonmicrophysical tendency [e.g., (δS w,i /δt) dyn ∆t dif ] and changes caused by diffusional growth/evaporation of liquid phase or deposition/sublimation of ice phase.
Using Eqs. (A20) -(A21), (A23) -(A24), definitions (A9) -(A17), expression for the water vapor mixing ratio q = 0.622(e/p), and dependence of the saturation vapor pressure over water e sw and ice e si on temperature, one can derive the following equations for S w and S i (KS96): Coefficients P w , P i , R w , and R w in (A25) − (A26) are given by where de sw /dT and de si /dT are any analytical formulai that express dependence of saturation pressure with respect to water e sw and ice e si on temperature.
If the microphysical time step ∆t dif is small enough, the coefficients (A27) − (A30) can be considered as constants, and the analytical solution of (A25) − (A26) during the time τ ∆t dif can be written as KS96: and where To account the fact that S w and S i are non-constant during time step the following iteration procedure is used. Expressions (A23)−(A24) and solution (A31)−(A32) permit us to calculate new water (k = 1) and ice (k > 1) particle masses m kj in jth bin: where 0 < τ w/i 1 are parameters, s is iteration number, and m (t 0 ) kj are given by (A22). It was found that effective stopping criterion of the iteration process (A36)-(A37) is where δ is a minimum mass increment permitted. If creterion (A38) is satisfied, we use m (s+1) kj and the method by Kovetz and Olund (1969), which conserves both concentration (A3), (A13) and mass (A5), (A14), to calculate new values of distribution functions f kj (t 0 + τ ) on regular mass grids. To derive expressions for the changes of LWC ∆q w and IWC ∆q i during timestep τ , we use (A22) and definition of hydrometeor content (3), which can be rewritten as Because of the fact that dJ = 1, we replace integral by summation and obtain Since ∆q w and ∆q i are known, mass and energy conservation laws are used to calculated new temperature T(t 0 + τ ) and water vapor mixing ration q(t 0 + τ ) at the end of microphysical time step.  (6) and (7), respectively, that influence ice crystal concentration increase due to the different ice initiation mechanisms.     Table 6. Average values of ice water content (IWC p ), effective radius (R ep ), concentration (N p ), and precipitation flux (P p ) for plates and ice water content (IWC d ), effective radius (R ed ), concentration (N d ), and precipitation flux (P d ) for dendrites in experiments with ice microphysics.  Table 7. Average values of liquid water content (LWC), effective radius (R ew ), and concentration (N w ) for liquid phase and ice water content (IWC), effective radius (R ei ), and concentration (N i ) for ice phase in experiments with ice microphysics.