Fractal geometry of aggregate snowflakes revealed by triple‐wavelength radar measurements

Radar reflectivity measurements from three different wavelengths are used to retrieve information about the shape of aggregate snowflakes in deep stratiform ice clouds. Dual‐wavelength ratios are calculated for different shape models and compared to observations at 3, 35, and 94 GHz. It is demonstrated that many scattering models, including spherical and spheroidal models, do not adequately describe the aggregate snowflakes that are observed. The observations are consistent with fractal aggregate geometries generated by a physically based aggregation model. It is demonstrated that the fractal dimension of large aggregates can be inferred directly from the radar data. Fractal dimensions close to 2 are retrieved, consistent with previous theoretical models and in situ observations.


Introduction
Most ice phase precipitation at the surface is composed of aggregates of ice crystals [Hobbs et al., 1974], and aggregation often dominates the formation of large particles in deep ice clouds aloft [Field and Heymsfield, 2003;Westbrook et al., 2007]. In spite of this, the scattering properties of these aggregates (which depend on their shape) are not well characterized at present, and this leads to uncertainty in remote sensing retrievals of ice cloud properties and snowfall rates. Similarly, information about the shape of aggregates is needed to predict their mass and fall speed as a function of size [e.g., Mitchell, 1996], which is key to the interpretation of in situ observations and the parameterization of the ice phase in numerical weather prediction and climate models. Idealized theoretical models of the ice crystal aggregation process have been developed [Westbrook et al., 2004;Maruyama and Fujiyoshi, 2005] which suggest that the aggregates are fractal in their geometry (with a power law scaling between mass and size); but the predictions from these idealized models need to be tested against real particles in real clouds.
In this letter we develop an idea recently proposed by Kneifel et al. [2011]. We made measurements of radar reflectivity Z at three frequencies. The ratio of Z measured at a pair of frequencies is a function of particle size and particle shape; Kneifel et al. [2011] showed that the relationship between two such ratios (derived from measurements at three frequencies) allows us to distinguish between different particle shapes present in a given sample of cloud and hence validate or disprove theoretical models of their structure and scattering properties. Note that we use "shape" here to mean the complete structure of the aggregate snowflakes, rather than some qualitative descriptor of their overall form.
The study by Kneifel et al. [2011] was purely theoretical, and there is a lack of suitable triple-wavelength data available. Leinonen et al. [2012] analyzed a selection of vertical profiles from airborne radar measurements in convective clouds at 13, 35, and 94 GHz from the Wakasa Bay field campaign in 2003 and concluded the data were inconsistent with spherical/spheroidal models of the snowflakes. Kulie et al. [2013] extended that analysis to incorporate more profiles and more case studies from the Wakasa Bay experiment, but they report numerous practical difficulties in the analysis of these data. In particular, the relative calibration between the three radars was difficult to assure, and in a number of the clouds there was substantial attenuation of the 94 GHz data due to supercooled water; as a result their data show a very high degree of scatter.
The new measurements of aggregate snowflakes presented in this letter were obtained using three sensitive ground-based radars as part of the Chilbolton Triple-Wavelength Snowflake Experiment at the Chilbolton Observatory in southern England during 2014. Data were collected at 3, 35, and 94 GHz (wavelengths of 97.5, 8.6, and 3.2 mm, respectively), and the ground-based experiment allows for very close colocation in STEIN ET AL.
©2014. The Authors. space and time, while integration over longer time periods to obtain more accurate reflectivity samples is possible.

Aggregate Models and Dual-Wavelength Ratios
The radar reflectivity at a given wavelength and frequency is where (D) is the radar cross section for a particle of maximum dimension D and N(D) the particle size distribution (PSD). The factor 0.93 is the dielectric factor for liquid water at 3 GHz; we choose to normalize our 35 and 94 GHz reflectivity measurements by the same factor such that for a cloud of Rayleigh scattering ice particles,Z 3 = Z 35 = Z 94 (since the magnitude of the dielectric factor for ice is essentially frequency independent).
The key quantities of interest in this paper are the dual-wavelength ratios DWR 3−35 = 10 log 10 (Z 3 ∕Z 35 ) and DWR 35−94 = 10 log 10 (Z 35 ∕Z 94 ). For ice particles that are small compared to the wavelength (the Rayleigh regime) the radar cross-section ∝ −4 m 2 , where m is the particle mass [see, e.g., Westbrook, 2014, and references therein]. For our 3 GHz measurements, ice particles in the atmosphere will always be in the Rayleigh regime. For the higher frequencies, where particles of maximum dimension comparable to (or greater than) ∕4 may exist, interference occurs between the scattered waves originating from the front and rear of snowflake (as viewed by the vertically pointing radar beam). This leads to an additional dependence on the ratio D∕ , and we will characterize this deviation from the Rayleigh limit by a factor f so that ∝ −4 m 2 f . From equation (1) it is thus apparent that in this regime dual-wavelength ratios depend on the mean size of the particles, and their shape, as encapsulated in the function f . Crucially for what follows, the relationship between DWR 3−35 and DWR 35−94 depends purely on the shape, and plotting this relationship will allow us to distinguish between various shape models [Kneifel et al., 2011].
We assumed an exponential PSD in the calculations shown here; almost identical results are found when the more complex Field et al. [2005] PSD shape is used instead. We set the upper limit of integration to D max = 20 mm since this encompasses the range of aggregate sizes typically observed in the atmosphere.
Dual-wavelength ratios were calculated for several different scattering models, including soft spheres [Hogan et al., 2000], horizontally oriented oblate spheroids of aspect ratio 0.6 [Matrosov, 1998;Hogan et al., 2012], and more complex aggregate models [Westbrook et al., 2004;Petty and Huang, 2010]. Westbrook et al. [2004] (from now W04) developed an idealized model of the aggregation process that generates thousands of synthetic aggregates. These aggregates were found to be fractal with a power law scaling between mass and size m ∝ D d f , and a fractal dimension d f of 2. Westbrook et al. [2006Westbrook et al. [ , 2008 calculated the radar cross sections for the W04 aggregates using the Rayleigh-Gans approximation (RGA) and computed the average f for their large ensemble of aggregate geometries, and we find that the following function provides a close approximation to their simulated data: where x = 4 r g ∕ and the radius of gyration r g ≈ 0.3D for aggregates [Westbrook et al., 2006[Westbrook et al., , 2008. Tyynelä et al. [2013] presented radar scattering computations of aggregates generated using a very similar model; they show that the RGA underestimates the radar cross section of the aggregates a little but that the fractional bias is approximately independent of frequency, and so the use of the RGA results should not affect our DWR calculations.

Petty and Huang
[2010] constructed four realistic aggregate snowflake geometries and simulated their scattering properties using the Discrete Dipole Approximation to compute (D) for various wavelengths. These four models are aggregates of dendritic crystals ("DA1, " "DA2, " and "DA3") and needle crystals ("NA"). The dimensions of the aggregates were simply scaled up or down to simulate aggregates of various sizes, and hence, m ∝ D 3 for these particles. Data points for were relatively sparsely sampled in D space, and we interpolated logarithmically between values and extrapolated a limiting Rayleigh scattering behavior to compute for the small particles in the PSD. Full details can be found in the supporting information.

10.1002/2014GL062170
Finally, we considered spheres and spheroids composed of homogeneous air-ice mixtures, using Mie theory to compute for spheres, and Hogan et al.'s [2012] method for oblate spheroids. A mass-size relationship was used to compute the effective density of these mixtures. Various realistic power law relationships (m = aD b ) were tried: we found that the results were insensitive to the choice of a, and over the range of b suggested by Mitchell [1996] (1.7-2.2), so only the calculations performed with the Brown and Francis [1995] relationship will be reported and referred to as "MBF" and "MBF-O" (for oblate spheroids). Kneifel et al. [2011] also noted a similar insensitivity to mass-size parameters in their calculations. We also performed calculations with a fixed snowflake density of 100 kg m −3 ; these results were appreciably different to the Brown and Francis results and will be referred to as "MFIX" and "MFIX-O. "

The Triple-Wavelength Experiment
The Chilbolton Triple-Wavelength Snowflake Experiment was specifically designed to exploit the information content from simultaneous radar measurements at three frequencies. The sampling of the three radars in time was synchronized to within 0.1 s, and full pulse-to-pulse power and phase measurements were recorded. Spectral processing was used, and ground clutter was minimized by removing stationary targets from the Doppler spectra. The 35 GHz and 3 GHz radars at Chilbolton have a matched beam width of 0.25 • ; the 94 GHz radar is situated next to the 35 GHz radar and has a beam width of 0.5 • . The cloud radars are situated immediately adjacent to one another; the 3 GHz radar with its much larger 25 m antenna is sited less than 50 m away.
In the data presented here, the radar measurements are integrated over 10 s so that with a typical wind speed at 5 km of about 20 m s −1 , a 200 m section of cloud is observed. The range resolution of the radars is not identical (varies from 60 to 75 m), and so the data were interpolated onto a common 75 m grid.
The pointing of the radar antennas was checked by an engineer before the experiment began; any slight differences in sample volumes due to mispointing would be negligible compared to the separation of the radars themselves and the length of cloud sampled over the 10 s integration period.
Five case studies have so far been collected in the experiment; in this letter we present two of these. These two cases were chosen because they had the largest reflectivities and dual-wavelength ratios, and hence the best chance of robustly distinguishing between the various aggregate models.
The 3 GHz radar is absolutely calibrated to within 0.5 dB using the method of Goddard et al. [1994]. For the cloud radars, we adjust the reflectivity data so that the profiles match the 3 GHz values at cloud top where particles are small and thus Rayleigh scattering at all three wavelengths [Hogan et al., 2000]. The excellent sensitivity of our 3 GHz instrument (≈ −20 dBZ at 10 km) as well as the higher-frequency radars permits this very accurate relative calibration, something which was not possible with the Wakasa Bay measurements [Leinonen et al., 2012].
Specifically, we assumed that measurements near cloud top with Z 3 < −10 dBZ should be from Rayleigh scatterers, i.e., DWR 3−35 and DWR 3−94 are 0 dB in this region. So we calculated the DWR for our data, computed the median DWR found in this region, and subtracted it from Z through the entire profile.
We also excluded some of the data points. In some areas of cloud there was strong variability of reflectivity from profile to profile, which could lead to an unreliable relative calibration. To avoid this, the sample standard deviation s was calculated (in decibels) from reflectivities measured over a 50 s window at each range gate. Measurements for which s > 3 dB were excluded from the calibration procedure. The contour of s = 3 dB is shown in Figures 1a-1c and 3a. Similarly, if there were only one or two measurements in a profile which fulfilled the Z 3 < −10 dBZ condition, we did not feel confident that the relative calibration would be reliable, and so we exclude such profiles also.
The standard error of the reflectivity measurements was calculated following Hogan et al. [2000] using the Doppler spectral widths to compute the number of independent pulses in each sample. The measurement error was less than ±0.5dB at 3 GHz, and an order or magnitude smaller at the higher frequencies.
Specific gaseous attenuation between cloud base and cloud top was estimated using model output over Chilbolton provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). Over the range of interest in this study, the maximum specific gaseous attenuation at 94 GHz above 3 km for case 1 and case 2 was 0.09 and 0.13 dB km −1 , respectively. We neglected these small losses in our analysis. Because of the large antenna, it is necessary to apply a near-field correction to the 3 GHz data at heights below ∼ 6 km. This correction factor was derived by comparing 3 GHz reflectivity profiles against those measured by the 35 GHz instrument (which has a much smaller antenna) in a number of Rayleigh scattering ice clouds. The magnitude of the correction was 1 dB at 5 km, rising to 3 dB at 3 km. We chose not to consider data below 3 km height, since our results would have become increasingly sensitive to the precision of this correction factor. Figures 1a-1c show radar reflectivity profiles measured at Chilbolton for a thick stratiform ice cloud observed on 17 April 2014, ahead of a large-scale frontal rainband. Cloud top is approximately 9 km; reflectivity values increase from ≈ −15 dBZ near cloud top, to as much as +25 dBZ near the cloud base indicating substantial ice particle growth as the particles fall through the cloud layer; finally, there is a shallow region at cloud base where the particles sublimate and the reflectivity decreases sharply. Temperature profiles from the ECMWF model show that cloud top was at ≈ −45 • C, while the cloud base temperature increased from ≈ −15 • C to −5 • C over the observation period. The 0 • C isotherm was at approximately 2.7 km (below cloud base). Liquid water path retrieved from colocated microwave radiometer measurements were less than 50 g m −2 , which implies negligible liquid water attenuation of the radar beams. Figure 1d shows an example vertical profile from the cloud. At cloud top, all three reflectivities are matched and track one another closely. As one moves down through the profile Z 94 becomes systematically smaller than Z 3 and Z 35 , indicating the presence of larger particles. Below 4.3 km, Z 35 becomes significantly smaller than Z 3 , indicating that the particles are now large enough that the 35 GHz reflectivities are no longer dominated by Rayleigh scattering. At 4 km, DWR 35−94 ≈ 8 dB and DWR 3−35 ≈ 6 dB. It is in this region where DWR 35−94 and DWR 3−35 are both significantly greater than 0 dB that we can apply our triple-wavelength analysis.

Case Study 1: 17 April 2014
As noted in section 3, the errors on the reflectivity data are very small. However, the colocation of the radars is not perfect despite their close proximity, and some small-scale inhomogeneity is always present between the sample volumes, leading to additional noise in the DWR data. As an indication of the magnitude of this STEIN ET AL. We calculated the DWR for all of the pixels inside the black contour on Figures 1a-1c and compared the measurements against theoretical curves from the aggregate models; the results are shown in Figure 2. To avoid cluttering the figure, data with reflectivities Z 3 below +5 dBZ are not shown since these data have DWR ≈ 0 dB and contain no shape information. We find that the measurements are in close agreement with the W04 fractal aggregate model. The sphere and spheroid models predict DWR 35−94 which are several decibels too large. The Petty and Huang DA3 and NA aggregate models also show poor agreement with the radar observations, again predicting DWR 35−94 values which are too high. The DA1 and DA2 aggregate models lie closer to the measurements; however, DWR 35−94 increases monotonically with DWR 3−35 in the measurements, while the DA1 and DA2 models increase at small DWR 3−35 before curving back again at large DWR 3−35 , which we do not observe.
In addition to aggregate models, Kneifel et al. [2011] considered the pristine crystal models provided in the Liu [2008] scattering database. We did initially perform this comparison as well; all but one of the models were inconsistent with our data. We found that their dendrite crystal model did give results compatible with the measurements; however, the high reflectivities and Doppler velocities of ≈ 1 m s −1 observed in our data convinced us that the particles were likely to be aggregates. Further evidence to support this assumption is provided in Figure 1e which shows the reflectivity-weighted mean particle diameter D Z Using an exponential PSD and the W04 model, D Z can be derived as a function of DWR, since both quantities are independent of particle concentration [see, for example, Matrosov, 1998]. Retrieval of D Z from DWR 35−94 using this function shows that the radar reflectivities (Z 3 ) were dominated by particles up to 10 mm in size, which is consistent with large aggregated particles, rather than single vapor-grown pristine crystals.
Interestingly, we observe that D Z increases sharply below 4.5 km height (corresponding to a temperature of ≈ −10 • C). This may suggest an increase in aggregation efficiency as the particles become warmer as suggested by Hosler et al. [1957].

Case Study 2: 22 May 2014
We present a second case observed on 22 May 2014, when thick stratiform cloud associated with a complex system of convective cells and heavy rain passed northward over Chilbolton. The 3 GHz reflectivities for this case are shown in Figure 3, with the DWR 3−35 versus DWR 35−94 data in Figure 3c. Cloud top was 8 km (≈ −40 • C), while the melting layer was at approximately 1.9 km.
Compared to the first case, higher values of reflectivity Z 3 (up to +30 dBZ) and DWR are observed. In terms of the relationship between DWR 35−94 and DWR 3−35 the results closely resemble those for the 17 April 2014 case: the majority of high DWR values lie close to the behavior predicted by the W04 model, although the center of the distribution appears to lie at slightly higher values of DWR 35−94 for a given DWR 3−35 . Unlike case 1 there was low-level cloud and rain present, and only 24 profiles had successful relative calibration. For the period prior to 1150 UTC, the depth over which −20 < Z 3 < −10 dBZ was regularly less than 180 m. For the period between 1150 UTC and 1200 UTC, both Z 35 and Z 94 were completely extinguished before Z 3 < −10 dBZ. Echoes below 3 km were not considered in the analysis as noted earlier.
Since we perform a relative calibration at cloud top, this attenuation should not affect our results unless there is significant attenuating liquid water mixed with the ice particles above 3 km. Such attenuation would lead to a value of Z 94 that is artificially high at lower altitudes and therefore artificially decrease DWR 35−94 for a given DWR 3−35 (which are affected much less by liquid water).
The reason for the slightly higher values of DWR 35−94 may be microphysical as Doppler velocities, shown in Figure 3b, were higher than case 1, with particles falling at ≈ 2 m s −1 in regions of higher dual-wavelength ratio, suggesting that the aggregates in this case are rimed.

Retrieving the Fractal Dimension
We note that both in Figures 2 and 3b, the value of DWR 35−94 appears to approach an asymptotic limiting value as DWR 3−35 (and hence the size of the particles) increases. Similar behavior was observed by Leinonen et al. [2012] (although they noted saturation values for DWR 35−94 of 7 dB and 13 dB for their two profiles).
The origin of this "saturation" effect may be explained by the fractal geometry of the aggregates. Sorensen [2001] shows that for aggregates which are comparable to ∕4 in dimension or larger, we can expect the function f to approach a power law scaling f (x) ∝ x −d f . For large particle diameters, we can approximate Under these assumptions, we obtain DWR → d f × 10 log 10 ( 1 ∕ 2 ) .
For 1 = 8.6 mm (35 GHz) and 2 = 3.2 mm (94 GHz) we obtain DWR 35−94 = d f × 4.3 dB. We observe approximately 8 dB for case study 1 and 9 dB for case study 2, and we may therefore retrieve an estimate of the fractal dimension of the snowflakes of d f = 1.9 and d f = 2.1, consistent with the theoretical model of W04. However, the data are inconsistent with the predictions of the Maruyama and Fujiyoshi [2005] aggregation model which found d f between 1.1 and 1.3.
The retrieved values of d f close to 2 are consistent with aircraft image analysis  and several mass-size power laws from the literature which have an exponent close to 2 [Brown and Francis, 1995;Mitchell, 1996;Heymsfield et al., 2010;Cotton et al., 2013].
We also note that theory restricts d f < 3, which implies that we should always expect DWR 35−94 < 13 dB if the particles are fractal.
This effect furthermore implies that there is a limit to the extent that D Z can be retrieved from DWR 35−94 alone if the particles are sufficiently large; for the W04 model, this corresponds to a D Z of approximately 10 mm; above this value increasing particle size leads to a negligible increase in DWR 35−94 .

Discussion and Conclusions
We have used ground-based measurements from radars at three different wavelengths to evaluate different models of aggregate geometry and retrieve the fractal dimension. We have shown that two different case studies had a similar relationship between DWR 3−35 and DWR 35−94 , consistent with the fractal model proposed by W04. The observation that DWR 35−94 approaches a limiting value is consistent with fractal geometry and allows the fractal dimension of the aggregates to be estimated from observations: d f ≈ 2 in our cases. The result that ice particles can be modeled as fractal snowflakes rather than oblate spheroids is important for radar retrievals of ice water content or snowfall rates. For instance, Hogan and Westbrook [2014] show that ice water content retrieved from 94 GHz radar reflectivities may be in error by as much as a factor of 5 approximating the particle by a soft spheroid model rather than representing the true (fractal) structure.
Much less scatter in the ground-based radar measurements is observed than previously reported from aircraft observations [Kulie et al., 2013], and we feel our data are the most convincing evidence thus far that Kneifel et al.'s [2011] idea of distinguishing between different snowflake shapes is feasible in practice. Unlike results from Leinonen et al. [2012] and Kulie et al. [2013], there is no indication that a mixture of different types of aggregates are present for either of the cases studied. Indeed, the similarity of the results in the two different clouds reported here, despite the rather different meteorological conditions, suggests that the variability between aggregate shape populations from cloud to cloud may be modest, which is an important result for radar retrieval algorithms and microphysics parameterization if true. The measurements in this letter form part of the Chilbolton Triple-Wavelength Snowflake Experiment, and we plan to collect a number of additional case studies in the coming months to quantify this further.
It would be useful to compare the observations to other aggregate scattering models. Hogan and Westbrook [2014] used W04's aggregates to derive an equation for the scattering behavior of the particles. They decomposed the distribution of mass along the radar beam into a mean shape (which determines the scattering at small size parameters) and smaller-scale deviations from that mean (which determines the scattering behavior at large size parameters). The power spectrum of the deviations had a power law form. The equation contains three fitted constants, with different parameters recommended for horizontal or random orientation. We found that their equation gave results compatible with our data, provided that their parameters for aggregates of plates and random orientation were assumed. Using the parameters for orientated aggregates led to DWR 35−94 values which were too high compared to the observations. Similarly, Leinonen et al. [2013] generated aggregate geometries and parameterized the mass distribution within the aggregate in terms of Gaussian functions. Again, a number of fitted constants were required to characterize the mass distribution functions. It seems likely that triple-wavelength measurements such as those reported here can be used to refine these parameterized models.
We have neglected attenuation due to the aggregate snowflakes themselves in our analysis: this could be a source of error. However, the RGA used by Westbrook et al. [2006] also allows an estimate of this attenuation from the retrieved particle size and concentrations (Tyynelä et al. [2013] shows that RGA estimates of extinction are accurate to ∼20%.). At 94 GHz, the maximum specific attenuation thus calculated for case 1 and case 2 was 0.05 and 0.14 dB km −1 , respectively. The two-way integrated attenuation due to snow over the depth of the cloud was nowhere greater than 0.12 and 0.66 dB, respectively. Given these small values, the attenuation due to snow should not significantly affect the results of this study. Likewise, although significant attenuation could occur due to rainfall and liquid clouds in case 2, the relative calibration method applied in our study effectively accounts for these losses.