X-band polarimetric radar observations of winter storms in northeastern Colorado on 20–21 February, 9 March, and 9 April 2013 are examined. These observations were taken by the Colorado State University–University of Chicago–Illinois State Water Survey (CSU-CHILL) radar during the Front Range Orographic Storms (FROST) project. The polarimetric radar moments of reflectivity factor at horizontal polarization ZH, differential reflectivity ZDR, and specific differential phase KDP exhibited a range of signatures at different times near the −15°C temperature level favored for dendritic ice crystal growth. In general, KDP was enhanced in these regions with ZDR decreasing and ZH increasing toward the ground, suggestive of aggregation (or riming). The largest ZDR values (~3.5–5.5 dB) were observed during periods of significant low-level upslope flow. Convective features observed when the upslope flow was weaker had the highest KDP (>1.5° km−1) and ZH (>20 dBZ) values. Electromagnetic scattering calculations using the generalized multiparticle Mie method were used to determine whether these radar signatures were consistent with dendrites. Particle size distributions (PSDs) of dendrites were retrieved for a variety of cases using these scattering calculations and the radar observations. The PSDs derived using stratiform precipitation observations were found to be reasonably consistent with previous PSD observations. PSDs derived where riming may have occurred likely had errors and deviated significantly from these previous PSD observations. These results suggest that this polarimetric radar signature may therefore be useful in identifying regions of rapidly collecting dendrites, after considering the effects of riming on the radar variables.
Winter storms often have severe impacts on society, both economically and through increased risks to life and property. However, accurate short-term forecasting of these events remains a challenge. One of the greatest areas of uncertainty in these forecasts comes from a poor understanding of the microphysical processes that contribute to the development of heavy snow.
Snow growth through aggregation occurs most rapidly when dendritic ice crystals, with an intricate structure of branches that more readily interlock, are present. Dendrites are favored to grow in temperatures near −15°C and ice supersaturations above 0.15 (e.g., Bailey and Hallett 2009). Aggregates that form rapidly through the collection of these ice crystals can increase the snowfall accumulation rate (Lamb and Verlinde 2011) and decrease visibility near the ground (Rasmussen et al. 1999). Thus, from an operational perspective, the growth of dendrites and their subsequent aggregation are important microphysical processes and should be identified in available observations.
Dual-polarization radar provides information related to the orientation, shape, and diversity of the sampled hydrometeors. The National Weather Service array of WSR-88D instruments has recently been upgraded to dual polarization, allowing for this valuable information to be used operationally (e.g., Doviak and Zrnić 1993; Zrnić and Ryzhkov 1999; Bringi and Chandrasekar 2001; Ryzhkov et al. 2005b; Kumjian 2013a,b,c, and references therein). These radars have been used previously to study snow growth; enhancements in specific differential phase KDP and differential reflectivity ZDR above the melting layer have been observed by Ryzhkov and Zrnić (1998), Trapp et al. (2001), Wolde and Vali (2001), Kennedy and Rutledge (2011), Andrić et al. (2013), Bechini et al. (2013), Schneebeli et al. (2013), Kumjian et al. (2014), and Griffin et al. (2014) and have been linked to the presence of dendrites and platelike ice crystals.
Kennedy and Rutledge (2011) found that ZDR and KDP maxima detected near the −15°C level by an S-band radar were correlated with increased snowfall rates at the ground. These observations, along with reflectivity factor at horizontal polarization ZH, were compared with electromagnetic scattering calculations using the -matrix method (e.g., Waterman 1971; Mishchenko 2000) for distributions of ice particles. The distributions used in their study were based on airplane ice particle size distribution (PSD) observations (Lo and Passarelli 1982). The observed KDP and ZH values were reproduced with these PSDs; however, the resulting ZDR signature was greater than that found in the observations.
Andrić et al. (2013) found similar ZDR and KDP signatures in S-band radar observations of winter storms. This study also presented scattering calculations, but used the Rayleigh approximation for spheroids and PSDs of plates, dendrites, and aggregates from a one-dimensional microphysical model. As in Kennedy and Rutledge (2011), the calculated radar variables failed to reproduce all of the observed polarimetric signatures. Large concentrations of horizontally oriented ice crystals produced through secondary ice generation, not present in their model, were proposed to explain the error in the simulated polarimetric variables.
Bechini et al. (2013) also observed enhancements in ZDR and KDP with X- and C-band radars around the −15°C level within stratiform precipitation, and attributed the signature to dendritic crystals. The scattering calculations used to support this conclusion were based on those performed by Kennedy and Rutledge (2011), with T-matrix calculations performed at both X- and C-band. It was found that KDP scales inversely proportional to wavelength and ZH is independent of wavelength, suggesting that these particles were electromagnetically small. No discussion of the correspondence between the scattering calculations and the radar observations was provided.
The lack of agreement in these studies between electromagnetic scattering calculations and ZH, ZDR, and KDP observations of mixtures of dendrites and aggregates necessitates both a better understanding of the variability of the radar signatures in regions of possible dendritic growth and more advanced modeling of ice particles. This study uses X-band polarimetric radar data from several winter storms in Colorado to examine a range of signatures in regions favorable for dendritic crystal growth. Advanced scattering calculations are also used to determine PSDs that match each of the observed radar signatures. Unfortunately, there are no in situ measurements available in these storms, so the physical validity of the resulting PSDs is assessed by comparison with previous observational studies of ice crystal PSDs.
A description of the data and interpretation of the polarimetric radar variables used in this study are introduced in the next section. An overview of the meteorological background for the events of study is presented in section 3. Section 4 provides a detailed analysis of the radar signatures found within the events and statistics characterizing these signatures. A description and analysis of the scattering calculations used in this study are found in section 5. Section 6 contains a sensitivity analysis of the retrieved PSDs to assumed parameters. The plausibility of these PSDs and their physical implications are discussed in section 7. The main conclusions of this study are provided in section 8.
2. Observational data and methods
For this study, data from the Front Range Orographic Storms project (FROST; Kumjian et al. 2014) in northeastern Colorado were used. Significant winter storm events on 20–21 February, 9 March, and 9 April 2013 are the focus of this paper. Data collected during these events include radar observations from the Colorado State University–University of Chicago–Illinois State Water Survey (CSU-CHILL) radar (Brunkow et al. 2000) and soundings launched at the Marshall Field Site (MFS), depicted graphically in Kumjian et al. (2014, their Fig. 1). The CHILL radar and the MFS are at elevations of 1432 and 1742 m MSL, respectively.
The CHILL radar has recently been upgraded to a dual-wavelength system with X- and S-band frequencies (Junyent et al. 2015). In FROST, only the X-band frequency was used because of its smaller beam width (0.3°) and greater sensitivity to KDP. CHILL data are available for 1900 UTC 20 February–1200 UTC 21 February, 1029–2202 UTC 9 March, and 0041–2233 UTC 9 April. The radar data contain a series of PPI and RHI scans performed every ~12 min. The PPI and RHI scans were mostly performed at a set of fixed elevation and azimuth angles, respectively. The azimuth angles of the RHI scans occasionally varied depending on the appearance of precipitation targets. A full description of this field project can be found in Kumjian et al. (2014).
Pristine ice crystals generally have mean horizontal orientations and scatter anisotropically, making polarimetric radar an attractive choice to study these particles. The dual-polarization variables most often used to identify oblate particles are ZDR and KDP. Differential reflectivity ZDR (dB) is related to the ratio of the reflectivity factors at horizontal and vertical polarization (Zh and Zυ, respectively; mm6 m−3). It provides information about the reflectivity-weighted aspect ratio of the sampled targets. Specific differential phase KDP is one-half of the range derivative of the propagation differential phase shift ΦDP between horizontally and vertically polarized electromagnetic waves; ΦDP accumulates at a rate depending on the size, number concentration, and dielectric constant of the hydrometeors (i.e., is related to mass) as well as their aspect ratio, producing a differential slowing of the orthogonally polarized waves.
For the processed CSU-CHILL data, KDP is calculated from ΦDP using the Wang and Chandrasekar (2009; hereinafter WC09) method, where an adaptive algorithm smooths ΦDP to a degree inversely related to the signal-to-noise ratio. Before any range derivatives are computed, however, the data are first flagged in regions of high noise using the dispersion of ΦDP for a sample of consecutive gates. These flagged regions are set to KDP = 0° km−1, which may introduce a large number of erroneous zero values that can skew statistical analyses and microphysical interpretations. To eliminate these points, an alternate method of estimating KDP (described in the appendix) is used for scans containing large ZDR values in regions of high noise. Data that use this new scheme are presented in later sections.
3. Description of events
Information about the synoptic-scale features present during these events is taken from the North American Regional Reanalysis (NARR; Mesinger et al. 2006). One feature common to all events was a 500-hPa trough progressing eastward through the central Rocky Mountains (Fig. 1). The motion and position of this trough and its associated surface features helped dictate whether low-level upslope flow was present in northeast Colorado.
a. 20–21 February 2013
The 20–21 February 2013 event was associated with a broad trough over the mountainous west. A closed 500-hPa low was present in southern Arizona at 2100 UTC 20 February (Fig. 1a). At this time, the 750-hPa (roughly 1 km above the CHILL radar, hereinafter above radar level, or ARL) winds were mainly easterly to southerly in northeastern Colorado. By 1200 UTC 21 February, the 500-hPa trough axis was centered over central New Mexico and the low-level flow had become more northerly (Fig. 1b). Cold-air advection was ongoing throughout the duration of the sampling period; the −15°C level decreased in height from 2.8 to 2.5 km ARL between 0000 and 1200 UTC 21 February, based on radiosonde measurements from Denver, Colorado (DNR). The sounding taken at the MFS at 0017 UTC indicated the −15°C level was at about 2.6 km ARL (Fig. 2a).
One of the notable features visible during this event was a narrow band of enhanced ZH (>25 dBZ) that approached the radar from the southeast between around 2300 UTC 20 February and 0200 UTC 21 February. An RHI scan taken at 0058 UTC (along the 181.8° azimuth, roughly perpendicular to the major axis of this band) shows enhancements in both ZDR and KDP (>2.0 dB and 1.5° km−1, respectively) between 2.5 and 4.0 km ARL and between 28 and 35 km range (Fig. 3). After this period, lighter, generally less-organized echoes were observed as weak upslope flow became the primary contributor to ascent.
b. 9 March 2013
At 1200 UTC 9 March, shortly after the onset of data collection, a closed 500-hPa low was centered over the Four Corners region with generally northeasterly flow at 750 hPa in northeastern Colorado (Fig. 1c). By 1800 UTC, the closed 500-hPa low was positioned over southern Colorado, resulting in a stronger northerly component to the 750-hPa wind near the CHILL radar (Fig. 1d). This general flow regime continued through the termination of data collection. During this event, low-level cold-air advection was occurring in northeastern Colorado. Mandated soundings taken at DNR indicate the −15°C level decreased in height from 3.1 to 2.6 km ARL between 1200 UTC 9 March and 0000 UTC 10 March. The sounding taken at the MFS at 1319 UTC (Fig. 2b) indicated the −15°C level was around 3 km ARL.
Between 1100 and 1400 UTC, large variations in echo-top height with range were visible in PPI and RHI scans. An example of a feature with significant vertical extent is found in an RHI taken at 1249 UTC, shown in Fig. 4. This time period was coincident with weaker, northeasterly 750-hPa flow. As these winds became stronger and more north-northeasterly, PPI and RHI scans show a more stratiform appearance to the precipitation, with less variation in the echo-top height seen with range. The PPI in Fig. 5, taken at an elevation angle of 2.8° at 2044 UTC 9 March, shows enhancements in ZDR and KDP between the range rings corresponding to −10° and −15°C. Given that the radar is located north of DNR, the temperature level heights may have been lower near the CHILL than those indicated by the sounding. Thus, the true range of the −15°C ring may have been positioned closer to the KDP and ZDR enhancements.
c. 9 April 2013
The 9 April event was associated with a closed 500-hPa low centered over southern Utah at 0300 UTC. The more northern position of this feature relative to the other events allowed for thunderstorms to develop within range of the CHILL radar. As the 500-hPa trough progressed eastward, the 750-hPa winds became more northeasterly in northern Colorado by 0900 UTC, shown in Fig. 1e. By 1800 UTC, these winds backed to north-northeasterly, enhancing the upslope flow and increasing the advection of cold air (Fig. 1f).
Following the thunderstorm observations, ZH > 20 dBZ was observed in PPI and RHI scans. These areas of precipitation mainly translated from southeast to northwest. An RHI scan (Fig. 6) taken at 0510 UTC through a section of one intense ZH band shows substantial KDP enhancements (>1.5° km−1) between 4 and 7 km ARL and 25–35 km range (Fig. 6d). Positive ZDR (~1.0–1.5 dB) was also found in this region (Fig. 6b).
From 1200 to 1600 UTC, initially light (~10 dBZ) and sparse ZH echoes gradually intensified to ~20 dBZ and increased in areal coverage. These echoes had ZDR > 1.5 dB through a layer generally below 2 km ARL. This coincided with backing, low-level upslope flow and cold-air advection, where 750-hPa temperatures were around −15°C in the vicinity of CHILL. By 1600 UTC, a more significant band of enhanced ZH was observed moving slowly toward the north-northwest. At around 2000 UTC, this band extended north of the radar with enhanced values of ZDR (>3.0 dB) visible in the PPI scans (not shown).
4. Analysis of radar signatures
To explore the variability of the radar observations taken near the dendritic growth zone, distributions of the radar moments with height from several RHI scans were examined. These distributions were constructed by first binning the radar gates by height. At each height bin, the radar observations were then binned by increments of ZH, ZDR, and KDP. Only data within the horizontal distances from the radar containing enhancements in the polarimetric variables were used in the binning procedure, isolating the signatures of interest.
The first case (hereinafter “case 1”) presented is from 1249 UTC 9 March. As shown earlier, an intense convective feature was visible south-southwest of the radar at this time (Fig. 4). Within this RHI, data were extracted at all heights between horizontal distances of 49–56 km from the radar, where the greatest enhancements in KDP and ZDR were present. Figure 7a shows the count of radar gates within each height and ZDR bin, with a peak in median ZDR of 1.1 dB at a height of 4.2 km ARL. The corresponding distributions for KDP and ZH are shown in Figs. 7b and 7c, respectively. The peak median KDP is 1.6° km−1 at about 3.7 km ARL; the median ZH increases from 19 to 23 dBZ between 4.2 and 3.7 km ARL. The sounding-indicated −15°C level at this time was about 3.3 km ARL, several hundred meters below the ZDR and KDP maxima. Below their respective peaks, both ZDR and KDP decrease significantly, with ZDR values of 0.5 dB and KDP near 0° km−1 at 1.5 km ARL; ZH decreases within this layer to 19 dBZ at 1.5 km ARL.
Enhancements in ZDR and KDP observed during the period of more stratiform, upslope precipitation at 1951 UTC 9 March (hereinafter “case 2”) are shown in Fig. 8. These data were extracted from the corresponding RHI at horizontal distances 10–20 km from the radar (Fig. 9). During this time, the height of the −15°C level was approximately 2.8 km ARL, as indicated by the 0000 UTC DNR sounding. A peak in median ZDR of 1.4 dB occurs at 2.6 km ARL, with values between 0.5 and 1.0 dB below 2.2 and above 2.8 km ARL. The peak in median KDP of 0.8° km−1 is found at 2.5 km ARL, coinciding with ZH = 13 dBZ at 2.6 km ARL. Median ZH increases steadily to around 23 dBZ at 1.4 km ARL.
The latter portion of the 9 April event was also mainly driven by low-level upslope flow. However, one unique signature visible in RHI scans was high (3.5–5.5 dB) ZDR collocated with ZH < 5 dBZ (Figs. 10a,b; hereinafter “case 3”). Because this RHI points nearly due south, the inbound radial velocities above 2.5 km ARL indicate a southerly component to the flow at midlevels; the outbound velocities below 2 km ARL indicate a northerly component to the low-level flow (Fig. 10c). The value of KDP is also modestly enhanced at ~0.3° km−1 (Fig. 10d). Note that the KDP estimates for case 3 use the alternate procedure described in the appendix.
The distributions with height of the radar variables for case 3 are shown in Fig. 11. These data were extracted from the RHI in Fig. 10 at ranges of 20–30 km from the radar. The peak in median ZDR is ~4.5 dB at 3.3 km ARL, although significant counts of gates with ZDR values >4.5 dB are found at this height, suggesting that the ZDR of pure ice crystals may be above 4.5 dB. The vertical distribution of KDP shows maxima around 0.7 and 2.9 km ARL; at the height where the median ZDR is greatest, the median KDP is still enhanced at ~0.3° km−1. There are nonzero counts of gates containing KDP between −1.0° and 2.0° km−1, suggesting a high amount of noise in the KDP estimation at this altitude. Note that these data would be set to 0° km−1 using the WC09 scheme.
Because the most recent MFS sounding was launched ~11 h prior to this RHI, the 0000 UTC sounding from DNR is used to better represent the temperature profile. This sounding indicates that nearly the entire profile is below −10°C. Between roughly 780 and 650 hPa, temperatures vary from −20° to −15°C. Above 650 hPa, the temperature increases to around −13°C, before it decreases back to −15°C at about 580 hPa. As a result, two layers may have been conducive to dendritic growth. The 750-hPa temperatures (approximately 1 km ARL) from the NARR model analysis for 1800 UTC 9 April were below −15°C near the radar (Fig. 1f). The northeasterly flow at this time suggests cold-air advection, with the low-level air mass sampled by the DNR sounding likely moving over the radar at an earlier time. The 0000 UTC DNR sounding therefore may be a fair representation of the low-level temperatures near CHILL when the RHI in Fig. 10 was taken at 1940 UTC.
Values of ZDR and ZH within the 1940 UTC RHI1 are generally negatively correlated for ZH > 6 dBZ (Fig. 12a). For ZH < 6 dBZ, the relationship becomes less obvious. The two relative frequency maxima below 0 dBZ suggest the presence of different crystal populations: the maxima for ZDR < 2 dB may represent more isometric crystals; the secondary peak for ZDR > 4 dB indicates nonspherical particles. Joint relative frequency analyses performed for cases 1 and 2 (not shown) also suggest an inverse relationship between ZDR and ZH, although this relationship was not as robust due to the smaller range in ZDR observed in these cases.
Across these three illustrative cases, ZH, ZDR, and KDP all exhibit large variability in regions with temperatures near −15°C. The ZH and ZDR tend to be negatively correlated, reflecting the strong weighting of ZDR by particle size: if aggregation dominates hydrometeor growth, the larger, more isotropically scattering aggregates will decrease ZDR substantially and increase the observed ZH. How KDP varies with ZDR and ZH is less well defined because KDP depends on both ice particle number concentration and aspect ratio and generally is not as affected by large aggregates. However, KDP increases as dendrite concentration increases; the increased dendrite population promotes enhanced aggregation and thus increased ZH. We see this relationship in the three cases, where the KDP maxima are positively correlated with the median ZH at the same height.
The high degree of variability in the radar signatures presented above suggests that the microphysical characteristics of these cases are distinct. For instance, this variability suggests that particles at various stages of growth (e.g., pristine crystals and aggregates) are mixed together in the same radar gate. Because in situ ice particle measurements are unavailable for these cases, we cannot be certain that dendrites were indeed responsible for the observed radar signatures. However, forward modeling, in which the radar signatures of dendritic ice crystals are simulated, can support or reject the hypothesis that these ice crystals are present. To simulate the populations of ice particles found in the dendritic growth zone, electromagnetic scattering calculations and PSDs of these crystals are needed.
5. Scattering calculations
Previous studies have used PSDs from microphysical models (Andrić et al. 2013) or observational studies (Kennedy and Rutledge 2011) to calculate the radar observables. However, these studies failed to fully reproduce the observed radar signatures, suggesting there were errors in modeling the ice crystal scattering properties or in the PSDs. The variability in the radar signatures near −15°C shown in section 4 suggests there is also variability in the ice crystal PSDs associated with these signatures.
Therefore, to account for PSD error and variability, we take a different approach in which PSDs for ice crystals are determined based on the radar observations, after first separating out the contributions to the polarimetric variables from aggregates. If these pristine ice PSDs are physically realistic, they may provide plausible descriptions of the ice particle size spectrum for certain polarimetric radar signatures in the dendritic growth zones. The first step needed to retrieve PSDs from the radar observations is to calculate scattering properties of individual crystals with varying size.
a. Description and monodispersed distribution results
To improve upon previous work that used simplified scattering calculations, we use electromagnetic scattering calculations for pristine ice particles from the Pennsylvania State University (PSU) ice crystal database (e.g., Botta et al. 2013; Lu et al. 2014), which employs the generalized multiparticle Mie (GMM; Xu 1995; Xu and Gustafson 2001) method for a variety of ice crystal morphologies. The GMM method involves composing a particle of nonoverlapping spheres and solving the scattering of each sphere analytically. To ensure this method’s accuracy, the mass distribution within the particle must be represented well by the collection of spheres. See Botta et al. (2013) and Lu et al. (2014) for further details.
Based on prior work and the observed temperatures in each case, hexagonal plates and dendrites were chosen to represent the likely hydrometeors sampled by the radar. Crystal aspect ratios are based on empirical formulas adapted from Auer and Veal (1970) (Table 1). To obtain a larger range of ice crystal aspect ratios, dendrite thicknesses were multiplied by 0.5 and plate thicknesses by 0.5 and 2.0. The database’s dendrites (plates) have maximum dimensions varying from 0.5 to 5.65 mm (from 0.10 to 3.27 mm) (Lu et al. 2014). We assume the dendrites are oriented according to a two-dimensional Gaussian canting-angle distribution (Ryzhkov et al. 2011) with a mean of 0° and standard deviation σ varying from 1° to 20°. Previous studies have assumed σ = 10°–15° for pristine ice crystals (e.g., Matrosov et al. 2005; Kennedy and Rutledge 2011).
For illustrative purposes, we assume the monodispersed ice crystals produce a total ZH = 3 dBZ (i.e., representative of case 3). The resulting calculations for dendrites with thicknesses of 1.0 and 0.5 times that determined by the empirical relationship (hereinafter d1.0 and d0.5, respectively) and plates with thickness 2 times that prescribed by the empirical relationship (hereinafter referred to as p2.0) are shown in Fig. 13. Note that ZDR decreases as σ increases because of the smaller total projection of the horizontally polarized signal onto the axes of the ice particles’ maximum dimensions. Since Zh (mm6 m−3) is approximately proportional to the particle equivalent diameter to the sixth power, the smallest crystals require much larger number concentrations to produce ZH = 3 dBZ and thus produce KDP values several orders of magnitude higher than the largest crystals.
Differential reflectivity varies more with size for the p2.0 crystals (Fig. 13c) than for the d1.0 (Fig. 13b) crystals. For example, ZDR of the d1.0 crystals with σ = 20° only changes from 4.6 to 4.8 dB; however, ZDR of the p2.0 crystals with σ = 20° increases from 2.7 dB for the 0.1-mm crystals to 5.2 dB for the 2.5-mm crystals. The difference in ZDR between the smallest and largest p2.0 crystals increases to about 4.2 dB for σ = 1°. The relationship between maximum dimension and ZDR for the plates is attributed to the particles becoming more oblate with increasing maximum dimension and retaining the same effective density. The dendrites also become more oblate with size; however, their branches become farther apart (Botta et al. 2013), lowering their effective density and thus counteracting the decrease in aspect ratio (Lu et al. 2014).
Such calculations for monodispersed ice crystals can only be used to simulate radar observations corresponding to very early stages of ice crystal growth (low ZH and large ZDR). In case 3 above, the peak median ZDR value of 4.5 dB was associated with KDP = 0.3° km−1 and ZH = 3 dBZ. The scattering calculations show that p2.0 plates with maximum dimensions of 0.6 mm, a number concentration ~104 m−3 and σ = 18° best match the observations, whereas σ closer to 10° results in similar KDP but ZDR around 5.5 dB. The distribution of ZDR values near 3.3 km ARL in case 3 (Fig. 11a) shows a significant number of points clustered around 5.0 dB, which may be more representative of the ZDR contribution solely from pristine ice crystals.
In most situations, however, differential depositional growth conditions may widen the PSD. Thus, monodispersed distributions are less physically realistic for a population of ice crystals in a more mature development stage. Further, aggregation may play an increasingly important role as crystals develop secondary habit characteristics. Indeed, the other cases above were associated with lower ZDR, implying that an additional hydrometeor population that decreases the total ZDR is present.
b. Distributions of particle sizes and mixtures with aggregates
To simplify the retrievals, we assume that only aggregates and pristine ice crystals are present in radar sampling volumes and thus the radar-observed ZH, ZDR, and KDP are each composed of contributions from both ice crystals and aggregates. Since the pristine ice PSD is desired, we must separate the contributions to the radar variables from aggregates and pristine crystals. This requires knowledge of ZH, ZDR, and KDP values associated with snow aggregates. Large aggregates make a negligible contribution to KDP (e.g., Thompson et al. 2014) because they tumble and have low effective density. Similarly, aggregates produce small (but nonzero) ZDR values. Midlatitude observations (e.g., Ryzhkov and Zrnić 1998; Ryzhkov et al. 2005a; Homeyer and Kumjian 2015) reveal aggregate ZDR values between 0.0 and 0.6 dB at X, C, and S bands. Scattering calculations of aggregates composed of stellar crystals in the PSU database (not shown) also fall within this range. Thus, we assume that any population of aggregates, by itself, has a ZDR value between 0.0 and 0.6 dB and KDP = 0° km−1, with a ZH value to be determined. The sensitivity of the PSD retrievals to these assumptions is presented in section 6.
Because we assume aggregates generate no KDP, all signals from these particles are contained in the ZH and ZDR observations. We can write Zh and Zdr in terms of the aggregate and ice crystal contributions as follows:
where subscripts dr, h, and υ refer to the radar observables in linear units (mm6 m−3). The superscripts T, A, and I refer to the “total” (observed), “aggregate,” and “ice crystal” values, respectively. Given an assumed value for (between 0.0 and 0.6 dB), an additional piece of information () is needed to solve for the two unknowns ( and ); comes from the scattering calculations. In general, depends on the PSD and therefore must be calculated iteratively, as described in subsequent sections. For the d1.0 crystals, however, shows little variation with size and thus is nearly independent of the PSD.
where varies from 0 to , depending on . For example, if = , ice crystals are completely responsible for , and thus the sampling volume is assumed to consist only of ice crystals. Figure 14 shows the retrieved (dBZ) as a function of observed and for assumed values of and . The retrieved can then be used as a constraint on the pristine ice PSD. The sensitivity of the PSD retrieval to and other parameters is presented in section 6. These PSDs are often modeled by generalized gamma distributions: N(D) = N0Dμ exp(−λD), where N0, μ, and λ are the intercept, shape, and slope parameters. We will primarily use the exponential distribution, the special case when μ = 0. Another commonly used version of the generalized gamma distribution is where μ = 2, hereinafter referred to as the gamma2 distribution.
Two parameters are needed to define these PSDs (λ and N0), so another radar constraint is required in addition to . Because we assume ≈ 0, it follows that ≈ . Thus, provides the second constraint on the PSD. Using scattering calculations from the ice crystal database, and are computed (Fig. 15) for a wide (N0, λ)-parameter space2 following the formulas in Ryzhkov et al. (2011). The intersection of the and curves gives the N0 and λ values needed to match the radar observations. For a constant λ, increasing N0 increases both and as larger concentrations of crystals at all sizes are present. Increasing λ with constant N0 leads to decreases in and as both the total number concentration and median particle size decrease. The generally small angles between the and curves in Fig. 15 imply that both PSD parameters are highly sensitive to these radar constraints; this sensitivity is explored further in section 6.
In general (i.e., for the p1.0, p2.0, and d0.5 crystals), ZDR varies significantly with crystal size and thus and the constraint are dependent on the PSD. Therefore, we use an iterative procedure to simultaneously estimate , N0, and λ (shown schematically in Fig. 16), which follows four steps:
Make a first guess for based on the mean ZDR value of the dendrites used in the PSD.
Calculate using the first guess for (step 1), observations of and , the assumed , and trivially substitute for .
Estimate N0 and λ as described above.
Use these N0 and λ values to update the PSD-dependent .
Steps 2–4 are repeated using the updated until changes between consecutive estimates of N0 and λ are sufficiently small (<0.1%). Only a few iterations are needed for N0 and λ to converge. Unlike previous studies (Kennedy and Rutledge 2011; Andrić et al. 2013; Bechini et al. 2013), we now have ice crystal PSDs and polarimetric radar observations that are fully consistent.
c. PSD results
Table 2 shows PSDs derived using this methodology, along with the observed radar variables from each case. Observations used to constrain the PSDs come from the median height profiles (Figs. 7, 8, and 11; described in section 4), at 3.7, 2.5, and 3.3 km ARL for cases 1–3, respectively. Observations from the 0058 UTC 21 February case are the median values (not shown) of ZH, ZDR, and KDP derived from the RHI in Fig. 3 at a height of 3.1 km ARL and at distances 28–36 km from the radar. The retrieved N0 values for the diverse radar signatures vary by about an order of magnitude. A large range in λ is evident as well, with a ~5.3 mm−1 difference between maximum and minimum values.
The shallowest PSD is retrieved for case 1 (cf. Fig. 4), where the greatest contribution to comes from the dendrites with maximum dimensions of ~1.5 mm (Fig. 17). For dendrites larger than 1.5 mm, contributions to decrease more slowly with size than the contributions to retrieved for the other cases. The greatest contribution to , however, comes from dendrites with maximum dimensions between 3 and 4 mm, due to their increased concentrations and the size dependence of ZH.
The retrieved PSD for case 2 (cf. Fig. 9) is plotted in Fig. 18. A relatively high number concentration of smaller particles is derived for this case, with λ = 3.29 mm−1. The greatest contributions to from this PSD come from ice crystals with maximum dimensions between ~0.5 and 1.5 mm. These crystals have moderate intrinsic KDP values, but are in sufficiently high concentration to produce much of the observed enhancement. At greater maximum dimensions, the contributions to decrease substantially. Dendrites with maximum dimensions between 1.2 and 1.6 mm contribute the most toward , with a slower decrease in the contributions to with size. This is due to ZH depending more strongly on particle size than KDP.
The largest λ was retrieved for the case with largest (case 3). As mentioned in section 5a, these radar observations correspond fairly well to a monodispersed PSD of plates having a maximum dimension of 0.5 mm. Only the p2.0 crystal type was used here because the lower ZDR values of these crystals require lesser to match . This reflects the dearth of aggregation ongoing with plate crystals near echo top. The relatively low and modest for these crystals resulted in λ = 7.13 mm−1 (Fig. 19). The sharp peak in the contribution centered on the 0.5-mm crystal size bin suggests this PSD behaves similarly to the monodispersed distribution for plates of this size.
6. PSD retrieval sensitivity analysis
As discussed earlier, the retrieved PSD parameters can be quite sensitive to the radar observations, especially for small λ. This sensitivity arises due to the similar slopes of the and curves in (N0, λ) space (cf. Fig. 15). To gauge the impact of observational uncertainty on the PSD parameters, 20% perturbations in and were imposed on the observations for each case (Table 3). Case 1, which had the highest and lowest , showed the greatest sensitivity to changes in , whereas case 3 (lowest and highest ) had the smallest percentage changes. The main impact of these ZDRT perturbations was to change the relative contributions to from and . Decreases in increase both λ and N0 in order to provide higher contributions to from smaller-sized ice crystals. Perturbing (and thus the constraint) has the greatest impact for PSDs with small λ due to the exponential dependence of the PSD on λ. In contrast, perturbations in had more consistent impacts across the cases, because the constraint depends on both and , whereas (≈) is a direct constraint. Increases in lead to increases in N0 and λ, owing to the greater contribution of smaller particles needed to produce larger with the same . Negative perturbations in had the inverse effect (Table 3).
Additional uncertainty originates from the assumed value,3 which can have a significant impact on the constraint. To quantify this sensitivity, PSDs were determined for each case using = 0.0 and = 0.6 dB (Table 4). Again, sensitivities were largest for case 1, which had the largest , smallest , and thus largest initial . The sensitivity amplifies as decreases, implying aggregates increasingly dominate the backscatter. For case 1 ( = 1.0 dB), increasing to 0.6 dB causes a large decrease in , leading to large increases in λ and N0 for the PSD to maintain = 1.6° km−1. In contrast, very minor changes (<2%) occur for case 3 due to the small contribution to by aggregates.
PSD parameters are also affected by the type of ice crystals selected to represent sizes < 0.5 mm (Table 5). Using p2.0 plates (which have lower ) causes minor decreases in N0 and λ, because these particles have small impacts on , and especially . In contrast, using d0.5 (instead of d1.0) crystals for particles > 0.5 mm results in large decreases in both PSD parameters (Table 5). For a given number concentration, these thinner dendrites have lower values. Therefore, the PSD needed to satisfy the constraint requires a greater contribution from larger particles, which leads to decreased λ and N0.
Finally, assuming a gamma2 PSD instead of an exponential PSD also changes the amount each crystal of a given size contributes to the radar observables (Table 6). Generally, gamma2 PSDs had much lower contributions from the largest (>2 mm) and smallest (<0.5 mm) particles to and because of the narrower PSD shape. Consequently, 0.5–2.0-mm dendrites provide larger contributions to and to compensate.
In contrast to previous studies (Kennedy and Rutledge 2011; Bechini et al. 2013; Andrić et al. 2013), PSDs, along with the contributions to ZH from aggregates, that match the observations of ZDR, ZH, and KDP have been retrieved. However, the physical validity of these ice PSDs and their correspondence to observations must be established to consider them reasonable models of the ice crystal size spectrum for each case. A composite of several recent aircraft probe studies is presented in Heymsfield et al. (2013, hereinafter H13; their Fig. 5). Note that H13 included aggregates in the mean PSD, whereas the PSDs retrieved in this study represent only pristine ice crystals.
Figure 20 depicts the correspondence between retrieved PSDs and the H13 observations. Because the PSD retrievals showed the greatest sensitivity to dendrite aspect ratio, we show the retrieved PSDs for cases 1 and 2 using d1.0 (left panels) and d0.5 crystals (right panels). For case 1, the retrieved PSD using d1.0 crystals falls nearly entirely within 2 standard deviations of the mean observed PSD (Fig. 20a). The PSD retrieved using d0.5 crystals for case 1 has nearly zero slope, falling almost completely outside 2 standard deviations of the observations (Fig. 20b). This near-zero slope suggests that, for case 1, using thinner d0.5 crystals fails to retrieve a physically realistic PSD associated with the radar observations.
In contrast, both PSDs retrieved for case 2 have smaller differences with respect to the observed mean PSD. The d1.0 PSD is within 1 standard deviation of the observed PSD for sizes between 1.0 and 3.0 mm. Because of its larger slope, however, the retrieved PSD is more than 2 standard deviations below the mean PSD for crystals >3.5 mm. Some of this discrepancy may be attributable to the inclusion of aggregates in the mean PSD (H13), increasing the number of larger particles relative to the retrieved PSDs that only account for pristine ice crystals. The PSD obtained for case 2 using the d0.5 crystals is slightly more than 2 standard deviations above the mean PSD between 1.8 and 3.2 mm, and within 2 standard deviations elsewhere. Comparing the retrieved PSDs for case 2 using the d1.0 and d0.5 crystals suggests that using a thickness factor between 0.5 and 1.0 could produce a PSD that falls completely within 2 standard deviations of the H13 composite mean of observed PSDs.
Within a given sampling volume, the true population of hydrometeors may include more than just pristine dendrites and aggregates. Observed ZH, ZDR, and KDP therefore could be different from what would be observed with pristine dendrites and aggregates, introducing errors in the PSD retrievals. For example, light riming occurring within the dendritic growth zone will tend to drive the dendrites’ aspect ratio toward unity while simultaneously increasing their density. These changes in particle properties will modify the scattering properties of the crystals, although changes in ZDR from these effects tend to be of opposite sign. Heavy riming and the subsequent development of graupel will decrease ZDR toward zero for lump graupel or negative for conical graupel (e.g., Aydin and Seliga 1984; Oue et al. 2015). The increased particle aspect ratios will lower ZDR, decreasing the constraint, and KDP. Thus, N0 and λ will be positively biased, reflecting a PSD of dendrites instead of a more appropriate PSD of rimed dendrites and/or graupel. Any ice multiplication processes active at −15°C would likely increase the observed KDP, as a larger population of horizontally oriented ice crystals would be present. The sensitivity analysis above shows that this would decrease the slope and intercept of the retrieved dendrite PSD. Thus, the PSD retrieval technique works best when aggregation and vapor depositional growth are the dominant particle growth mechanisms, and temperatures and supersaturations are favorable for dendritic growth; erroneous PSDs may be retrieved when riming and ice multiplication are strongly contributing to the observed radar variables.
Despite these potential errors, the retrieved PSDs may be able to explain some of the variability in the radar observations. In case 3, the derived PSD implies a highly sloped exponential distribution, with ice crystals clustered within a narrow range of sizes. This suggests crystals in early stages of depositional growth, producing a narrow distribution of small crystals yet to exhibit secondary habit characteristics, and thus yet to exhibit any substantial degree of aggregation. Case 2 had a smaller PSD slope but a more prominent presence of aggregates. The retrieval suggests pristine plates and dendrites were present in sufficiently large concentrations to enhance the KDP, but aggregates contributed enough to the observed ZH to decrease the observed ZDR from that of pristine dendrites. In case 1, aggregates provide a much larger contribution to the observed ZH. The retrieved PSD is shallower than the previous two cases, suggesting that many of the smallest crystals had been collected by the larger aggregates. However, the concentration of larger dendrites (>1 mm) is much greater than the other cases, leading to strongly enhanced KDP. The convective appearance to the radar echoes in this case, along with the higher observed ZH, suggests a more mature stage of development in the cloud that could have allowed for vigorous growth of dendrites via vapor deposition and thus more efficient aggregation. As mentioned earlier, any riming occurring within this region of the storm could introduce more uncertainty in this PSD estimate.
This study presents an analysis of X-band, polarimetric radar observations collected during several winter storm events in northern Colorado, with particular focus on the signatures found within regions favored for dendritic growth. A large variability in the observed ZH, ZDR, and KDP was found, partially explained by the differing meteorological conditions occurring between and during each event. Periods of significant low-level upslope flow tended to be associated with more stratiform precipitation and, at times, substantial ZDR enhancements. During these periods of weaker upward motion, snow growth may be dominated by aggregation and vapor deposition. When stronger upward motion occurs and the cloud depth increases, perhaps in the presence of weak convective instability, riming, aggregation, and vapor deposition may all contribute to snow growth. These cases often exhibited the largest ZH and smallest ZDR values (indicative of large and/or numerous aggregates or rimed particles), yet large KDP enhancements, suggesting the rapid growth and/or sizeable population of dendritic crystals. The presence of plentiful, large dendrites increases the rate and efficiency of aggregation, and therefore could explain the observed correlation between KDP and ZH in regions favored for dendritic growth.
Ultimately, in situ measurements of ice particle types and size distributions coincident with radar observations within dendritic growth zones are needed to validate whether dendrites and aggregates are indeed responsible for the observed radar signatures. In the absence of such measurements, however, we can try to better quantify plausible descriptions of the ice particles and PSDs responsible for these radar signatures using scattering calculations. We performed such calculations using the sophisticated GMM technique that allows us to more realistically account for the complex structures of platelike and dendritic ice crystals, an improvement over simplified methods and particle shapes used in previous studies (e.g., Kennedy and Rutledge 2011; Andrić et al. 2013). PSD slope (λ) and intercept (N0) parameters for exponential and gamma2 distributions were retrieved from the radar observations for a variety of cases representing different stages of snow growth. The radar-retrieved PSDs are consistent with what is expected for these different stages based on microphysical arguments. That the observed ZDR is positively correlated with λ (Kennedy and Rutledge 2011) provides some evidence that these PSDs capture physically realistic processes such as aggregation and spectral broadening with vapor growth. However, the possible occurrence of riming during the convective case (case 1), not accounted for in the retrieval method, makes it difficult to draw more general conclusions.
The retrieved PSDs are most sensitive to the ice crystal aspect ratios used in the scattering calculations. During a period of stratiform snow, PSDs were retrieved from observations using thickness factors of 0.5 and 1.0 applied to an empirical size–dimensional relationship for dendrites. The results fell mainly within two standard deviations of the composite mean in situ PSD measurements presented in Heymsfield et al. (2013). A PSD retrieved using an intermediate thickness factor would likely reduce the differences with respect to the mean PSD over all sizes. Size–dimensional relationships that depend on environmental growth conditions (e.g., Harrington et al. 2013; Sulia et al. 2013, 2014) may reduce this uncertainty in the PSD retrievals. This also underscores the need for more in situ observations of ice crystals to better constrain axis ratio relationships for dendritic and platelike crystals.
When processes other than vapor depositional growth and aggregation occur within the dendritic growth zone, the retrieved PSD becomes less reliable. For example, riming corrupts the retrieval process by modifying the shape and density of existing dendrites and thus the observed ZH, ZDR, and KDP. Heavy riming will make the pristine ice crystals appear more spherical, thereby reducing ZDR and KDP. To correct for this process, additional quantitative information about the degree of riming is needed. This demonstrates the importance of considering all of the potential microphysical processes that may be occurring when attempting to analyze polarimetric radar data; neglecting to do so may lead to significant errors in interpreting radar observations of the dendritic growth zone.
Again, in situ measurements are needed to validate or refute the retrieval technique presented herein. Our results do, however, provide the first self-consistent set of polarimetric radar observations and plausible descriptions of PSDs in various stages of snow growth. The radar-retrieved PSDs are at least consistent with the behavior expected from microphysical arguments. If verified by in situ measurements, the relationships between observed ZH, ZDR, and KDP and retrieved PSDs can provide information about the relative efficiencies of dendritic growth and aggregation, and thus be useful in operationally identifying regions of locally heavy snowfall and ultimately improving numerical model parameterizations of snow growth processes.
The authors thank Patrick C. Kennedy, Steve Rutledge, and Francesc Junyent (CSU) for help coordinating, collecting, and processing the CSU-CHILL data used in this study. We thank the three anonymous reviewers for their comments that significantly improved this manuscript. Funding for this work (RS and MK) comes from NSF Grant AGS-1143948; YL is funded through NSF Grant AGS-1228180. We thank Hans Verlinde, Eugene Clothiaux, Kultegin Aydin, and Jerry Harrington (all at PSU) for useful discussions about this work, which is part of the first author’s M.S. thesis at The Pennsylvania State University.
Alternate KDP Estimation Procedure
In particularly noisy ΦDP regions, often where ZH is relatively low, the WC09 KDP estimation procedure sets these values to 0° km−1. Substantial ZDR enhancements due to the presence of pristine ice crystals at these locations may therefore be incorrectly associated with KDP = 0° km−1. To better quantify the KDP associated with these high-ZDR, low-ZH signatures, an alternate KDP estimation method is introduced. Because this procedure was applied in radar scans exclusively sampling ice particles, no folding of ΦDP occurred. Also, the backscatter differential phase shift was considered small given the high probability that most particles in these low-ZH echoes are electromagnetically small.
The first step in this method is removing data containing primarily noise. A gate was considered noise when it had a normalized coherent power (NCP) value <0.4. NCP is a measure of the predictability of the phase between pulses (Dixon and Hubbert 2012), where values near 1.0 are considered mostly signal and those near zero mostly noise. After masking the high-noise gates, a low-pass range filter was applied to the ΦDP field. In this study, we used the Lanczos filter (Duchon 1979). This filter is computed using the Fourier transform of a rectangular function in frequency space and a factor of “sigma,” (set to 1.0 in this study) which minimizes the Gibbs oscillations that occur with Fourier representations of discrete data. Once the filter is calculated using 2n + 1 gates (with n corresponding to the number of Fourier modes used to construct the filter), the resulting weights can be applied to the range profile of ΦDP. At each gate j, a weighted average of the ΦDP values from gate j − n to j + n is calculated using the weights determined by the filter. This filtered value of ΦDP at gate j can only be determined if none of the gates from j − n to j + n have been flagged as noise.
The low-pass cutoff (the width of the rectangular response function) used in this case is the frequency of oscillations in the KDP field estimated using the WC09 method, corresponding to approximately 1 peak every 53 gates (1 peak per 4.7- km range). Thirty-one gates (n = 15; 1.3- km range) were used to create the Lanczos filter weights. Gates close to the radar and near echo top with less than 15 ΦDP values at greater and/or lesser ranges were removed. A linear regression centered on each gate was used to calculate the slope of the filtered ΦDP profile. While this scheme decreases the overall number of KDP estimates, estimates are now provided for some of the gates previously flagged as noise (and thus set to zero) in the WC09 algorithm. The value of n = 15 is less than the number needed by Duchon (1979) to produce a unit response midway between zero and the cutoff frequency, leading to some contribution to the smoothed ΦDP profile from the higher-frequency range oscillations. However, the linear regression effectively smooths ΦDP further as its range derivative is computed, damping the higher-frequency oscillations somewhat.
The gates for which estimates of KDP were originally provided (i.e., regions of low noise) by the WC09 scheme are well correlated with those using this estimation procedure. A comparison of the vertical profiles of KDP (as in Fig. 11b) estimated using the alternate and WC09 methods is shown in Fig. A1. Below 2.8 km ARL, the median KDP values are nearly identical, with slightly narrower widths and similar shapes to the distributions seen in the new KDP estimation technique.
Because dendrites in the database are available only for maximum dimensions of >0.5 mm, plates are used for sizes between 0.1 and 0.5 mm. Plates are the primary habit preceding dendritic growth (e.g., Yokoyama 1993; Lamb and Verlinde 2011; Pruppacher and Klett 1997), so using plates for the smallest sizes is physically realistic. The exact size at which plates begin to show characteristics of dendrites depends on the environmental conditions.