## Abstract

Two algorithms are presented for ice crystal classification using 95-GHz polarimetric radar observables and air temperature (*T*). Both are based on a fuzzy logic scheme. Ice crystals are classified as columnar crystals (CC), planar crystals (PC), mixtures of PC and small- to medium-sized aggregates and/or lightly to moderately rimed PC (PSAR), medium- to large-sized aggregates of PC, or densely rimed PC, or graupel-like snow or small lumpy graupel (PLARG), and graupel larger than about 2 mm (G). The 1D algorithm makes use of *Z*_{h}, *Z*_{dr}, LDR_{hv}, and *T,* while the 2D algorithm incorporates the three radar observables in pairs, (*Z*_{dr}, *Z*_{h}), (LDR_{hv}, *Z*_{h}), and (*Z*_{dr}, LDR_{hv}), plus the temperature *T.* The range of values for each observable or pair of observables is derived from extensive modeling studies conducted earlier. The algorithms are tested using side-looking radar measurements from an aircraft, which was also equipped with particle probes producing simultaneous and nearly collocated shadow images of cloud ice crystals. The classification results from both algorithms agreed very well with the particle images. The two algorithms were in agreement by 89% in one case and 97% in the remaining three cases considered here. The most effective observable in the 1D algorithm was *Z*_{dr}, and in the 2D algorithm the pair (*Z*_{dr}, *Z*_{h}). LDR_{hv} had negligible effect in the 1D classification algorithm for the cases considered here. The temperature *T* was mainly effective in separating columnar crystals from the rest. The advantage of the 2D algorithm over the 1D algorithm was that it significantly reduced the dependence on *T* in two out of the four cases.

## 1. Introduction

The classification of ice crystals in clouds using radars is expected to be useful for cloud microphysical studies. Knowing the crystal type in the radar resolution volume will also be helpful in improving the quantitative estimation of cloud ice water content (IWC) with radar (Aydin and Tang 1997a). Polarimetric radars at millimeter-wave frequencies, such as 95 GHz, have the potential for classifying ice crystals in clouds. The short wavelength (3.2 mm), compared to S- or C-band meteorological radar wavelengths (10 and 5 cm, respectively), increases the radar's sensitivity to small ice crystals (diameters less than 2 mm) that may be present in many clouds (e.g., see Aydin and Walsh 1999, and the references therein). The size of a millimeter-wave radar can be made small enough to be mounted on aircraft for close range observations, which is important at these higher frequencies because of the increased atmospheric attenuation (Lhermitte 1987; Pazmany et al. 1994; Galloway et al. 1997).

The sensitivity of polarimetric radar observables to size, shape, orientation, and dielectric constant (which is affected by the density of the ice crystals and the amount of liquid water in the case of melting particles) of ice crystals make them suitable for use in bulk classification algorithms. The utility of S-band polarimetric radars for distinguishing ice and water phase hydrometeors has been extensively demonstrated (Seliga and Bringi 1978; Hall et al. 1984; Bringi et al. 1984; Aydin et al. 1986; Carey and Rutledge 1998; Nanni et al. 2000; Kennedy and Detwiler 2003). Reinking et al. (2002) show the potential use of elevation angle dependence of depolarization ratio (for slant linear and elliptical polarizations) at 35 GHz for differentiating drizzle droplets, planar crystals, columnar crystals, and irregular ice particles. More recently, fuzzy classification techniques have been proposed for identifying a wide range of hydrometeor types using S-band radars (Vivekanandan et al. 1999; Liu and Chandrasekar 2000; Straka et al. 2000; Zrnic et al. 2001). This paper presents and illustrates two algorithms for cloud ice crystal classification based on 95-GHz polarimetric radar measurements. In situ verification of the results from these algorithms is performed using particle probe images obtained simultaneously with the radar data.

## 2. Classification algorithms

The algorithms presented here for cloud radars are similar to the “fuzzy logic” techniques that are being used for hydrometeor classification schemes with S-band meteorological radars (Vivekanandan et al. 1999; Liu and Chandrasekar 2000; Straka et al. 2000; Zrnic et al. 2001). They employ 95-GHz polarimetric radar signatures of various crystal types, which are established mainly based on computational studies utilizing the finite difference time domain method (Kunz and Luebbers 1993; Taflove and Hagness 2000) for columnar and planar crystals and the *T*-matrix method (Waterman 1969; Mishchenko et al. 1999) for aggregates and graupel. The one- and two-dimensional (1D and 2D) algorithms are described below.

The radar observables used in these algorithms are the effective reflectivity factor at horizontal polarization *Z*_{h} in dB*Z* units, the differential reflectivity *Z*_{dr} = *Z*_{h} − *Z*_{υ} in decibel units, and the linear depolarization ratio LDR_{hv} = *Z*_{hv} − *Z*_{υ} in dB units, where *Z*_{hv} is the cross-polarized reflectivity factor (in dB*Z* units) for vertically polarized transmission and horizontally polarized reception. The effective reflectivity factors in this paper are based on the dielectric constant of water. There are other polarimetric observables such as the copolarized correlation coefficient (*ρ*_{hv}) and the specific differential phase (*K*_{dp}) that could be helpful in the classification schemes but are not incorporated in this study because the radar data used in demonstrating the algorithms did not have these observables. Preliminary simulations show that |*ρ*_{hv}| near side incidence might be useful in the classification algorithms for differentiating columnar crystals from other crystals types without the aid of the air temperature. On the other hand, the utility of *K*_{dp} in the classification process is not clear at this time because its measurement requires long paths in the cloud, which may not always be available, and its range resolution is not as good as the other observables used in this study. Finally, the air temperature *T,* which may be inferred from the altitude of the measurements, is also used as a parameter in the classification algorithms.

The five ice crystal classes used for classification purposes are given in Table 1. These should be modified and/or increased as more information becomes available on the radar signatures of different ice crystals. Given the radar signatures for the crystal types discussed in sections 2a and 2b, it appears that these classes are presently the most feasible groupings resolvable with the three radar observables (*Z*_{h}, *Z*_{dr}, and LDR_{hv}) plus the temperature *T.*

Note that the fuzzy logic hydrometeor classification algorithms developed for S-band frequency radars also classify ice crystals, but in different and more general categories than the ones used here. For example, Liu and Chandrasekar (2000) consider the following classes in addition to various categories of rain and hail: low-density dry ice, high-density dry ice, wet ice crystals, dry graupel, and wet graupel. On the other hand, Zrnic et al. (2001), similar to the one proposed by Vivekanandan et al. (1999), classify ice crystals as horizontally oriented ice crystals, vertically oriented ice crystals, dry snow, wet snow, and graupel and/or small hail. None of these algorithms differentiate columnar or planar crystals or their aggregates and rimed forms.

### a. The 1D algorithm

The 1D-classification algorithm is based on associating a weighting function corresponding to each parameter (radar observables and temperature) for each ice crystal class. A simple trapezoidal weighting function as shown in Fig. 1 is used for *Z*_{h}, *Z*_{dr}, and LDR_{hv}. For some crystal types the temperature (*T*) weighting function requires two such trapezoids. A weighting function is unity over most of the range of values expected for a given parameter. The weighting function drops linearly from a maximum value of 1 to a minimum of 0, and remains 0 outside the expected range. The trapezoidal-weighting-function transition points are given in Table 2 for the radar observables and in Table 3 for *T.* These values may be fine tuned as more experimental observations are made with in situ confirmation. The range of values for the radar observables corresponding to each ice crystal type are obtained from model simulations (Aydin and Tang 1997a,b; Walsh 1998; Aydin and Walsh 1999; Aydin et al. 2000; Aydin and Mazlum 2001; Aydin and Somaya 2001), which are generally in good agreement with the experimental observations of Wolde and Vali (2001). The temperatures at which different ice crystal types can exist are available in the literature (Ono 1970; Takahashi et al. 1991; Pruppacher and Klett 1997, their section 2.2); this parameter is mainly useful in differentiating columnar crystals from aggregates and rimed planar crystals.

The basic fuzzy logic classification algorithm can be described as follows (Zrnic et al. 2001). Let the class of ice crystals be denoted by the integer *j*, which ranges from 1 to 5 in this study. The ice crystals of a resolution volume are assigned to the class that produces the largest value for the sum *S*_{j}:

Here *W*_{j}(*Y*_{i}) is the weighting function for the *j*th class of ice crystals corresponding to the parameter *Y*_{i} that represents the three radar observables *Z*_{h}, *Z*_{dr}, and LDR_{hv}, and the temperature *T.* If additional parameters are used, then the summation would be increased to include those as well. *A*_{i}, which is the weight factor (or weighting function multiplier) for the parameter *Y*_{i}, allows for the possibility of emphasizing some parameters over others, if there is such a need. For example, initially all parameters have *A*_{i} equal to unity. However, in cases where there is ambiguity between planar crystals and another class (i.e., when they have the same sum *S*_{j}), then the weight factor of *Z*_{dr} is doubled. The ambiguity is resolved when one of these classes produces a larger sum *S*_{j}. This approach works because the PC class (planar crystals) has *Z*_{dr} values that do not overlap with the other classes. A similar approach is used to resolve ambiguities between the G class (graupel larger than about 2 mm) and the rest. Ambiguities between the PSAR class (mixtures of PC and small- to medium-sized aggregates and/or lightly to moderately rimed PC) and the PLARG class (medium- to large-sized aggregates of PC and/or densely rimed PC, or graupel like snow, or small lumpy graupel) are not possible to resolve with the four parameters that are used here. In such cases these locations are tagged and indicated as ambiguous regions. The temperature parameter is very effective in minimizing ambiguities between the CC class (columnar crystals) and the other classes, specifically the PSAR class. Note that the spatial variability of ice crystal classes might also be incorporated in resolving some of the ambiguities by requiring some amount of spatial coverage and continuity.

### b. The 2D algorithm

The 2D algorithm is similar to the 1D version with the exception that it utilizes two-dimensional weighting functions for the radar observables, which are paired up as (*Z*_{dr}, *Z*_{h}), (LDR_{hv}, *Z*_{h}), and (*Z*_{dr}, LDR_{hv}). Figures 2 and 3 show these functions for the columnar and planar crystal types. The vertices of the polygonal weighting functions are given in the figure captions. A weighting function has a value of zero outside the polygon and transitions linearly to a value of unity inside the polygon. This transition occurs over the intervals 0.5, 1, and 2 dB for *Z*_{dr}, LDR_{hv}, and *Z*_{h}, respectively, similar to the 1D case. For the sloped lines in the polygons, if the slope is greater than 1, then the taper range for the parameter of the horizontal axis is used; if the slope is less than (or equal to) 1, then the taper range for the vertical axis parameter is used. The remaining crystal types have rectangular-shaped regions in the two dimensional plots, which are derived from their corresponding 1D weighting functions; for example, consider the (*Z*_{dr}, *Z*_{h}) weighting function for the PSAR class, the vertices would be (1.5, −16), (1.5, 9), (5.5, −16), and (5.5, 9). This is a result of the large uncertainties in the shapes, fall behavior, and densities of nonpristine crystals, which affect their radar signatures yielding the broad rectangular shaped clustering patterns. The temperature (*T*) uses the same 1D weighting function in this algorithm as well. Initially all three pairs of parameters and *T* are given the same weight factor (*A*_{i}) of unity. Similar to the 1D case, when there is ambiguity between the planar crystal class and another class, the weight factor for the pair (*Z*_{dr}, *Z*_{h}) is doubled, which usually leads to the resolution of the ambiguity by increasing the sum *S*_{j} for one of the classes. The same approach is used when there is a ambiguity involving the G class (graupel larger than about 2 mm).

## 3. Testing the algorithms

The classification algorithms described above are applied to unique datasets obtained during the WYICE'97 winter–spring campaign (February–April 1997) near Laramie, Wyoming. A 95-GHz radar system (details available online at http://www-das.uwyo.edu/wcr/;) on board the Wyoming King Air aircraft obtained side-looking dual-polarization measurements (*Z*_{h}, *Z*_{dr}, and LDR_{hv}) while flying through the cloud at a speed of about 100 m s^{−1}. Particle imaging probes (2D-P and 2D-C) attached to the wings of the aircraft provided shadow images of ice crystals simultaneously with the radar measurements. The radar range resolution was 30 m (0.2-*μ*s pulse duration) and the first sampling volume that could be used started at 120-m radial distance from the radar. The first two radial sampling volumes were considered to be almost collocated with the particle probe images obtained under the aircraft wings. The near uniformity of the radar observables over several sampling volumes along range supported this assumption for the cases considered here.

### a. Radar data

The radar data consisted of range profiles obtained in about 87 ms using a pulse sequence of VVHH (0.2-*μ*s pulses 50-*μ*s apart) every 1 ms, leading to approximately 58 independent samples per profile. [The antenna diameter is 0.3 m and the aircraft speed is 100 m s^{−1}; therefore, one has an independent sample every (0.3/2)/100 = 1.5 ms, and in 87 ms there are 87/1.5 = 58 independent samples.] A running average was also performed on the radar data over two range gates and seven profiles. The averaging process led to 2 × 7 × 58 = 812 independent samples and corresponded to an effective resolution volume of dimensions approximately (60 m) × (60 m) × (1.5 m) at 120-m range. The standard deviations of the radar observables after averaging were about 0.15 dB for *Z*_{h}, 0.21 dB for LDR_{hv}, and 0.03–0.11 dB (for zero-lag copolarized correlation coefficient ranging from 0.85 to 0.99) for *Z*_{dr} (Sachidananda and Zrnic 1985; Pazmany et al. 1994; Bringi and Chandrasekar 2002, their section 6.5). The absolute value of *Z*_{h} is considered to be accurate to within ±3 dB and the minimum detectable LDR_{hv} is −22 dB (Wolde and Vali 2001). It should be noted that increasing the minimum detectable value of LDR_{hv} to −17 dB did not cause any noticeable change in the results presented here. The uncertainty in the absolute value of *Z*_{h} had some effect and will be discussed at the end of this section. The bias in *Z*_{dr} was determined to be about −1.3 dB based on vertical pointing measurements and was corrected. The bias in LDR_{hv} is estimated to be less than 1 dB; a 1-dB bias had negligible effect on the results presented here.

### b. Results

Figure 4a shows the measured radar observables *Z*_{h}, *Z*_{dr}, and LDR_{hv} looking to the side of the aircraft on 18 February 1997. The vertical axis is the time in seconds, which can be converted to distance in the along-track (flight) direction by multiplying with the aircraft speed *υ* ≈ 100 m s^{−1}. The altitude and air temperature during this flight are between 4274 and 4307 m and −14.4° and −13.4°C. The aircraft roll is generally within ±6° and does not exceed ±10°, while the pitch is confined to less than 5°.

Figures 4b and 4c show the ice crystal classification results from the 1D and 2D algorithms using all four parameters (*Z*_{h}, *Z*_{dr}, LDR_{hv}, and *T*). Overall they look similar and there is an 89% agreement in the first range gate results, which compare very well with the particle images, as seen in the sample set of Fig. 5. The images are from the 2D-P and 2D-C probes, the first and second strips in each set, respectively. The vertical dimension of the 2D-P probe is 8 times that of the 2D-C probe; that is why the 2D-P probe images appear to be much smaller.

The algorithms produce results consistent with the particle probe images. The particle images in Fig. 5 corresponding to the time intervals 1.5–3, 23–25, and 82–84 s are predominantly stellar crystals; hexagonal plates are also visible in the 1.5–3-s interval, as well as some rimed crystals. Both algorithms (Figs. 4b and 4c) indicate the PC class (light blue) for these time intervals, except near 3 s, where the PSAR class (yellow) appears, which is consistent with the rimed crystal images. In the time intervals 10–12.5 and 47–48 s the images show rimed stellar crystals and their aggregates, which are in agreement with the PSAR class produced by the algorithms. The time interval 61–62 s shows larger aggregates of rimed stellar crystals, which is in agreement with the PLARG class (red) indicated by the algorithms over a few pixels.

Note that the ice crystal classes shown in Figs. 4b and 4c have sums *S*_{j} [Eq. (1)] that are at least 25% above the next highest sum corresponding to another class for the same location (in some cases two classes may tie with the second highest sum). This 25% difference is significant, indicating a high level of confidence in the classification results.

Datasets from three other days in 1997, 11 February (beginning at 2250:06 UTC), 24 March 24 (beginning at 0011:58 UTC), and 2 April (beginning at 2000:48 UTC), also produced very good comparisons between the classification results and the particle probe images. Both algorithms indicate the following classes: PLARG for 11 February and for 24 March, and CC for 2 April, which agree very well with the particle images shown in Fig. 6. February 11th has densely rimed crystals and graupel-like snow, 24 March is mainly composed of large rimed aggregates of stellar crystals, and 2 April contains columnar crystals.

### c. Effectiveness of parameters

It is of interest to determine the effectiveness of each parameter, or pair of parameters, in the algorithms. The datasets chosen for illustrating the classification algorithms were selected such that the aircraft pitch and roll angles were within ±5° and ±10°, respectively. This is because the polarimetric radar signatures presented here are for side incidence, which can also be used for angles close to side incidence. (Larger deviations will require modifications in the weighting functions for the radar observables.) Tables 4 and 5 summarize these results for the 1D and 2D algorithms, where using all the parameters are compared with using fewer parameters. For the 18 February 1997 dataset and the 1D algorithm, *Z*_{dr} is the most effective parameter, *Z*_{h} and LDR_{hv} contribute the least, while *T* influences about 26% of the data (i.e., when only three parameters, *Z*_{h}, *Z*_{dr}, and LDR_{hv}, are used, the agreement with using all four parameters is 74%). If *Z*_{dr} is removed from the algorithm, then the remaining parameters produce sums (*S*_{j}) that are the same for more than one class and a decision cannot be made (the ambiguity cannot be resolved by increasing the weight factor of a parameter). This leads to a very low percentage of agreement. If only one parameter is used in making the decision (i.e., three parameters are removed), then *Z*_{dr} alone produces 73% agreement, and each of the other parameters produce no agreement at all. For the 18 February 1997 dataset and the 2D algorithm, the (*Z*_{dr}, *Z*_{h}) pair is the most effective among all the parameters. When this pair is removed from the algorithm, the agreement reduces to 58%. When only *T* is removed, the agreement is 83%, which is a significant improvement (9%) over the 1D algorithm without *T.* Using the (*Z*_{dr}, *Z*_{h}) pair alone in the 2D algorithm produces 77% agreement. Similarly, (*Z*_{dr}, LDR_{hv}) or (LDR_{hv}, *Z*_{h}), or *T* alone, produce 42%, 2%, and 0% agreement, respectively. It is worth noting that 0% agreement (i.e., no agreement) does not necessarily result from different classes being indicated by the algorithms, but rather due to ambiguous results when fewer parameters are used.

For the 11 February and 24 March 1997 datasets (both were composed of the PLARG class) it is clear that in the 2D classification algorithm the determining pair is (*Z*_{dr}, *Z*_{h}) and the other pairs of parameters, as well as the temperature, have a negligible effect on the classification results. For the 24 March dataset in the 1D algorithm *Z*_{dr} alone determines the ice crystal class and the other parameters do not contribute at all. On the other hand, the 1D algorithm applied to the 11 February case shows a dependence on both temperature and *Z*_{dr}. This is in contrast to the 2D algorithm that shows no dependence on temperature. Finally, for the 2 April dataset (predominantly composed of columnar crystals) it is seen that the temperature is the determining parameter in both algorithms assuming that the radar is not looking in a cloud where the G class (graupel larger than about 2 mm) can exist. The agreement between the 1D and 2D algorithms for these three cases were above 97%.

Among all the radar observables considered here the largest uncertainty is in *Z*_{h}, which has an estimated accuracy of ±3 dB. This could cause changes in the results of the classification algorithms shown in Fig. 4. Table 6 shows the comparison with biases of +3 and −3 dB in *Z*_{h}. In each case the *Z*_{h} + 3 dB results have 90% or better agreement with the results shown in Fig. 4, except for the 2D algorithm on 18 February, which has an 87% agreement. The 3-dB increase in *Z*_{h} changed some pixels classified as PSAR to PLARG for the 11 February, 18 February, and 24 March cases and from CC to PSAR for the 2 April case. The *Z*_{h} − 3 dB results show at least 92% agreement in all cases except 11 February, which has agreements of 57% and 46% for the 1D and 2D algorithms. For 11 February a 3-dB decrease in *Z*_{h} changed nearly half of the pixels classified as PLARG to PSAR.

## 4. Conclusions

Two algorithms were presented for ice crystal classification using 95-GHz polarimetric radar observables and air temperature *T.* The parameters for the 1D algorithm were *Z*_{h}, *Z*_{dr}, LDR_{hv}, and *T.* The 2D algorithm used pairs of radar parameters, (*Z*_{dr}, *Z*_{h}), (LDR_{hv}, *Z*_{h}), and (*Z*_{dr}, LDR_{hv}), and *T.* The algorithms were tested on four datasets from experiments conducted in Wyoming in 1997. An aircraft carrying a 95-GHz polarimetric radar and particle imaging probes flew through different clouds, providing side-looking radar measurements simultaneously with particle shadow images. Ice crystal classes derived from the algorithms compared very well with the in situ shadow images of ice crystals. Both algorithms produced similar results with 89% agreement in one case and above 97% agreement in three cases. One advantage of the 2D algorithm was that it reduced the dependence on air temperature *T* in two out of the four cases presented here. In all four cases LDR_{hv} had a negligible contribution. Further testing will be needed on more datasets to determine its utility at side incidence. It should be noted that LDR_{hv} (or LDR_{vh}) can discriminate columnar crystals from other crystal types at vertical incidence (looking straight up or down) or, if available, from the elevation angle dependence of this parameter. The most useful radar observable in the 1D algorithm was *Z*_{dr}, and in the 2D algorithm the pair (*Z*_{dr}, *Z*_{h}) was most effective.

The air temperature *T* is mainly helpful in discriminating columnar crystals from the rest. However, this parameter cannot take into account circumstances involving sedimentation in stratiform clouds and advection in convective clouds, where the crystals may be moved from their original location to a region with a different temperature. This could lead to errors. Radar parameters are not affected by such relocation unless different crystal types are mixed together. The advantage of a radar-based classification algorithm is that the entire cloud (or a good portion of it) could be scanned in a short time interval and the spatial variability of ice crystal classes could be mapped. This would be very useful in tracking the sedimentation, advection, riming, and aggregation processes over space and time. Tracking of the riming process may provide an indirect way of identifying regions of supercooled droplets.

## Acknowledgments

A major part of this research was supported by the National Science Foundation under Grant ATM-9813512. The authors would like to thank Dr. R. D. Kelly of the University of Wyoming for providing the radar and particle probe data, and Samuel Haimov for helpful discussions on the radar's calibration.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{DR}measurement considerations for fast scan capability radar.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Dr. Kultegin Aydin, Department of Electrical Engineering, 203 Electrical Engineering East, The Pennsylvania State University, University Park, PA 16802. Email: k-aydin@psu.edu