## Abstract

A Bayesian optimal estimation retrieval is used to determine probability density functions of snow microphysical parameters from ground-based observations taken during four snowfall events in southern Ontario, Canada. The retrieved variables include the parameters of power laws describing particle mass and horizontally projected area. The results reveal nontrivial correlations between mass and area parameters that were not apparent in prior studies. The observations provide information mainly about the mass coefficient , somewhat less information about the mass exponent and the projected area coefficient , and minimal information about the projected area exponent . The expected values for retrieved mass power-law parameters = 0.003 28 and = 2.25 are consistent with those from several prior studies that looked at the mass of aggregate-like particles and precipitating ice aloft as functions of maximum particle dimension. Differences from other studies appear related to differences in the dimensions used to define particle size. The retrieval allows the analysis of relatively large volumes of continuous observations, greatly enhancing sampling relative to single-particle analyses. The retrieved properties are used to constrain 94-GHz (W band) radar scattering properties for a variety of snow particle shapes. Synthetic reflectivities calculated using these scattering properties and observed particle size distributions show that a branched, spatial aggregate-like particle produces good agreement with coincident observed W-band reflectivities. Uncertainties in the synthetic reflectivities, estimated by applying a simple error-propagation model, are substantial and are dominated by the uncertainties in and .

## 1. Introduction

Millimeter-wavelength (e.g., Ka and W band) radars have use for the remote sensing of precipitation, especially from space-based platforms owing to their sensitivity and relatively compact form (Mead et al. 1994). *CloudSat* (Stephens et al. 2008), for example, operates a nadir-viewing W-band profiling radar that observes latitudes from 82°N to 82°S, making it an important platform for sensing high-latitude precipitation. Scattering by precipitation-sized ice particles at these wavelengths, however, is sensitive to particle shape as well as mass, and accounting for the variation of snow microphysical properties and shapes is a substantial difficulty for retrieving snowfall with these instruments. Such retrieval problems require a priori assumptions about the microphysical and scattering properties of observed snowfall. To adequately characterize uncertainties in snowfall retrievals, the a priori information must characterize the variability of these properties as well as their expected values.

For snowfall at the surface, microphysical properties have typically been evaluated using time-consuming analyses of individual particles (e.g., Mitchell et al. 1990; Heymsfield and Westbrook 2010). The limited sampling in such an approach makes it difficult to characterize the expected values and variability of these properties. As an alternative, several recent field campaigns, including ground validation studies associated with satellite missions (Hudak et al. 2006; L’Ecuyer et al. 2010; Hudak et al. 2012) as well as a stand-alone study (Löhnert et al. 2011), have made intensive observations of snowfall properties using nearly collocated radars, disdrometers, and precipitation gauges. Although the locations sampled by these campaigns are currently limited, such observations offer the potential to characterize snowfall properties in varied locations and meteorological regimes.

Wood et al. (2014, hereinafter W14) introduced an optimal estimation (OE; Rodgers 2000) retrieval for analyzing such observations. The retrieval produces Gaussian probability density functions (PDFs) representing the expected values and uncertainties for parameters describing the variation of particle mass *m* and horizontally projected area with particle size. Such PDFs have value as explicit a priori constraints for remote sensing retrievals, and may also be used to characterize uncertainties in scattering properties caused by uncertainties in microphysical properties. Determining mass and parameters simultaneously ensures their physical consistency, and the automation of most of the relevant measurements ensures substantial sampling during snowfall events.

This work characterizes snow microphysical properties via the W14 retrieval, and then applies those properties to investigate the modeling of radar scattering properties of snow particles. Observations of snowfall rate, snow particle size distribution, size-resolved fall speeds, and 9.35-GHz (X band) radar reflectivity are used to estimate the parameters of power laws describing particle mass and as functions of particle size. The retrieval is applied to four midlatitude snowfall events from the Canadian *CloudSat*/*CALIPSO* Validation Project (C3VP; Hudak et al. 2006). The retrieved PDFs of snow microphysical parameters are then applied to constrain particle models used to evaluate W-band radar scattering properties via the discrete dipole approximation (Draine and Flatau 1994), and the performance of these models is evaluated using ground validation observations.

## 2. C3VP events

The snow events used in this work were observed at the Meteorological Service of Canada’s Centre for Atmospheric Research Experiments (CARE) at Egbert, Ontario, Canada, approximately 80 km north of Toronto. Observations used as inputs to the retrieval consist of 9.35-GHz radar reflectivities from the McGill University vertically pointing X-band radar (VertiX; Fabry and Zawadzki 1995), snowfall rates from a Vaisala, Inc., FD12P (Vaisala Oyj 2002) scaled to match snow accumulations from a double fence intercomparison reference (DFIR; Goodison et al. 1998), size-resolved fall speeds from Colorado State University’s 2D video disdrometer (2DVD; Kruger and Krajewski 2002; Huang et al. 2010), and size distributions from the National Aeronautics and Space Administration’s Snowflake Video Imager (SVI; Newman et al. 2009). The characteristics of these observations are described more completely in W14.

Four snowfall events were selected because of completeness of the required observations and because they represent a modest range of snowfall conditions. Observations by personnel on the ground at CARE (R. T. Austin et al. 2007, unpublished manuscript) along with daily operations logs (CIRA 2013) from the 10th Cloud Layer Experiment (CLEX-10), which operated jointly with C3VP, provide characteristics of three of the events. Synoptic event 1 (SYN1; 6 December 2006) involved a weak low passing northeastward over Ontario that produced snowfall at CARE mainly between about 1200 and 1530 UTC. Aircraft observations near CARE showed liquid phase near cloud top with mixed phase and ice below. Snowfall at CARE was described as light and dry early in the day transitioning to moderate wet snow later. VertiX echo-top heights were about 4 km above ground level (AGL) during the precipitation and SVI size distributions showed tails extending to 4–8 mm (Fig. 1). Temperatures during the most significant snowfall were near freezing. The 24-h accumulation for the event was 3.2 mm liquid water equivalent (LWE).

Lake-effect event 1 (LE1; 7 December 2006) consisted of lake-effect snow squalls that resulted from the cold air mass and northwesterly winds that followed the system of the previous day. CARE received snowfall over most of the day, with a 24-h accumulation of 10.2 mm LWE. Temperatures were near freezing early in the day and decreased with time, reaching 255 K at the day’s end (Fig. 2). VertiX echo-top heights were shallower than the previous day, varying from 1 to 3 km AGL. SVI size distributions were similar to those of the previous day but more variable over time. A period after 2100 UTC showed high concentrations of particles smaller than 2 mm and was associated with the coldest temperatures of the day.

Lake-effect event 2 (LE2; 27–28 January 2007) was a second lake-effect snow event that resulted as a warm front near CARE shifted to the south during the evening of 27 January and cold northwesterly winds entered the area. Snow fell mainly between 0100 and 0400 UTC, at temperatures between 267 and 270 K. Snowfall rates at CARE were initially light, but increased rapidly as a heavy snowband lingered over the site (Fig. 3). Large snowflakes, near 10 mm in diameter, fell during periods of heavy snow and visibility was near zero. SVI size distributions showed particles with sizes up to 10 mm early in the event. Total accumulations for the day were 4.6 mm LWE.

Synoptic event 2 (SYN2; 14 February 2007) occurred between intensive observing periods. Although C3VP observer reports of the conditions on the surface and aloft are lacking, the event has been recognized as a massive synoptic winter storm, producing extensive snowfall with substantial accumulations and societal impacts over northeastern United States and southeastern Canada (Grumm and Stuart 2007). At CARE, this system was significantly deeper than the other three events, with VertiX echo-top heights extending to about 6 km AGL (Fig. 4). Observations from precipitation gauges show that snowfall occurred throughout most of the day and produced accumulations of 8.3 mm LWE. This event was also the coldest of the four, with temperatures ranging from 256 to 261 K during the snowfall and with SVI size distributions that were narrow in comparison with those of the earlier events. Characteristics of the events are summarized in Table 1.

The causes for the large DFIR/FD12P accumulation ratios for events LE2 and SYN2 (Table 1) are not clear. The FD12P uses both an optical sensor and a heated capacitive sensor to estimate precipitation rates for snow and the capacitive sensor is subject to undercatch (Vaisala Oyj 2002). For event LE2, however, winds observed at the CARE meteorological tower were light at 0.5–2.5 m s^{−1} during the snowfall, suggesting that undercatch should not have been significant. Winds were stronger for event SYN2, generally below 5 m s^{−1}, but ranging as high as 6 m s^{−1} for short periods. Additionally, the size distributions for event SYN2 suggest high concentrations of small particles were present (Fig. 4). For this event, undercatch by the capacitive sensor may have been significant and biased the unscaled FD12P precipitation rates low.

Observations were averaged over independent 5-min samples and uncertainties estimated as described in Wood et al. (2013) and W14. Ground clutter caused the nearest usable VertiX reflectivities to be 488 m above the ground. From the reflectivities at 488 m, reflectivities at the surface and the corresponding uncertainties were estimated by considering vertical reflectivity gradients and the likely time separation between precipitation features observed aloft and their appearance at the surface, as described by W14. The resulting samples were then screened to ensure temperatures less than 273 K and wind speeds less than 5 m s^{−1} at the surface, giving 375 distinct samples from the four events: 33 from event SYN1, 173 from LE1, 43 from LE2, and 126 from SYN2.

The VertiX observations are susceptible to attenuation by wet snow accumulations on the radar’s conically shaped radome. Wet snowfall may have occurred for event SYN1 from roughly 1200 to 1500 UTC, and for event LE1 from 0000 to 0400 UTC, during periods of near-freezing temperatures. The VertiX was monitored and cleared of snow periodically, but the times at which this was done are not known. Comparisons of time series of observations over the CARE site by Environment Canada’s King City radar against coincident VertiX observations suggest the VertiX experienced about 2 dB of attenuation by the end of the SYN1 period. For the LE1 event, the time series comparisons suggest about 4 dB of attenuation, but it appears the radome was cleared after 0400 UTC, prior to the remaining 20 h of LE1 observations. Bias corrections for the VertiX observations were determined using these King City radar comparisons, as described in W14, and provide partial compensation for this attenuation.

## 3. The snow microphysics retrieval

The retrieval, described and evaluated by W14 using synthetic observations, assumes that particle mass *m* and horizontally projected area can be related to particle size using power laws (e.g., Mitchell 1996):

where is the particle maximum dimension, the distance between the two farthest-separated points on the particle’s surface, and is simply the area of the projection of the particle’s shape onto the plane normal to the particle’s downward motion. For a given sample of the observations, the retrieval gives the joint distribution of , , , and , as well as one additional variable . This additional variable relates the particle dimension observed by a particular disdrometer to , compensating for the varying ability of disdrometers to measure the true dimensions of particles of different shapes (Wood et al. 2013):

with the measured size distributions transformed as

A single value of is applied over all particle sizes in any observed distribution. Together, these variables form the state vector , whose PDF is to be retrieved,

where and have been log transformed because of their large range and the need for their distributions to approximate the Gaussian shape.

To find the retrieved estimate of the state, the cost function

is minimized with respect to using Newtonian iteration (Rodgers 2000). The observation vector contains 9.35-GHz radar reflectivity, snowfall rate, and three measures derived from fall speeds observed in size bins centered at values of 4, 2, and 1 mm (W14). The forward model approximates the true physical relation between and , and may require other influence parameters , where the tilde indicates these parameters may be known imperfectly. In particular, this forward model, described more completely in W14, takes trial values of the elements of the state vector , then uses the observed particle size distribution (PSD) to simulate X-band reflectivity, snowfall rate, and the three fall speed measures for comparison against the observed values in . The covariance matrix represents the combined measurement and forward model uncertainties as a multivariate Gaussian distribution. Prior information about the state vector is represented with the vector of expected values and their covariances (Table 2), also as a multivariate Gaussian distribution. W14 provides an assessment of the a priori state and uncertainties for the retrieval. At convergence, a statistic is calculated by applying in (6), and a value near , the number of observations, suggests correct convergence.

Selecting Gaussian distributions acknowledges the limited information available about the actual forms of these distributions, and also leads to a reasonably tractable form for the retrieval. For distributions with specified widths, Gaussian distributions maximize the entropy of the distribution (Rodgers 2000; Shannon and Weaver 1949) and impose the least constraint on the retrieval. Selecting other PSD forms without evidence of their appropriateness would introduce unjustified constraints.

The retrieval also estimates the error covariance matrix for the retrieved state vector ,

where only the upper triangular portion is shown, are variances in retrieved parameters, and are the covariances between them.

Together, and describe a five-dimensional Gaussian PDF that represents the retrieved state. Values for , , , and in this work are in cgs units, with in centimeters and and taking appropriate units to give mass in grams and in square centimeters. The results for each event were composited by Monte Carlo sampling from each retrieved PDF within the event, then evaluating the mean values and covariance matrix for the pooled sample. Ten realizations were performed and the resulting means and covariance matrices averaged to obtain final values. Retrievals were performed in two different configurations, representing two different fall speed forward models [MH05, from Mitchell and Heymsfield (2005), and HW10, from Heymsfield and Westbrook (2010), as described in W14]. The results presented here are for the MH05 configuration except as noted, but results for the HW10 configuration were similar.

## 4. Snow microphysics retrieval results

Considering the mass parameters, retrieved values of range mainly from 1.75 to 2.75 and from about 0.001 to 0.01, with and positively correlated (Fig. 5). The pooled means and standard deviations, taken as the square roots of diagonal elements of the pooled covariance matrices, are indicated by the blue error bars. The black triangles in Fig. 5a show parameter values compiled by Mitchell (1996) for a range of particle shapes (Table 3). The retrieved values from this study are similar to those from Mitchell for several shapes associated with large, irregularly shaped particles: densely rimed dendrites (R2b), aggregates (S3 and S1a), rosettes (C2a), and side planes (S1), where the terms in parentheses are the corresponding habit designations from Magono and Lee (1966). For all events, the pooled standard deviations are substantially smaller than the a priori standard deviations (gray bars), suggesting the mass parameters are well constrained by the observations, particularly .

Differences between events are apparent in the mass parameters, with values for the two synoptic events (SYN1 and SYN2) larger than those for the lake-effect events (LE1 and LE2). The strongest contrasts are between events SYN1 and LE2. The mass parameter pooled means for these two events differ by almost a full standard deviation in both and , perhaps capturing differences in the meteorological conditions producing the snowfall. In spite of the differences, the reflectivity values for the two events are similar, as indicated by the sizes of the plotted points in Fig. 5. In contrast, while event SYN2 was also a synoptic-scale event and has parameters similar to event SYN1, it is characterized by much smaller reflectivities and temperatures than is event SYN1. The narrow PSDs observed for event SYN2 (Fig. 4) suggest weak aggregation, in contrast to the broader observed PSDs and the large, sticky particles reported by ground personnel at CARE for event SYN1; however, an examination of the SVI images shows substantial numbers of irregular, aggregate-like particles with *D* = 1 to 5 mm for event SYN2. In contrast to the other events, event LE1 has a more mixed range of temperatures and reflectivities, and has retrieved parameters falling between those of the two synoptic events SYN1 and SYN2 and the other lake-effect snow event LE2.

Unlike the mass parameters, the area parameters and are similar among the events and differ only slightly from the a priori estimates (Fig. 6). The pooled standard deviations are also similar to the a priori values. Overall, these results suggest the observations provide some information about the area parameters but constrain them only weakly. Again, values from Mitchell (1996) that represent larger, irregular particles are shown for comparison (Fig. 6a, black triangles). The two points nearest the expected values for the retrieved state correspond to densely rimed dendrites (R2b), aggregates (S3 and S1a), rosettes (C2a), and side planes (S1) (Table 3), the same shapes that best matched the results for the mass parameters.

For , smaller values are associated with narrower distributions (Fig. 7), where values in this plot were obtained using weighted orthogonal distance regression (Boggs et al. 1992) to fit the observed SVI size distributions to the log transformed form of the negative exponential size distribution, with weights based on uncertainties calculated following Wood et al. (2013). There is also a trend toward smaller values of with decreasing temperatures, with the smallest values of associated with the coldest event, SYN2. The pooled standard deviations are somewhat smaller than the a priori value, suggesting the observations weakly constrain this parameter. The pooled mean values are smaller than the a priori and vary somewhat among the events. For the SVI, reductions in are associated with reductions in particle aspect ratios (Wood et al. 2013). The trend toward smaller with narrower distributions and colder temperatures would be consistent with a trend from aggregate to more pristine particles, but the differences among the a priori and event-specific values are small compared to the standard deviations.

A pooled estimate was also produced for all four events combined, following the process described above. Because the number of retrievals per event varied considerably, weights were applied to the sampled points to balance the importance of the events, giving as results

with the covariance matrix

As anticipated from Figs. 5–7, the variances and covariances associated with the mass parameters have decreased compared to the a priori, while those associated with the area parameters and have changed little (Table 4). The posterior state for the HW10 configuration is essentially unchanged from that for the MH05 configuration. Values for the diagonal elements of the averaging kernel matrix , a performance metric described in W14 for this retrieval, show that is strongly constrained by the observations ( values from 0.9 to 1.1), and are influenced by both the observations and the a priori ( values of about 0.3), and and primarily reproduce the a priori ( values near 0.1).

One of the key questions to be answered by this analysis is whether covariances exist that were not present in the a priori. From (9), it can be seen that covariances that mix mass and area parameters {e.g., } are nonzero, whereas they were zero in the a priori. The corresponding correlation matrix (Fig. 8) shows nontrivial correlations between and (=0.20), between and (=0.09), and between and (=0.14). The first two, in particular, likely arise because of the dependence of fall speed on the ratio of mass to (Mitchell and Heymsfield 2005).

While the expected values of and from this study are similar to estimates from Mitchell, they are considerably different from the often-used values described by Brown and Francis (1995) ( = 1.9, = 0.002 94 in cgs units), which were taken from the results of Locatelli and Hobbs (1974). Locatelli and Hobbs reported these values for unrimed aggregates of bullets, columns and side planes, and also for aggregates of densely rimed dendrites or radiating assemblages of dendrites. Heymsfield et al. (2010), using aircraft observations from six field campaigns, have suggested that values of near 2.1 are more appropriate for cloud and precipitating ice. They found corresponding values of (cgs units) of 0.003 59 for warm-topped nonconvective clouds, 0.005 74 for cold-topped nonconvective clouds, and 0.006 30 for convectively generated clouds. The warm-topped cloud cases included observations from C3VP. The mass–dimension relations from Mitchell (1996) and Heymsfield et al. (2010) are all based on measurements of maximum particle dimension, while Locatelli and Hobbs (1974) used the diameter of an equal-area circle.

Several other studies have suggested values of that are substantially larger than the value of 0.003 11 from this study. Brandes et al. (2007) found = 0.008 90 and = 2.1 for snow along Colorado’s Front Range; however, their particle size was an equivalent volume diameter obtained from 2DVD observations and likely substantially smaller than . Further, their mass–dimension relation was given as a function of the size distribution median volume diameter, rather than actual particle size. Muramoto et al. (1995) presented a density relation that, based on their definition of particle volume, can be converted to a mass–dimension relation with = 0.009 87 and = 2.594. For particle size, they used the width of the particle as observed by a side-viewing camera with no correction for viewing geometry, which also would underestimate (Wood et al. 2013). Similar to Brandes et al. (2007), they presented their density relation as a function of the mean particle size of the observations, rather than as a function of actual particle size. Magono and Nakamura (1965) gave a density relation for observations of wet and dry snow. Using only their dry snow observations and converting their densities to masses using their definition of particle volume, a best-fit mass–dimension relation gives = 0.009 07 and = 1.82. Constraining to 2.25 gives = 0.007 22. They collected the snow particles on a flat surface, measured the longest horizontal dimension and the dimension normal to that, and then used the geometric mean of those dimensions as the particle size. Again, this underestimates .

These results suggest that differences in how particle size is measured can have a substantial impact on estimates of . Using as defined in (3) and letting be the value of determined when particle size is given by , the value associated with some other measure of particle size can be evaluated using

Estimating as 0.8 and using the expected value of of 2.25, will be a factor of 1.7 larger than . A significant part of the differences in mass–dimension coefficients may therefore be explained by these differences in the treatment of particle size. Differences may also be due to the use of independent variables other than (distribution median volume diameter or mean particle size).

## 5. Application to radar retrieval of snowfall

A primary goal of the preceding analysis is to provide observational constraints on mass and shape that can be used for modeling snow particle scattering properties and their uncertainties for use in millimeter-wavelength radar retrievals of snowfall rate. Using these observational constraints helps insure a physically consistent relationship between the particle microphysical and scattering properties that influence the relationships between radar reflectivity and snowfall rate that are assumed implicitly or explicitly within such retrievals.

Because of shape sensitivity, the Rayleigh or Mie approaches commonly used to model scattering properties at lower frequencies are not applicable for other than small particles (Schneider and Stephens 1995; Liu and Illingworth 1997). Scattering properties of larger snow particles may be simulated with Mie theory using spheres composed of a mixture of air and ice (the soft sphere approximation); however, comparisons against less approximate methods have shown the inability of such models to reproduce backscattering properties across multiple frequencies (Petty and Huang 2010). Several more complex techniques may be used at millimeter wavelengths (see, e.g., Bohren and Singham 1991 for an overview). The discrete dipole approximation (DDA; Draine 1988; Draine and Flatau 1994), the method used for this work, allows for arbitrary geometry by replacing the continuous particle with an array of discrete dipoles on a spatial lattice.

The relationships for particle mass and given by (8) do not fully constrain particle shape, so additional assumptions about shape are required before DDA calculations can be made. For millimeter wavelengths, it is likely that the scattering properties are not strongly sensitive to the fine structure of the spatial distribution of mass but will be sensitive to the larger-scale structure (Matrosov 2007). The objective of an assumed particle shape should be to reasonably capture the gross features of the spatial distribution of mass, rather than replicate particular fine features. The method for constraining DDA particle models using an assumed shape and the retrieved mass and relationships is described in the appendix. Here we describe the various shape assumptions examined in this work, then present and assess the DDA results.

### a. Shape assumptions

Using constrained DDA particle models ( appendix), scattering calculations were made for a variety of shapes intended to represent the features of more pristine, planar particles and of larger, irregular aggregates. For planar particles, a branched platelike shape with six branches, designated as shape SPp, was used (Fig. 9, upper left). With this planar shape, the horizontally projected area may be altered by changing the width of the branches. For narrow widths, the shape is much like a stellar crystal (P1d; Magono and Lee 1966) and has a small horizontally projected area. As the branch width increases, the shape approaches that of a crystal with broad branches (P1c) and at the limit of maximum branch width, is a hexagonal plate (P1a). For purposes of comparison, calculations were also done for hexagonal plates (shape HPp) that met the mass constraints from the snow microphysics retrieval but not the constraints on horizontally projected area.

Spatial particles were represented with clusters of thick hexagonal branches and with scalene ellipsoids (Fig. 9). These spatial shapes are considered to be simplistic, somewhat abstract representations of aggregate particles. Three different configurations were examined: B6pf, with six branches all lying in the horizontal plane (Fig. 9, upper middle); B8pr-30, with eight branches, six of which intersected the horizontal plane at angles of 30° (Fig. 9, lower left); and B8pf-45 with eight branches, six of which intersected the horizontal plane at angles of 45° (Fig. 9, lower right). The orientation of the branches controls the aspect ratio of the particle and for a given aspect ratio, the horizontally projected area may be altered by changing the branch thickness. For shape B6pf, both aspect ratio and horizontally projected area get smaller as the branch thickness is reduced. Shape B8pr-30 has an aspect ratio of about 0.5, and B8pr-45 has an aspect ratio near 0.70.

The scalene ellipsoid Ep (Fig. 9, upper right) is a simple shape that also has the ability to meet the constraints on horizontally projected area and aspect ratio, but whose shape is fundamentally different than the branched particles. Given a desired size , the length of the major horizontal axis is set to this size and the length of the vertical axis is set to . This choice provides the same aspect ratio as the B8pr-30 shape. The length of the minor horizontal axis is adjusted to match the required horizontally projected area. In the event the required length is less than , the minor horizontal axis is set to and the required area is achieved by placing continuous porosities at random locations extending vertically through the particle.

### b. DDA results

Scattering calculations were performed for particle sizes from = 0.025 mm to = 18 mm using the discrete dipole approximation scattering software (DDSCAT), version 7.1 (Draine and Flatau 2010). Particles were assumed to be oriented randomly with their longest dimension lying nominally in the horizontal plane (i.e., rotations in the horizontal plane were sampled uniformly), and were illuminated with a vertically incident, linearly polarized plane wave. Canting angles were applied, having values ranging over ±10° but sampled uniformly in the cosine of the angle following DDSCAT’s standard method. This approach is an approximation based roughly on the canting angle distributions found by Matrosov et al. (2005b) for pristine dendritic particles. DDSCAT provides the differential backscatter cross section normalized by :

where is solid angle, indicates the derivative is evaluated in the backscattering direction, and is the equivalent volume radius,

with the solid ice density of 0.917 g cm^{−3}. For simulating radar reflectivities, the backscatter cross section is calculated from as

Figure 10 shows the resulting for the planar SPp shape. For comparison, for Rayleigh spheres, Mie spheres, and the HPp shape are shown also. Since the DDA dipoles are taken to be solid ice, as are the Rayleigh and Mie spheres, particles with the same contain the same mass. For below about 0.7 mm, for the SPp and HPp shapes are similar. At these small sizes, the constraint on horizontally projected area gives large area ratios, causing the SPp particles to be similar in shape to the HPp particles.

In the Rayleigh size range, for these platelike particles exceed those for spheres, consistent with predictions by models for Rayleigh backscatter by oblate spheroids (Atlas et al. 1953). At the point where the HPp cross sections fall below those for Rayleigh spheres, the size parameter based on ,

where is the radar wavelength, has a value of about 1.4. At larger than about 0.7 mm, for the HPp shape fall below those for the SPp shape since HPp particles have the same mass as the SPp particles but larger horizontally projected areas. Consequently, HPp particles are thinner than the SPp, while the SPp are thicker and have much of their mass concentrated near their centers. Additionally, for larger than about 0.73 mm, the HPp particles have an insufficient number of dipoles to ensure that the entire plate area is occupied by dipoles. As a result, there are porosities extending through the plate. These through-porosities occupy 20%–30% of the plate area. These factors appear to be sufficient to cause the HPp cross sections to fall below those of the SPp particles. Neither shape shows evidence of the resonance at = 0.8 mm that is apparent in the results for the Mie spheres.

Figure 11 shows for the spatial particles. As before, values for Rayleigh and Mie spheres are shown for comparison. The for the more compact shape, B6pf, are similar to those for spheres for sizes up to = 0.35 mm. For the more spatially extended particles B8pr-30, B8pr-45, and Ep, falls below those for Rayleigh spheres at = 0.15 to 0.2 mm, corresponding to size parameters of 0.3–0.4. For the branched particles, as the aspect ratio of the particle increases, decreases, although for some sizes values for B8pr-30 and B8pr-45 are almost equal. The B8pr-30 and B8pr-45 shapes show Mie-like resonances, albeit with much smaller amplitude than the Mie sphere resonance at = 0.8 mm, while the B6pf shape shows none. The Ep shape shows markedly smaller than the branched particles over most of the size range, and exhibits strong resonance features.

Over most of the shown size range, for the two extremes of the branched particles (B6pf versus B8pr-45) differ by at least an order of magnitude, indicating that additional information is needed to adequately constrain the scattering properties for these spatial particles. Korolev and Isaac (2003) suggest that, for particles smaller than = 1.0 mm, aspect ratios should be no larger than 0.6–0.8. Magono and Nakamura (1965), using photographs of snow particles taken in elevation view, found that for particles with < 10 mm, aspect ratios were near 1.0, and that for larger particles, the horizontal dimension was substantially larger than the vertical dimension. Matrosov et al. (2005a) found that scattering models based on particles with aspect ratios of 0.6 gave better agreement with aircraft-observed dual-frequency radar ratios than did models based on particles with aspect ratios of 1.0. These results suggest that the B8pr-30 or B8pr-45 shapes should be more representative of the backscattering properties of true snow particles, especially in larger sizes, than is the B6pf shape.

### c. Assessments

The modeled scattering properties were assessed using observed PSDs and collocated W-band radar reflectivities. PSDs were obtained from the SVI, which operated nearly continuously at CARE during the 2006/07 C3VP observing season. Observations at 1-min resolution were averaged using distinct 5-min samples as was done for the snow microphysics retrievals. These distributions, based on the Feret diameter, were converted to distributions on using = 0.78, the value obtained from the snow microphysics retrieval.

The observed radar reflectivities were provided by the airborne cloud radar (ACR) (Sadowy 1999), a 95-GHz profiling radar deployed on the ground at CARE during C3VP. The ACR pointed vertically, and the range bin nearest the surface was centered at 197 m AGL. Comparisons of the reflectivities in this bin versus reflectivities in the adjacent bin above suggest this lowest bin was not substantially affected by ground clutter for reflectivities above about −15 dB*Z*_{e}, so observations from this lowest bin were used and assessments were limited to cases with reflectivity greater than −15 dB*Z*_{e}. Reflectivities in linear units at about 2.8-s resolution were averaged in time using 5-min samples, consistent with the treatment of the SVI observations. Although a formal calibration of the ACR was not performed immediately prior to C3VP, a previous intercomparison between the ACR and the University of Massachusetts Cloud Profiling Radar System showed average differences of 0.3 dB*Z*e (Sekelsky et al. 1999). Here, calibration errors are neglected, and it is noted that any biases in the ACR calibration could affect the results presented below.

#### 1) Reflectivity comparisons

Data were screened for air temperatures colder than 273 K as observed on the 10-m meteorology tower at CARE. After matching observations, valid cases were obtained for twelve distinct snow events occurring over 13 days from 2 December 2006 to 26 February 2007. Attenuation in the observed ACR reflectivities is expected to be negligible because of the short range of 197 m to the selected range bin and because the radome was regularly cleared of accumulated snow by operators attending the radar during snow events (R. T. Austin et al. 2007, unpublished manuscript). For these comparisons, then, reflectivities were simulated from the modeled scattering properties assuming negligible attenuation. The unattenuated equivalent reflectivity factor for non-Rayleigh scatterers is given by

where , and is the complex refractive index of liquid water. DDA backscatter efficiencies (13) were interpolated to the SVI sizes, then converted to backscatter cross sections that were used with the SVI PSDs in a discrete form of (15) to obtain simulated reflectivities.

Of the five shapes considered, the B8pr-30 shape provided the best agreement with the observed reflectivities, with a bias of −0.03 dB and root-mean-square (RMS) difference of 5.4 dB over all of the PSDs that were well sampled by the disdrometer (Fig. 12). The more compact shapes B6pf and SPp substantially overestimated reflectivities over most of the observed reflectivity range and produced considerably greater scatter versus the ACR observations. The less compact B8pr-45 and Ep underestimated reflectivities, with no improvement of RMS differences versus that for the B8pr-30 shape. Obviously, the use of the SPp shape over the full size range observed by the SVI represents an unrealistic and severe extrapolation. The largest observed for this type of pristine shape (e.g., branched, stellar, or dendritic crystals) are typically a few millimeters (Heymsfield and Kajikawa 1987; Mitchell 1996), while the particle models extend to = 18 mm. Nonetheless, the result serves to illustrate the magnitude of errors that may be caused by the ill-considered use of such models.

The difference in bias between the B8pr-30 and the B8pr-45 shapes corresponds with differences in the vertical aspect ratios: the vertical aspect ratio for the B8pr-30 particle is about 0.5 while that for the B8pr-45 is near 0.7, giving a particle that is more extended along the direction of radar beam propagation than is the B8pr-30 shape. Additionally, for a given particle size, the branches of the B8pr-45 particle are likely somewhat wider than those of the B8pr-30 particle. This increase in width is necessary for the different shapes to have equal . Wider branches would cause the B8pr-45 particle to have somewhat larger volume than the B8pr-30 particle. Since for a given size, the particles have the same mass, the dipoles in the B8pr-45 shape will not be as closely packed as in the B8pr-30 shape. In contrast, the Ep particle has the same vertical aspect ratio as the B8pr-30, but for most the volume of the Ep particle substantially exceeds that of the B8pr-30, leading to a much less dense arrangement of dipoles. The resulting reflectivity bias for the Ep particle is similar to that of the B8pr-45, but the backscatter cross sections are quite different for the two shapes (Fig. 11).

Clearly, combinations of the shapes used in this study could be found that give biases similar to that of the best-fit shape. For example, one can imagine a combination of the highly reflective SPp shape at small sizes with the less reflective B8pr-45 shape at large sizes. Since these shapes are constructed using the same and relations, any such mixed-shape cases would give the same snowfall rates and snow water contents as the single-shape cases. In this scenario, the reflectivity bias of any mixed-shape case will simply fall between the reflectivity biases of each of the component single shapes.

The comparisons shown in Fig. 12 also highlight the possible effects of the limited sample volume of the SVI compared to the radar. The blue points are cases for which the SVI detected fewer than 100 particles over the 5-min sample with those particles distributed in five or fewer size bins, suggesting the size distributions may have been poorly sampled. An examination of the ACR operator’s log (R. T. Austin et al. 2007, unpublished manuscript) showed that many of these cases were associated with the initiation or termination of snowfall at the surface, or with low-level stratocumulus lacking precipitation.

#### 2) Uncertainties in modeled reflectivity

Random placement of dipoles on the lattice ( appendix) causes small variations in the scattering properties. To evaluate these variations, four distinct realizations of the random dipole arrays for the B8pr-30 shape were constructed and used to calculate . At small particle sizes, the were largely insensitive to the dipole locations as shown by the small fractional uncertainties (Fig. 13), consistent with the expectation that scattering properties are primarily sensitive to the mass or equivalent volume diameter in the Rayleigh regime. At larger sizes, fractional uncertainties were about 0.05, with values for particular sizes as large as 0.15 and as small as 0.015.

To evaluate the influence of random dipole locations on modeled reflectivities, the replicate realizations of the B8pr-30 scattering properties were used along with the SVI PSDs to calculate radar reflectivity per (15). Evaluations were limited to cases for which the SVI size distribution was well sampled (more than 100 particles observed or particles distributed over more than five size bins). As expected, since these random variations are uncorrelated over the range of particle sizes, the resulting uncertainties in the modeled radar reflectivity are negligible (Fig. 14), with typical fractional uncertainties in *Z*_{e} of 0.02.

The covariance matrix in (9) represents uncertainties in the microphysical parameters. These uncertainties contribute to uncertainties in scattering properties and thence to uncertainties in modeled reflectivities. For retrieval schemes using millimeter-wavelength radar observations, knowledge of these reflectivity uncertainties is critical. The reflectivity error variance can be estimated using simple uncertainty propagation as

where is the Jacobian of the modeled reflectivity with respect to the microphysical parameters

Evaluating Jacobian terms numerically requires that reflectivities be calculated separately with perturbed and unperturbed scattering properties. Perturbed scattering properties are obtained by perturbing each microphysical parameter then constructing new dipole arrays for each particle size, from which DDA calculations provide the perturbed scattering properties. Perturbations to the microphysical parameters must be large enough that the resulting change in is distinguishable from the uncertainty resulting from random dipole placement. If is the fractional uncertainty in , a reasonable approach (Dennis and Schnabel 1983) is to take

where is the parameter of interest. Normally is determined by the precision of the calculation of ; however, the use of random dipole locations for the DDA calculations of the scattering properties introduces additional uncertainties. From the results above, the fractional uncertainty for in linear units is typically less than 0.02, suggesting that , and perturbations of 15% were used.

Jacobians were evaluated using the PSDs from the SVI dataset. Reflectivity increases with increasing but decreases with increasing (Figs. 15a,b). Increasing causes mass to increase over the entire size distribution, while increasing causes masses to decrease for < 1 cm and increase for > 1 cm,

Negative values for indicate that the reflectivities are dominated by contributions from particles smaller than 1 cm. As particle size distributions broaden (indicated by decreasing in Fig. 15), the influence of large particles becomes more significant and the derivative becomes less negative. Conversely, reflectivity decreases with increasing but increases with increasing (Figs. 15c,d). Increasing causes horizontally projected areas to increase over the entire size distribution, while increasing causes to decrease for < 1 cm and increase for > 1 cm, parallel to the behavior for . For a given mass and shape, increasing produces a less compact particle, which tends to decrease reflectivity.

Uncertainties in modeled reflectivities range from 5 to 15 dB and also depend markedly on the slope of the size distribution, with narrow distributions having much smaller uncertainties than broad distributions (Fig. 16), consistent with the behavior of the Jacobian for and . It is potentially useful to know how significantly the uncertainties in individual parameters , , , and contribute to uncertainties in the reflectivity. In the evaluation of via (16), however, the contribution due to uncertainty in a single parameter cannot be isolated because of the presence of covariances in . In an approximate sense, estimates of these contributions can be obtained from the products of the derivatives and the parameter uncertainties (Table 5). Uncertainties for the parameters were estimated as the square roots of the variances shown in (9) and derivatives were estimated as the simple means of the values shown in Fig. 15. The results suggest that uncertainties in reflectivity are dominated by contributions from and .

## 6. Conclusions

Descriptions of snow microphysical properties, consisting of expected values, uncertainties, and correlations of microphysical parameters, are essential constraints for the remote sensing of snowfall. Although estimates of parameters for mass– and area–dimension relationships for various habits are available from a number of prior studies, these studies largely lacked information needed to estimate the distributions of these parameters. Additionally, in many cases estimates of both mass and area parameters were not available from the same sample of particles, thus there was little guarantee of consistency between the mass and area representations. To address these issues, a retrieval algorithm was developed that estimates the PDFs for these microphysical properties using multisensor observations of snowfall, designed around the observations available from a highly instrumented ground site used for a snowfall remote sensing field validation campaign: Rayleigh-regime radar reflectivity, snowfall rate, particle fall speeds, and size distributions.

The algorithm was applied to observations from four snowfall events involving both synoptic front and lake-effect processes. The measurements were found to primarily determine , and to a lesser degree and . Relatively little information was provided for and . The results showed nontrivial correlations between and , between and , and between and . These correlations, likely the result of the dependence of fall speed on the mass-to- ratio, are new results that were not apparent in the a priori estimates of these parameters because of the limitations of previous analyses. These correlations form off-diagonal elements in the covariance matrix describing the multidimensional PDF for these parameters, and may have significant effects on Bayesian retrievals that incorporate these PDFs as a priori information. The retrieved microphysical properties are not inconsistent with observations from other sources; however, it is essential to ensure that the parameters are determined using similar definitions of particle size.

The results from these analyses provide information essential to developing physically consistent representations of snow particle microphysical and scattering properties needed by snowfall retrievals using millimeter-wavelength observations. Given the complexity of snow crystals and the highly variable shape of aggregates, the goal for such representations is to be sufficiently realistic for retrieval purposes. Particle models were constructed using the retrieved microphysical properties, and so should be consistent with observed X-band reflectivities, fall speeds, and snowfall rates. While the models did not incorporate explicit representations of pristine shapes, models with reasonable aspect ratios and inhomogeneous structure were able to reproduce observed W-band reflectivities.

Characterizing the microphysical properties in the form of PDFs allows such retrievals to better quantify forward model uncertainties, which in turn allows the uncertainties in retrieved snowfall rates to be better quantified. Using particle models constructed with the retrieved microphysical properties, we found that the uncertainties in the retrieved microphysical parameters, especially those associated with particle mass, cause substantial uncertainties in modeled reflectivities. The estimated uncertainties due to uncertainties in (=8.0 dB) and (=6.3 dB) are large compared to observational uncertainties that might be expected for W-band radars like the ACR (Sekelsky et al. 1999), and also large compared to differences among the aggregate-like Ep, B8pr-30, and B8pr-45 particles (≈3 dB). Retrievals that account for the dependence of the reflectivity uncertainties should offer improved performance over those that do not.

The approach used to estimate these uncertainties assumes a linear error-propagation model and is likely a crude approximation for large uncertainties. Even so, in the applied context of retrievals that minimize differences between modeled and observed millimeter-wavelength radar reflectivities, these reflectivity model uncertainties are likely significant. As demonstrated by W14, however, improvements to ground-based observing systems may better constrain microphysical properties and reduce uncertainties. Additionally, extending these analyses to a more diverse range of snowfall regimes may reveal regime-dependent variations in microphysical properties and suggest approaches to further reduce their uncertainties. Several recent field experiments, including the Light Precipitation Validation Experiment (LPVEx) and GPM Cold Season Precipitation Experiment (GCPEx), have examined snowfall in different meteorological regimes and/or used enhanced suites of instrumentation. Although the OE method developed here was targeted specifically to observations from C3VP, it is readily adaptable to observations from these and similar future experiments.

## Acknowledgments

Parts of this research by NBW and TSL were performed at the University of Wisconsin–Madison and at Colorado State University for the Jet Propulsion Laboratory, California Institute of Technology, sponsored by the National Aeronautics and Space Administration. AJH acknowledges support from the JPL *CloudSat* Office and NASA Global Precipitation Measurement program Contract NN-X13AH73G. Computing resources for discrete dipole modeling were provided by the Community Computing Facility of the National Center for Atmospheric Research. Thanks are given to G.-J. Huang of Colorado State University, F. Fabry of McGill University, and L. Bliven of NASA Goddard Space Flight Center for making their C3VP datasets available and sharing their expertise. We appreciate the efforts of three anonymous reviewers who provided helpful feedback on the paper.

### APPENDIX

#### Constrained Discrete Dipole Modeling Method

To construct a dipole array for a particle with maximum dimension , first particle mass *m* and horizontally projected area are determined using the power laws and values of , , , and from the snow microphysics retrieval. The nondimensional area ratio can then be calculated as

From an assortment of nondimensional particle shapes, an instance of the shapes can be selected that matches the area ratio to within a small error, taken to be 1%.

A three-dimensional cubic lattice with lattice spacing *d* is defined inside the instance (Fig. A1, top). Possible values for *d* are obtained by dividing the maximum dimension by integer numbers of dipoles. For accuracy when calculating , the resulting lattice spacing should comply with

(Draine and Flatau 1994), where *n* is the complex refractive index of the dipole material and is the wavelength of the incident radiation. Any spacings that do not comply are discarded. Values of *d* were limited to no more than 73 μm, suitable for solid ice dipoles at 250 K and frequencies up to 183 GHz based on the ice refractive index data of Warren (1984).

Given *d*, the mass of a single dipole is and, knowing the required particle mass *m*, the required number of dipoles is

If exceeds 10^{4}, a value again based on accuracy criteria (Draine and Flatau 1994), and if the lattice has enough nodes on which to place the dipoles, the lattice spacing is acceptable. In practice, a range of *d* will produce acceptable lattices, so a *d* is chosen that makes the number of dipoles as small as possible without falling below 10^{4}.

Dipoles are placed randomly on the lattice nodes in such a way that and are maintained at the required values (Fig. A1, middle and bottom). The resulting particle, rather than being solid ice, contains porosities. The porosities are normally not permitted to extend uninterrupted vertically through the particle in order to preserve . This approach gives a particle model that matches the desired shape and that also matches the required , , and *m*. The random placement of dipoles means that the modeled scattering properties will have some amount of random variation. Scattering properties at millimeter wavelengths, however, should not be strongly sensitive to the fine elements of particle structure (Matrosov 2007). Additionally, since radar reflectivity is obtained by integrating over the particle size spectrum, the effects of random variations in scattering properties on radar reflectivity should be further reduced.

## REFERENCES

*Numerical Methods for Unconstrained Optimization and Nonlinear Equations.*Prentice-Hall, 378 pp.

*Proc. 22nd Conf. on Weather Analysis and Forecasting/18th Conf. on Numerical Weather Prediction,*Park City, UT, Amer. Meteor. Soc., P2.49. [Available online at http://ams.confex.com/ams/pdfpapers/123823.pdf.]

*CloudSat*validation project.

*Proc. Fourth European Conf. on Radar in Meteorology and Hydrology,*Barcelona, Spain, ERAD,

*Proc. Meteorological Satellite Conf.,*Sopot, Poland, EUMETSAT, P61_S4_08. [Available online at http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_CONF_P61_S4_08_HUDAK_V&RevisionSelectionMethod=LatestReleased&Rendition=Web.]

*2010 Fall Meeting,*San Francisco, CA, Amer. Geophys. Union, Abstract A13K-02. [Available online at http://abstractsearch.agu.org/meetings/2010/FM/sections/A/sessions/A13K/abstracts/A13K-02.html.]

*Inverse Methods for Atmospheric Sounding.*World Scientific Publishing, 240 pp.

*Proc. Ninth ARM Science Team Meeting,*San Antonio, TX, U.S. Department of Energy,

*The Mathematical Theory of Communication.*University of Illinois Press, 117 pp.

## Footnotes

The National Center for Atmospheric Research is sponsored by the National Science Foundation.