A bulk parametrization of melting snowflakes with explicit liquid water fraction for the COSMO model

A new snow melting parametrization is presented for the non-hydrostatic limited-area COSMO (“consortium for small-scale modelling”) model. In contrast to the standard cloud microphysics of the COSMO model, which instantaneously transfers the meltwater from the snow to the rain category, the new scheme explicitly considers the liquid water fraction of the melting snowflakes. These semi-melted hydrometeors have characteristics (e.g., shape and fall speed) that differ from those of dry snow and rain droplets. Where possible, theoretical considerations and results from vertical wind tunnel laboratory experiments of melting snowflakes are used as the basis for characterising the melting snow as a function of its liquid water fraction. These characteristics include the capacitance, the ventilation coefficient, and the terminal fall speed. For the bulk parametrization, a critical diameter is introduced. It is assumed that particles smaller than this diameter, which increases during the melting process, have completely melted, i.e., they are converted to the rain category. The values of the bulk integrals are calculated with a finite difference method and approximately represented by polynomial functions, which allows an efficient implementation of the parametrization. Two case studies of (wet) snowfall in Germany are presented to illustrate the potential of the new snow melting parametrization. It is shown that the new scheme (i) produces wet snow instead of rain in some regions with surface temperatures slightly above the freezing point, (ii) simulates realistic atmospheric melting layers with a gradual transition from dry snow to melting snow to rain, and (iii) leads to a slower snow melting process. In the future, it will be important to thoroughly validate the scheme, also with radar data and to further explore its potential for improved surface precipitation forecasts for various meteorological conditions.


Introduction
Accurate prediction of surface snowfall is a particularly challenging aspect of numerical weather prediction (NWP).It requires capturing the general dynamics leading to cloud and precipitation formation, and a meaningful representation of microphysical processes by cloud parametrizations.Particularly difficult are surface precipitation predictions in situations with near-surface temperatures slightly above 0 • C, when snowflakes start to melt before reaching the ground (e.g., Frick and Wernli, 2012).In such situations, slight errors in the vertical temperature and humidity structure and/or in the parametrization of the melting process can lead to a misforecast of the surface precipitation type (rain instead of snow or vice versa).Thériault et al. (2006) highlighted the great importance of the vertical temperature structure and Matsuo and Sasyo (1981) the impact of humidity for the resulting precipitation type.High values of relative humidity lead to an accelerated melting process while for low values the cooling due to sublimation of snow delays the onset of melting.The melting itself also influences the thermodynamics of the atmosphere due to the associated latent cooling (e.g., Raga et al., 1991).Lin and Stewart (1986) described perturbations due to cooling by melting snow that might lead to mesoscale thermal circulations.Additionally, Szyrmer and Zawadzki (1999) identified horizontal variabilities of atmospheric properties in the melting layer that are able to initiate convective cells.
These aspects illustrate the need for developing cloud microphysical parametrizations that address the issue of realistically representing the formation, sedimentation, and melting of snow (e.g., Seifert and Beheng, 2006;Woods et al., 2007;Thompson et al., 2008;Szyrmer and Zawadzki, 2010).
Published by Copernicus Publications on behalf of the European Geosciences Union.

C. Frick et al.: A bulk parametrization of melting snowflakes
The present study focuses on the numerical implementation of the melting process of snowflakes.In general, frozen hydrometeors start melting when their surface temperature exceeds 0 • C (Mitra et al., 1990).The melting process is not instantaneous and extends over a layer of several hundred metres for completely converting a snowflake into a raindrop (Mitra et al., 1990).The region in which this transition occurs is called the melting layer and has been intensively investigated (e.g., Bocchieri, 1980;Czys et al., 1996;Rauber et al., 2001).This layer consists of partially melted particles, which have specific characteristics that differ from the ones of dry snow and rain, e.g., in terms of shape and fall velocity.In radar observations the partially melted hydrometers of the melting layer lead to a band of enhanced radar reflectivity, the so-called "bright band".This phenomenon has been well described and investigated (e.g., Austin and Bemis, 1950;Yokoyama and Tanaka, 1984;Fabry and Zawadzki, 1995;Braun and Houze, 1995;Fabry and Szyrmer, 1999).
The melting of single snowflakes has been investigated in several laboratory studies (e.g., Matsuo and Sasyo, 1981;Mitra et al., 1990) and field experiments (e.g., Knight, 1979;Fujiyoshi, 1986;Barthazy et al., 1998).The findings of these studies provide the background for the numerical description of the melting process and its implementation into bulk microphysical parametrizations of NWP models.For the usage in operational NWP models the process is typically highly simplified, for instance by converting the meltwater generated within one time step instantaneously into rain.As a consequence, most schemes cannot represent partially melted snowflakes, which might affect the interaction between the melting process and the vertical profiles of temperature and humidity, and lead to erroneous fall speeds and an accelerated melting process.In a recent case study of a poorly predicted, high-impact wet snowfall event in northwestern Germany in November 2005, Frick and Wernli (2012) illustrated the difficulty in simulating such events with conventional snow melting schemes.Also when predicting the vertical temperature profile fairly realistically, the melting process in the considered NWP models was too fast leading to a prediction of surface rain instead of wet snow.Additionally, such a simplified treatment of the melting process in NWP models impedes a direct comparison between model data and radar observations in situations with a bright band.
The objective of this study is to develop a new melting parametrization that is able to represent wet snowflakes for the COSMO model (for more information about the COSMO model see, e.g., Doms and Schättler, 2002, or visit www.cosmo-model.org).If successfully implemented into the COSMO model, such a scheme might lead to improved snowfall predictions under melting conditions and the model data can be more easily compared and validated with radar observations of the "bright band".In Sect.2, the development of the new parametrization is described, starting with the melting behaviour of a single snowflake (Sect.2.2).The bulk approach for considering an ensemble of snowflakes is presented in Sect.2.3 and its implementation in the COSMO model version 4.14 in Sect.2.4.The source code of the presented parametrization is available from the authors.In Sect.3, some first results are presented to illustrate the differences in surface precipitation between simulations using the new and the standard melting scheme.Note that these results mainly provide an impression of the new scheme's potential; however a thorough validation of the new scheme is beyond the scope of this study and the subject of future investigation.

Modelling concept
The melting of snowflakes can be described based on the energy budget of an individual melting snowflake.It follows that the mass change of the corresponding meltwater above the temperature T 0 = 273.15K is given by (e.g., Mitra et al., 1990;Pruppacher and Klett, 1997) with T = T − T 0 where T is the temperature of the air, and q = e/T − e s,w (T 0 )/T 0 where e gives the vapour pressure and e s,w the saturation vapour pressure over water.L F is the latent heat of melting, L v the latent heat of evaporation, l h the diffusivity of heat, D v the diffusivity of water vapour, and R v the gas constant of water vapour.The left term in the brackets represents the contribution to melting due to diffusional heat transport and the right term accounts for evaporation.The ventilation coefficient F v and the capacitance C s depend on the snowflake characteristics while the remaining part of the equation is a function of the ambient air conditions.
Most bulk microphysical parametrizations make the simplyfing assumption that C s and F v only depend on the characteristic diameter (e.g., maximum dimension) D s of the snowflake.In this case, Eq. ( 1) can be written as where G(T , e) is a function that combines the thermodynamic, i.e., environmental parameters.For a bulk microphysical parametrization this rate equation has to be integrated for an ensemble of snowflakes.Usually it is assumed that the number density distribution of snowflakes can be described by an exponential or Gamma distribution, e.g., where λ is the slope and N 0 the so-called intercept parameter of the distribution.First, the diameter D has to be specified.Most measurements show an exponential or Gamma distribution w.r.t.equivalent diameter D eq (Gunn and Marshall, 1958;Brandes et al., 2007), but others (e.g., Field et al 3) with geometric diameter D s which simplifies the formulation of the collection rates.We follow the latter approach and use D s in Eq. (3).Second, the exponent µ has to be considered.For the calculation of mass changes, an exponential distribution, i.e., µ = 0, is a reasonable assumption for the number density distribution because the massweighted part of the spectrum is more dominant.Therefore, we apply an exponential distribution (µ = 0) in geometric diameter D s for the number density distribution of snowflakes.Large snowflakes are mostly aggregates of crystals and therefore it is well founded to expect a m s ∼ D 2 s relation (Westbrook et al., 2004a, b).More specifically we assume m s = αD 2 s with α = 0.069 kg m −2 (Wilson and Ballard, 1999;Field et al., 2005).The snow water content of snowflakes L s , and the mixing ratio q s are then given by with ρ for the density of air.The resulting mass change due to melting for an ensemble of snowflakes is calculated from The integral can be solved analytically and yields the mass change of the ensemble of snowflakes.In almost all bulk microphysical parametrizations the meltwater produced by Eq. ( 5) is instantaneously converted into rain (Lin et al., 1983;Rutledge and Hobbs, 1984;Cotton et al., 1986;Ferrier, 1994;Mölders et al., 1995;Reisner et al., 1998;Morrison et al., 2005;Seifert and Beheng, 2006;Lim and Hong, 2010).In these schemes, there is no internal mixing of liquid and solid water within the snow category, i.e., no separate prognostic variable for the meltwater on snowflakes.Meyers et al. (1997) discussed various attempts to fix the inherent problems of bulk melting schemes.This simplification leads to errors in the prediction of snowfall, the thermodynamics of the melting layer, and might lead to misforecasts of surface snowfall (Frick and Wernli, 2012), as already described in Sect. 1.The conceptional idea of the new melting scheme is to allow the representation of partially melted snowflakes, which leads to a splitting of the mass of an individual snowflake m s into a liquid part m w and a solid part m i , i.e., the melting snowflake has a liquid water fraction The melting snowflakes are represented by an internal mixture of liquid water and ice instead of having just externally mixed rain and snow categories.
Of course, the formulation of a model using internally mixed melting snowflakes and the incorporation in a bulk scheme is not straightforward.To our knowledge the only bulk parametrization including such an approach is the scheme of Szyrmer and Zawadzki (1999, SZ99 hereafter).For spectral bin models, where the implementation of a prognostic meltwater is considerably easier, Phillips et al. (2007) describe such a prognostic melting scheme.Our new scheme follows SZ99 with the main modification that we use the bulk meltwater mixing ratio instead of the critical diameter as the additional prognostic variable, which helps to ensure mass conservation and makes the implementation in a threedimensional model more consistent.
The following subsections describe the parametrizations of the properties of individual melting snowflakes, the formulation of the bulk scheme, and the implementation in the COSMO model.

Individual melting snowflakes
Following Mitra et al. (1990, M90 hereafter), who performed extensive laboratory studies on the melting of snowflakes in a vertical wind tunnel, the melting process of a single snowflake with included meltwater can be completely described.
The mass change of the meltwater of an individual snowflake above the critical surface temperature can be written as following Eq.(2) of M90.As for a snowflake without meltwater, cf.our Eq.( 2), G(T , e) only depends on the ambient air conditions, while capacitance C m and ventilation coefficient F v of the melting snowflake depend on the maximum dimension D s of the mass equivalent dry snowflake and the liquid water fraction .Based on their laboratory results M90 provide parametrizations of these parameters with the additional dependency on .
For the calculation of the capacitance M90 applied the approximation for an oblate spheroid.The axis ratio is assumed to be 0.3 for a dry dendritic crystal, and 1.0 for a raindrop.The axis ratio for melting snowflakes is approximated by a linear interpolation, i.e., a( ) = 0.3 + 0.7 (8) and the capacitance is then given by Pruppacher and Klett (1997, p. 547, Eq. 13-78),  and c) ventilation coefficient of a melting snowflake as a function of the equivalent diameter for various values of .Rain and snow correspond to = 1 and = 0, respectively.d) Empirical function Ψ( ) that describes the transition between the fall speed of a snow flake and a rain drop.Shown are the approximate relationship Eq. ( 15) and the original data (grey dots) of Mitra et al. (1990).
with C m (D s ,0) = C s and C m (D s ,1) = C r .D m is the maximum dimension of the melting snowflake, which can be calculated as follows assuming an oblate spheroid shape of the melting snowflake (see above) and in agreement with Eq. (8) of M90.Here ρ m is the density of the melting snowflake.As suggested by M90 we interpolate ρ m (D s , ) between the density of liquid wa-ter, ρ w = 1000kgm −3 , and the density of the dry snowflake ρ s (D s , ): For the density of a dry snowflake with the axis ratio of the melting snowflake it follows from the assumption of the oblate spheroid shape that  15) and the original data (grey dots) of Mitra et al. (1990).
with C m (D s , 0) = C s and C m (D s , 1) = C r .D m is the maximum dimension of the melting snowflake, which can be calculated as follows assuming an oblate spheroid shape of the melting snowflake (see above) and in agreement with Eq. (8) of M90.Here ρ m is the density of the melting snowflake.As suggested by M90, we interpolate ρ m (D s , ) between the density of liquid water, ρ w = 1000 kg m −3 , and the density of the dry snowflake ρ s (D s , ): For the density of a dry snowflake with the axis ratio of the melting snowflake it follows from the assumption of the oblate spheroid shape that but only until a maximum value of ρ s = 500 kg m −3 because higher densities are not reasonable for snowflakes.The empirical correction factor α cap ( ) in Eq. ( 9) is about 0.8 for dry snowflakes and for melting snowflakes M90 again suggest a linear interpolation, i.e., The resulting C m is shown in Fig. 1a for various values of .With increasing diameter, the capacitance of dry snow increases faster than the one of rain and the capacitance of partially melted particles follows the interpolation between both.M90 have shown that such a simple model for the capacitance is sufficient to achieve a good agreement with the experimental data.
Another important result of M90 is that the terminal fall velocity of a melting snowflake can be parameterized by where v s and v r are the terminal fall velocities of the mass equivalent dry snowflake and raindrop, which we calculate following Khvorostyanov and Curry (2005, KC05 hereafter).Both depend on the corresponding maximum dimensions which are in fact functions of the mass equivalent diameter D eq of a liquid sphere.For the calculation of the maximum dimension of the raindrop from D eq we follow Khvorostyanov and Curry (2002).For dry snowflakes D s is calculated by using m s = (π/6) ρ w D 3 eq and D s = (m s /α) 1/2 with α = 0.069 and, in addition, a cross sectional area of A = 0.45 π/4 D 2 s is assumed following Field et al. (2008).The empirical function ( ) is given by Fig. 2 of M90.An approximation for ( ) derived from the measurements is with α = 0.246, as shown in Fig. 1d.The resulting fall velocity of a melting snowflake is presented in Fig. 1b.As indicated by the empirical function ( ), the fall velocity of a melting snowflake shows a slow linear increase for small , but increases rapidly to the one of a raindrop for larger than 0.7.The resulting sedimentation velocity of a melting snowflake is necessary to describe the sedimentation process, but also affects the melting process due to the ventilation coefficient.
For the determination of the ventilation coefficient of a melting snowflake our definition of the length scale in the Reynolds number deviates from M90.Instead we are consistent with KC05 and others in using the maximum dimension of the particle.For large snowflakes or raindrops the drag coefficient becomes approximately constant.In this large particle limit the Reynolds number is only a function of mass, i.e., with ν a the kinematic viscosity of air, and therefore the Reynolds number of a melting snowflake N Re (D s , ) does not change strongly during the melting process.This motivates us to calculate the Reynolds number of melting snowflakes by linear interpolation w.r.t. between the mass equivalent dry snowflake and raindrop, instead of trying to come up with a more complicated model for the particle density or the geometry of a melting snowflake, which would be necessary to calculate N Re (D s , ) explicitly.Having an approximation for the Reynolds number of the particle the ventilation coefficient is calculated using the empirical formula of Hall and Pruppacher (1976) F v (D s , ) = 1.00 + 0.14 X 2 (D s , ), for X < 1 0.86 + 0.28 X(D s , ), else with Re (D s , ) and the Schmidt number N Sc = 0.64 consistent with M90.The resulting ventilation coefficient is shown in Fig. 1c as a function of the equivalent diameter.The ventilation coefficients for rain and snow are nearly equal, meaning that a dry snowflake has nearly the same ventilation coefficient as the mass equivalent raindrop.With this parametrization all terms of Eq. ( 7) are specified and therefore the melting process of an individual snowflake is completely described.

Bulk parametrization
In a scheme with a prognostic liquid water fraction, the mass m s of an individual snowflake is decomposed into an ice part m i and a liquid part m w .Therefore, we define for the bulk content of ice and meltwater of the wet snowflakes.Note that m i and m w are in general functions of D s since snowflakes of different sizes have different liquid water fractions.Here we assume that snowflakes smaller than D * have already melted completely and are no longer snowflakes, but raindrops, as discussed later.In addition, the size distribution of melting snowflakes f m (D s , ) becomes a function of the liquid water fraction.Using Eq. ( 7) we find the rate equation From the previous section, we already know C m (D s , ) and F v (D s , ), i.e., Eqs. ( 9) and (17).To evaluate the integral on the r.h.Using a bulk approach we predict only the bulk quantities L s,i and L s,w , not the size-resolved itself.Following SZ99, we assume They showed that approximately ∼ D − ˜ eq with ˜ = 1.3, cf.their Eq.( 18), and with D 3 eq ∼ D 2 s it follows = 2/3 ˜ = 0.87.This approach takes into account that smaller snowflakes melt faster than larger particles (e.g., Willis and Heymsfield, 1989).
To specify the size distribution of melting snowflakes we follow, again, SZ99 and use that in steady state the number flux at a certain melted diameter is constant with height due to mass conservation.Because D s is the maximum dimension of the dry mass equivalent snowflake, not the maximum dimension of the melting snowflake, D s is in fact a mass coordinate and number flux conservation can be written as In contrast to SZ99 we do not attempt to match both the snow distribution above and the raindrop distribution below the melting layer, because for size distributions of dry snow, which are exponential in geometric diameter rather than melted diameter, such a matching is not exactly possible using the flux conservation alone.Therefore, we decide to match the distribution of melting snow to the distribution of dry snow only.For the size distribution of dry snow we use, for simplicity, the local value of the snow water content, i.e., λ is a function of L s = L s,i + L s,w .Note that by doing so we chose to build a local scheme.Alternatively, one could try to estimate the snow content of the unmelted distribution above.In a time-dependent 3-D framework both approaches are, of course, approximations.
The critical diameter D * can be calculated iteratively from the local values of L s,i and L s,w for a given set of parameters N 0 and λ, i.e., from L s = L s,i +L s,w and N 0 one can calculate λ, and then D * iteratively from Eq. ( 18).This cumbersome inversion problem is maybe one reason why prognostic melting schemes have not become very popular so far.
Figure 2 illustrates the described size distributions for a snow mixing ratio of 3 g kg −1 and a value of 0.2 for the bulk liquid water fraction.The decrease of meltwater from smaller to larger snowflakes and the behaviour of the corresponding ice part illustrate well the effect of the size resolved liquid water fraction (D s ).Vertical lines mark the evolution of D * within one time step and therefore the conversion from snow to rain.The critical diameter of the melting snow size distribution D * represents the size of the smallest predicted snowflakes.Smaller snowflakes have already melted completely due to the fact that small snowflakes melt faster compared to larger ones.Therefore, D * increases during melting Using a bulk approach we predict only the bulk quantities L s,i and L s,w , not the size-resolved itself.Following SZ99 we assume They showed that approximately ∼ D − κ eq with κ = 1.3, cf.their Eq.( 18), and with D 3 eq ∼ D 2 s it follows κ = 2/3 κ = 0.87.This approach takes into account that smaller snowflakes melt faster than larger particles (e.g., Willis and Heymsfield, 1989).
To specify the size distribution of melting snowflakes we follow, again, SZ99 and use that in steady state the number flux at a certain melted diameter is constant with height due to mass conservation.Because D s is the maximum dimension of the dry mass equivalent snowflake, not the maximum dimension of the melting snowflake, D s is in fact a mass coordinate and number flux conservation can be written as In contrast to SZ99 we do not attempt to match both the snow distribution above and the rain drop distribution below the melting layer, because for size distributions of dry snow, which are exponential in geometric diameter rather than melted diameter, such a matching is not exactly possible using the flux conservation alone.Therefore we decide to match the distribution of melting snow to the distribution of dry snow only.For the size distribution of dry snow we use, for simplicity, the local value of the snow water content, i.e., λ is a function of L s = L s,i + L s,w .Note that by doing so we chose to build a local scheme.Alternatively, one could try to estimate the snow content of the unmelted distribution above.In a time-dependent 3-d framework both approaches are, of course, approximations.
The critical diameter D * can be calculated iteratively from the local values of L s,i and L s,w for a given set of parameters N 0 and λ, i.e., from L s = L s,i + L s,w and N 0 one can calculate λ, and then D * iteratively from Eq. ( 18).This cumbersome inversion problem is maybe one reason why prognostic melting schemes have not become very popular so far.
Figure 2 illustrates the described size distributions for a snow mixing ratio of 3 g kg −1 and a value of 0.2 for the bulk liquid water fraction.The decrease of meltwater from smaller to larger snowflakes and the behavior of the corresponding ice part illustrate well the effect of the size resolved liquid water fraction (D s ).Vertical lines mark the evolution of D * within one time step and therefore the conversion from snow to rain.The critical diameter of the melting snow size distribution D * represents the size of the smallest predicted snowflakes.Smaller snowflakes have already melted completely due to the fact that small snowflakes melt faster compared to larger ones.Therefore, D * increases during melting and this shift yields the conversion rate from snow to rain: (23) How we evaluate this integral will be discussed at the end of the next section.
Having specified the fall speeds of melting snow flakes v m (D s , ) in the previous section the calculation of the bulk sedimentation velocities is straightforward: with m s = m i + m w = αD 2 s .Note that the ice part and the meltwater of the melting snowflake ensemble have different sedimentation velocities, namely vs,i and vs,w .This is similar as for number and mass in a two-moment bulk microphysical scheme.In principle, this leads to a gravitational sorting with, e.g., wetter snowflakes falling faster than drier Fig. 2. Size distributions obtained with the new melting parametrization using a bulk liquid water fraction of 0.2 and a snow mixing ratio of 3 g kg −1 .The two values of D * (indicated by the vertical dashed lines) illustrate a shift of the critical diameter due to the conversion from snow to rain.and this shift yields the conversion rate from snow to rain: How we evaluate this integral will be discussed at the end of the next section.
Having specified the fall speeds of melting snowflakes v m (D s , ) in the previous section the calculation of the bulk sedimentation velocities is straightforward: Note that the ice part and the meltwater of the melting snowflake ensemble have different sedimentation velocities, namely vs,i and vs,w .This is similar as for number and mass in a two-moment bulk microphysical scheme.In principle, this leads to a gravitational sorting with, e.g., wetter snowflakes falling faster than drier snowflakes, but as we will see in the next section, this effect is small and probably negligible.

Implementation in the COSMO model
In the two-category mixed-phase cloud microphysical scheme which is currently used operationally in the COSMO model the parametrization of the melting process of snow (see also Sect.2.1) is based on an approximate form of the rate equation (Doms and Schättler, 2002) , where q v is the water vapour mixing ratio.Integrating the rate equation of a single melting snowflake over an exponential size distribution, as it is performed in the COSMO model microphysics and using v s = v s 0 D 1/2 as sedimentation velocity yields: with coefficients given by (29) Note that the COSMO model uses the snow mixing ratio q s = L s /ρ, instead of snow content L s .
In the new scheme, whose source code is available from the authors, the equation for the melting rate corresponding to Eq. ( 27) is formulated in terms of the ice mixing ratio q s,i of the snowflakes and reads with The coefficient M s melt , in the following called melting integral, includes the capacitance and the ventilation coefficient and is only a function of the total snow mixing ratio q s and the bulk liquid water fraction L = q s,w /q s .That makes it possible to pre-calculate this over the entire range of values of the two variables.The normalization factor N s melt can be interpreted as the average capacitance of the dry snowflakes which, in the old formulation, is incorporated in the parameter c s melt .To calculate the integrals in Eqs. ( 32), ( 24) and ( 25) we apply the finite difference method of Berry and Reinhardt (1974) using a logarithmic mass coordinate with 550 bins and mass doubling every 14 grid points.
In many microphysics schemes such coefficients are then stored in look-up tables (e.g., Walko et al., 1995;Meyers et al., 1997), which can be an efficient way of implementation.Alternatively, one can approximate the results of the numerical integration by an analytic, e.g., polynomial, function.We have chosen to follow the latter approach for two reasons, first, to avoid memory bottlenecks when accessing large look-up tables and, second, because functional fits can be easily published which makes it possible for other groups to reproduce our results, or even apply this parametrization in other microphysical schemes.To approximate, e.g., the melting integral, we chose Pade-type rational functions: with i + j ≤ n.To calculate the coefficients p ij and q ij a nonlinear least-square fit is performed using a Levenberg-Marquardt method (Press et al., 1992).For an efficient approximation we define the normalized snow coordinate ξ = 0.2 log (q s ) + 3 (35) where ξ varies between 0 and 1.For values of q s below 10 −3 g kg −1 we assume that ξ = 0, for q s > 100 g kg −1 we set ξ = 1.This coordinate is distributed uniformly over the logarithmic scale of q s .Figure 3 shows the resulting approximations with n = 3 for the melting integral and the sedimentation velocities as a function of ξ and L. The corresponding coefficients of Eq. ( 34) are presented in Table 1 and Table 2.
Figure 3a shows the results for the melting integral which decreases with an increasing bulk liquid water fraction, meaning that melting is most efficient for dry snowflakes.Figures 3c and 3d show the bulk sedimentation velocities of ice and meltwater.Both increase with ξ and L. The sedimentation velocity of meltwater is slightly larger than the one of ice but the difference is very small.According to the limitations of our assumptions and the underlying measurements the presented parametrization is only reasonable until the bulk liquid water fraction exceeds a value of 0.85 during the melting process.Beyond this threshold the partially melted snow will be converted completely into rain within the time step.
www.geosci-model-dev.net/6/1925/2013/To calculate the conversion from snow to rain the meltwater within one time step is given by and is normalized with the ice mass mixing ratio of snow yielding the normalized melted water Note that ζ ∈ [0, 1].Given the snow ice mass mixing ratio at the current time step q (t) s,i the value at the new time step, t + t, is and D (t+ t) * as well as q (t+ t) s,w can in principle be calculated by increasing D * until the ice mixing ratio of the distribution equals q (t+ t) s,i .The conversion term is then given by with and here q melt is an intermediate value of the meltwater mixing ratio, i.e., after melting but before the conversion to rain.Similar to the melting integral, this conversion coefficient M s conv , which is a function of ξ , L and ζ , can be pre-calculated numerically.Or in other words, by defining the normalized melt water ζ and calculating M s conv for all ζ ∈ [0, 1] we avoid the expensive and cumbersome calculation of D (t+ t) * and q (t+ t) s,w during the runtime of the host model.Because of the normalization of ζ the dependency on ξ (or q s ) is relatively weak and can be eliminated by averaging over ξ .The resulting parametrization which is a function of normalized meltwater ζ and bulk liquid water fraction L is shown in Fig. 3b and the corresponding coefficients of Eq. ( 34) are presented in Table 1.The conversion coefficient increases with the normalized melted water ζ , i.e., more meltwater is available for the conversion, and has the desirable property that M s conv = 1 for ζ = 1, i.e., in this case all where it refreezes completely.For melting snow in an area of temperatures below −2 • C the meltwater q s,w is therefore ascribed to the ice part q s,i within one time step.
The described parameterization is implemented into the COSMO model version 4.14.The new prognostic variable q s,w is included using a new generalized tracer implementation in the COSMO model (Roches and Fuhrer, 2012).For both mixing ratios, q s,w and q s,i , the same advection scheme is applied.As initial and boundary conditions for the new prognostic meltwater we use q s,w = 0 since this variable is not yet implemented in the data assimilation system and the large-scale models that provide the boundary conditions.

First results
As a first application of the new parameterization scheme, hindcast simulations are performed of a cold-season precipitation event in Central Europe in 2010, which featured a near surface snow melting layer, and of the high-impact wet snowfall event in northwestern Germany in November 2005 described in Frick and Wernli (2012).For both cases, COSMO simulations were run with the standard and the new snow melting scheme, respectively, with a horizontal resolution of 7 km and 40 vertical layers using ECMWF analyses as initial and boundary data.The simulation domain covers Central Europe, approximately from 14 • W to 22 • E and from 44 • N to 61 • N. In the following, a brief comparison of the simula- s with α = 0.069 for the geometry of snow, = 0.87 for the decay of the meltwater in the snow distribution, and N 0 = 8 × 10 6 m −4 for the intercept parameter of the distribution of dry snow aloft f D = N 0 exp(−λD s ).Shown are the analytic approximations using rational functions.The melting integral and the conversion coefficient are dimensionless quantities.
ice melts within the time step and all meltwater is immediately transferred to rain.
In the atmosphere, lifted melting layers occur frequently, meaning that melting snow reaches an area with temperatures where it refreezes completely.For melting snow in an area of temperatures below −2 • C the meltwater q s,w is therefore ascribed to the ice part q s,i within one time step.
The described parametrization is implemented into the COSMO model version 4.14.The new prognostic variable q s,w is included using a new generalized tracer implementation in the COSMO model (Roches and Fuhrer, 2012).For both mixing ratios, q s,w and q s,i , the same advection scheme is applied.As initial and boundary conditions for the new prognostic meltwater we use q s,w = 0 since this variable is not yet implemented in the data assimilation system and the large-scale models that provide the boundary conditions.

First results
As a first application of the new parametrization scheme, hindcast simulations are performed of a cold-season precipitation event in Central Europe in 2010, which featured a near surface snow melting layer, and of the high-impact wet snowfall event in northwestern Germany in November 2005 described in Frick and Wernli (2012).For both cases, COSMO simulations were run with the standard and the new snow melting scheme, respectively, with a horizontal resolution of 7 km and 40 vertical layers using ECMWF analyses as initial and boundary data.The simulation domain covers Central Europe, approximately from 14 • W to 22 • E and from 44 • N to 61 • N. In the following, a brief comparison of the simulations with the two schemes is presented for both cases.The aim of the comparison is to illustrate that the new scheme (i) produces wet snow instead of rain in some regions with surface temperatures slightly above 0  Shown is sea-level pressure (black contours, every 2 hPa) and potential temperature at 850 hPa (colors, in K).The black dot marks Dresden and the black rectangle the area shown in Figure 5.
tions with the two schemes is presented for both cases.The aim of the comparison is to illustrate that the new scheme (i) produces wet snow instead of rain in some regions with surface temperatures slightly above 0 • C, (ii) is able to simulate an atmospheric melting layer with a gradual transition from dry snow to melting snow to rain, (iii) leads to a decelerated snow melting process, and (iv) can provide novel and useful forecast information in situations of wet surface snowfall.An in-depth validation of the new snow melting scheme, based upon re-forecasts of a large set of situations, is beyond the scope of this study.

Simulation of a snow melting layer on 16 November 2010
The first event on 16 November 2010 featured snowfall in eastern Germany in a region with surface temperatures in the range of ±5 • C. The synoptic situation over Central Europe at 12 UTC 16 November 2010 was characterized by a pronounced upper-level trough, extending from the North Sea into the Mediterranean (not shown).The surface pressure field features a persistent high-pressure system over southern Scandinavia and an intense cyclone in the Gulf of Genoa (Fig. 4).Potential temperature at 850 hPa shows relatively cold air over France and Germany and a prominent frontal zone extending along the border between Germany and the Czech Republic and across Poland.It is along this frontal zone that the snow melting event occurred.Surface observations in Dresden (black dot in Fig. 4) indicate a 2-m temperature of 5.1 • C and 98% relative humidity with continous (moderate) rain at 12 UTC and light drizzle at 14 UTC.Additionally, radar observations in Dresden (not shown) provide evidence for precipitation in the frontal area and show an enhanced radar reflectivity indicating the likely occurrence of melting hydrometeors.
Results shown here are from a simulation started at 12 UTC 16 November 2010, i.e., just two hours prior to the time of the analysis of the snow melting event in eastern Germany.Such a short simulation has been chosen in order to minimize the feedback of the modified microphysical scheme on the dynamics of the event, which would render a comparison of the performance of the two schemes more difficult.We note however that earlier simulations started at 00 and 06 UTC 16 November 2010 reveal qualitatively similar results.At 14 UTC (Fig. 5), the horizontal section on model level 8 (about 920 hPa) shows a west-east oriented band of cold air with temperatures below 0 • C (see green contour) in eastern Germany.Both simulations produce snow in this area, but the simulation with the standard melting scheme produces less snow in the region to the north of the cold band (i.e., in an area with temperatures larger than 0 • C) compared to the simulation with the new melting scheme (compare Figs. 5a,b).The increased snow in this region mainly consists of partially melted snow as indicated by the snow meltwater mixing ratio shown by the red contours in Fig. 5b.
Vertical sections across the snow band, i.e., along the black lines shown in Fig. 5, clearly show the decelerated snow melting process due to the new parameterization (Fig. 6).The green line represents the 0 • C isoline, which marks approximately the top of the melting layer.For the bottom of the melting layer the isoline of 1.3 • C wet bulb temperature (dark green lines) appears to be a reasonable indicator.With the new scheme, snow penetrates to slightly lower levels (slightly below the 1.3 • C contour of wet bulb temperature) and, consistently, the rain contours are shifted to lower altitudes (compare Figs. 6a,b).The simulated meltwater in the simulation with the new melting scheme (Fig. 7a) indicates a qualitatively realistic representation of the snow melting process.The simulated melting layer has a vertical extension of about 40 hPa and is approximately sandwiched between the contours of 0 • C temperature and 1.3 • C wet bulb temperature.Of course, such a melting layer characterized by a transition from dry snow to wet snow to rain cannot be represented by the standard snow melting scheme.The vertical resolution of the presented COSMO simulations provides about 10 levels in the lowest 100 hPa.Therefore, the investigated melting layer extends over approximately 4 model levels.
The peak values of the meltwater mixing ratio exceed 0.04 g kg −1 at about the 940 hPa level (Fig. 7a).In contrast, the maximum of the meltwater fraction (black lines) occurs approximately 20 hPa below this level.This vertical shift can be understood as follows: consider snow falling in a layer with temperatures above the freezing level.As snow starts to melt the meltwater mixing ratio increases and reaches its maximum at a certain distance below the freezing level.Initially the snow mixing ratio decreases only slightly (i.e., the rain mixing ratio increases only slightly) due to the finite time it takes till the onset of the snow to rain conversion.Later, Shown is sea-level pressure (black contours, every 2 hPa) and potential temperature at 850 hPa (colours, in K).The black dot marks Dresden and the black rectangle the area shown in Fig. 5. from dry snow to melting snow to rain, (iii) leads to a decelerated snow melting process, and (iv) can provide novel and useful forecast information in situations of wet surface snowfall.An in-depth validation of the new snow melting scheme, based upon re-forecasts of a large set of situations, is beyond the scope of this study.

Simulation of a snow melting layer on 16 November 2010
The first event on 16 November 2010 featured snowfall in eastern Germany in a region with surface temperatures in the range of ±5 • C. The synoptic situation over Central Europe at 12:00 UTC, 16 November 2010 was characterised by a pronounced upper-level trough, extending from the North Sea into the Mediterranean (not shown).The surface pressure field features a persistent high-pressure system over southern Scandinavia and an intense cyclone in the Gulf of Genoa (Fig. 4).Potential temperature at 850 hPa shows relatively cold air over France and Germany and a prominent frontal zone extending along the border between Germany and the Czech Republic and across Poland.It is along this frontal zone that the snow melting event occurred.Surface observations in Dresden (black dot in Fig. 4) indicate a 2 m temperature of 5.1 • C and 98 % relative humidity with continuous (moderate) rain at 12:00 UTC and light drizzle at 14:00 UTC.Additionally, radar observations in Dresden (not shown) provide evidence for precipitation in the frontal area and show an enhanced radar reflectivity indicating the likely occurrence of melting hydrometeors.
Results shown here are from a simulation started at 12:00 UTC, 16 November 2010, i.e., just two hours prior to the time of the analysis of the snow melting event in eastern Germany.Such a short simulation has been chosen in order to minimise the feedback of the modified microphysical scheme on the dynamics of the event, which would render a comparison of the performance of the two schemes more difficult.We note however that earlier simulations started at 00:00 and 06:00 UTC, 16 November 2010 reveal qualitatively similar results.At 14:00 UTC (Fig. 5), the horizontal section on model level 8 (about 920 hPa) shows a west-east oriented band of cold air with temperatures below 0 • C (see green contour) in eastern Germany.Both simulations produce snow in this area, but the simulation with the standard melting scheme produces less snow in the region to the north of the cold band (i.e., in an area with temperatures larger than 0 • C) compared to the simulation with the new melting scheme (compare Fig. 5a and b).The increased snow in this region mainly consists of partially melted snow as indicated by the snow meltwater mixing ratio shown by the red contours in Fig. 5b.
Vertical sections across the snow band, i.e., along the black lines shown in Fig. 5, clearly show the decelerated snow melting process due to the new parametrization (Fig. 6).The green line represents the 0 • C isoline, which marks approximately the top of the melting layer.For the bottom of the melting layer the isoline of 1.3 • C wet bulb temperature (dark green lines) appears to be a reasonable indicator.With the new scheme, snow penetrates to slightly lower levels (slightly below the 1.3 • C contour of wet bulb temperature) and, consistently, the rain contours are shifted to lower altitudes (compare Fig. 6a and b).The simulated meltwater in the simulation with the new melting scheme (Fig. 7a) indicates a qualitatively realistic representation of the snow melting process.The simulated melting layer has a vertical extension of about 40 hPa and is approximately sandwiched between the contours of 0 • C temperature and 1.3 • C wet bulb temperature.Of course, such a melting layer characterised by a transition from dry snow to wet snow to rain cannot be represented by the standard snow melting scheme.The vertical resolution of the presented COSMO simulations provides about 10 levels in the lowest 100 hPa.Therefore, the investigated melting layer extends over approximately 4 model levels.
The peak values of the meltwater mixing ratio exceed 0.04 g kg −1 at about the 940 hPa level (Fig. 7a).In contrast, the maximum of the meltwater fraction (black lines) occurs approximately 20 hPa below this level.This vertical shift can be understood as follows: consider snow falling in a layer with temperatures above the freezing level.As snow starts to melt the meltwater mixing ratio increases and reaches its maximum at a certain distance below the freezing level.Initially the snow mixing ratio decreases only slightly (i.e., the rain mixing ratio increases only slightly) due to the finite time it takes till the onset of the snow to rain conversion.Later, i.e., at levels below the level of maximum snow meltwater mixing ratio, the snow to rain conversion becomes more efficient and both the snow and meltwater mixing ratios start Shown are snow mixing ratio (colors, in g kg −1 ), rain mixing ratio (black lines, in g kg −1 ), the 0 • C isoline (green), and the isoline of 1.3 • C wet bulb temperature (dark green) for the simulations with (a) the standard and (b) the new snow melting scheme.
i.e., at levels below the level of maximum snow meltwater mixing ratio, the snow to rain conversion becomes more efficient and both the snow and meltwater mixing ratios start to decrease, while the liquid water fraction of the snow still increases due to ongoing melting.

Simulation of wet surface snowfall on 25 November 2005
COSMO hindcast simulations for the wet snowfall event in northwestern Germany in November 2005 have been described by Frick and Wernli (2012).They showed that in the area affected by wet snow short-range COSMO hindcasts with the standard snow melting scheme produce a large amount of surface precipitation in liquid form (about 50%) even when capturing the low-tropospheric temperature structure fairly accurately.Here one of these hindcast simulations, initialized at 12 UTC 24 November, is repeated with the new snow melting scheme.Figure 8 shows the snow mixing ratio on the lowest model level at 03 UTC 25 November for both simulations with the standard and the new melting scheme, respectively.It is at this time that wet snowfall started at the monitoring station in Essen (Frick and Wernli, 2012).The overall distribution is fairly similar (compare Figs. 8a,b), but the domain integrated snow fraction, i.e., the ratio from surface snowfall to total precipitation, is increased by a few percents.In a fairly large area in Belgium, the Netherlands and western Germany, in a band of about 100 km to the north of the 0 • C isotherm, snow reaches the surface with a significant meltwater mixing ratio (see red contours in Fig. 8b), i.e., as wet snow.This is also revealed by the vertical section across the area with strongest (wet) surface snowfall (Fig. 7b), which shows the formation of a layer with non-zero Shown are snow mixing ratio (colors, in g kg −1 ), rain mixing ratio (black lines, in g kg −1 ), the 0 • C isoline (green), and the isoline of 1.3 • C wet bulb temperature (dark green) for the simulations with (a) the standard and (b) the new snow melting scheme.
i.e., at levels below the level of maximum snow meltwater mixing ratio, the snow to rain conversion becomes more efficient and both the snow and meltwater mixing ratios start to decrease, while the liquid water fraction of the snow still increases due to ongoing melting.

Simulation of wet surface snowfall on 25 November 2005
COSMO hindcast simulations for the wet snowfall event in northwestern Germany in November 2005 have been described by Frick and Wernli (2012).They showed that in the area affected by wet snow short-range COSMO hindcasts with the standard snow melting scheme produce a large amount of surface precipitation in liquid form (about 50%) even when capturing the low-tropospheric temperature structure fairly accurately.Here one of these hindcast simulations, initialized at 12 UTC 24 November, is repeated with the new snow melting scheme.Figure 8 shows the snow mixing ratio on the lowest model level at 03 UTC 25 November for both simulations with the standard and the new melting scheme, respectively.It is at this time that wet snowfall started at the monitoring station in Essen (Frick and Wernli, 2012).The overall distribution is fairly similar (compare Figs. 8a,b), but the domain integrated snow fraction, i.e., the ratio from surface snowfall to total precipitation, is increased by a few percents.In a fairly large area in Belgium, the Netherlands and western Germany, in a band of about 100 km to the north of the 0 • C isotherm, snow reaches the surface with a significant meltwater mixing ratio (see red contours in Fig. 8b), i.e., as wet snow.This is also revealed by the vertical section across the area with strongest (wet) surface snowfall (Fig. 7b), which shows the formation of a layer with non-zero Fig. 6.Vertical section across the snowfall region in eastern Germany at 14:00 UTC, 16 November 2010 (along the black line shown in Fig. 5).Shown are snow mixing ratio (colours, in g kg −1 ), rain mixing ratio (black lines, in g kg −1 ), the 0 • C temperature isoline (green), and the isoline of 1.3 • C wet bulb temperature (dark green) for the simulations with (a) the standard and (b) the new snow melting scheme.
to decrease, while the liquid water fraction of the snow still increases due to ongoing melting.

Simulation of wet surface snowfall on 25 November 2005
COSMO hindcast simulations for the wet snowfall event in northwestern Germany in November 2005 have been described by Frick and Wernli (2012).They showed that in the area affected by wet snow short-range COSMO hindcasts with the standard snow melting scheme produce a large amount of surface precipitation in liquid form (about 50 %) even when capturing the low-tropospheric temperature structure fairly accurately.Here one of these hindcast simulations, initialised at 12:00 UTC, 24 November, is repeated with the new snow melting scheme.Figure 8 shows the snow mixing ratio on the lowest model level at 03:00 UTC, 25 November for both simulations with the standard and the new melting scheme, respectively.It is at this time that wet snowfall started at the monitoring station in Essen (Frick and Wernli, 2012).The overall distribution is fairly similar (compare Fig. 8a and b), but the domain integrated snow fraction, i.e., the ratio from surface snowfall to total precipitation, is increased by a few percents.In a fairly large area in Belgium, the Netherlands and western Germany, in a band of about 100 km to the north of the 0 • C isotherm, snow reaches the surface with a significant meltwater mixing ratio (see red contours in Fig. 8b), i.e., as wet snow.This is also revealed by the vertical section across the area with strongest (wet) surface snowfall (Fig. 7b), which shows the formation of a layer with non-zero meltwater mixing ratio in the simulation with the new snow melting scheme, extending from the freezing level down to the surface.

Statistical analysis of the snow fraction
To quantify the effect of the new snow melting parametrization in a statistical way, a 24 h simulation with hourly output meltwater mixing ratio in the simulation with the new snow melting scheme, extending from the freezing level down to the surface.

Statistical analysis of the snow fraction
To quantify the effect of the new snow melting parameterization in a statistical way, a 24-hour simulation with hourly output has been performed for the snowfall event discussed in section 3.1, initialized at 00 UTC 16 November 2010.This simulation has been chosen in order to capture the entire snowfall that occurred during this day.The simulated precipitation is investigated in the main snowfall region indicated by the dashed box in Fig. 5.For each gridpont in this area associated with a snow mixing ratio larger than zero and a temperature above 0 • C the ratio of snow to total precipitation (snow fraction, q s /(q s + q r )) and the distance to the top of the melting layer (i.e., the 0 • C isosurface) are calculated.
The resulting relations for the standard and the new melting scheme are presented in Fig. 9.As expected both schemes show a decreasing snow fraction with increasing distance below the top of the melting layer.However, the slope of the relation is much steeper for the standard scheme, which indicates a more rapid conversion from snow to rain.Note that the relatively large scatter of the relation is due to the variability of humidity, precipitation intensity, and advection of precipitation at the different grid points.20 hPa below the top of the melting level the standard scheme reduces the snow fraction to typically less than 20%, while with the new Fig. 7. Vertical cross sections from simulations with the new snow melting scheme of the meltwater mixing ratio (colours, in g kg −1 ).Panel (a) shows the same section as Fig. 6 at 14:00 UTC, 16 November 2010.Also shown are the liquid water fraction (black contours for values of 0.2, 0.3 and 0.4), the 0 • C isoline of temperature (green contour), and the isoline of 1.3 • C wet bulb temperature (dark green contour).Panel (b) shows a section across the wet snowfall area (see black line in Fig. 8) at 03:00 UTC, 25 November 2005.meltwater mixing ratio in the simulation with the new snow melting scheme, extending from the freezing level down to the surface.

Statistical analysis of the snow fraction
To quantify the effect of the new snow melting parameterization in a statistical way, a 24-hour simulation with hourly output has been performed for the snowfall event discussed in section 3.1, initialized at 00 UTC 16 November 2010.This simulation has been chosen in order to capture the entire snowfall that occurred during this day.The simulated precipitation is investigated in the main snowfall region indicated by the dashed box in Fig. 5.For each gridpont in this area associated with a snow mixing ratio larger than zero and a temperature above 0 • C the ratio of snow to total precipitation (snow fraction, q s /(q s + q r )) and the distance to the top of the melting layer (i.e., the 0 • C isosurface) are calculated.
The resulting relations for the standard and the new melting scheme are presented in Fig. 9.As expected both schemes show a decreasing snow fraction with increasing distance below the top of the melting layer.However, the slope of the relation is much steeper for the standard scheme, which indicates a more rapid conversion from snow to rain.Note that the relatively large scatter of the relation is due to the variability of humidity, precipitation intensity, and advection of precipitation at the different grid points.20 hPa below the top of the melting level the standard scheme reduces the snow fraction to typically less than 20%, while with the new has been performed for the snowfall event discussed in Sect.3.1, initialized at 00:00 UTC, 16 November 2010.This simulation has been chosen in order to capture the entire snowfall that occurred during this day.The simulated precipitation is investigated in the main snowfall region indicated by the dashed box in Fig. 5.For each gridpont in this area associated with a snow mixing ratio larger than zero and a temperature above 0 • C the ratio of snow to total precipitation (snow fraction, q s / (q s + q r )) and the distance to the top of the melting layer (i.e., the 0 • C isosurface) are calculated.
The resulting relations for the standard and the new melting scheme are presented in Fig. 9.As expected both schemes show a decreasing snow fraction with increasing distance below the top of the melting layer.However, the slope of the relation is much steeper for the standard scheme, which indicates a more rapid conversion from snow to rain.Note that the relatively large scatter of the relation is due to the variability of humidity, precipitation intensity, and advection of precipitation at the different grid points.20 hPa below the top of the melting level the standard scheme reduces the snow fraction to typically less than 20 %, while with the new scheme still about half of the hydrometeors correspond to (partially melted) snow.In regions where this altitude corresponds to the surface, the two schemes produce a strikingly different surface precipitation signal with hardly any snow (standard scheme) versus much larger amounts of potentially heavy wet snow (new scheme).
With the standard scheme, virtually no snow exists at levels below about 30 hPa below the top of the melting level, whereas with the new scheme, the snow mixing ratio still amounts to about 30 % at this altitude.Below this level, i.e., in regions where the liquid water fraction attains the largest values, the snow to rain transition appears to be particularly slow with the new melting scheme.This is due to the fact that scheme still about half of the hydrometeors correspond to (partially melted) snow.In regions where this altitude corresponds to the surface, the two schemes produce a strikingly different surface precipitation signal with hardly any snow (standard scheme) versus much larger amounts of potentially heavy wet snow (new scheme).
With the standard scheme, virtually no snow exists at levels below about 30 hPa below the top of the melting level, whereas with the new scheme, the snow mixing ratio still amounts to about 30% at this altitude.Below this level, i.e., in regions where the liquid water fraction attains the largest values, the snow to rain transition appears to be particularly slow with the new melting scheme.This is due to the fact that the intercept parameter D * is already that large and the snow mixing ratio q s that small due to snow to rain conversion that the ensemble mainly consists of few large snowflakes that melt relatively slowly.
In general, the melting rate is reduced by the new melting scheme.This finding is also supported by the evolution of the ice fraction (q s,i /(q s + q r ), not shown) which has a slightly more compact density distribution but a slope comparable to the one of the snow fraction.
This statistical analysis indicates that the melting process is decelerated by the new parameterization and snow is able to reach much farther below the top of the melting layer (up to 70 hPa compared to 30 hPa with the standard scheme).

Conclusions
A new snowflake melting scheme has been introduced and implemented into the COSMO model version 4.14.The source code is available from the authors 1 and it is intended 1 Co-authors: axel.seifert@dwd.de,heini.wernli@env.ethz.ch to include the new scheme in later versions of the official COSMO code.The new melting parameterization allows the explicit prediction of partially melted snow (i.e., wet snow) using a prognostic meltwater mixing ratio.The formulations of capacitance, ventilation coefficient and fall velocity of a single melting snowflake follow mainly the empirical and theoretical findings of M90.For an ensemble of snowflakes, as described by the COSMO model, flux conservation is used to determine the size distribution of melting snowflakes.A size resolved liquid water fraction is defined following SZ99 and the meltwater and ice part of snow are treated as seperate parameters.The parameterization is implemented by using approximations in form of rational functions for the hydrometeor fall velocities, the melting integral, and the conversion rate from snow to rain.
First case studies of snowfall events associated with a melting layer illustrate the potential of the new snow melting scheme.Hindcast simulations with the new scheme of a frontal snowfall event in eastern Germany in November 2010 produce a melting layer with semi-melted snow particles in the lower troposphere with a vertical extension of about 40 hPa.The meltwater fraction attains peak values of about 0.3 at the bottom of the melting layer.In contrast, the peak of the meltwater mixing ratio is reached at a slightly higher altitude.Qualitatively, this simulated melting layer corresponds nicely to the "bright band" as observed by radar.At the surface, wet snow is predicted in fairly large areas with surface temperatures above 0 • C and, compared to the simulation with the standard scheme, the total area of surface snowfall is substantially larger.This indicates that with the new scheme, the snow melting process in the atmosphere is decelerated and the snow particles reach to lower levels compared to the standard scheme.The second case study of the high-impact wet snowfall event in northwestern Ger- the intercept parameter D * is already that large and the snow mixing ratio q s that small due to snow to rain conversion that the ensemble mainly consists of few large snowflakes that melt relatively slowly.
In general, the melting rate is reduced by the new melting scheme.This finding is also supported by the evolution of the ice fraction (q s,i / (q s + q r ), not shown) which has a slightly more compact density distribution but a slope comparable to the one of the snow fraction.
This statistical analysis indicates that the melting process is decelerated by the new parametrization and snow is able to reach much farther below the top of the melting layer (up to 70 hPa compared to 30 hPa with the standard scheme).

Conclusions
A new snowflake melting scheme has been introduced and implemented into the COSMO model version 4.14.The source code is available from the authors 1 and it is intended to include the new scheme in later versions of the official COSMO code.The new melting parametrization allows the explicit prediction of partially melted snow (i.e., wet snow) using a prognostic meltwater mixing ratio.The formulations of capacitance, ventilation coefficient and fall velocity of a single melting snowflake follow mainly the empirical and theoretical findings of M90.For an ensemble of snowflakes, as described by the COSMO model, flux conservation is used to determine the size distribution of melting snowflakes.A size resolved liquid water fraction is defined following SZ99 and the meltwater and ice part of snow are treated as separate parameters.The parametrization is implemented by using approximations in form of rational functions for the hydrometeor fall velocities, the melting integral, and the conversion rate from snow to rain.
1 Co-authors: axel.seifert@dwd.de, heini.wernli@env.ethz.chFirst case studies of snowfall events associated with a melting layer illustrate the potential of the new snow melting scheme.Hindcast simulations with the new scheme of a frontal snowfall event in eastern Germany in November 2010 produce a melting layer with semi-melted snow particles in the lower troposphere with a vertical extension of about 40 hPa.The meltwater fraction attains peak values of about 0.3 at the bottom of the melting layer.In contrast, the peak of the meltwater mixing ratio is reached at a slightly higher altitude.Qualitatively, this simulated melting layer corresponds nicely to the "bright band" as observed by radar.At the surface, wet snow is predicted in fairly large areas with surface temperatures above 0 • C and, compared to the simulation with the standard scheme, the total area of surface snowfall is substantially larger.This indicates that with the new scheme, the snow melting process in the atmosphere is decelerated and the snow particles reach to lower levels compared to the standard scheme.The second case study of the high-impact wet snowfall event in northwestern Germany in November 2005 (Frick and Wernli, 2012) shows a clear qualitative improvement of the simulation of the surface precipitation type when using the new snow melting scheme.With the new scheme, semi-melted snow is simulated at the surface in regions where wet snowfall was observed and simulations with the standard scheme produced surface rainfall.
These first results are promising and indicate that the new scheme can lead to improved forecasts of surface precipitation in critical situations with near-0 • C surface temperatures.Longer-term simulations and a more in-depth validation are required to quantify the potential improvement when using the new scheme.Such a validation should include observations from surface stations (precipitation type and amount) and also from (dual-polarimetric) radar.A comparison with radar data would allow an assessment of accuracy for the vertical structure of the simulated melting layer, which is sensitive to the assumptions made in the snow melting scheme.

Fig. 1 .
Fig. 1. a) Capacitance, b) fall velocity,and c) ventilation coefficient of a melting snowflake as a function of the equivalent diameter for various values of .Rain and snow correspond to = 1 and = 0, respectively.d) Empirical function Ψ( ) that describes the transition between the fall speed of a snow flake and a rain drop.Shown are the approximate relationship Eq. (15) and the original data (grey dots) ofMitra et al. (1990).

Fig. 1 .
Fig. 1.(a) Capacitance, (b) fall velocity, and (c) ventilation coefficient of a melting snowflake as a function of the equivalent diameter for various values of .Rain and snow correspond to = 1 and = 0, respectively.(d) Empirical function ( ) that describes the transition between the fall speed of a snowflake and a raindrop.Shown are the approximate relationship Eq. (15) and the original data (grey dots) ofMitra et al. (1990).
s. we have to further specify 1. the liquid water fraction of individual snowflakes (D s ); C. Frick et al.: A bulk parametrization of melting snowflakes 3. the threshold diameter D * of the distribution.

6
Frick et al.: A bulk parameterization of melting snowflakes 1. the liquid water fraction of individual snow flakes (D s ) 2. the size distribution of melting snowflakes f m (D s , ) 3. the threshold diameter D * of the distribution.

Fig. 2 .
Fig.2.Size distributions obtained with the new melting parameterization using a bulk liquid water fraction of 0.2 and a snow mixing ratio of 3 g kg −1 .The two values of D * (indicated by the vertical dashed lines) illustrate a shift of the critical diameter due to the conversion from snow to rain.

Fig. 3 .
Fig. 3. Parameterizations of (a) melting, (b) conversion, and the mass-weighted sedimentation velocities of (c) ice and (d) liquid mass.Here we have assumed ms = αD 2s with α = 0.069 for the geometry of snow, κ = 0.87 for the decay of the meltwater in the snow distribution, and N0 = 8 × 10 6 m −4 for the intercept parameter of the distribution of dry snow aloft fD = N0 exp(−λDs).Shown are the analytic approximations using rational functions.The melting integral and the conversion coefficient are dimensionless quantities.

Fig. 3 .
Fig. 3. Parametrizations of (a) melting, (b) conversion, and the mass-weighted sedimentation velocities of (c) ice and (d) liquid mass.Here we have assumed m s = αD 2s with α = 0.069 for the geometry of snow, = 0.87 for the decay of the meltwater in the snow distribution, and N 0 = 8 × 10 6 m −4 for the intercept parameter of the distribution of dry snow aloft f D = N 0 exp(−λD s ).Shown are the analytic approximations using rational functions.The melting integral and the conversion coefficient are dimensionless quantities.

Fig. 4 .
Fig. 4. ECMWF analysis fields at 12 UTC 16 November 2010.Shown is sea-level pressure (black contours, every 2 hPa) and potential temperature at 850 hPa (colors, in K).The black dot marks Dresden and the black rectangle the area shown in Figure5.

Fig. 4 .
Fig. 4. ECMWF analysis fields at 12:00 UTC, 16 November 2010.Shown is sea-level pressure (black contours, every 2 hPa) and potential temperature at 850 hPa (colours, in K).The black dot marks Dresden and the black rectangle the area shown in Fig.5.

Fig. 5 .Fig. 6 .
Fig. 5. Snow mixing ratio (colors, in g kg −1 ), meltwater mixing ratio (red contours, in g kg −1 , only in panel (b)), and the 0 • C isoline (green contour) on model level 8 at 14 UTC 16 November 2010 for the simulations with (a) the standard and (b) the new snow melting scheme.The dashed box marks the main area of interest while the black line indicates the location of the vertical cross sections shown in Figs. 6 and 7a.a) standard melting scheme b) new melting scheme

Fig. 5 .Fig. 5 .Fig. 6 .
Fig. 5. Snow mixing ratio (colours, in g kg −1 ), meltwater mixing ratio (red contours, in g kg −1 , only in panel b), and the 0 • C isoline (green contour) on model level 8 at 14:00 UTC, 16 November 2010 for the simulations with (a) the standard and (b) the new snow melting scheme.The dashed box marks the main area of interest while the black line indicates the location of the vertical cross sections shown in Figs. 6 and 7a.

Fig. 7 .Fig. 8 .
Fig. 7. Vertical cross sections from simulations with the new snow melting scheme of the meltwater mixing ratio (colors, in g kg −1 ).Panel (a) shows the same section as Fig. 6 at 14 UTC 16 November 2010.Also shown are the liquid water fraction (black contours for values of 0.2, 0.3 and 0.4), the 0 • C isoline of temperature (green contour), and the isoline of 1.3 • C wet bulb temperature (dark green contour).Panel (b) shows a section across the wet snowfall area (see black line in Fig. 8) at 03 UTC 25 November 2005.a) standard melting scheme b) new melting scheme

Fig. 7 .Fig. 8 .
Fig. 7. Vertical cross sections from simulations with the new snow melting scheme of the meltwater mixing ratio (colors, in g kg −1 ).Panel (a) shows the same section as Fig. 6 at 14 UTC 16 November 2010.Also shown are the liquid water fraction (black contours for values of 0.2, 0.3 and 0.4), the 0 • C isoline of temperature (green contour), and the isoline of 1.3 • C wet bulb temperature (dark green contour).Panel (b) shows a section across the wet snowfall area (see black line in Fig. 8) at 03 UTC 25 November 2005.a) standard melting scheme b) new melting scheme

Fig. 8 .
Fig. 8. Snow mixing ratio (colours, in g −1 ), meltwater mixing ratio (red contours for 2, 6, 8, 10 mg kg −1 , only in panel b), and the 0 • C isoline (green contour) on model level 1 at 03:00 UTC, 25 November 2005 for the simulations with (a) the standard and (b) the new snow melting scheme.The black line indicates the location of the vertical cross section shown in Fig. 7b.

FrickFig. 9 .
Fig. 9.The frequency distribution of the ratio of snow to overall precipitation (snow fraction) for all gridpoints with snow and temperatures above 0 • C as a function of the vertical distance below the top of the melting layer, (a) for the standard and (b) the new snow melting scheme.Results are from a November 2010 hindcast simulation (see text).

Fig. 9 .
Fig. 9.The frequency distribution of the ratio of snow to overall precipitation (snow fraction) for all gridpoints with snow and temperatures above 0 • C as a function of the vertical distance below the top of the melting layer, (a) for the standard and (b) the new snow melting scheme.Results are from a November 2010 hindcast simulation (see text).

Table 1 .
The coefficients used within the Pade-type rational function (Eq.34) for the melting integral and the conversion coefficient.

Table 2 .
The coefficients used within the Pade-type rational function (Eq.34) for the two sedimentation velocities.