Diabatic cooling from hydrometeor phase changes in the stratiform melting layer is of great interest to both operational forecasters and modelers for its societal and dynamical consequences. Attempts to estimate the melting-layer cooling rate typically rely on either the budgeting of hydrometeor content estimated from reflectivity Z or model-generated lookup tables scaled by the magnitude of Z in the bright band. Recent advances have been made in developing methods to observe the unique polarimetric characteristics of melting snow and the additional microphysical information they may contain. However, to date no work has looked at the thermodynamic information available from the polarimetric radar brightband signature. In this study, a one-dimensional spectral bin model of melting snow and a coupled polarimetric operator are used to study the relation between the polarimetric radar bright band and the melting-layer cooling rate. Simulations using a fixed particle size distribution (PSD) and variable environmental conditions show that the height and thickness of the bright band and the maximum brightband Z and specific differential phase shift are all sensitive to the ambient environment, while the differential reflectivity is relatively insensitive. Additional simulations of 2700 PSDs based on in situ observations above the melting layer indicate that the maximum Z, , and within the melting layer are poorly correlated with the maximum cooling rate while is strongly correlated. Finally, model simulations suggest that, in addition to riming, concurrent changes in aggregation and precipitation intensity and the associated cooling may plausibly cause observed sagging brightband signatures.
Hydrometeors undergoing a phase transition exchange latent heat with the environment and are a fundamental driver of atmospheric motion across a range of spatial and temporal scales. In stratiform precipitation, the exchange of latent heat as snow and graupel melt cools the environment. This melting-induced cooling may deepen isothermal layers (e.g., Stewart et al. 1984), which can suddenly alter the surface precipitation type and the attendant societal impacts (e.g., Wexler 1955; Bosart and Sanders 1991; Kain et al. 2000). It can also have dynamical consequences and generate downdrafts and gravity waves (e.g., Szeto et al. 1988), turbulence and convection (e.g., Findeisen 1940; Stewart et al. 1984), mesoscale wind perturbations (e.g., Atlas et al. 1969; Heymsfield 1979), and mesoscale fronts (e.g., Heffernan and Marwitz 1996). These impacts can amplify over time and even alter downstream cyclogenesis (e.g., Stewart and Macpherson 1989). While the magnitude of cooling in stratiform precipitation is typically smaller than in deep moist convection, it occurs within a relatively narrow layer and over a much larger spatial area.
Owing to these pronounced impacts, possessing information about the cooling rate within the melting layer is highly desirable. The lack of in situ thermodynamic observations in the melting layer has led to the development of a number of methods to retrieve thermodynamic information from remote sensing data. Of these data sources, weather radars are a natural source of data to exploit as their measured signal is composed of backscattering from the hydrometeors responsible for the cooling and features a pronounced signature associated with the melting layer at weather radar wavelengths. This signature, first identified by a layer of enhanced equivalent radar reflectivity factor Z called the “bright line” by its founders (now known as the “bright band”; e.g., Byers and Coons 1947), was correctly understood to be associated with the melting layer soon after its discovery (e.g., Cunningham 1947; Austin and Bemis 1950). To a first-order approximation, the enhancement in Z is due to the rapid increase in the relative permittivity of melting snowflakes before they collapse into smaller raindrops and decrease their concentration through flux divergence. Leary and Houze (1979) estimated hydrometeor content above and below the melting layer from Z and, through a budgeting approach, obtained an estimate of the melting-layer cooling rate. However, this approach relies on estimating hydrometeor content from Z using a fixed relationship, which may be prone to appreciable and unknown errors (e.g., Li and Srivastava 2001; Carlin et al. 2016). More recently, lookup tables of vertical profiles of latent heating and cooling generated from convection-allowing model ensembles have become more prevalent (e.g., Shige et al. 2004, 2007; Tao et al. 2006). However, Shige et al. (2007) note the sensitivity of the lookup tables to the simulations used to derive them, and the use of single-moment bulk microphysics may preclude the proper modeling of the melting process and its associated profile of cooling rates (e.g., Phillips et al. 2007), resulting in a monotonic relation between melting layer Z and the magnitude of melting-layer cooling that is overly simplistic.
Despite these uncertainties, the radar bright band also exhibits distinct polarimetric characteristics that nevertheless suggest it is a valuable source of microphysical information. As with Z, the polarimetric characteristics of the bright band were also observed relatively early on, with observations revealing local linear depolarization ratio (Browne and Robinson 1952; Humphries 1974; Anderson 1974), circular depolarization ratio (Humphries and Barge 1979), and differential reflectivity (Bringi et al. 1981) maxima below the height of the Z maximum noted in subsequent decades. Illingworth and Caylor (1989) and Zrnić et al. (1993) noted a pronounced reduction in the copolar correlation coefficient in the melting layer, attributed both to hydrometeor diversity (e.g., Balakrishnan and Zrnić 1990) and backscatter differential phase δ caused by large wet aggregates (e.g., Zrnić et al. 1993; Trömel et al. 2014). This δ can result in an oscillatory profile of specific differential phase shift through the melting layer, which is otherwise positive with a maximum near the height of the Z maximum (Zrnić et al. 1993). These conspicuous polarimetric signatures play a primary role in most modern melting-layer detection algorithms (e.g., Giangrande et al. 2008; Wolfensberger et al. 2016), but have yet to be used to make inferences about melting-layer thermodynamics.
Quasi-vertical profiles (QVPs; Ryzhkov et al. 2016) have emerged as a convenient way to study the evolution of polarimetric radar signatures. QVPs display a time series of azimuthally averaged radar data from a high-elevation scan in a time–height format. This averaging greatly increases the number of samples, which reduces errors and improves data quality compared to typical plan position indicator scans and allows for the detection of subtle polarimetric signals that may have otherwise been obscured by noise (e.g., Griffin et al. 2018). However, the averaging procedure obscures heterogeneities within the averaging domain and works best in more homogeneous precipitation, such as stratiform rain or snow. QVPs have been used to study the characteristics of the melting layer (e.g., Trömel et al. 2017) and its transient features, such as so-called sagging bright bands (e.g., Kumjian et al. 2016).
Models are frequently used to gain insight into the melting layer and its impacts (e.g., Tao et al. 1995). Owing to the deficiencies of bulk schemes in representing the melting process, these models often employ spectral bin microphysics to explicitly simulate mixed-phase particles and variable terminal velocity and water fraction across the particle size distribution (PSD). Early modeling studies focused primarily on the representation of the physical processes of a melting snowflake (e.g., Matsuo and Sasyo 1981; Mitra et al. 1990). Because of the intimate connection between the melting layer and the radar bright band, however, many studies include electromagnetic scattering components of varying complexity to compare profiles of observed Z in bright bands to modeled Z in one-dimensional (e.g., Yokoyama and Tanaka 1984; Klaassen 1988; Hardaker et al. 1995; Fabry and Szyrmer 1999; Gray et al. 2001; Olson et al. 2001; Zawadzki et al. 2005) and three-dimensional (e.g., Phillips et al. 2007; Planche et al. 2014; Iguchi et al. 2014) models. A lesser number of studies have simulated polarimetric variables in the bright band (e.g., Russchenberg and Ligthart 1996; D’Amico et al. 1998; Giangrande 2007; Trömel et al. 2014). These studies have sought to better understand the physical processes responsible for the radar bright band and evaluate the performance of microphysics schemes by their ability to reproduce realistic brightband signatures, but none have investigated the thermodynamic information of the polarimetric brightband signature.
In this study, a one-dimensional spectral bin model of melting snow and a coupled polarimetric operator are used to investigate the relationship between the polarimetric characteristics of the melting layer and the latent cooling rate within it. A description of the model is provided in section 2. The sensitivity of the bright band to the environment is examined in section 3, while section 4 delves into the potential for polarimetric thermodynamic retrievals in the melting layer. Finally, section 5 uses the model to investigate the potential causes of sagging bright bands and is followed by a brief summary and discussion of future work in section 6.
2. Description of one-dimensional melting snow model
a. Microphysical model
The one-dimensional melting snow model (1D-MS) is a one-dimensional Lagrangian spectral bin model of melting snow aggregates similar to those used in previous studies (e.g., Szyrmer and Zawadzki 1999; Zawadzki et al. 2005; Giangrande 2007; Grim et al. 2009; Trömel et al. 2014) with some modifications. (A full list of variables and their units is included in Table A1 in appendix A.) PSDs containing up to 80 size bins are tracked downward from the model top as particles sublimate, melt, and evaporate. Interactions between bins (e.g., aggregation) are excluded, with one snowflake corresponding to one raindrop and the particle diameter of each bin changing with height (i.e., in the absence of evaporation/sublimation, the bins are mass conserving). The density of snow varies with the riming factor , which varies from 1 for unrimed snow to 5 for heavily rimed snow (Zawadzki et al. 2005), and inversely with equivolume snowflake diameter according to (Brandes et al. 2007)
and is capped at a maximum value of 0.5 g cm−3 (Zawadzki et al. 2005). Initial bin masses at the top of the model are equal to those of raindrops ranging from 0.05 to 7.95 mm in diameter in 0.1-mm increments, with the corresponding found according to
An example of how particle diameters change with height is shown in Fig. 1a. Through the assumption of flux conservation, the PSD at any height z is found according to
where D is the equivolume diameter of the particle at any stage of melting and () is the particle (dry snow) terminal velocity. While any PSD form can be specified, gamma distributions of the form
where , , and , the intercept, shape, and slope parameters for snow, respectively, are used. Following Szyrmer and Zawadzki (1999), is a function of meltwater fraction [defined as ()−1, where () is the mass of water (ice) in the particle] and the terminal velocity of a raindrop of equivalent mass is given by Brandes et al. (2002) with a correction for local air density (Foote and du Toit 1969):
where and .
The transfer of heat by radiation, curvature and solute effects, and the collection of cloud droplets are neglected as they are small in comparison to the primary microphysical processes (Szyrmer and Zawadzki 1999). The particle’s temperature is assumed to be homogeneous, with no sensible heat stored in the particle during melting (i.e., all available heat goes toward the phase change). The local environment and particle composition and temperature determine whether sublimation, melting, or evaporation occurs.
If = 0.0 and the environment is subsaturated with respect to ice, sublimation occurs and cools the particle. Following Pruppacher and Klett (1997), the equilibrium particle temperature is found according to
where is the difference in vapor density between the environment and the particle surface with respect to ice at temperature , given by
The thermodynamic terms in these and subsequent equations are further defined in appendix B. Equation (6) is solved iteratively to within a 0.01°C threshold. Because is relatively slow compared to the model grid spacing (dz = 10 m), is assumed to be equal to for dry snow. Following Rogers and Yau (1989), the rate of change of ice mass due to sublimation is found according to
The total exchange of mass by a particle within a model layer is found by considering the particle residence time, given by where w is the local vertical velocity.
Cooling due to sublimation can result in 0°C when the particle reaches the 0°C level, which can delay the onset of melting by up to a few hundred meters (Matsuo and Sasyo 1981). Melting commences once reaches 0°C, which is assumed to remain constant for the duration of the melting process. The rate of ice loss due to melting is equal to
where the heat available for melting is given by
and where the vapor density deficit is with respect to water. If , cooling due to evaporation consumes more heat than is available from conduction, resulting in no melting and the consideration of sublimation instead. The decrease in ice mass due to melting is converted to meltwater (i.e., ). Refreezing is not considered.
Evaporation is treated in an analogous manner to sublimation and occurs if in an environment subsaturated with respect to water, with the rate of change of water mass due to evaporation given by
However, unlike snow, the of raindrops cannot be assumed to be equal to because of the larger with respect to dz and the thermal relaxation time for raindrops (about a few seconds for large drops; e.g., Caplan 1966; Tardif and Rasmussen 2010). Once melting is complete, the incremental change in for raindrops within a model layer is given by (Tardif and Rasmussen 2010)
with capped at [found using Eq. (6) but replacing and with and , respectively]. An example comparison of and for a variety of raindrop sizes is shown in Fig. 1b. Sublimation delays melting by 130 m. Once melting begins, the different sized particles remain at 0°C and take different distances to completely melt, ranging from 190 m for a particle with a ≈ 0.5 mm to 660 m for a particle with ≈ 4.0 mm. Because of variations in mass and the local environmental temperature where melting ceases, different size raindrops warm at different rates; smaller raindrops reach quickly once melting is complete, while the largest drops never quite reach their .
The density of the snow core is assumed to increase during melting. By modifying Eq. (1) for volume instead of diameter, the volume of snow can be expressed as
where on the right-hand side of Eq. (13) comes from the preceding model layer. From this, the new snow core density is found according to
Once reaches , it is kept constant regardless of .
The rate of change of the environmental temperature at each height is found following Grim et al. (2009):
Similarly, the rate of change of the environmental water vapor mixing ratio is found according to
where is the density of dry air at height z. Steady-state conditions are assumed to exist for 60 s at a time.
b. Electromagnetic scattering model
The values of Z, , and are calculated from the 1D-MS output using the Ryzhkov et al. (2011) polarimetric operator. This operator treats particles as homogeneous oblate spheroids. Although more detailed microphysical and electromagnetic scattering models of melting aggregates are beginning to be explored (e.g., Botta et al. 2010; Johnson et al. 2016; Leinonen and von Lerber 2018), oblate spheroids greatly simplify the treatment of these processes and have been shown to be an acceptable approximation for hydrometeors at longer weather radar wavelengths (e.g., X, C, and S bands; Matrosov et al. 1996; Hogan et al. 2012) and reproduce observed polarimetric signatures fairly well (e.g., Ryzhkov et al. 2011). Scattering amplitudes are calculated using PyTMatrix (Leinonen 2014), a Python package containing an implementation of the Mishchenko (2000) T-matrix code. The relative permittivity of ice and air are considered constant while the relative permittivity of water varies with temperature following Ray (1972). Past studies have used a wide variety of approaches to calculate the effective relative permittivity of melting snowflakes (e.g., Yokoyama and Tanaka 1984; Hardaker et al. 1995; Klaassen 1988; Meneghini and Liao 1996; Russchenberg and Ligthart 1996; D’Amico et al. 1998; Fabry and Szyrmer 1999; Battaglia et al. 2003), which remains a source of uncertainty in radar-modeling studies of melting snow. The widely used Maxwell-Garnett mixing formula (Maxwell Garnett 1904), which treats mixtures as randomly distributed spherical inclusions within a matrix, provides different results depending on the topology of the mixture. For a three-phase mixture, there are 12 unique combinations of inclusions and matrices. However, not all of these mixtures are equally plausible. For example, it is known that the water in melting snow aggregates is initially sucked inward toward the vertices of the ice structure (e.g., Knight 1979; Mitra et al. 1990). As a result, mixing formulas that assume water as the dominant medium typically result in unrealistically large brightband signatures that conflict with available observations (e.g., Fabry and Zawadzki 1995). Sensitivity tests with the 1D-MS (not shown) confirm these biases. Alternative approaches, such as the Polder–van Santen mixing formula (Polder and van Santen 1946), exist that attempt to minimize biasing too heavily toward that of the assumed dominant media. Many studies (e.g., Meneghini and Liao 1996; Matrosov 2008; Ryzhkov et al. 2011) suggest using a more complex weighted Maxwell-Garnett mixing formula approach that varies between a matrix of dry snow (itself an air matrix with ice inclusions) with water inclusions () and a matrix of water with dry snow inclusions () based on the volumetric water fraction (, defined as ). Here we follow an approach similar to Ryzhkov et al. (2011) in which
to simulate the initial dry snow shell followed by a quick transition to water as the dominant dielectric medium.
The particle axis ratio varies linearly with between that of raindrop of equivalent mass (given by Brandes et al. 2005) and that of dry snow aggregates , which is typically assumed to be constant across the size spectrum (e.g., Vivekanandan et al. 1994; Ryzhkov et al. 2011; Thompson et al. 2014). Observational evidence suggests = 0.6 is typical for unrimed aggregates (Korolev and Isaac 2003; Matrosov et al. 2005; Hogan et al. 2012; Garrett et al. 2015). Based on observations of increasingly spherical particles with riming (Garrett et al. 2015), is derived as a linear function of according to
Particles are assumed to have a mean canting angle of 0° (i.e., no preferred orientation), with a standard deviation that varies linearly with during melting between that of rain (10°; Ryzhkov et al. 2002; Huang et al. 2008) and dry snow (40°; Hendry et al. 1987; Garrett et al. 2015).
Figure 2 shows output from the 1D-MS compared to a range-defined QVP (Tobin and Kumjian 2017) from the Vance Air Force Base, Oklahoma (KVNX), WSR-88D radar in the stratiform portion of a mature mesoscale convective system (MCS) on 20 May 2011. The 1D-MS was initialized using temperature, humidity, and particle data collected during the Midlatitude Continental Convective Clouds Experiment (MC3E; Petersen and Jensen 2012) by the University of North Dakota (UND) Cessna Citation II in north-central Oklahoma for the same event (Wu and McFarquhar 2016; Xue et al. 2017). The temperature and humidity profiles were taken from the aircraft’s measurements during ascent and verified against the 1200 UTC sounding from the nearby Atmospheric Radiation Measurement (ARM; Stokes and Schwartz 1994) Southern Great Plains facility. Joint mean observed PSDs were constructed for 10-s intervals from the two-dimensional cloud probe (2D-C; Knollenberg 1981) below D = 1.0 mm and from the high-volume precipitation spectrometer (HVPS; Lawson et al. 1993) at and above D = 1.0 mm and fit to an incomplete gamma distribution following Wu (2017). An adjustment to the particle bins by a factor of (which remained near the assumed value of 0.6) is required to account for the discrepancy between the particle maximum diameter that characterizes the observed PSDs and the equivolume diameter used in the 1D-MS. The 1D-MS was run using each of the fit PSDs from between 1332:24 and 1339:47 UTC (i.e., while the aircraft was between 0° and −2°C) to account for the possible range of PSDs observed above the melting layer, with the maximum particle size found according to the relation proposed in Heymsfield et al. (2013),
where and are in centimeters and per centimeter, respectively, and where a correction factor of is again needed. Snow was assumed to be unrimed. Note that the treatment of all particles as snow aggregates is a simplifying assumption. A subjective review of particle imager data from both the 2D-C and HVPS, along with the results of the Holroyd (1987) habit classification algorithm, suggest that hexagonal plates and irregular particles were dominant at smaller sizes and aggregates were dominant at larger sizes. Accordingly, this assumption appears justifiable, particularly for evaluating Z and , which are dominated by the largest particles. At the time of the aircraft’s flight, the system was already prohibitively far from the radar to create a range-defined QVP with useful vertical resolution. To achieve more accurate and representative vertical profiles of the polarimetric variables, the range-defined QVP profiles are shown for between 1027:47 and 1132:51 UTC as the bulk of the stratiform region passed over KVNX. Although the aircraft sampled the stratiform region 2–3 h later, the MCS was mature and seemingly steady-state during this time period, justifying this comparison. A bias correction of −0.2537 dB was applied to the range-defined QVP using the bias estimation method described in Zittel et al. (2014) and Richardson et al. (2017). This estimate was found by taking the weighted mean of the median values from each of the three bias estimation approaches from KVNX during an 8-day window centered on the event. For both the 1D-MS results and the range-defined QVP profiles, the distribution of profiles is shown in shading with the median profile denoted. To make a direct comparison to observations, the 1D-MS results have been modified to account for attenuation, δ, and beam broadening. Because the range-defined QVP combines data from all elevation angles, it is not straightforward to simulate the impact of beam broadening on the resultant profile. Here, an elevation angle of 8.0° is used as it is an intermediate elevation angle and results in moderate beam-broadening effects.
Overall, there is modest agreement between the observed and simulated profiles, although some notable differences are evident. The peak magnitudes of Z and within the bright band are well captured, along with the magnitude of Z below the melting layer. The general shape of the Z and profiles is also well represented, with realistic rates of increase during melting and a maximum at a lower height than the Z maximum. The 600-m depth of the intrinsic Z bright band (prior to accounting to beam broadening, not shown) is also in line with observations, with Fabry and Zawadzki (1995) suggesting a mean brightband depth of ≈500 m for rain Z between 30 and 35 dBZ, and the peak simulated values of all three variables are in agreement with those observed by Wolfensberger et al. (2016). However, a number of differences between the simulated and observed profiles are present. Some of these differences are methodological: the melting-layer profile in QVPs is artificially broadened by the method used to remove the contribution of δ (Griffin et al. 2018), which involves linearly interpolating the profile of each radial through the melting layer prior to the azimuthal averaging. In addition, the input PSDs used in the 1D-MS already take aggregation into account and represent PSDs immediately above 0°C (≈3.6 km AGL), while the observed profile of Z shows ongoing aggregation toward the melting layer. Other discrepancies, however, likely point to processes currently excluded from the 1D-MS. For example, spontaneous breakup of the largest raindrops could account for the observed decrease of Z and toward the surface, while the exclusion of aggregation within the melting-layer is likely responsible for the underpredicted depth of the bright band and the magnitude of the values within and beneath it. Beyond these considerations, a number of uncertainties remain that complicate this direct comparison, including the representativeness of the assumed environment and initial PSDs (in both space and time), the steady-state assumption for the range-defined QVP, the effect of azimuthal averaging if the melting-layer height and depth vary in space, and the exclusion of dense and very oblate pristine ice crystals from the 1D-MS, which likely explains, at least in part, the discrepancies in above the melting layer. Keeping these various assumptions and sources of uncertainty in mind, these results suggest that the 1D-MS is capable of simulating the melting layer and its general polarimetric characteristics.
3. Impact of environment on brightband characteristics
Before investigating the retrieval of latent cooling rates within the melting layer, it is worthwhile to investigate how the melting layer and bright band are modified by factors other than the PSD. Figure 3 shows the impact of varying the environmental lapse rate and relative humidity (RH) profile on the melting layer and bright band for both unrimed aggregates ( = 1.0) and heavily rimed aggregates ( = 4.0) at S band (λ = 11.0 cm). In situ observations from near the melting layer of a wide range of MCSs are used to constrain the parameter space. These observations suggest lapse rates of 6°C km−1 are typical (e.g., Willis and Heymsfield 1989; McFarquhar et al. 2007; Smith et al. 2009), with McFarquhar et al. (2007) observing ranging from approximately 4° to 8°C km−1. These same field campaigns suggest that RH profiles typically exhibit negative gradients with respect to temperature below the 0°C level, with values ranging from −1% to −5% °C−1 and with gradients from −3% to −4% °C−1 typical. While RH values near the 0°C level are typically quite moist (e.g., 90% or higher; Willis and Heymsfield 1989; McFarquhar et al. 2007; Smith et al. 2009), these values can occasionally drop to as low as 60% (McFarquhar et al. 2007; Heymsfield et al. 2015). Here, lapse rates are varied from 4.0° to 8.0°C km−1 and RH at the 0°C level is varied from 60% to 100% with a RH gradient of −3% °C−1 to capture a wide range of conditions observed near the melting layer in MCSs. For all environments, the sample input PSD was characterized by = 1720 m−3 mm−(1+μs), = 0.34 mm−1, = −1.22, and = 20.0 mm, and the input PSD for heavily rimed snow was that which resulted in the equivalent distribution of rain at the surface. The qualitative conclusions of Fig. 3 were consistent for other PSDs examined (not shown).
Figures 3a and 3b show the distance below the environmental 0°C level at which melting begins. Melting begins at the same height for both = 1.0 and = 4.0 because is assumed equal to for dry snow and thus is solely a function of the environmental T and RH. Consistent with previous findings (e.g., Matsuo and Sasyo 1981; Rasmussen and Pruppacher 1982; Heymsfield et al. 2015), melting is delayed for colder and drier environments, with increased sensitivity to RH in colder environments. For environments near saturation, melting begins within the first 100 m below the environmental 0°C level. For an RH at 0°C of only 60%, however, melting can be delayed by ≈600 m when = 8.0°C km−1 to over 1200 m when = 4.0°C km−1.
The depth of the melting layer as a function of and RH is shown in Figs. 3c and 3d. The melting layer is defined as the layer where melting has begun but ice is still present. For both = 1.0 and = 4.0, the depth of the melting layer increases as decreases because of a smaller sensible heat flux and slower melting. Counterintuitively, the environmental RH has little effect on the melting-layer depth. This limited sensitivity has two causes: 1) drier environments increase the distance below the 0°C level at which melting begins, resulting in melting occurring in warmer layers with larger sensible heat fluxes that offset the impact of the decreased RH, and 2) drier environments increase sublimation above the melting layer, decreasing the amount of ice entering the melting layer and thus the distance needed for complete melting. For a given environment, the depth of the melting layer when = 4.0 is also larger by roughly 50% than when = 1.0 because of faster particle fall speeds. These faster fall speeds result in reduced sensitivity to the local RH at a given height and less mass lost to sublimation above the melting layer.
The environment can also modulate the magnitude of the radar bright band. Figure 3e,f shows the maximum Z within the bright band as a function of the environment. The maximum Z is quite sensitive to the RH because of the evaporation of meltwater and the aforementioned sublimation above the melting layer, which leads to smaller and fewer particles entering the melting layer. In contrast, has relatively less impact on Z because of its sensitivity to the largest particles whose increases appreciably shortly after melting begins. Many of these same sensitivities are seen for maximum (Figs. 3i,j), which decreases in drier environments because of both the evaporation of meltwater and a decrease in ice mass entering the melting layer because of sublimation. In contrast with Z and , the maximum (Figs. 3g,h) is relatively insensitive to the environment. This insensitivity is due to being unaffected by concentration, the assumption of constant and during sublimation (which may not be the case; Sulia and Harrington 2011), and the weighting of toward the largest particles, which experience relatively little sublimational mass loss.
Finally, Figs. 3k and 3l show the sensitivity of the maximum cooling rate (i.e., the minimum ) within the melting layer to the environment. As with the bright band, both the magnitude of the maximum cooling rate and its sensitivity to the environment are decreased for heavily rimed particles because of faster and the greater depth over which melting is distributed (e.g., peaking at −17.5 K h−1 for = 1.0 and −8.7 K h−1 for = 4.0). However, in contrast to the maximum Z and , which are positively correlated with both and RH, the maximum cooling rate generally increases for increasing and decreasing RH. For environments at or near saturation, larger lapse rates increase the heat conducted into particles and the efficiency of melting, resulting in larger cooling rates in a given layer. For drier environments, the impact of the environment depends on the balance between the changing environmental temperature and vapor pressure deficit (itself an interplay between the RH and environmental temperature). For the same subsaturated RH, the conduction term exceeds the evaporation term [term 1 and term 2 in Eq. (10), respectively] more quickly for larger . This results in melting beginning sooner (which decreases the amount of ice lost to sublimation prior to melting) and to larger values and subsequently larger cooling rates. The exception to this pattern is in very dry environments, where sublimational mass loss above the melting layer is so large that maximum cooling rates begin to decrease for all values of . However, even in these very dry environments, the interplay between the conduction and evaporation terms results in larger values for larger and thus larger cooling rates.
4. Relation between brightband characteristics and cooling rate
To investigate the connection between the polarimetric variables and the cooling rate within the melting layer, the 1D-MS was initialized using a range of PSDs based on in situ observations of stratiform precipitation regions (e.g., Heymsfield et al. 2002; McFarquhar et al. 2007; Neumann 2016). These PSDs were characterized by three-parameter gamma distributions and intended to capture a realistic range of PSD variability above the melting layer. However, the three parameters of the gamma distribution are not necessarily independent, with numerous relations between parameters reported in the literature (see Table 2 of Patade et al. 2015). We primarily utilize the parameter relations of McFarquhar et al. (2007) in this work, although the approach of McFarquhar et al. (2015) using ellipsoids characterizing equally realizable combinations of gamma PSD parameters likely better captures the true range of uncertainty of in situ PSDs and may be explored in future work.
The slope parameter ranges from approximately 0.2 to 1.5 mm−1 in the region immediately above the melting layer (e.g., McFarquhar et al. 2007; Heymsfield et al. 2002, 2013, 2015; Patade et al. 2015; Neumann 2016; Matrosov and Heymsfield 2017). The in situ observations used to derive these values of were characterized by particle maximum diameter. Therefore, the values of used to initialize the 1D-MS require the aforementioned adjustment factor of and ranged from approximately 0.25 to 1.8 mm−1.
The slope parameter and intercept parameter (which can span many orders of magnitude) can be constrained by their relation to . McFarquhar et al. (2007) provide such relations, given by
where is per centimeter and is in cm−(4+μs). To not reduce the PSDs to single-moment distributions and retain multiple degrees of freedom in the PSD parameter space, ten samples of and are randomly selected from a Gaussian distribution centered on the mean and [found from Eqs. (21) and (22)] for each value of . The standard deviation of each parameter distribution was subjectively determined to be 0.2 from Figs. 18 and 19 of McFarquhar et al. (2007). The conversion from units of cm−(4+μs) in Eq. (22) to units of m−3 mm−(1+μs) used in the 1D-MS is dependent. Thus, every drawn from the distribution is converted to units of m−3 mm−(1+μs) using all 10 sampled values, resulting in 2700 total PSDs. Comparisons of and retrieved using this sampling approach with those reported in McFarquhar et al. (2007) are shown in Figs. 4 and 5, respectively. There is overall good agreement between the sampled values and those reported in McFarquhar et al. (2007), providing confidence that this sampling approach successfully reproduces the observed variability of the in situ observations. The maximum particle diameter is found according to Eq. (20).
In addition to varying the initial PSDs, both the environmental and RH profiles were varied as the bright band can be sensitive to the environment (see section 3). To incorporate the effects of variable environmental conditions on the thermodynamic retrieval efficacy, simulations using all 2700 PSDs were performed for the range of environments described in section 3 (i.e., of 4.0°–8.0°C km−1 every 1.0°C km−1 and RH gradients from −4.5% to 0.0% °C−1 every 1.5% °C−1, with the RH at 0°C assumed to be 100%) for a total of 54 000 simulations. The same set of input PSDs were used with each environment for consistency. However, additional experiments were performed in which every environment used its own sampling of 2700 distributions that reached similar conclusions (not shown), providing confidence that the sampling procedure is not introducing biases or appreciable sampling error into the analyses. The ranges of PSD parameters, , and RH profiles used are summarized in Table 1. All of the simulations were for unrimed snow ( = 1.0) with scattering calculations performed at S (λ = 11.0 cm), C (λ = 5.45 cm), and X (λ = 3.2 cm) bands.
Figure 6 shows the maximum Z, (defined as the change in Z at the 0°C level to its maximum within the melting layer), , and compared to the maximum cooling rate within the melting layer for all 2700 modeled distributions at S, C, and X bands. For brevity, only the environment characterized by = 6.0°C km−1 and an RH gradient of −3.0% °C−1 is shown; the qualitative conclusions remained the same for all environments. Ordinary least squares regressions were performed in logarithmic space for all variables to combat the heteroscedasticity of the data.
The maximum brightband Z is essentially uncorrelated with the maximum cooling rate in the melting layer, with values near 0.1, maximum cooling rates that span an order of magnitude for a given value of Z, and root-mean-square errors (RMSEs) of ≈6.7 K h−1 at all wavelengths (Figs. 6a–c). RMSEs are larger for warmer and drier environments (not shown). The lack of correlation between Z and the cooling rate results from differences in which particles dominantly contribute to each variable: the bulk of the cooling is due to the smallest particles that are present in high concentrations, while Z in the Rayleigh regime is a strong function of larger particles that are more sparse and contribute minimally to the cooling rate. The effects of aggregation may exacerbate this lack of correlation by further enhancing Z while only minimally affecting the cooling rate.
Similarly, is weakly correlated with the maximum cooling rate, with comparable RMSEs to Z and values that remain below ≈0.2 (Figs. 6d–f). However, unlike with Z, and the maximum cooling rate are inversely related. This seemingly counterintuitive relation results from PSDs with many small particles (and thus large cooling rates) generally having larger and thus smaller , which leads to smaller .
Almost no relation exists between the maximum brightband and the maximum cooling rate in the melting layer (Figs. 6g–i). As with Z, is sensitive to the largest particles in the volume that are not responsible for the bulk of the cooling and may be further enhanced by aggregation. Unlike with Z, however, is insensitive to particle concentration (which is linearly related to the cooling rate). These factors result in high RMSEs and low values, with unable to provide appreciable information about the maximum cooling rate within the melting layer.
In stark contrast, the maximum exhibits a striking linear relation to the maximum cooling rate, with values of 0.96 and an approximate fivefold reduction in RMSEs compared to Z, , and (Figs. 6j–l). For all environments and radar wavelengths, the exponent in the regression equation remains near 0.9, with the coefficient a function of the environment that increases for increasing and RH gradients and radar wavelength. Similar to Z, RMSEs are smaller for moister environments; the RMSE is reduced by 50% for an RH gradient of 0.0% °C−1 compared to that shown here.
The cause of the pronounced correlation between the maximum and maximum cooling rate is examined in Fig. 7, which shows the contribution of each particle size bin to the cooling rate, Z, and by height and the moment M of the PSD each variable is proportional to locally (found by taking the gradient of the logarithm of each variable with respect to the logarithm of D). The cooling rate is dominated by contributions from smaller size bins (peaking around 1.0–1.2 mm) and reaches its maximum shortly after melting begins at approximately 150 m below the 0°C level (Fig. 7a). Additionally, the cooling rate due to melting is proportional to ≈M1.3 of the PSD at all sizes and heights (Fig. 7d). Neglecting the variability of , the cooling rate due to melting is equal to the particle concentration ( of the PSD) multiplied by the rate of ice loss due to melting. The melting rate equation [Eq. (9)] includes two PSD-dependent terms: the particle capacitance c, which is approximately proportional to following Eq. (B1), and the ventilation coefficients , which (assuming ) are proportional to . Thus, the cooling rate due to melting is roughly constant and proportional to ≈M1.3 of the PSD (Fig. 7d).
The evolution of Z and is more complex. In contrast to the cooling rate from melting, the contributions to Z peak in much larger size bins and at a lower height (490 m below the 0°C level; Fig. 7b). As melting begins, Z is proportional to ≈M4 because of the inverse relationship between D and [Eq. (1)] and becomes proportional to ≈M6 once melting is complete (Fig. 7d). However, during melting the forward- and backscattering amplitudes display nonmonotonic behavior with respect to D owing to the variable across the size spectrum and increasing rapidly once exceeds 0.1. This nonmonotonic behavior is further enhanced by the stabilization of particle orientations during melting. These factors result in a wide range of moments of proportionality across the size spectrum during melting, with a weighted mean of at the height of maximum Z. These differences in maxima height, bin contributions, and moments of proportionality help to explain the poor relationship between maximum Z and cooling rate. Conversely, while is proportional to in dry snow (Bukovčić et al. 2018) and close to – in rain (Sachidananda and Zrnić 1986), the weighted mean moment of proportionality during melting ranges between and and is at the height of the maximum, nearly the same moment of proportionality as the maximum cooling rate due to melting (Fig. 7d). This is due to the sensitivity of to both particle mass and axis ratio (Gorgucci et al. 2002), which at a particle’s peak contribution during melting are proportional to ≈M1.4 and ≈M−0.3, respectively. This fact, along with having dominant contributions from similar size bins as the peak cooling due to melting, helps to explain their pronounced relationship within the melting layer.
These conclusions appear to be robust in the face of uncertainty in determining the effective relative permittivity of melting snowflakes. Table 2 compares the RMSE of retrievals performed using three Polder–van Santen mixing formula topology combinations and the percent change with respect to the default weighted Maxwell-Garnett mixing formula used thus far. While the magnitude of the variables changes slightly depending on the mixing formula, the qualitative conclusions for each variable’s retrieval efficacy remain unchanged, with RMSEs varying by ≈7% or less.
However, a cursory investigation reveals that the relation between the maximum and the maximum cooling rate is sensitive to . Figure 8 demonstrates the impact of the assumed on the correlation between the maximum and the maximum cooling rate in the melting layer. While there is ample evidence of the nonsphericity of ice and snow particles, it is instructive to compare against the extreme case of spherical particles. In comparison to the assumption of = 0.6 for unrimed snow aggregates, the correlation is strongly degraded if = 1.0, with an of 0.22 and an increase in RMSE of an order of magnitude over Fig. 8a and approaching that of Z (Fig. 6a). In addition, the values of are smaller by an order of magnitude despite simulating the same PSDs, and in many cases are so low as to be operationally impractical to use. When = 0.6, small particles responsible for the bulk of cooling contribute appreciably to as they melt from oblate snow crystals into nearly spherical raindrops. In contrast, sparse large particles with minimal contributions to do not undergo a large change in axis ratio while melting, beginning and ending as oblate particles. This situation is reversed if snow particles are assumed to be initially spherical, with small particles remaining nearly spherical and contributing little to during the melting process and large particles dominating the contribution to , resulting in much smaller values and a degraded relationship with the cooling rate. Suffice it to say that, subject to the model assumptions described in section 2 and the feasibility of its estimation, appears to be an attractive choice for retrieving the cooling rate in the melting layer and warrants further research.
5. Application of the 1D-MS: Investigation of sagging brightband signatures
One of the most pronounced features to be studied using QVPs is the “sagging” brightband signature, a temporary and often sudden downward excursion of the bright band (Ryzhkov et al. 2016; Kumjian et al. 2016; Erlingis et al. 2018). Kumjian et al. (2016) examined sagging bright bands in detail and identified several common characteristics, including an increase in the Z and brightband maxima and depth, a decrease in the height of the brightband maxima, an increase in Z and decrease in above the bright band, and an enhancement in both Z and below the bright band. An example QVP of a sagging bright band from the Saint Louis, Missouri (KLSX), WSR-88D radar is shown in Fig. 9 from the trailing stratiform portion of a strong MCS. In agreement with past observations, there is an enhancement of Z and during the period of sagging and a lowering of the height of the brightband maxima by approximately 400 m. In the sagging region is also enhanced, with values exceeding 0.4° km−1, along with a distinct lowering of the field suggesting a delay in melting. However, there are some discrepancies, including the lack of a clear increase and decrease in Z and above the bright band, respectively. Despite these observations, the cause of the sagging brightband signature is not fully understood. Using in situ thermodynamic and particle data during a sagging brightband event in stratiform frontal rain, Kumjian et al. (2016) found strong evidence supportive of riming, and an earlier version of the 1D-MS successfully reproduced the observed brightband signature when the riming factor and precipitation intensity were covaried. However, the observed signatures of Z and are also consistent with aggregation, with in situ observations from a second case in the trailing stratiform region of a classical MCS containing aggregates and environmental conditions unsupportive of riming. Kumjian et al. (2016) propose and discount, but do not simulate, a number of other potential causes of brightband sagging, including the role of isothermal layers. Thus, it is worth investigating these potential causes of brightband sagging using the latest 1D-MS.
Four possible causes of sagging bright bands were investigated in isolation to study their impacts on the resultant brightband signature: changes in aggregation, precipitation intensity, relative humidity, and vertical velocity. Following Kumjian et al. (2016), the parameters of interest were modulated using the Gaussian-based profile (σ = 12 min) shown in Fig. 10 for a total of 150 min to facilitate comparisons with Fig. 9. For all of the causes examined, environmental temperature and moisture feedbacks are turned off; the latent heating rates shown are instantaneous but do not affect the evolution of the temperature field. The 1D-MS lacks a number of processes that, in reality, act in concert with latent heating and cooling to modify the environmental temperature and moisture field including vertical motion (e.g., adiabatic warming), turbulent mixing, and advection. These processes can act to make the growth of isothermal layers self-limiting and may even offset the effect of diabatic cooling entirely (Kain et al. 2000). However, these factors are not straightforward to include in the 1D-MS. By turning off environmental feedbacks, a quasi balance between these effects is approximated and allows the instantaneous effect of the simulated changes to be evident. Unless otherwise stated, = 6.0°C km−1, the RH gradient is −3.0% °C−1, and = 1.0 is assumed.
To simulate the effects of aggregation above the melting layer, is decreased while conserving ice water content (IWC), given by (Boudala et al. 2006)
where is the complete gamma function and where α = 9.361 × 10−5 g mm−β and β = 2.08 are the parameters in the implicit 1D-MS mass–diameter relation for snow (of the form ) and fall roughly within the range of values reported in the literature (e.g., Mitchell et al. 1990; Szyrmer and Zawadzki 2010). The IWC is held constant at 1.0 g m−3 while is varied from a background value of 1.3 mm−1 to a perturbed value of 0.3 mm−1. Equation (23) is used to find the corresponding using the found from Eq. (21). The resultant brightband signature bears a striking resemblance to observations of sagging bright bands reported in the literature, with an increase in Z and maxima and a lowering of their height of occurrence, an increase in the depth of the bright band, an increase in Z and decrease in above the bright band, and an increase in Z and below the bright band during periods of enhanced aggregation (between 60 and 90 min; Figs. 11a,b). The height of the maximum sags by 330 m, in good agreement with the reported sagging of 351 m in Kumjian et al. (2016). However, the height at which melting begins is unaffected (Fig. 11d) and is decreased during periods of sagging (Fig. 11c), in contrast to Fig. 9c. Despite these discrepancies, it is noteworthy that the simulated effects of aggregation alone reproduce almost all of the observed features of sagging bright bands.
To simulate the effects of changing precipitation intensity, is held constant at 1.0 mm−1 while the IWC is perturbed from 1.0 to 2.0 g m−3. In contrast to aggregation, varying the precipitation intensity does not result in a visual sagging of the bright band (Figs. 11e–g) or changes in above and below the bright band. However, both the increase in the magnitude of (Fig. 11g) and the thickness of the bright band are simulated.
To simulate the effect of infiltrating dry air, such as might occur with a rear-to-front inflow jet (e.g., Houze 1994), the relative humidity is varied from a constant 100% to 70%. The height of the maximum and the onset of melting descend by 380 m (Fig. 11j) and 440 m (Fig. 11l), respectively. However, sublimation above the melting layer decreases Z (Fig. 11i) and (Fig. 11k) and thins the bright band, in contrast to what is observed.
Finally, to simulate the effect of varying vertical velocity (assumed to be 0.0 m s−1 until now), w is varied from 0.0 to −1.5 m s−1 for the background PSD used in the precipitation intensity simulation. This is at the extreme end of what is observed, with most studies indicating mesoscale downdrafts in stratiform precipitation not typically exceeding 1 m s−1 (e.g., Yang and Houze 1995a,b; Schuur and Rutledge 2000; Kumjian et al. 2016). Even with a downdraft of 1.5 m s−1, the height of the maximum only descends by 80 m (Fig. 11), along with little to no change in the magnitude of Z, , , and brightband thickness (Figs. 11m–o). Thus, the impact of vertical velocity on particle fall speeds alone cannot explain sagging bright bands.
Although these analyses are instructive, it is likely that in reality many of these variations are coupled. For example, it stands to reason that more intense precipitation rates may enhance aggregation, or that downdrafts near the melting layer may form from cooling due to sublimation in dry air aloft. While no one physical process examined can fully explain every observed sagging brightband characteristic, some combination of processes may be able to. An example of a potentially more realistic modulation is shown in Figs. 11q–t. Here, both aggregation and precipitation intensity are covaried as described previously, and environmental temperature and moisture feedbacks are now included. With all three of these factors, all of the observed brightband characteristics during periods of sagging are simulated. The height at which melting begins descends (Fig. 11t), with increases in brightband depth and Z (Fig. 11a), (Fig. 11r), and (Fig. 11s) within the bright band, along with the corresponding changes in Z and above and below the melting layer. Because of the aforementioned exclusion of restorative processes that prevent the indefinite deepening of isothermal layers, the return of the bright band to its presagging level is not simulated. While the goal of these simulations is not to disprove the contribution of riming to sagging bright bands, these results suggest that sagging due to other causes—namely, enhanced precipitation and aggregation and the resultant cooling of the atmosphere—cannot be entirely ruled out, and it is possible that there are multiple pathways that result in sagging brightband signatures.
6. Summary and discussion
In this study, a one-dimensional model of melting snow is used to study the polarimetric characteristics of the radar bright band and their connection to the environment and microphysical processes within the melting layer. The brightband thickness and the magnitude of Z and are shown to be sensitive to the environmental lapse rate and relative humidity profile. Based on simulations of a wide variety of snow PSDs, Z appears to hold little value for retrieving the cooling rate within the melting layer. However, for a given environment, the maximum in the bright band appears to be robustly correlated with the maximum cooling rate. This promising result suggests that, given knowledge of the environment (e.g., from a short-term numerical weather prediction model) and an accurate estimate of in the melting layer, the evolution of the maximum cooling rate in the melting layer may be retrieved. Future work should investigate the potential for time–height estimates of latent heating of real-data cases using the 1D-MS initialized with QVP-derived estimates of mean volume diameter and total number concentration (Ryzhkov et al. 2018), as done for evaporation in Xie et al. (2016, their Fig. 14). Until now, information derived from polarimetric radar data has been limited to microphysical quantities (e.g., quantitative precipitation estimation) and qualitative applications (e.g., tornado debris detection). The novel approach presented here for retrieving thermodynamic information from dual-polarization radar data suggests the potential for a shift toward a more holistic approach for using polarimetric radar data for both observational retrievals and for assimilation into numerical weather prediction models.
Although these results pertaining to the maximum cooling rate are promising, a number of aspects remain to be explored. The spatial distribution and height of the maximum cooling rate was not investigated but may be of importance (e.g., Szyrmer and Zawadzki 1999), and the impact of assimilating improved retrievals of melting-layer cooling rate into numerical weather prediction models should be investigated. There are also a number of improvements to the 1D-MS that should be considered for their impact on these results, particularly the inclusion of aggregation, collision–coalescence, and breakup processes, as well as some consideration of buoyancy and vertical motion (e.g., as done in Srivastava 1985, 1987). From a broader perspective, the adoption of QVPs and the large amount of information they contain begs for a comprehensive one-dimensional model coupled to a polarimetric operator that expands upon the 1D-MS and includes multiple ice habits and options for all microphysical processes within the column (e.g., initial ice generation, deposition, secondary ice generation, breakup, etc.).
The work presented here stemmed from the author’s dissertation research. The authors thank Amanda Murphy and John Krause for processing the quasi-vertical profile radar data, Greg McFarquhar and Wei Wu for their helpful discussions and providing us the in situ MC3E data, and Darrel Kingfield for his bias calculations, as well as the anonymous reviewers for their constructive comments that improved this manuscript. Funding was provided by NOAA/Office of Oceanic and Atmospheric Research under NOAA–University of Oklahoma Cooperative Agreement NA11OAR4320072, U.S. Department of Commerce. Additional funding was provided by the U.S. Department of Energy Atmospheric System Research grant DE-SC0014295. Part of the computing for this project was performed at the OU Supercomputing Center for Education and Research (OSCER) at the University of Oklahoma (OU).
Table of Variables
1D-MS Thermodynamic Terms
The particle capacitance c reflects the impact of the particle’s shape on its rate of mass transfer. For spherical particles, c is the particle radius. For oblate spheroids, c is given by (McDonald 1963)
where is the axis ratio of the particle at any stage of melting and the ellipticity of the spheroid x is equal to .
The ventilation coefficients for vapor and heat account for the removal of vapor and heat from the falling particle because of air motion. Following the empirical relations of Hall and Pruppacher (1976), the ventilation coefficient for vapor can be defined by
where , is the Schmidt number given by
and is the Reynolds number given by
where is the particle characteristic length, which for an oblate spheroid can be expressed as (Pruppacher and Klett 1997)
is the diffusivity of water vapor (Hall and Pruppacher 1976),
and is the kinematic viscosity of air,
where the dynamic viscosity of air is given by
The ventilation coefficient for heat follows a similar form as Eq. (B2), except uses the Prandtl number instead of the Schmidt number (i.e., ), where is given by
where the thermal diffusivity of air is defined as
and the thermal conductivity of air is given by
The latent heat of sublimation and fusion are assumed to be constant, while the latent heat of vaporization is found following Kumjian and Ryzhkov (2010):