## Abstract

To obtain the full description of the dynamical and microphysical finescale structures required for the computation of the radar-derived brightband parameters, a numerical model has been developed. A bulk microphysics module was introduced into a nonhydrostatic, fully compressible dynamic framework. A microphysical parameterization scheme, with five water categories (vapor, cloud water, snow, melting snow, and rain), describes the interactions related to the evolution of the melting layer (melting and diffusional exchanges of mass of each hydrometeor category). Dynamic, thermodynamic, and microphysical processes are fully coupled. The main characteristics of the bulk parameterization scheme for melting of snow are the following: 1) wet snow is described by its water content and by an additional prognostic variable, namely, the diameter of the smallest snowflake not yet completely melted; 2) the fall velocity of the melting snowflakes is based on the laboratory observations; and 3) a size-dependent ventilation coefficient of the melting particles is used. With this new formulation of the melting process some approximate analytical relations between variables that characterize the melting layer are obtained. The results of simulations show that the nonuniformity of the snow content causes horizontal variability of various atmospheric properties within the melting layer, which leads to the generation of convective cells. The impact of the induced finescale dynamics on the microphysical structure within the melting zone is analyzed.

## 1. Introduction

Most midlatitude rainfall is initiated through the ice formation process followed by melting, and it falls as widespread stratiform rain. The layer of the atmosphere where the transition from the solid to the liquid phase takes place is a region of important discontinuity of the radar-measured properties. When the transition takes place at a well-defined height it appears as a layer of enhanced reflectivity that is called bright band.

Since the beginning of quantitative radar measurements, the bright band has been considered a source of error in precipitation estimates, and the need to remove its effects has been recognized repeatedly (Smith 1986;Klaassen 1988; Joss and Waldvogel 1990; Fabry and Zawadzki 1995; Hardaker et al. 1995). The effect of melting particles on electromagnetic wave propagation, especially on the signal attenuation and on the increasing scattered power, is important in microwave communication engineering (Battan 1973; Bellon et al. 1997). Brightband characteristics are related to the type of precipitation particles and are suggestive of the processes involved in precipitation formation (Atlas 1957;Fabry and Zawadzki 1995; Huggel et al. 1996). For example, contrary to riming, snowflake aggregation results in the bright band’s increased prominence. The dynamic effects of the melting-induced air cooling may be important on the microscale, convective scale, and mesoscale.

The study of the complex physical phenomena involving the coupling of particle melting with the dynamical and thermodynamical processes within the bright band benefits from different approaches. Laboratory studies (Knight 1979; Matsuo and Sasyo 1981a;Fujiyoshi 1986; Oraltay and Hallett 1989; Mitra et al. 1990; Fujiyoshi and Muramoto 1996) give valuable information about the melting process at the microscale. Some aircraft penetration data (Stewart et al. 1984, 1996; Willis and Heymsfield 1989; Raga et al. 1991) and measurements on the ground (Ohtake 1969, 1970; Yokoyama et al. 1985; Waldvogel et al. 1993) within and in the vicinity of the phase-transition layer provide direct but limited samples of a highly variable process. Radar observations (Takeda and Fujiyoshi 1978; Yokoyama et al. 1984; Caylor and Illingworth 1989; Fabry and Zawadzki 1995) provide quantitative data in time and space with a good resolution. To help their interpretation, radar data can be combined with in situ measurements (Stewart et al. 1984; Zrníc et al. 1993; Hagen et al. 1993; Barthazy et al. 1996). Doppler capabilities can be exploited to gain insight on particle size distributions from observations of the Doppler spectra (Drummond et al. 1996). The radio acoustic sounding system capability can be used to derive temperature profiles in the melting layer. However, the full understanding of the physical processes involved cannot be obtained by observation only. For example, the vertical distribution of the fraction of melted water in the precipitation particles cannot be derived unambiguously from observations. Numerical models, on the other hand, are capable of a detailed description of the melting layer if they incorporate sufficient physical constraints;those obtained from radar observations can be particularly helpful.

Two distinct modeling avenues have been explored in the past. In the first, a one-dimensional vertical simulation using quasi-stochastic equations allows the analysis of microphysical details of the evolution of an assumed initial snowflake distribution during the precipitation fall through the melting zone. Thermodynamic conditions of the environment may be made to interact with the microphysical processes. The features of the melting hydrometeors predicted by these models have been verified by in situ measurements (Matsuo and Sasyo 1981b; Willis and Heymsfield 1989), as well as by a comparison of the model-generated reflectivity profiles with either observational data or profiles produced by other similar models (Ekpenyong and Srivastava 1970;Yokoyama and Tanaka 1984; Yokoyama and Ishizaka 1985; Klaassen 1988; Hardaker et al. 1995). Due to the inability of the 1D model to account for horizontal advection, the effects arising from the horizontal variability of the conditions controlling the phase change cannot be adequately represented.

In the second approach, a two- or three-dimensional dynamical model is coupled with a microphysical module. This representation is capable of capturing the distributions of both dynamical and microphysical fields in the region of interest. However, in this type of model, with a larger simulated domain, a full particle size spectrum is computationally prohibitive. The microphysical processes are thus simplified into a bulk parameterization. Typically, these numerical models have been used to investigate the effects of melting particles’ removal of heat from the atmosphere on mesoscale circulation (Szeto et al. 1988a,b; Barth and Parsons 1996). A feedback mechanism between melting microphysics and the orographic flow was examined with a two-dimensional model by Marwitz and Waight (1986) and Wei and Marwitz (1996).

In addition, some theoretical models (Moore and Stewart 1985; Stewart and Lin 1992) and field experiments (Atlas et al. 1969, Stewart and Macpherson 1989;Stewart et al. 1996) highlight the nonnegligible influence of melting of ice-phase precipitation on the mesoscale dynamics of the atmosphere. A general review of the effects of the melting of snowflakes in the atmosphere can be found in Heffernan and Marwitz (1996).

The spatial variations of precipitation causing the nonuniform cooling by melting and consequent pressure perturbations are held responsible for the small mesoscale wind perturbations within the melting layer (Atlas et al. 1969). Harrold and Browning (1967), Heymsfield (1979), and Carbone (1982) have reported similar observations. The development of organized mesoscale circulation in response to a nonuniform melting-induced cooling has been shown using a simple model by Lin and Stewart (1986) and with a two-dimensional model by Szeto et al. (1988a,b).

In turn, the melting process may be affected by air dynamics, including the generated small-scale perturbations embedded in the larger-scale flow: the saturation conditions associated with the vertical motion determine the latent heat accompanying the vapor diffusion on precipitating particles. This heat is available for melting in addition to sensible heat transferred from the ambient air. Thus, the rate of melting depends critically on the condensation/evaporation rate. Matsuo and Sasyo (1981b) have showed theoretically the importance of the relative humidity of air for the melting process. Using a kinematic microphysical model, Gedzelman and Arnold (1993) pointed out the effects of air relative humidity on the form of precipitation on the ground. Donaldson and Stewart (1993) included in their microphysical model the diffusional growth of partially melted particles and showed that the trajectories of these particles vary depending on the degree of environmental saturation. Hanesiak and Stewart (1995) point out that the saturation condition aloft can alter precipitation type at the surface and that heating linked with adiabatic descent considerably contracts diabatic cooling.

Additional feedback mechanism may come from the fact that the melting zone of the precipitating stratiform region may be associated with the production of cloud liquid water (Stewart et al. 1984; Lin and Stewart 1986;Houze 1993) and consequent precipitation enhancement. This liquid water within the melting zone may be generated by local convection (as well as uplift associated with a larger-scale forcing) together with the melting-induced air cooling.

The main objective of this study is to develop a model as complete as possible of the stratiform precipitation melting region. We need a detailed description of microphysical structures for the computation of the distribution of the radar-derived parameters so that a quantitative comparison between model results and observations can be made. This requires a realistic partition of particles into rain and melting snow, the profiles of the particle size distributions, and the profiles of the melted fraction in each melting particle. The simulated radar measurements, together with the details of the developed electromagnetic scattering computations, are presented in a companion paper (Fabry and Szyrmer 1999, hereafter FS).

Through the analysis of the model results presented here we search to better understand the causes and consequences of the small-scale variability in microphysical and dynamic properties. The model, with its finescale processes, bridges the gap between the one-dimensional modeling, with detailed microphysics but without interacting dynamics, and the mesoscale models, with realistic dynamics but with a rather crude description of the melting process and/or incomplete microphysics. It incorporates within a dynamic framework a new parameterization scheme of melting process. The scheme also allows an analytical description of the process of melting of the entire population of the precipitating snowflakes. Realistic simulations are obtained here by introducing radar-observed snow at the upper boundary of the model. Air circulation is induced by the horizontal nonuniformity of the snow that produces a spatially variable melting rate and the consequent heat exchange between precipitation and ambient air.

This paper is organized as follows. The next section describes the parameterization scheme for melting snow and some theoretical analytical solutions derived from this scheme. The general features of the model are described in section 3. The results of the numerical simulations and comparisons with radar data are given in section 4. Concluding remarks close the presentation.

## 2. Parametric description of the melting of snow

In this section we will describe the main assumptions used to derive our parameterization of the melting process and give the expressions relating particle density, fall velocity, and ventilation coefficient to their size and degree of melting. From this, we will describe the melting process of individual particles and of a population of particles.

### a. Main assumptions

The first assumption is that melting snow will be considered as wet ice inclusions in air with most of the melted water flowing toward the center of the particle. Thus, the volume of the melting particle is determined by its solid component. When in the course of melting, the volume of water equals that of ice, that is, the ice structure is completely embedded within water, the particle is transferred from the melting snow category into the rain category. When the ice skeleton is completely surrounded by liquid, the remaining ice is considered to be unimportant for the heat exchange with the environment and for the modeled vertical structure of radar reflectivity.

Persistence of the ice frame during most of the melting and the absence of a water shell around the large melting particles, all but in the final stage, was observed in field work as well as in wind tunnel laboratory experiments. The melted water is held in the snowflake interior while the ice protrudes out (Knight 1979; Fujiyoshi 1986), or water penetrates into the ice through small gaps and pores (Matsuo and Sasyo 1981a; Mitra et al. 1990). Therefore, observations indicate that the melting particle volume is determined by the solid ice structure during most of the melting. However, we are aware that the melting behavior of different individual ice crystals and snowflakes may present differences depending on snow type and on environmental conditions (Oratlay and Hallett 1989). Our description probably applies better to aggregates than to single crystals.

Second, a one-to-one correspondence between the snowflake falling across the 0°C level and the raindrop into which it melts is assumed. This implies that aggregation and breakup in the melting layer are neglected.

Despite the many field observations, theoretical works and developed models of the microphysical processes in the melting layer, the question of the importance of aggregation and breakup still remains unresolved and deserves a detailed discussion here. In their laboratory study, Mitra et al. (1990) rarely observed spontaneous breakup of melting dendrite aggregates; breakup seems to occur only in snowflakes with a strongly asymmetric mass distribution. The same conclusion can be drawn from other studies (Knight 1979;Fujiyoshi 1986). The one-to-one correspondence between snow particles and raindrops is suggested by the field observations of Du Toit (1967) and Ohtake (1969). However, other in situ observations show the larger particles in the middle portion of the melting. It is possible that these contradictory findings are due to the intermittent nature of the turbulence in the melting layer that could induce localized aggregation or breakup. We believe that, *on the average,* the aggregation of snowflakes that may occur in the melting layer is not increased by the melting process. In fact, less aggregation should be expected in the melting layer than in the kilometer just above it where snow becomes sticky and collection efficiency is already close to unity. These effects are well illustrated by the observations of Fabry and Zawadzki (1995). In addition, results of computations (not shown here) of the collision rate between particles falling through the melting making use of the quasi-stochastic collection equation (with collision efficiency of 1) indicate that only in the zone just below the top of the melting layer is aggregation comparable to the one in snow just above, and then falls rapidly thereafter. In any case, long-term radar observations of weak and moderately intense bright bands by Fabry and Zawadzki (1995) indicate that the combined effect of aggregation and breakup, although certainly present, on the average accounts for less than 1 dB of change in reflectivity from snow above to rain below the melting layer. This implies that the flux of the sixth moment of the particle size distribution is conserved. Mass flux is also conserved and, thus, at least two parameters of the particle size distribution are the same above and below the melting layer. Therefore, we must conclude that either aggregation and breakup are negligible or they compensate each other closely. The latter is quite unlikely. The situation may be different in bright bands of extreme intensity and depth.

### b. Density of snow and of melting snow, and fraction of melted water

Atmospheric ice occurs in a variety of habits and covers a wide range of particle size with various degrees of riming and aggregation. Describing a generalized snow particle is difficult and we must rely on average expressions for particles, especially the dendrite- and plate-type snowflakes that are most common (Jiusto and Bosworth 1971).

Snowflake bulk density *ρ*_{s}, defined as the particle’s mass *m* divided by the volume of a sphere of diameter, *D*_{s}, follows a power-law relationship:

with *y*_{s} < 0 according to a generally observed decrease of the snowflake density with size.

The diameter of the melted snowflakes, *D*_{w}, is related to that of the dry snowflakes, *D*_{s}, by mass conservation and (1):

As the snowflake melts, its density varies from the value *ρ*_{s} to a value approaching that of water, *ρ*_{w}. According to our assumptions, we consider the change in mass of a melting particle via diffusional growth to be negligible. Therefore, the melted mass *m*_{w} and the density of melting particles *ρ*_{m} are given by

where *V*_{s}, *V*_{m}, and *V*_{r} denote the volume of dry, melting, and completely melted snowflakes, respectively.

The degree of melting for an individual particle is described by the melted fraction *f* ≡ *m*_{w}/*m.* Its maximum, *f*_{max}, is attained when the water covers the ice entirely, that is, *V*_{m} = *V*_{w}; then the particle is considered to be a raindrop. Using (2), it can be found that

### c. Fall velocities

Fall speed power-law relationships for various types of natural snow can be found in the literature in terms of maximum or mean dimension, melted diameter, or snowflake mass (among others, in Langleben 1954; Magono and Nakamura 1965; Jiusto and Bosworth 1971; Locatelli and Hobbs 1974; Kajikawa 1989). On one hand, the numerical values of the constants vary rather considerably: crystals and their aggregates can exhibit different modes of falling. On the other hand, it is well known that most snowflakes fall at about 1–1.5 m s^{−1}. This narrow range of average fall velocity suggests a rather weak dependence on parameters such as shape, density, and size of particles. We have adopted the snowflake fall velocity relation in terms of the melted diameter, as proposed by Langleben (1954) for mixture of dendrites and aggregates of plates and given by

The Gunn and Kinzer’s (1949) terminal velocity data for raindrops can be approximated by a power-law function suggested by Sekhon and Srivastava (1971):

Thus, the values of the exponent in the fall speed relation for raindrops and snowflakes, described by their melted diameter, are very close; this allows for some simplifications in the description of the size distribution of melting snow particles as we will see below.

To account for the increase in particle fall speed with height, a factor, (*ρ*_{0}/*ρ*)^{0.5}, is applied to (4) and (5).

During melting, the fall velocity of melting snowflakes increases as water accumulation increases their bulk density. The fall velocities of melting snowflakes, *U*_{m}, have been studied experimentally by Mitra et al. (1990). They define the following function of *f*: *y*(*f*) = (*U*_{m} − *U*_{s})/(*U*_{r} − *U*_{s}). From a given *y*(*f*), the value of *U*_{m} can be computed as a function of particle size and its melted fraction. Moreover, since the exponents in the fall speed–size relationship for rain and snow particles [(4) and (5)] are very close, *U*_{s}(*D*_{w}) can be approximated by (*α*/*a*)*U*_{r}(*D*_{w}), and *U*_{m}(*D*_{w}, *f*) can be computed from the approximated relation,

where

The following simple form of *g*(*f*) assures a good fit to the measurements of Mitra et al. (1990):

with *α*/*a* = 4.6 and *c*_{g} = 0.5[(*α*/*a*) − 1] = 1.8.

Figure 1 shows *U*_{m} as a function of the melted mass fraction, as predicted by (6) [with *U*_{r}(*D*_{w}) given by (5)], for some selected *D*_{w}. The values derived from the points of the curve of Mitra et al. (1990) are also shown. We can see that the simplified expression(6) with (7) provides a very good approximation to experimental data.

### d. Ventilation coefficients

The ventilation coefficients depend on particle size and on physical characteristics of the particle. Thus, the value of these coefficients changes constantly during melting, and the description of the microphysics of melting may be sensitive to the ventilation coefficients for the melting snowflakes, *F*_{m}.

According to Pruppacher and Klett (1997) the ventilation coefficients for rain and snow can be expressed in terms of the Reynolds number Re as

for Re greater than 2.5.

To the best of our knowledge, there are no studies of the ventilation coefficient of melting snowflakes. In the parameterization of melting it is customary to adopt some average value for *F*_{m}. However, given this coefficient’s strong dependence upon particle size, we prefer to adopt an approximate expression for *F*_{m} that is a function of its melted diameter *D*_{w} and of the degree of melting described by *D*_{m}. The chosen approximated form, following Papagheorghe (1996), is given by

The numerical constants *A*_{m} and *B*_{m} have been chosen to fit to the values of *F*_{r} and *F*_{s} expressed in terms of melted diameter assuming a constant snow density of 0.1 g cm^{−3}. Figure 2 shows the ventilation factor for melting particles at different stages of melting together with the ventilation coefficient of raindrops and snowflakes as described above. Given the uncertainties in snowfall velocity, shape, and ventilation coefficient in general, (8) seems to be an adequate approximation, providing realistic values of the ventilation coefficient for melting snowflakes with mean bulk density; it also allows a concise description of the melting process within the constraints of bulk parameterization.

### e. Melting process of an individual snowflake

The rate at which the ice structure melts can be calculated through the condition of the balance of heat at the particle surface. We neglect the following heat transfers.

The transfer of heat by conduction through the possible water layer on the particle surface. Since most of the melted water seems to be trapped inside the particle, heat passes directly from the surrounding air to the unmelted ice frame.

The storage of sensible heat within the particle. We assume that until the melting is completed the temperature of the mixture of ice and water within the particle remains 0°C.

The transfer of heat by radiation that is extremely small compared to other transfers.

Therefore, the latent heat used by the melting particle to melt the ice has to be balanced by the heat transferred through conduction and by the latent heat released/removed at the particle surface by condensation/evaporation:

with

and

where *δT* and *δρ*_{υ} describe the temperature and vapor density differences between the air and the surface of the melting particle (assumed at 0°C and at water-saturated conditions). The substitution of the approximated form of *F*_{m}, given by (8), in (10), and (11), makes the expressions independent of the degree of melting. Introducing (10) and (11) in (9) leads to the following expression for the rate of mass melting:

The term *δQ,* a measure of heat available for melting per unit time, is defined as

The two terms correspond to the conduction of heat from the warmer surrounding air and to the latent heat released onto the melting particles during water vapor condensation.

If the air is saturated, *δρ*_{υ} can be approximated for 0°C < *T* < 4°C by a linear function of *δT*:

with *η* = 0.340 g m^{−3} K^{−1} for a temperature of 0.5°C and *de*_{s}/*dT* approximated by a constant equal to 45.8 Pa K^{−1} (from the Clausius–Clapeyron equation for this temperature). For temperatures lower than 1.5°C, the approximation is within 3%; the error becomes greater than 10% for *T* > 4°C. With (14) at saturation *δQ* becomes the following function of *δT* (equal to air temperature in °C):

where *q*_{sat} ≡ *K* + *D*_{υ}*L*_{υ}*η* ≈ 4.75 · 10^{3} J cm^{−1} s^{−1} K^{−1}. The expressions used to evaluate *η* and *q*_{sat} in (15) are given in appendix B. It is interesting to note that at saturation, the ratio of heat transferred by conduction to that due to condensation of water vapor equals approximately *K*/*D*_{υ}*L*_{υ}*η* ≈ 0.5.

The rate of melting as given by (12) is an increasing function of the initial size of the particle, independent of the degree of melting. The liquid mass fraction of the particle can be obtained by integrating (12) over time from the onset of melting as

where [*δQ*Δ*t*] describes the heat exchanged from the beginning of the melting along the trajectory of precipitation.

Let us define a size *D*_{w} = *d*_{w} such that *f*(*d*_{w}) = *f*_{max}(*d*_{w}). At any given time, particles of this size are the smallest melting particles. Smaller particles have attained *f*_{max} earlier and consequently were transferred to the rain category. Thus, at any time, the melting snow size distribution is truncated at the diameter *d*_{w}.

Replacing *f*_{max} in (16) by (3), for *D*_{w} = *d*_{w}, the following relationship between the particle size and the heat required for its complete melting is obtained:

where

For the average density of *ρ*_{s} = 0.1 g cm^{−3}, *ζ* = 0.9. The approximation *ζ* ≈ 1 is good for larger particles and reasonable for the smaller ones, depending on the assumed constants *x*_{s} and *y*_{s} and on the assumed mass distribution within the snowflakes.

### f. Condensation effects on melting particles

With fallen distance the temperature difference between the air and the melting particles increases. Therefore, the difference of water vapor densities between the saturated vapor at the surface of the melting snowflake at 0°C and the vapor in the warmer environment, *δρ*_{υ}, also increases significantly. For example, the water-saturated air at +2.5°C gives a supersaturation of 20% on snowflakes composed of both solid and liquid at 0°C. It has been shown above that in the saturated air about a half of the heat used to melt comes from the heat released in the condensational growth of the melting hydrometeors. Thus, the heat due to condensation may have an important impact on the melting process, and its role has been already pointed out by Matsuo and Sasyo (1981b), Donaldson and Stewart (1993), and Hanesiak and Stewart (1995).

If the generated excess of water vapor does not exceed the saturation value, the smaller and completely liquid raindrops (and eventually cloud droplets if present) evaporate, and the condensation on the partially melted particles can be viewed as a mass transfer from the smaller raindrops to the larger melting snowflakes, similarly to the Bergeron process. (However, in this case, the condensational heat associated with mass transferred from rain is balanced by evaporative cooling.) The rate of the mass increase caused by the condensation on an isolated particle can be found from (9), (10), and (12) as

Thus, after melting is completed, *m*_{cond} = 0.068*m*; that is, the increase of mass during melting due to condensation of water vapor, *m*_{cond}, is less than 7% of the original mass of the particle, *m*; in the subsaturated environment, it may be much smaller.

It can be easily demonstrated that the error introduced by neglecting *m*_{cond} causes the calculated melted fraction to be underestimated by less than 6%, and that the error is at maximum just after the onset of melting.

We can also evaluate the increase of the reflectivity factor due to the condensation on the melting particles as (1.068)^{2} = 1.14 times, that is, 0.57 dB (neglecting the increase of fall velocity), which is in close agreement with the “educated guess” of 0.5 dB given by Fabry and Zawadzki (1995). Moreover, due to the removal of the water vapor by the melting particles, the evaporation of rain increases, so the gain of mass below the melting layer could be lower than calculated from the condensational growth. In contrast, as far as the reflectivity factor within the melting layer is concerned, the growth by condensation of the larger melting particles dominates the evaporation of the smaller raindrops.

It can be concluded that, in general, whereas the rate of the mass change via the vapor transfer during melting, as well as other related effects, is very small, the thermal effect due to diffusional growth significantly affects the melting process. An accurate description of the melting layer must take into account the condensation heat, and since the vertical motion is closely related to the air saturation conditions, its impact on the melting process in the atmosphere may be important.

### g. Size distributions of precipitation particles

For raindrops below the melting layer we assume the inverse exponential Marshall–Palmer distribution with the intercept parameter *N*_{or} = 8 × 10^{4} m^{−3} cm^{−1}. The slope of the distribution, *λ,* changes with the precipitation rate. At steady state, the assumption of no aggregation and of no breakup results in the conservation of the flux of the number concentration of particles for each size interval. Thus, the distribution of snowflakes in terms of melted diameter, *N*_{s}(*D*_{w}), just above the melting layer, is related to the raindrop distribution, *N*_{r}(*D*_{w}), below this layer, by

In the same way, the size spectrum of the mixed-phase particles within the melting region, written in terms of melted diameters, can be expressed as

For the population of partially melted particles, *N*_{m}(*D*_{w}) can be rewritten using (18) with *ζ* = 1 as

### h. Formulation of the parameterization of the melting process

The total mass content of the whole population of melting snow particles, truncated at *d*_{w}, is given by

containing the mass of melted water given by

obtained with the aid of (18). For a given slope of the distribution, *λ,* and using (22), the integrals on the right-hand side of the last two equations can be expressed in terms of the truncated moments of the distribution *N*_{r}(*D*_{w}), *F*^{(x)}_{r}(*d*_{w}), where *x* describes the moment’s order. The average melted fraction, 〈*f*〉 ≡ *M*_{w}/*M*_{m}, consequently can be expressed as

Thus, for a given *λ* and for an assumed snow density, the melted mass fraction can be calculated from the varying in time and space the cutoff diameter *d*_{w}.

The total rate of melting, obtained by integrating of (12) over all the particles in the melting snow category, can be written with aid of (23) and (24) as

where *δQ* is a measure of the rate of heat available locally, while *d*_{w} is directly related to the measure of the total heat used for melting from the beginning of the process by (17).

It should be noted that the melting rate describes here the mass transfer between two phases within the melting particle category. The rate at which the mass is transferred from the melting snow category to the rain category, −*dM*_{m}/*dt* = *dM*_{r}/*dt,* is determined by the increment of the minimal diameter that can be obtained from (17):

where the value of *δQ* is assumed constant during the time step Δ*τ.*

The mass-weighted terminal velocity of the melting snow, 〈*U*_{m}〉, is calculated using the assumption of mass flux conservation:

## 3. Numerical model and experimental design

Melting layer processes are simulated using a dynamic model in which dynamics, thermodynamics, and microphysics are fully coupled. The dynamic framework is the same as that of the Canadian Mesoscale Community Compressible Model (MC2) described by Tanguay et al. (1990). Briefly, the integration of the fully elastic set of equations (three equations for momentum and equations for pressure and temperature deviations) uses a semi-implicit marching algorithm designed to effectively slow down the high-frequency acoustic waves. It is combined with a three-time-level semi-Lagrangian advection scheme with cubic spatial interpolation. The microphysical parameterization that has been adapted to this advection scheme introduces only a small numerical dispersion and a tolerable amount of diffusion (Ostiguy and Laprise 1990).

### a. General features of the microphysics scheme

The representation of microphysical interactions uses a bulk physical representation that here includes five water categories: vapor, cloud water (with a monodisperse distribution), rain, snow, and melting snow. Each category is characterized by its mixing ratio as a prognostic variable. However, melting snow is described by an additional advected variable, the minimum diameter of melting snowflakes, *d*_{w}, directly related to the total heat exchange and used to diagnose the degree of melting of the mixed phase precipitation. The coupling between dynamics and microphysics is done through the conservation equation for each microphysical variable.

The source terms result from the following active microphysical processes: (a) snow melting and (b) diffusional mass exchanges of all the modeled forms of condensed water. These mass transfers are determined by requiring that the flux of sensible heat between the particle and its surrounding air be balanced by the flux of the latent heat associated with the phase change, in accordance with the description of the condensation/evaporation of rain and the deposition/sublimation of snow given in Zawadzki et al. (1993). The changes in temperature accompanying all these processes act as sources for the dynamic fields. Since the generated cloud liquid water quantity is rather small (only concentrated in the more vigorous updrafts with mixing ratio smaller than 0.15 g kg^{−1}), its accretion by the melting particles or raindrops may be neglected in our simulations: the contribution of collectional growth to the precipitation rate increase is negligibly small and the effect on the melting process of the sensible heat transferred from accreted liquid droplets is also not significant (in the upper part of the melting region, the quantity of transferred heat is extremely small due to the temperature of the droplets close to 0°C; in the lower part, only the larger melting particles, with the very great mass compared to that of accreted cloud water, are present).

The mass transfer via the condensation process on melting snow and the cooling of air due to the heat removed by conduction from the surrounding air by the melting particles population are calculated from the melting rate by the following equations:

obtained by integrating (19), and (10) and (12). The melting rate is calculated from (26) where the values of *δQ* and 〈*f*〉 have been obtained with the aid of (13) and (25).

The mass transferred from the melting snow to rain during one time step is obtained as the difference between the value of *M*_{m} with the cutoff diameter *d*_{w}, and that calculated from (23) with *d*_{w} increased by Δ*d*_{w} found from (27).

All microphysical transformations are calculated after computation of the prognostic dynamic and thermodynamic equations. Microphysical feedback on dynamic/thermodynamic fields takes place through phase-change-induced heat transfer between water reservoirs and the surrounding air. Hydrometeor loading and air humidity are taken into account in the buoyancy term.

### b. Domain of computation and spatial boundary conditions

The model has been used in its two-dimensional version with the basic domain covering 11 km in the horizontal and 4 km in the vertical, and the bottom level at ground. Grid spacing was chosen to be 30 and 45 m in the vertical and horizontal, respectively. This corresponds to the radar data resolutions; these data were used to initialize the snow in the model and to compare with the model results.

A staggered grid is employed in the model: vertical velocities are calculated at the center of the top and bottom faces of the grid box, horizontal velocities at the center of the sides of the box, and the other variables, such as pressure and temperature deviations as well as all the microphysical fields, are computed at the center of the box.

Subgrid turbulence is included in the model using a first-order closure scheme based on the parameterization proposed by Klemp and Wilhelmson (1978).

At the bottom and top boundaries, the vertical motion is assumed to vanish with the thermodynamic conditions remaining constant throughout the numerical integration. An absorbing layer is used at the top of the domain to avoid artificial modification of dynamics. The 1-km zones at both lateral sides of the domain are not taken into account since the used periodic horizontal boundaries may have an impact on the results in this region. Several experiments have been performed in order to test model sensitivity to changes in the size of domain (doubled horizontal extension of the domain; lifted the upper boundary of the model) and spatial/temporal resolution, with negligible differences in the results.

The model is integrated with a time step of 2 s. Since the snow input in the upper layer is maintained constant in time, a quasi-equilibrium situation with the slow evolution of the model-generated fields develops at about the time that the precipitation reaches the ground (after ∼20 min of simulated time). The model results shown here are obtained after 150 additional time steps. The stationary snow distribution at the top of the model, together with the cyclic horizontal boundary conditions, makes the much longer simulations not relevant. Such lateral conditions imply regularity in the evolution of the different fields in each periodic domain that is not characteristic of atmospheric melting layers. However, with the generated structures much smaller compared to the domain dimension and during the limited time frame of the numerical experiments, the effects of the boundary conditions are negligible (as demonstrated in sensitivity studies).

### c. Initial conditions

Initially, horizontal homogeneity is assumed for all fields, except precipitating snow. The atmosphere is considered at rest. The environmental temperature is initialized using a pseudoadiabatic lapse rate profile with a ground temperature of 13.5°C and with the melting level located at 2.4 km. Pressure is assumed to be initially hydrostatic. The humidity profile corresponds to ice saturation in the presence of snow precipitation in the upper part of the domain. Below this layer, the assumed initial saturation is 90% with respect to water or ice, depending on temperature. The cloud water mixing ratio is set equal to zero throughout the domain.

In order to have realistic space variability of snow, a radar-observed snowfall rate field is used to initiate the microphysical evolution. It is obtained from a vertically pointing X-band radar, in which the timescale was transformed into a space scale using storm velocity. (Appendix C summarizes the main characteristics of the used X-band radar.) The reflectivity is converted to a snowfall rate using a standard *S–Z* relationship. The radar-observed snow is introduced at each time step, in the top layer of the model down to 150 m above the zero isotherm level. These data were complemented by UHF profiler data.

Two 8-min episodes of the bright band associated with stratiform melting conditions observed on 4 June 1996 will be studied. The first numerical experiment (C I) simulates the conditions of light precipitation that occurred between 1920 and 1928, with the reflectivity of rainfall below the bright band of 25 dB*Z* (∼1 mm h^{−1}); the second (C II) refers to a moderate precipitation of 35 dB*Z* (∼5 mm h^{−1}) observed about 20 min earlier, from 1856 to 1904. The presented analysis of the model outputs described in this paper is supported by the results obtained in the second experiment. However, in FS, both cases are presented.

## 4. Description of the model results

### a. Comparison between radar observations and model outputs

The best verification of a numerical experiment is a comparison of the predicted to the observed fields. However, the detailed description of the kinematic and microphysical properties as given by the model is beyond our observational capabilities. We verify the radar measurements of reflectivity factor, the reflectivity-weighted velocity observed in and below the bright band during the event, and the average reflectivity profiles described in Fabry and Zawadzki (1995).

Reflectivity computations making use of the backscattering model for melting snow are described in detail in FS. These electromagnetic scattering computations are made for the hydrometeors in each of the domain’s grid points to derive the following X-band and UHF radar-observable variables: the reflectivity factor, attenuation (total scattering and extinction), and reflectivity-weighted velocity taking into account the simulated vertical air velocity. The quantities describing the hydrometeor population required for electromagnetic calculations (the minimal diameter and melted fraction of the melting snowflakes and the slope of the size distribution for each category) are obtained from dynamic/microphysic model outputs.

The calculated X-band reflectivity pattern using model 5, described in FS, and the snow density given by *ρ*_{s}(*D*_{s}) = 0.015*D*^{−1.0}_{s}, where *D*_{s} and *ρ*_{s} are given in centimeters and grams per cubic centimeters for the second simulated case is compared to the observed reflectivity fields in Fig. 4.

The important horizontal gradient of the simulated reflectivity located 4 km from the left lateral boundary is the most evident difference with the observed pattern. The simulated rapid decrease of reflectivity is related to the strong horizontal gradient of snow mass above the melting layer (as can be seen in Fig. 7a). This gradient is reinforced in the immediate vicinity of the 0°C isotherm by the horizontal advection associated with the adjacent areas of vigorous downdraft and updraft predicted by the model in this region (see Fig. 5), the effect of which may be overestimated by the limitation of 2D modeling. The smoothing inherent in the model output is responsible for the lack of small-scale variability present in the radar data. Otherwise, the resemblance between the two reflectivity fields is rather good.

### b. Induced dynamical structure and related horizontal variability of various features within the melting layer

Figure 5 shows the model-simulated near-equilibrium velocity fields driven by the melting of precipitating snow. The model-generated dynamical structure shows high-frequency oscillations of the vertical wind (and consistent with it, the horizontal wind). The stronger downdrafts are associated with the larger snow contents that are also linked to the strongest convergence–divergence zones associated with the convection. The velocity pattern follows closely, and reinforces, the very small inhomogeneity in the snow precipitation content (as seen in Fig. 7). The peaks of the vertical velocity are associated with regions of the most important horizontal gradients of the snow content. In general, the updraft cores are located slightly higher than the downdraft cores. The level of the domain-averaged position of downdraft maxima is slightly beneath the zone of the maximal cooling (see Fig. 6e). The typical magnitude of the simulated perturbation of the vertical velocity is a few tens of centimeters per second in this case. The horizontal velocity field shows a strong shear associated with the melting layer.

Horizontal variability of some of the microphysical features is shown in Fig. 6. It can be noted that maxima of the condensation rate (Fig. 6b) correspond to minima of the melting rate (Fig. 6a). Contrary to the melting rate—mainly related to the mass concentration—the condensation rate is determined by the saturation level (driven by the vertical motion). The combination of the two processes, one heat releasing and the other heat absorbing, reinforces the horizontal structure of the resulting cooling rate as shown in Fig. 6c.

The comparison of the net cooling rate (Fig. 6e) with the cooling rate caused by melting and the rate of temperature change related to the diffusion process on raindrops (Fig. 6d), points out the important role of the adiabatic heating/cooling in decreasing the difference in the cooling rate of ascending and descending air parcels. The relatively smaller spatial variability of the net rate of cooling assures that the dynamical, thermodynamic, and microphysical structure predicted by the model changes very slowly in time.

The horizontal variability of some simulated microphysical characteristics within the melting region at a level located between the peaks of the horizontally averaged melting rate and reflectivity factor are presented in Fig. 7, together with snow content above the melting layer and rain content below this layer.

### c. Vertical profiles of various microphysical characteristics through the melting layer

Figure 8 presents the horizontally averaged profiles of various microphysical characteristics within the zone of melting as predicted by the model for the simulated case. The profiles of the simulated radar reflectivity factor and the temperature lapse rate are also shown. In order to give an idea of the variability of the presented characteristics, the profiles at the extreme values of the brightband thickness are traced on all the graphs. Each of these curves represent the horizontal averages over three adjacent grid points (i.e., over a distance of 135 m that corresponds to a 6-s-sequence of radar observation).

Figure 8a shows the computed profile of the average melting rate. As it can be seen in (26), this rate is given by the 1.7th moment of the particle size distribution; thus, a rather significant melting rate can be expected in the temperature layer not very far from 0°C as a consequence of the melting of the important mass of ice contained in the smaller, but more numerous, particles. In contrast, the term *δQ* that controls the heat available for melting increases with the distance below the level of 0°C. The resulting profile shows that the melting rate reaches its maximum at a level rather close to the 0°C isotherm. Below the brightband peak, the decrease of melting precipitation mass causes a decrease of the melted mass rate, even if the term *δQ* becomes more important due to the large value of *δT.*

In Fig. 8b the profile of *d*_{w} is shown. The very slow increase of *d*_{w} just below the 0°C isotherm is associated with the very small lapse rate in this layer. It can be noted that *d*_{w} is larger in locations where the vertical extent of the melting zone is smaller, that is, in regions of ascending air. This can be explained by the greater amount of latent heat of condensation available there for melting. At the level at which *d*_{w} reaches the median volume diameter *D*_{o}, the mass of the completely melted snowflakes is half of the initial mass of snow. The averaged value of *D*_{o} is obtained from the rain content below the melting layer. In the studied cases, the levels corresponding to *d*_{w} = *D*_{o} are located below the level of the maximal melting rate, and slightly above the reflectivity peak (Fig. 8e).

Besides *d*_{w}, the mass-weighted fraction of the melted water, 〈*f*〉, gives a bulk description of the stage of the melting process [the profiles of 〈*f*〉 and of *d*_{w} are related by (25)]. Figure 8c shows the vertical profile of the horizontally averaged fraction 〈*f*〉. The pronounced increase of 〈*f*〉 near the onset of melting can be noted. However, when only the larger and fast falling melting particles remain (with 〈*f*〉 between 0.6 and 0.7), the rate of change with height of 〈*f*〉 becomes slower.

The changes with height of the water contents of melting snow and of rain (*M*_{m} and *M*_{r}, respectively) are shown in Fig. 8d. The profiles show large variations in the mass concentration across the melting zone, that is, in the transition zone between low fall-velocity snowflakes and high fall-velocity raindrops. The radar reflectivity profile calculated for the simulated populations of rain and melting snow particles can be seen in Fig. 8e.

The global effect of heat exchange on ambient air temperature during melting and diffusional growth processes is shown in Fig. 8f where the initial and modified temperature profiles are plotted. The major feature of the modified temperature profile is the formation of a nearly isothermal layer with a very small lapse rate starting in the vicinity of 0°C down to the level of the maximal melting rate. It shows that the melting-induced cooling plays a determining role in establishing the final temperature profile. The reflectivity peak is located slightly below the bottom of the quasi-isothermal layer. These results are consistent with the field observations of Willis and Heymsfield (1989) who noted a distance of approximately 100 m between the levels of maximum reflectivity and the base of the isothermal layer. The relatively small thickness of the predicted layer, with a weak vertical temperature gradient, can be explained by a short time of simulation (20 min) starting with a pseudoadiabatic temperature profile. Probably, for the same reason, the level of maximum reflectivity corresponds to the temperature of 1°–1.5°C. Longer periods of melting will deepen the isothermal layer [the observations of Stewart et al. (1984) give a value of 2°C at peak reflectivity in the presence of rainfall of 1.5–2.0 mm h^{−1}]. In the presented case, melting is nearly complete at 4°C.

It is obvious that the distribution of cooling throughout the melting layer is not uniform. The vertical distribution of the heating/cooling rates associated with different diabatic processes can be examined in Fig. 9, where the horizontally averaged profiles of the rate of temperature change due to melting, condensational growth, and rain evaporation are shown. The net predicted cooling rate is also presented in the Fig. 9. Its peak value is in general agreement with the estimations of the cooling associated with the stratiform precipitation region within a midlatitude squall line of Braun and Houze (1995).

The difference between the profile of the cooling due to melting and condensation, and the profile of the net cooling (Fig. 9) is mainly caused by the adiabatic heating/cooling induced by the velocity field. It can be concluded that the dynamical structure within the melting zone may play an important role in the resulting modification of the temperature profile.

Without any background descending motion of air in the domain, the simulated evaporative cooling by rain is relatively small. However, it may become more important below the level of the peak in the cooling-by-melting rate, in a more subsaturated environment.

## 5. Discussion

The bright band is a result of many processes with several feedback mechanisms leading to the generation of finescale circulation. To obtain a complete description of the small-scale structures, and an accurate representation of their role, a numerical model has been developed that includes more complex interactions compared to previous studies. The melting process is described here by a new parameterization scheme suitable for melting layer modeling in stratiform conditions.

Some interesting analytical conclusions were drawn from our parameterization equations. Specially, the role of the exchanged diffusional latent heat has been pointed out. Thus the sign of the diffusional processes depends on the state of saturation when the following conditions are met.

when

*ρ*_{υ}≥*ρ*_{υs}(*T*) condensation takes place on all hydrometeors. At saturation, the contribution to melting of the heat released during the condensation on wet snow is equal to that of the sensible heat extracted from the ambient air.when

*ρ*_{υs}(*T*) >*ρ*_{υ}>*ρ*_{υs}(0°C) condensation takes place on wet snow, with latent heat release, while rain and cloud evaporate consuming some of the heat.when

*ρ*_{υ}<*ρ*_{υs}(0°C) all hydrometeors evaporate and heat is extracted by the process. All the energy for melting is taken from the air sensible heat.

Even if the heat exchange during the difussional processes is important, the mass transferred between vapor and hydrometeors during their fusion is very small.

Besides the importance of the diffusional heat exchange, there are other feedbacks operating within the melting layer whose impact on the local characteristics of the microphysical/dynamical processes is analyzed. All the modeled feedback loops are shown schematically in Fig. 10 except for horizontal and vertical advection (the pressure perturbations are implied in the buoyancy term) at any grid point within the snow/rain transition region. We now attempt a qualitative summary of the feedbacks involved in melting-induced small-scale circulation. As has been shown, inhomogeneities in the horizontal distribution of snow above the melting layer cause horizontal nonuniformity of the melting rate and, hence, induce variability in negative temperature perturbation; this is followed by the development of horizontal gradients of buoyancy. The gradients are strengthened by the loading of precipitation hydrometeors, but with a smaller contribution than that of melting-induced cooling (as has been also shown by Moore and Stewart 1985). The resulting contrasts in buoyancy initiate vertical air motions: downdrafts associated with the larger precipitation content and updrafts with the smaller mass flux becoming positively buoyant with respect to surrounding parcels of air.

Since the relative humidity tends to decrease (increase) in ragions of descent (ascent), the initiated dynamical structure may be invigorated by the reinforcement of the buoyancy fluctuations caused by the two following mechanisms directly related to the local degree of saturation. First, the latent heat associated with vapor diffusional processes causes more intense air cooling in the subsaturated downdrafts, compared to the regions of ascending air. Second, supersaturation can be attained in the updraft branches of circulation, leading to the generation of the cloud liquid water associated with the release of latent heat and consequent local warming. The very localized production of cloud liquid water is predicted by the model to be in the most intense updrafts only.

In order to prevent an excessive buoyancy contrast, opposing mechanisms develop that help to maintain the circulation in a nearly steady state. These are mainly the adiabatic heating during descent and the cooling during ascent. The vertical advection transporting the warmer air upward in the regions of smaller cooling-by-melting followed by the horizontal advection constitutes an additional action modulating horizontal gradients and reducing the intensity of convection.

The two-dimensional aspect of the model used here, in which the water content gradients and wind convergence along the second horizontal direction are neglected, certainly affects the results. This may tend to increase the intensity of the convective motion if adding a component in the third dimension reinforces the gradient of melting snow content. On the other hand, we suspect that, under conditions in which the downdraft and updraft motions are not constrained to form in the same plane and thus the air is allowed to spread in all directions, the intensity of convection may be reduced. However, due to inherent nonlinearity, three-dimensional experiments are needed to simulate the airflow more precisely and to evaluate the effects of the constraints of two-dimensionality.

The melting-induced small-scale dynamical perturbations, controlled by the horizontal precipitation inhomogeneities, introduce only small variability around the averaged values of the different microphysical fields. Particularly, the calculated radar reflectivity profile is only slightly sensitive to the dynamical structures; the only critical parameter for the finestructure of the radar reflectivity seems to be the precipitation rate. Yet, this finescale dynamics may alter slightly the local cloud and precipitation fields and influence larger-scale circulation by exercing internal drag, as has been mentioned by Szeto et al. (1988a).

Given the great variability of snow parameters (shape, density, fall speed, etc.) our microphysics parameterization should be considered as representing some average behavior and the model-predicted melting layer characteristics are expected values ensuing from a highly variable process. However, as far as the comparison of the calculated fields to the radar observations is concerned, the greatest uncertainty lies in the assumed morphology of the melting particles and subsequent electromagnetic modeling as described in the companion paper by FS. We believe that it is no coincidence that the best agreement of the radar-observed variables with the corresponding model-derived variables corresponds to electromagnetic computations using snowflake morphology consistent with the one assumed in our microphysical parameterization. The agreement between model and observations supports both models.

## Acknowledgments

This work was partially supported by a contract with the Canadian Communication Research Centre. We would like to thank Aldo Bellon for a thorough reading of the manuscript. We would like to express our thanks also to the anonymous reviewers for helpful comments and critical questions.

## REFERENCES

### APPENDIX A

#### List of Symbols

### APPENDIX B

#### Expressions Used in Derivation of Eq. (15)

##### Expressions

Clausius–Clapeyron equation

Empirical relation for *e*_{s}

with coefficients:

*b*_{0} = 6.110 455

*b*_{1} = 0.444 235 1

*b*_{2} = 1.430 209 9 × 10^{−2}

*b*_{3} = 2.645 470 8 × 10^{−4}

*b*_{4} = 3.035 709 8 × 10^{−6}

*b*_{5} = 2.097 226 8 × 10^{−8}

*b*_{6} = 6.048 759 4 × 10^{−11}

*b*_{7} = −1.469 687 × 10^{−13}

The error associated with the relation for *e*_{s} is lower than 0.02% for the temperature above 0°C.

### APPENDIX C

#### Characteristics of the McGill X-Band VPR Radar Used in This Work

Wavelength 3.2 cm

Peak power 25 kW

Beamwidth 2°

Pulse length used 0.25 *μ*s

Height resolution 30 m

Time resolution 2 s

Doppler capability None

## Footnotes

*Corresponding author address:* Dr. Isztar Zawadzki, Dept. of Atmospheric and Oceanic Sciences, McGill University, 805 Sherbrooke St. West, Montréal, PQ H3A 2K6, Canada.