Parameterizing the competition between homogeneous and heterogeneous freezing in cirrus cloud formation – monodisperse ice nuclei

We present a parameterization of cirrus cloud formation that computes the ice crystal number and size distribution under the presence of homogeneous and heterogeneous freezing. The parameterization is very simple to apply and is derived from the analytical solution of the cloud parcel equations, assuming that the ice nuclei population 5 is monodisperse and chemically homogeneous. In addition to the ice distribution, an analytical expression is provided for the limiting ice nuclei number concentration that suppresses ice formation from homogeneous freezing. The parameterization is evaluated against a detailed numerical parcel model, and reproduces numerical simulations over a wide range of conditions with an average error of 6 ± 33%. This study presents an analytical parameterization of cirrus formation that addresses the issues identiﬁed and explicitly considers the competition between homogeneous and heterogeneous nucleation. We develop a novel method to account for the growth of the heterogeneously frozen ice crystals and incorporate their e ﬀ ect within the physically-based homogeneous nucleation framework of Barahona and Nenes (2008, hereinafter BN08). The new parameterization explicitly resolves the dependency of N c on dynamical conditions of cloud formation (i.e. temperature, pressure, and updraft ve-locity), aerosol number, size, the characteristics of IN and their e ﬀ ect on ice saturation ratio. To clearly present


Introduction
Cirrus clouds are key components of climate and can have a major impact on its radiative balance (Liou, 1986;Hartmann et al., 2001;Lohmann and Diehl, 2006). They can be affected by anthropogenic activities such as aircraft emissions in the upper troposphere (Seinfeld, 1998;Penner et al., 1999;Minnis, 2004;Stuber et al., 2006;Kärcher 15 et al., 2007), and from long range transport of pollution and dust (Sassen et al., 2003;Fridlind et al., 2004). When studying ice-cloud-climate interactions, the key cloud properties that need to be computed are ice crystal concentration distribution, and ice water content (e.g. Lin et al., 2002); in large-scale models such parameters are either prescribed or computed with a parameterization. Due to the limited understanding of the Heterogeneous freezing encompasses a group of nucleation "modes" (deposition, contact, immersion, and condensation) in which ice formation is mediated by insoluble solids (Pruppacher and Klett, 1997) or surfactant layers  that lower the energy barrier for the formation of an ice germ. Macroscopically, this can be expressed as a decrease in the freezing saturation threshold (or an increase in the 15 freezing temperature) with respect the value for homogeneous freezing. Aerosol particles that contribute such surfaces and undergo heterogeneous freezing are termed ice nuclei (IN) (Pruppacher and Klett, 1997); the nature of IN is very diverse and not well understood. They can be composed of mineral dust (DeMott et al., 2003b;Sassen et al., 2003), sulfates (Abbatt et al., 2006), crustal material (Zuberi et al., 2002;Hung 20 et al., 2003), carbonaceous material (Beaver et al., 2006;Cozic et al., 2006;Zobrist et al., 2006), metallic particles (Al, Fe, Ti, Cr, Zn, and Ca) (Chen et al., 1998;DeMott et al., 2003a), and biological particles (e.g. Möhler et al., 2007). Sources of IN have been associated with long range transport of dust (DeMott et al., 2003b;Sassen et al., 2003), pollution (Fridlind et al., 2004), and direct aircraft emissions (Penner et al., 1999 Haag et al., 2003;Archuleta et al., 2005;Field et al., 2006;Prenni et al., 2007). Field campaign data generally show that ice crystals collected at conditions favorable for homogeneous freezing contain a moderate fraction of insoluble material (e.g. DeMott et al., 2003a;Prenni et al., 2007), which suggests the potential for interaction between homogeneous and heterogeneous freezing during cirrus formation. There- 5 fore, a comprehensive understanding of ice formation requires a combined theory of homogeneous and heterogeneous freezing (e.g. Vali, 1994;Kärcher and Lohmann, 2003;Khvorostyanov and Curry, 2004;Zobrist et al., 2008). Modeling studies however suggest (Kärcher and Lohmann, 2003;Khvorostyanov and Curry, 2005;Liu and Penner, 2005;Monier et al., 2006) that at atmospherically relevant conditions, a sin- 10 gle mechanism (i.e. homogeneous or heterogeneous) dominates the freezing process. This is largely because IN tend to freeze first, depleting water vapor, and inhibiting homogenous nucleation before it occurs (DeMott et al., 1997); only if IN concentration is low enough can homogeneous freezing occur. Thus, IN and supercooled droplets may be considered separate populations that interact through the gas phase, but undergo 15 heterogeneous or homogeneous freezing, respectively, at different stages during cloud formation (DeMott et al., 1997). This distinction may become less clear for heterogeneous freezing in the immersion mode for low concentration of immersed solid, or, if the solid is not an efficient ice nuclei (i.e. S h is high, close to the value for homogeneous freezing) (Khvorostyanov and Curry, 2004;Marcolli et al., 2007); we will assume for the 20 rest of this study that we are far from such conditions, given that their importance for the atmosphere is not clear.
The presence of IN may drastically reduce the crystal concentration compared to levels expected for homogeneous freezing (DeMott et al., 1994(DeMott et al., , 1998, creating a potential anthropogenic indirect effect of IN on climate (DeMott et al., 1997;Kärcher and 25 Lohmann, 2003). Therefore, parameterizations of cirrus formation for use in climate models should explicitly resolve the interaction between homogeneous and heterogeneous freezing. Using a simplified approach to describe the growth of ice particles and box model simulations, Gierens (2003)  This study presents an analytical parameterization of cirrus formation that addresses the issues identified above and explicitly considers the competition between homogeneous and heterogeneous nucleation. We develop a novel method to account for the growth of the heterogeneously frozen ice crystals and incorporate their effect within the physically-based homogeneous nucleation framework of Barahona and Nenes (2008, 20 hereinafter BN08). The new parameterization explicitly resolves the dependency of N c on dynamical conditions of cloud formation (i.e. temperature, pressure, and updraft velocity), aerosol number, size, the characteristics of IN and their effect on ice saturation ratio. To clearly present the approach used to unravel the interactions between IN and supercooled droplets, and develop the details of the parameterization, we assume IN 25 are monodisperse and chemically uniform. The extension of the parameterization to polydisperse IN (in size and chemical composition) will be presented in a companion study. Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version Interactive Discussion

Parameterization development
The parameterization is based on the ascending parcel conceptual framework. Contained within the cooling parcel is a population of supercooled droplets externally mixed with IN that compete for water vapor when their characteristic freezing threshold is met. Heterogeneous and homogeneous freezing occur sequentially as the parcel water va-5 por saturation ratio over ice, S i , increases during the cloud formation process. Homogeneous freezing is restricted to liquid droplets and described using the formulation of BN08. The IN population can be characterized by its number concentration, N IN , and freezing threshold, S h . Upon freezing, ice grows on the IN, depleting and redistributing the available water vapor between the homogeneously and heterogeneously 10 frozen crystals (i.e. IN-frozen crystals) (e.g. Lin et al., 2002;Ren and Mackenzie, 2005;Kärcher et al., 2006). Water vapor competition also affects the size distribution of homogeneously frozen crystals. To account for this competition effect, the change in homogeneous freezing probability due to the presence of IN-frozen crystals (i.e. the amount of water vapor condensed upon IN), and the size of the IN-frozen crystals at 15 the homogenous freezing threshold, should be calculated. Both effects can be directly related to the impact of IN-frozen crystals on the rate of change of S i (Barahona and Nenes, 2008). A method to represent such interactions is presented below. According to the homogeneous nucleation framework of BN08, the fraction of frozen droplets for a pure homogeneous pulse is given by where symbols are defined in Appendix A and B. Equation (1) was obtained under the assumption that homogeneously nucleated ice crystals do not substantially impact S i , up to the point of maximum ice saturation ratio (in reality, S i is affected, and is accounted for by adjusting the effective growth parameterΓ and setting S max equal to Introduction as all other terms in Eqs. (1) and (2) where N o is the liquid droplet number concentration in the parcel. The first term in Eq. (4) is the contribution from homogeneously frozen droplets (Barahona and Nenes, 2008), while the second term is from heterogeneous freezing.
where symbols are defined in Appendix B. Equation (5) expresses the rate of change of parcel ice saturation ratio, which increases from cooling (first term of right hand side), and decreases from deposition on the growing IN (second term of right hand side).

10
Equation (6) represents the rate of water vapor deposition on IN-frozen crystals. Equation (7) expresses the rate of change of crystal size due to water vapor condensation. Equation (5), evaluated at S max providesṠ max,IN . If the initial conditions of the parcel are known, Eqs. (5) to (7) can be numerically solved, as closed solution of Eqs. (5)  very stiff as S i →S max , so even numerical integration of the system becomes difficult.
The behavior outlined above can be described in terms of the following dimensionless variables, where,Ṡ h =αV S h , is the rate of change of S i at the instant when IN freeze. x represents a normalization of S i by the difference between homogeneous and heterogeneous freezing thresholds, and y represents a normalized d S i d t with respect to its value at S h . The relation between x, y, S i , and 15 x vs. y is linear and the value of y at x=1 (termed "y max ") is equal to − Increasing N IN is reflected by increase in y max (y max 1 and y max 2 , for green and purple curves, respectively), in these cases the change in the sign of y (Fig. 1, bottom panel) is related to the onset of significant depletion of water vapor, which impacts is represented as y max =1 (red curves); at this point S i passes by a maximum at S max 20 so that there are at least two values of y for each value of x, one in the increasing branch of S i and the second as S i decreases (dashed red curves). Thus, at x=y=1 the functional dependency between x and y ceases to be univocal and lim Based on the above analysis an approximate mathematical relation can be obtained between x and y. Since x is less than 1, together with dx d y y=1 =0, a Taylor series expansion around y=1 gives, where K is a positive constant, and O(y, y 2 , ...) are higher order terms that are neglected (the consequences of this assumption are discussed in Sect. 3). Since to zeroth order, the curvature around the point y=1 is constant (equal to −K ), x is related to y as, x=−K y 2 +ay+b (being a and b integration constants). Since x=0 when y=0, and, dx d y =0, for y=1, then b=0, K = a 2 , and Ax = y(y − 2) where A=− 1 a is an arbitrary constant. As we desire x and y to be normalized, we select A such that x=1 when y=y max , (Fig. 1), or A=y max (y max −2). Substitution into Eq. (11) and rearranging gives 15 Equation (12)  , must be calculated (Eq. 5); however, due to the coupling between Eqs. (6) and (7) (7) can be used for the latter, as follows. First, rearranging Eq. (7) and integrating from S h to S max we obtain  8) and (9), and using Eq. (12) substituting Eqs. (14a) and (14b) into Eq. (13), gives integrating Eq. (15a), and after substituting A=y max (y max −2) and rearranging, gives substituting Eq. (15b) into Eq. (13), and solving for D max,IN gives,  (7), at the point of maximum ice saturation ratio in the parcel; y max can be computed after combining Eqs. (6) and (7) and using Eq. (14b),  (Gierens, 2003). Figure 1 shows that this condition is uniquely determined by y max =1; the limiting size of the ice crystals, D lim , at the maximum water vapor deposition rate can then be computed from Eq. (16) for ∆S * h calculated for y max =1,  (7) into Eq. (19), N lim can be obtained as

Applying the parameterization
The application of the parameterization is outlined in Fig. 2 (5) is used then to calculate N c from combined 5 heterogeneous and homogeneous freezing.

Results and discussion
The parameterization was evaluated against the numerical parcel model developed in BN08 modified to include the effect of IN (by allowing the deposition of water vapor onto the IN when S i exceeds S h ). The performance of the parameterization for pure homogeneous freezing conditions was discussed in BN08, so the evaluation here is focused on the competition between homogeneous and heterogeneous freezing. The conditions used for evaluation covered a wide range of T , V , α d , S h and N IN (Table 1) for a total of 800 simulations; initial p was estimated from hydrostatic equilibrium at T o using a dry adiabatic lapse rate. For simplicity, α d was assumed to be the same for the 15 homogeneously frozen and IN-frozen ice crystal populations (although this assumption can easily be relaxed).   Figure 3 shows that at high V and low T , large concentrations of IN (N lim ∼1 cm −3 ) are required to impact homogeneous freezing, which is consistent with existing studies (e.g. DeMott et al., 1997); given that high N IN are rare in the upper troposphere (e.g. Gayet et al., 2004;Prenni et al., 2007), our results are consistent with current understanding that homogeneous nucleation is the 5 primary freezing mechanism at these conditions. Overall, the parameterization reproduces very well the parcel model results (Fig. 3), and clearly captures the effect of N IN on N c . When calculating the error between parcel model and parameterization, care should be taken to exclude points where homogeneous nucleation is prevented, as both the parameterization and the parcel model  N IN →0, Fig. 3); which is about 1±10% for all simulations of Table 1. This suggests that most of the discrepancy between parameterization and parcel model results occurs at low N c (N IN close to N lim ), therefore, most of the error in the parameterization should result from the calculation of N lim . This is 5 supported by the simulations in Fig. 3, which show that the shape of the parameterized N c vs. N IN follows well the parcel model results, and discrepancy between the parameterization and parcel model mainly originates from slight "shifting" in N lim . Potential sources of error in the calculation of N lim can be identified by analyzing Eq. (20). First, S max has been assumed to remain constant during the homogeneous pulse. Although this is a good approximation for computing N c from homogeneous nucleation, ∆S h , which is dependent on S max (Eq. 18), may be in error and bias the calculation of D lim (hence N lim ). This is shown in Fig. 4, which presents D lim (dotted lines) and N lim (solid lines) as a function of S h , for T =215 K (S max =1.51, blue lines) and T =230 K (S max =1.46, black lines). At low S h (∆S h >0.1), D lim decreases almost 15 linearly with S h . For ∆S h <0.1, the dependency of D lim on S h changes so that D lim →0 when S h →S max . As a result, for ∆S h >0.1, N lim is almost linear in S h (Eq. 22) and insensitive to small errors in S max . As S h →S max , N lim increases steeply and becomes very sensitive to errors in S max . This is strongly evident for T =230 K when S h >1.4 (Fig. 4). At these conditions, ∆S h <0.06, the parameterization becomes very sensitive 20 to the value of S max (Fig. 3 for S h =1.4 and T =230 K, i.e. black lines and circles), and tends to overestimate N lim because the parameterized S max is below S max reached in the parcel model simulations. This effect is more pronounced at high V (Fig. 3) as S max tend to reach high values at such conditions. The overestimation in N lim at S h =1.4, however, is less pronounced at T =215 K (Fig. 3 (10), O(y, y 2 , ...). Because of this, Eq. (12), will not predict a change in the sign of y (e.g. green line, Fig. 1 Some bias may also be introduced when S h is low. At these conditions, recently frozen crystals are exposed to low S h therefore grow slowly, so even at N IN =N lim there is a substantial time lag before they can impact ice saturation ratio. Since Eq. (12) does not account for this, D lim is slightly overpredicted and N lim underpredicted (e.g. Fig. 4). This "growth lag" effect is even stronger at low T and low α d since condensation (growth) rates are further decreased. This is seen in the slight underprediction of the parameterization shown in Fig. 3 at S h =1.1 and T =215 K (upper panels, blue lines and 15 circles).

Evaluation against other parameterizations
Finally we compared Eq. (20) to the expressions provided for N lim by Gierens (2003, G03) and Liu and Penner (2005, LP05). LP05 (Fig. 5, green circles) is based on detailed parcel model simulations and resolves the dependency of N lim on T and V . G03 20 (Fig. 5, stars) uses a semi-analytical approach, where box model simulations are employed to fit the time scale of growth of the ice crystals; it resolves explicitly the dependency of N lim on T , V , S h and p. Figure 5 presents the results of the three expressions as functions of T (left panel, V =0.1 m s −1 ) and V (right panel, T =215 K), and S h =1.3 (blue) and 1.4 (black). Equation (20) is presented in Fig. 5 for α d =0.1 (dashed lines) 25 and 1 (solid lines). At a given V (Fig. 5, left panel), N lim decreases almost exponentially with increasing T . The difference between N lim predicted by G03 and Eq. (20)  Interactive Discussion the largest at T =200 K, but is minimal at T =230 K. The convergence at high T is likely from the large crystal size attained, which implies that the mass transfer coefficient is consistent with the approximation taken in G03 of neglecting kinetic effects. LP05 predicts N lim about one order of magnitude lower than Eq. (20), potentially because they included both deposition and immersion freezing in their simulations; deposition 5 IN would act as a second population (with different S h and N IN ) that would compete for water vapor with the existing crystals lowering N lim . Equation (20) cannot currently represent such a case, but can be expanded to include it and will be the subject of a future study. At fixed T , G03 and Eq. (20) predict a N lim ∼V 3/2 dependency (Fig. 5, right panel). LP05 shows a higher order dependency, likely due to the effect of deposition 10 IN depleting water vapor.

Summary and conclusions
We developed an analytical parameterization for cirrus formation that explicitly considers the competition between homogeneous and heterogeneous nucleation. Using functional analysis of the crystal growth and ice saturation ratio evolution equations for 15 a cloud parcel, an expression relating d S i d t to S i was found, and used to calculate the limiting size of the ice crystals, D lim , and limiting concentration of IN that would prevent homogeneous nucleation, N lim . These two parameters were employed to obtain a simple, yet accurate, expression to account for the reduction in d S i d t due to the presence of previously frozen crystals which, when incorporated within the physically-based homo-20 geneous nucleation framework of Barahona and Nenes (2008), allows the calculation of the number of crystals produced when heterogeneous and homogeneous freezing mechanisms compete for water vapor during cirrus cloud formation.
When evaluated over a wide range of conditions, the new parameterization reproduced the results of the detailed numerical solution of the parcel model equations,  (7). The approach presented here overcomes common difficulties in including the competition between homogeneous and heterogeneous modes in large scale models. In this study we assumed a single type, monodisperse, chemically uniform, IN population. Multiple IN populations and freezing modes that reflect cirrus formation for atmospheric relevant aerosol is the focus of a future study. Even though much work still remains to be done on the prediction of IN concentrations and freezing thresholds, this work offers a computationally efficient and rigorous approach to include such knowledge in atmospheric and cloud models, once available.

Appendix A
A summary of parameters developed for the BN08 parameterization is provided here.

10
The effective growth parameter in Eq.
(1) is defined as where γ=1.33, and D c,s max is the maximum size reached by the homogeneously nucleated ice crystals, given by Interactive Discussion with T in K. S max is obtained by solving J(S max )=10 16 s −1 m −3 (Koop et al., 2000), using the analytical fit of Ren and Mackenzie (2005).