Articles | Volume 13, issue 9
Research article
08 Sep 2020
Research article |  | 08 Sep 2020

Retrieval of lower-order moments of the drop size distribution using CSU-CHILL X-band polarimetric radar: a case study

Viswanathan Bringi, Kumar Vijay Mishra, Merhala Thurai, Patrick C. Kennedy, and Timothy H. Raupach

The lower-order moments of the drop size distribution (DSD) have generally been considered difficult to retrieve accurately from polarimetric radar data because these data are related to higher-order moments. For example, the 4.6th moment is associated with a specific differential phase and the 6th moment with reflectivity and ratio of high-order moments with differential reflectivity. Thus, conventionally, the emphasis has been to estimate rain rate (3.67th moment) or parameters of the exponential or gamma distribution for the DSD. Many double-moment “bulk” microphysical schemes predict the total number concentration (the 0th moment of the DSD, or M0) and the mixing ratio (or equivalently, the 3rd moment M3). Thus, it is difficult to compare the model outputs directly with polarimetric radar observations or, given the model outputs, forward model the radar observables. This article describes the use of double-moment normalization of DSDs and the resulting stable intrinsic shape that can be fitted by the generalized gamma (G-G) distribution. The two reference moments are M3 and M6, which are shown to be retrievable using the X-band radar reflectivity, differential reflectivity, and specific attenuation (from the iterative correction of measured reflectivity Zh using the total Φdp constraint, i.e., the iterative ZPHI method). Along with the climatological shape parameters of the G-G fit to the scaled/normalized DSDs, the lower-order moments are then retrieved more accurately than possible hitherto. The importance of measuring the complete DSD from 0.1 mm onwards is emphasized using, in our case, an optical array probe with 50µm resolution collocated with a two-dimensional video disdrometer with about 170µm resolution. This avoids small drop truncation and hence the accurate calculation of lower-order moments. A case study of a complex multi-cell storm which traversed an instrumented site near the CSU-CHILL radar is described for which the moments were retrieved from radar and compared with directly computed moments from the complete spectrum measurements using the aforementioned two disdrometers. Our detailed validation analysis of the radar-retrieved moments showed relative bias of the moments M0 through M2 was <15 % in magnitude, with Pearson’s correlation coefficient >0.9. Both radar measurement and parameterization errors were estimated rigorously. We show that the temporal variation of the radar-retrieved mass-weighted mean diameter with M0 resulted in coherent “time tracks” that can potentially lead to studies of precipitation evolution that have not been possible so far.

1 Introduction

The principal application of polarimetric radar has historically been directed towards more accurate estimation of rain rate (R) that is driven largely by the operational agencies for hydrological applications. It now strongly appears that, as a major step forward, the operational algorithm for the US Weather Surveillance Radar – 1988 Doppler (WSR-88D) network will be based on specific attenuation because, among other advantages, it is linearly related to rain rate at S-band (Ryzhkov et al.2014; Cocks et al.2019; Wang et al.2019). This method has also been evaluated quite extensively at X-band by Diederich et al. (2015), where the specific attenuation (Ah) is much larger than at S-band but not linear with R. The development of R(Ah) algorithms rests on a large body of work since the early 1990s and is related to attenuation correction using the differential propagation phase as a constraint (Bringi et al.1990; Smyth and Illingworth1998; Testud et al.2000; Bringi and Chandrasekar2001, and references therein).

The retrieval of drop size distribution (DSD) parameters has also been a strong impetus for radar polarimetry. In this context, there exists a large body of literature that is based mainly on the unnormalized (Ulbrich1983) or normalized (Illingworth and Blackman2002; Testud et al.2001) gamma model. This model is parameterized by a set of three quantities, namely {N0,μ,Λ} or {Nw,μ,Dm}, where N0 and Nw are “intercept” parameters, μ is the shape factor, Λ is the “slope”, Dm is the ratio of the 4th to 3rd moments of the DSD N(D), and D is the diameter of the raindrop (Ryzhkov and Zrnić2019, and references therein). The gamma model has also been used in the double-moment “bulk” microphysical schemes that predict the mass mixing ratio (or equivalently the 3rd moment M3; for moment order k, we write Mk and the total concentration of drops (or M0) (e.g., Meyers et al.1997). The lower-order moments of the DSD (M0 through M3+b, where b is the exponent of the fall-speed-D power law), are important in describing various microphysical processes such as collisional (break-up and coalescence), evaporation, and sedimentation (e.g., Milbrandt and Yau2005). However, radar polarimetry has not been focused on these lower-order moment retrievals because the radar observables (horizontal) reflectivity Zh, differential reflectivity Zdr, and specific differential phase Kdp are related to the higher-order moments such as M6, the ratio M7M6, and M4.5, respectively.

Defining a scaled diameter x=D/Dm, the normalized DSD is a function of x as h(x)=N(D)/Nw. The observation of Testud et al. (2001) regarding the “remarkable” stability of the shape of h(x) using measured DSDs was a significant advance because they did not impose an a priori form for h(x). Apart from the shape “stability” of h(x), they also showed a large compression in the “scatter” of h(x) compared to N(D). While they did not refer to their normalization as double-moment using M3 and M4 as the reference moments, Lee et al. (2004, henceforth, L04) generalized the scaling/normalization framework by introducing any two reference moments Mi and Mj of any order i, j>0. As per this framework, in a compact notation, N(D)=N0h(x) with a different x=D/Dm, where N0=Mi(j+1)/(j-i)Mj(i+1)/(i-j) and Dm is the ratio of (Mj/Mi)1/(j-i). In essence, the variance of the DSDs due to different rain types and intensities is largely controlled by the variability in N0 and Dm and much less so by h(x). Further, any moment Mk can be expressed as power laws of Mi, Mj, in which the coefficient is the kth moment of h(x) and the two exponents are predetermined by i and j. L04 also recognized that if h(x) is assumed to follow the generalized gamma (G-G) model with two shape parameters, then it could fit most naturally occurring DSD shapes. We refer the reader to Stacy (1962) for the expressions of the probability density function (pdf) of the G-G and its moments. The G-G form has been applied to model cloud droplet spectra, ice crystal and snow spectra (Delanoë et al.2014), as well as raindrop spectra (Raupach and Berne2017a; Thurai and Bringi2018). The three-moment normalization for this model is provided in Szyrmer et al. (2005). The generalization to K-moment normalization scheme given in Morrison et al. (2019) does not specify any particular form for h(x) other than that its moments should be finite. In agreement with Szyrmer et al. (2005), they found that three-moment normalization was sufficient to “compress” the scatter of h(x). Further, it minimized the errors in the estimation of the other moments expressed as power laws of the reference moments. For remote-sensing applications (cloud and drizzle), they found that the set of moments {M2,M3,M6} was one possible choice mentioning lidar backscatter (M2), microwave attenuation (M3), and radar reflectivity (M6). While the combination of M3 and M6 was not optimal for estimating the lower-order moments (in particular, M0), it was better than using M6 alone.

Recently, using the double-moment approach of L04, Raupach and Berne (2017a, RBa) showed that measured DSDs in stratiform rain with h(x) expressed in the G-G form have shape factors that are sufficiently “invariant” for practical use across different rain climatologies if the reference moments are chosen carefully. Their result essentially validated the “remarkable” stability conclusion of h(x) by Testud et al. (2001) which was based on limited data in oceanic rain. RBa speculated that the transition (i.e., between convective and stratiform rain) and convective rain DSDs would also have a sufficiently “invariant” h(x), but they did not have a large enough database to make such a conclusion.

The polarimetric (X-band) radar retrieval of moments using reference moments M3 and M6 and an “invariant” h(x) of the G-G form were first described in Raupach and Berne (2017b, RBb). Their retrieval of M6 was based on Zh, while M3 was retrieved in a two-step procedure using Zdr and Kdp. We discuss their results in detail later in this paper. Here, it suffices to mention that their measured DSDs were based on a network of Parsivel disdrometers, which did not have the resolution to measure the shape of h(x) for D<0.75 mm or so (as shown later by Raupach et al.2019). Additionally, with “noisy” Zdr and Kdp radar data (for Zh<37 dBZ), their validation of the moments (M0 through M7) using radar measurements was not conclusive but was sufficient to demonstrate that their approach gave results similar to other methods based on the normalized gamma model using “more classical” radar retrievals of {Nw,μ,Dm} (Gorgucci et al.2008; Kalogiros et al.2012).

Whereas RBb used measured DSDs and the polarimetric radar forward operator to derive the retrieval algorithms for M3 and M6, there has been a reverse moment-based polarimetric forward operator (Kumjian et al.2019). This reverse approach employs a very large database of measured and bin-resolved one-dimensional (1-D) model output DSDs to build a look-up table that maps the various moment pairs to the expected values of Zh, Zdr, and Kdp along with their standard deviations. Their application was to determine the moment pairs that could potentially be prognosed in numerical microphysical schemes of rain would be “optimally” constrained by polarimetric radar measurements. They found that the pair {M6,M9} was optimal in terms of lowest variability in {Zh,Zdr,Kdp} but that the pair {M3,M6} was suboptimal but still “useable”. Thus, RBb’s choice of {M3,M6} as the two reference moments was validated by Kumjian et al. (2019). The rationale through which M9 (whose sampling error would be very large using available disdrometers) entered the moment pair is not entirely clear. It could be because of correlating Zdr with absolute moments (M0 through M9) as opposed to the more physically based ratio of moments such as Dm=(M6/M3)1/3 in RBb or M7M6 as in Jameson (1983).

This work further develops on RBb but, instead of Kdp, we employ specific attenuation (Ah) given its operational use in estimating R. The iterative correction of measured reflectivity Zh using the total Φdp constraint, i.e., the iterative ZPHI algorithm, which is a variant of Testud et al. (2000), is used here to estimate Ah (Bringi et al.2001; Park et al.2005a, b). The reference moment M3, which is proportional to rain water content (W), is retrieved using a modification of Jameson (1993) by fitting AhW as a smoothed cubic spline with Dm. The prior step is the retrieval of Dm from Zdr and then the retrieval of Dm from Dm. This multi-step procedure was found to minimize the parameterization errors (also referred to as algorithm errors) in the estimation of M3. As in RBb, the reference moment M6 is derived as power-law fit to Zh. The other major difference with RBb is the use of collocated optical array probe (50 µm resolution) and two-dimensional video disdrometer (2DVD) inside a double-fence wind shield (Thurai et al.2017a). The “complete” DSD was measured from 0.1 mm onwards thus avoiding truncation at the small drop end. This leads to more accurate estimates of the lower-order moments as well as more accurate h(x). The methodology of G-G fits to h(x) are described in Thurai and Bringi (2018) and Raupach et al. (2019). The use of a very narrow (0.33) beam at X-band with high gain and a short vertical distance from radar pixel to the instrumented site were additional factors that differed from RBb. We also show coherent “time tracks” in the Dm versus M0, Dm versus W, and Dm versus M6 planes, where all variables are based on radar retrievals.

The rest of the article is organized as follows. In the next section, we briefly discuss the surface instrumentation (disdrometers). The CSU-CHILL radar and its use in characterizing the multi-cell storm complex (the case study) as well as data extracted over the instrumented site are given in Sect. 3. The retrieval of the two reference moments {M3,M6} follows in Sect. 4. Several different ways of validating the moment retrievals are presented in Sect. 5. We follow this with a discussion in Sect. 6 and summarize the case study in Sect. 7. We cannot draw firm conclusions from just one case study even though the analysis is quite detailed. Rather this work may be considered a proof-of-concept that will require further validation to be undertaken in the future. It is difficult to find radar data with revisit times <90 s over an instrumented site unless a dedicated experiment is proposed and funded. In our case, the event of 23 May 2015 was largely a target of opportunity as one of the co-authors (PCK) had the foresight to collect data on this day without considering that it would lead to a detailed case study of moment retrievals. An Appendix provides procedures for estimating the radar measurement error contribution to the variances of, firstly, the reference moments {M3,M6}, and then the variances of the other moments. The estimates of the variances of the ratio of correlated variables of the form XpYq are derived using a Taylor expansion to second order. The parameterization error variances are estimated for {M3,M6} and summed with the radar measurement error variances to yield the total error variance for each moment retrieved.

Throughout this paper, we use “H” as a subscript for reflectivity ZH to denote units of dBZ at horizontal polarization. The lower case “h” in Zh means units of mm6 m−3. The same applies to ZDR (in dB) or Zdr (ratio). The functions Var(⋅) and () yield the variance and mean of their arguments, respectively. We use Cov(X,Y) for the covariance between the variables X and Y. A set is denoted by curly brackets {,,}. The notation 𝔼{} is used for the statistical expectation; 〈⋅〉 for the average of its argument; Γ(⋅) for the gamma function; and 𝕀m{⋅} for the imaginary part of its complex argument.

2 Surface instrumentation

The principal surface-based instruments used in this study are the MPS (or Meteorological Particle Spectrometer, manufactured by Droplet Measurement Technologies) and third-generation 2DVD, both located within a 2∕3-scale Double Fence Intercomparison Reference (DFIR; Rasmussen et al.2012) wind shield. As reported in Notaroš et al. (2016), the 2∕3-scale DFIR was effective in reducing the ambient wind speeds by nearly a factor of 3 based on data from outside and inside the fence. An OTT Pluvio rain gage was also available for rain rate and rain accumulation comparison with the disdrometers.

Our retrieval algorithms (see Sect. 4) of the reference moments M3 and M6 were based on scattering simulations from the combined DSD data from two sites, namely Greeley, Colorado (GXY), and Huntsville, Alabama (HSV). The same disdrometer and wind shield configuration were deployed at both locations. However, the case study in this paper concerns the event of 23 May 2015 captured in Greeley, which also has a coincident CHILL X-band radar. Huntsville has a very different climate from Greeley, and its altitude is 212 m mean sea level (m.s.l.) as compared with 1.4 km m.s.l. for Greeley. According to the Köppen–Trewartha climate classification system (Trewartha and Horn1980), Greeley has a semiarid-type climate, whereas Huntsville has a humid subtropical-type climate (Belda et al.2014).

The MPS is an optical array probe (OAP) that uses the technique introduced by Knollenberg (1970, 1976, 1981) and measures drop diameter in the range from 0.05 to 3.1 mm (but the upper end of the usable range is limited to 1.5 mm due to reduced sampling volume). A 64-element photo-diode array is illuminated with a 660 nm collimated laser beam. Droplets passing through the laser cast a shadow on the array, and the decrease in light intensity on the diodes is monitored with the signal processing electronics. A two-dimensional image is captured by recording the light level of each diode during the period that the array is shadowed. The limitations and uncertainties associated with OAP measurements have been well documented (Korolev et al.1991, 1998; Baumgardner et al.2017). The sizing and fall speed errors primarily depend on the digitization error (±25µ). The fall speed accuracy according to the manufacturer (DMT) is <10 % for 0.25 mm and <1 % for sizes greater than 1 mm, limited primarily by the accuracy in droplet sizing. To calculate N(D), the measured fall speed is not used. Rather a cubic polynomial fit from the manufacturer (DMT) is employed. Details of the calculation of N(D) are given in the Appendix of Thurai et al. (2017a) and updated in Raupach et al. (2019).

The third-generation 2DVD is described in detail by Bernauer et al. (2015). Its operational characteristics are similar to earlier generations. In particular, the accuracy of size and fall speed measurement has been well documented (e.g., Schönhuber et al.2007, 2008; Thurai et al.2007, 2009; Huang et al.2008). Considering the horizontal pixel resolution of about 170µm and other factors, the effective sizing range is D>0.7 mm (Thurai et al.2017a). The fall velocity accuracy is determined primarily by the accuracy of calibrating the distance between the two orthogonal light “sheets” or planes and is <5 % for fall velocity <10 m s−1 (Schönhuber et al.2008). A comparison of fall speeds between the MPS and 2DVD has been reported by Bringi et al. (2018) from both the GXY and HSV sites, with excellent agreement. The only fall velocity threshold used for the 2DVD is the lower limit set at 0.5 m s−1 in accordance with the manufacturer guidelines for rain measurements. The instrument is designed to prevent drops from entering the housing where the cameras are positioned. Without going into details, it suffices to mention that small drops can enter via slits that allow the light to illuminate the cameras, or drops can hang on the slits. Both of these effects cause spurious images that the matching software cannot reject (Larsen and Schönhuber2018). Thus, caution is necessary when using the 2DVD fall speeds for sizes <0.6 mm (about 3–4 pixels).

Figure 1(a) Conceptual illustration of the complete DSD comprising the drizzle mode, the “shoulder” region, and the precipitation mode. The incomplete DSD is due to drop truncation by instruments that cannot measure the concentration of small drops. (b) An example of measured 3 min averaged DSD (R≈60 mm h−1) using a collocated optical array probe with a 2DVD showing the separate measurements (note the high resolution of the MPS and the drop truncation of the 2DVD). The composite or compete spectrum is obtained by using the MPS for D≤0.75 mm and 2DVD for D>0.75 mm. The dashed blue line is the G-G fit (with parameters μ=-0.3, c=6; see Eq. (1) for details) to the complete spectrum. Data from the 23 May 2015 case study at 20:45 UTC. (c) Same data points as panel (b) but with the standard gamma (black) and exponential (red) fits.


In our application, we utilize the MPS for measurement of small drops with 0.1D<0.75 mm and the 2DVD for larger-sized drops (see Raupach et al.2019 for the rationale). This is termed here the “complete” size spectrum, and 2928 3 min averaged spectra were available from the two sites with a minimum rain rate of 0.1 mm h−1 and a maximum of 286 mm h−1. More details of the rainfall types, measurement time periods, comparison with gages, and related analyses are available in Thurai et al. (2017a) and Raupach et al. (2019). Figure 1a illustrates the “complete” DSD with the “drizzle” mode defined by a peak in N(D) occurring when D<0.5 mm (Abel and Boutle2012). The “shoulder” is the diameter range where the N(D) either remains steady or falls off more “slowly” (generally found under equilibrium conditions: McFarquhar2004). The precipitation range is used here for larger-sized drops after the “shoulder”, if any. These ranges are used here only to illustrate Fig. 1. The “incomplete” spectra, in which small drops are not measured accurately due to resolution, sensitivity, or other issues (2DVD or Parsivel; Park et al.2017), frequently show the convex down shape at the small drop end. Here, we only use the complete N(D) by compositing the data from MPS and 2DVD. An example is shown in Fig. 1b, which illustrates the main features of the “complete” DSD during the time that peak rain rate (3 min averaged R of 60 mm h−1) was occurring at the instrumented site during the 23 May 2015 event. We refer the reader to the following link for further details: (last access: 27 August 2020). The shape is equilibrium-like but with a single “drizzle” mode, a well-defined shoulder, and faster (than exponential) fall-off in the tail (Low and List1982; McFarquhar2004; Straub et al.2010). These features cannot be captured by either standard gamma or exponential fits, as shown in Fig. 1c.

3 CSU-CHILL radar

The CSU-CHILL (Colorado State University-University of Chicago/Illinois State Water Survey) radar is described in Brunkow et al. (2000) and Bringi et al. (2011a). Details on the conversion to a dual-wavelength system and the current radar specifications are given in Junyent et al. (2015) (see the condensed version in Table 1).

Suffice to state here that the X-band polarimetric mode is “simultaneous transmit and receive” or SHV and the 3 dB beam width is very narrow at 0.33 with a gain of 55 dB. There are three separate feed or orthogonal mode transducers (OMTs) available: (a) an S-band feed (beam width in the far field is 1) whose performance is described in Bringi et al. (2011a), (b) an S–X-band dual-wavelength feed that was used in the data described herein, and (c) an X-band feed. For the 23 May 2015 event, our retrievals of the reference moments {M3,M6} are based on the X-band polarimetric data {Zh,Zdr,Φdp} only, where Φdp is differential phase shift. Only X-band data were used even though S-band data were available simultaneously. Our choice for using the X-band data was due to the very high resolution provided by the 0.33 beam and the larger range of values for the X-band specific attenuation for a given rain rate (relative to S-band). The case study convective event was a complex of multiple cells with strong azimuthal and elevation gradients across the “echo cores” which generally precludes accurate dual-wavelength estimation of range-resolved specific attenuation. The narrow X-band beam also allows a lower elevation angle (1.5) to be used before clutter contaminates the signal. The instrumented site was located at Easton, which is 13 km SSE of the radar (along the 171.25 azimuth). Details of the terrain variation between the radar and the Easton site are given in Kennedy et al. (2018).

Table 1Technical specifications of the CSU-CHILL X-band channel with the dual-offset Gregorian antenna.

Download Print Version | Download XLSX

3.1 Brief description of storm characteristics from radar

The synoptic environment on 23 May 2015 was conducive to thunderstorm development in northeastern Colorado. A low at the 500 hPa level was analyzed over Utah at 12:00 UTC. This system was forecast to move eastward and promote upward motion within the moist air mass that was in place over the eastern plains of Colorado. In recognition of this situation, the Storm Prediction Center (SPC) Convective Outlook valid for the afternoon hours included a slight risk of severe thunderstorm development over northeastern Colorado. Persistent low cloud coverage ended up limiting surface heating within ∼50 km of CSU-CHILL, reducing thunderstorm intensity. The SPC storm reports did not contain any severe category hail (diameter of 2.54 cm or larger) or surface wind speeds of 25 m s−1 or more. Volunteer weather observers reported several instances of small (0.64 cm or less) hail mostly in the Rocky Mountain foothills ∼60 km southwest of the radar. Afternoon surface temperatures were ∼14C in the CSU-CHILL/Greeley area. The 0C level in the Denver late afternoon sounding was at ∼3.5 km m.s.l. (2.1 km above ground level – a.g.l.).

Figure 2 shows a low elevation angle (1.5) plan position indicator (PPI) scan of S-band ZH at 20:45 UTC which was the time of peak rainfall at the instrumented site (also referred to as Easton) identified by the + marker in Fig. 2. The main echo feature is the near N–S orientation of multiple 55 dBZ cores extending from the Easton site to nearly 50 km to the south. The rainfall over the site lasted for over 90 min, and PPI scans at a fixed 1.5 elevation angle were repeated every 90 s. This good time resolution enabled the validation of the moment retrievals which otherwise would not have been possible (for example, with WSR-88D scan cycle times of around 5 min).

Figure 2Low elevation angle (1.5) PPI sweep of S-band reflectivity (ZH) at 20:45 UTC. The “+” marks the location of the instrumented site (MPS and 2DVD).

The general echo movement near Easton was estimated at 10 m s−1 towards the radar on average from the south. After the peak echo of 55 dBZ traversed the instrumented site, another cell produced very heavy rain at the radar site with no evidence of graupel/hail (visual observations by one of the co-authors, PCK). One volunteer observer located 15 km east of the Easton instrumentation site reported 0.64 cm hail mixed with heavy rain between 20:30 and 20:45 UTC. The CSU-CHILL radar data showed that this small hail was generated by an isolated, higher-reflectivity cell that was separated from the storms that crossed the instrumented site. We do not believe that hail occurred at the Easton site during the analyzed time period, as also confirmed by 2DVD fall speed observations.

There was no range–height indicator (RHI) scan at 20:45 UTC. Therefore, the vertical echo structure could not be determined at this time close to the peak rainfall over Easton, but RHI scans about 18 min earlier performed along the 171 azimuth are shown in Fig. 3. The top panel shows measured reflectivity at the X-band (uncorrected for attenuation), while the measured S-band reflectivity is plotted in the bottom panel. Several strong cells (>55 dBZ) are noted south of Easton at ranges of 27 and 32 km; the cell at 32 km shows significant attenuation. However, there is no significant attenuation at the 13 km range where the instrumented site is located. The 10 dBZ echo top reaches 8 km a.g.l.

Figure 3RHI sweep of ZH along the 171 azimuth at 20:27 UTC about 18 min before peak echo descended on the Easton instrumented site at a range of 13 km. The “+” marks are at 2 km intervals. (Top) X-band measured (uncorrected) reflectivity. The range profiles of radar data along the dashed line are shown in Fig. 4. (Bottom) S-band measured reflectivity.

Figure 4(a) Range profiles of measured and attenuation-corrected ZH at X-band at 20:27 UTC at an azimuth angle of 171 (elevation angle of 2). The Easton instrumented site is located at a range of 13 km, (b) measured and attenuation-corrected ZDR at X-band, (c) measured and filtered Φdp, and (d) specific attenuation (Ah).


Figure 5(a) PPI of measured X-band reflectivity at 20:43 UTC (elevation angle is 1.5). The range profiles in the panels below are along the red line (radial) to the instrumented site denoted by the “+” marker. (b)(e) As in Fig. 4, panels (a)(d), except at 20:43 UTC.

3.2 Attenuation correction

As mentioned earlier, we use the X-band radar for quantitative moment retrievals. It is apparent that the strong cells will cause attenuation, so the X-band measured Zh and Zdr have to be corrected for attenuation and differential attenuation, respectively. The method used herein is exactly the same as described in Mishra et al. (2016). To correct the measured Zh, we apply an iterative version of the ZPHI method, which uses a Φdp constraint (Testud et al.2000; Bringi et al.2001) that was originally developed at C-band but later extended to X-band by Park et al. (2005a, b). In short, the coefficient α in the linear relation between the specific attenuation at H polarization (Ah) and specific differential phase (Kdp) is determined by minimizing a cost function based on least squares (we refer to Bringi and Chandrasekar2001, for details), whereas the standard ZPHI method uses a fixed a priori value for α (Testud et al.2000). In addition, a power law of the form Ah=b2Zhb1 is assumed where b1=0.78 and b2 are constants (Park et al.2005a). The method gives the estimate of Ah at each resolution volume in the selected range interval (here 0–40 km). The upshot of using Ah instead of Kdp is that the former closely follows the variations in Zh without the smoothing needed for estimating the latter but with all the advantages of Kdp such as immunity to calibration offsets and partial beam blockage (Ryzhkov et al.2014). Note that using Ah for retrieval of W is restricted to precipitation comprising pure rain. In contrast, using Kdp (as in RBb) in pure rain entails spatial (range) smoothing which, in compact convective rain cells, “distorts” the spatial representation of the rain rate profile depicted by ZH. In our multi-step retrieval procedure, it is reasonable not to mix different smoothing scales for the radar observables. There are many variants of the attenuation-correction method at X-band, as elucidated, for example, by Anagnostou et al. (2004) and Gorgucci and Chandrasekar (2005). Here, the iterative filtering method of Hubbert and Bringi (1995) is used to separate the backscatter differential phase from the propagation phase. In essence, the estimate of Ah may be considered a by-product of attenuation correction of the measured Zh using the differential propagation phase over the selected path interval as a constraint.

The correction of the measured Zdr for differential attenuation is based on an extension of the method proposed by Smyth and Illingworth (1998) for C-band, which is described in Bringi and Chandrasekar (2001) as a “combined Φdp-Zdr” constraint. The extension to X-band is described in Park et al. (2005b), which is used herein with some modifications implemented for the CSU-CHILL radar. In brief, the Ah determined by the Φdp constraint is scaled by a factor ν and the measured Zdr is corrected for differential attenuation (Adp=νAh) such that a desired value is reached at the end of the beam. The desired value is the intrinsic or “true” Zdr at the end of the beam, which is estimated from the corrected Zh using a mean ZhZdr relation based on scattering simulations that use measured DSDs from several locations that encompass a wide variety of rain types. This sets a constraint for Zdr at the end of the beam (generally Zdr≈0 dB because of light rain at the end of the beam or because of ice particles above the 0C level). By the end of the beam, we mean the last range gate where the hydrometer echoes are detected. Range profiles of measured and corrected ZH and ZDR, the measured and filtered Φdp (which is used as constraints from 0 to 40 km), and Ah are shown in the four panels of Fig. 4 at 20:27 UTC along the radial to the instrumented site located at Easton. The ZH profiles show that very minor attenuation correction is needed at this time, while the ZDR is corrected by 2 dB at the end of the ray. The change in differential phase, i.e., ΔΦdp, is also small at 25. Consistent with these values, the Ah peak is 1.5 dB km−1, coinciding with the ZH peak at 25 km. At the Easton location (13 km range) the Ah is negligible.

Figure 5a shows the PPI (at an elevation angle of 1.5) of the measured X-band ZH at 20:43 UTC, at which time the peak 55 dBZ echo traversed over the instrumented site. The X-band reflectivity in Fig. 5a can be compared with S-band data in Fig. 2. The line of cells organized south of the radar causes significant attenuation of the X-band signal power. This is clear in the range profile in Fig. 5b, where the attenuation has increased dramatically with ZH corrected by 35 dB and ZDR by 9 dB. The ΔΦdp now increases by around 150. Assuming a nominal α of 0.25 km−1, the path-integrated attenuation would be 37.5 dB. The Ah values have increased with peaks of 3 dB km−1. At the Easton site, Ah≈1.5 dB km−1. From a moment-retrieval viewpoint, significant attenuation correction only begins beyond the Easton site (13 km range), so that the errors due to such correction will not be significant in this case. This situation persists after 20:43 UTC until the end of the analysis period (21:35 UTC or so).

3.3 Time series of radar measurements and DSD-based simulations

A “necessary” condition for accurate radar retrievals of DSD moments is that the time series of corrected Zh, Zdr, Kdp, and Ah extracted over the resolution volumes (or pixels) surrounding the Easton site agree “reasonably” well with the same observables simulated using measured DSDs and a scattering model (what is generally referred to as the forward radar model/operator). The criterion of “reasonable” agreement is difficult to quantify but is elucidated in Thurai et al. (2012) using error variance separation. The radar data were extracted around a polar area defined by a range interval ±0.18 km centered at the range (13 km) to Easton and ±0.2 in azimuth angle for a total of 15 pixels surrounding the Easton site. The height of the pixels at an elevation angle of 1.5 at the 13 km range is 340 m a.g.l. The radar data from each pixel are plotted as a time series in Fig. 6 which shows the pixel-to-pixel variability. A Lee (1980) filter (henceforth Lee filter) used to reduce speckle in images is adapted here to filter the pixel-to-pixel variability with a sliding window of ±11 (weighted) points; the filtered values are shown in Fig. 6 interpolated in time to that of the disdrometer. The “effective” time resolution after Lee filtering is 2.5 min. The radar time series were shifted backward in time by 60 s, as is common when matching the peak in Zh (e.g., May et al.1999). A more general analysis of the error characterization of radar–gage comparison is given in Anagnostou et al. (1999). However, such an analysis is not needed herein because of the narrow antenna beam and short range to the instrumented site. Thurai et al. (2012) applied the Lee filter to time-series data versus range filtering applied to range gates along a fixed ray profile and showed that they were nearly equivalent. The Lee-filtered values of the radar data show the time evolution of the main echo passage over the Easton site.

Figure 6Time series of X-band radar data compared with simulations based on measured complete DSDs and the scattering model described in the text. (a) Corrected ZH from radar showing pixel-to-pixel variations which have been filtered using Lee (1980). (b–d) Same but, respectively, for corrected ZDR, Kdp, and Ah.


The composite 3 min averaged DSDs (an example was shown in Fig. 1b) were used to simulate radar observables as a time series using the T-matrix scattering code (Barber and Yeh1975; Bringi and Seliga1977). The time resolution of 3 min corresponds to a spatial scale of 1.8 km (using an echo movement speed of 10 m s−1), which is less than the echo cell sizes estimated to be around 2–3 km. For an estimate of the decorrelation time of radar-retrieved D0, we refer the reader to Thurai et al. (2012), who studied stratiform rain with embedded weak convection using 4 s samples; the 1∕e-folding time was around 200 s, where e is Napier's constant. For a highly convective case of our present case study, the decorrelation time would be substantially smaller but probably similar to our radar sampling of 90 s. Further, in Bringi et al. (2015), the decorrelation distance in a highly convective squall line event was estimated to be 3.5 km for radar-retrieved R; the time resolution obtained using a radar sampling time of 40 s was 2.5 min. In our case, the echo speed of 10 m s−1 together with the “effective” Lee-filtered radar resolution of 2.5 min and the disdrometer 3 min averaging corresponds to spatial distances in the range 1.5–1.8 km. This is well within the estimated 1∕e-decorrelation distance of 3.5 km in Bringi et al. (2015).

While the DSD data were available at much higher time resolution, the choice of 3 min averaged DSD is a compromise between smaller DSD sample sizes when integration times are, say, 1 min, versus poorer representativeness of the spatial scales for longer time integrations (e.g., 5 min). The radar update time was around 90 s, which is short enough not to introduce excessive temporal representativeness errors. Note that in Fig. 6a, the DSD-based simulated ZH is about 5 dB lower than that measured by radar for the time period 20:00–20:30 UTC. The measured ZH was 18–25 dBZ, implying very low rain rates (∼0.5 mmh−1) and a low number density of drops sampled by the disdrometers. In addition, it follows from the RHI taken at 20:27 UTC in Fig. 3 that the cells are moderately slanted from NNW aloft, where generating cells might have formed at 5 km a.g.l. to SSE at the surface. Given the unsteady conditions in this complex of echoes, it is not surprising that the disdrometer-based ZH calculation is biased low by around 5 dB relative to low radar ZH values of 18–25 dBZ. These problems are mitigated when heavier rain rates traverse the site about 15 min later.

The scattering model is based on the mean shapes from the 80 m fall bridge experiment described in Thurai et al. (2007) and Gaussian canting angle distribution with mean = 0 and standard deviation σ=7 (from Huang et al.2008). The dielectric constant of water at a wavelength of 3 cm and a temperature of 8C (Ray, 1972) was used. The time series of the simulated radar observables are shown in Fig. 6 marked as “DSD”. The visual agreement between simulations and the Lee-filtered mean radar values are qualitatively quite good except for a small underestimation of simulated Zdr relative to radar measurements by around 0.5 dB at 20:45 UTC. The discrepancy at 21:40 UTC noted in Fig. 6a and d is because of heavy rain on radome (observed by PCK). Overall, the good agreement between corrected radar measurements and the DSD-based forward simulations shows good calibration of Zh and Zdr. The radar-retrieved specific attenuation closely follows the Zh due to the AhZh power-law assumption in the ZPHI method (only the fixed exponent is relevant), whereas the Kdp does not, as expected.

4 The methodology of radar retrieval of the DSD moments

As mentioned in Sect. 1, the methodology we used here follows RBb except for the use of specific attenuation rather than Kdp in the retrieval of M3 (however, both methods use Zdr in a multi-step retrieval described below). There are several advantages to this approach. First, “noisy” Ah is strictly positive, as opposed to Kdp, which can be “noisy” with both positive- and negative-valued fluctuations in measurements. This is an issue because horizontal orientation of raindrops is usually assumed for simulation of Kdp from DSD measurements, meaning that all simulated Kdp values used to train retrieval algorithms are positive. Second, the smoothing in range necessary for Kdp is not needed for Ah, which closely follows the spatial variability in Zh. The basis for the retrieval methodology lies in the double-moment normalization of L04. This method is explained in detail in RBb and, hence, we only summarize it in the next subsection.

4.1 Overview

Moment retrievals from polarimetric radar data are a relatively recent application of the scaling/normalization of the DSD. There are several aspects in this scaling as described by L04, namely that there is substantial reduction in the scatter in h(x) from the un-normalized scatter, implying that most of the variability of the DSD can be attributed to variability in N0 and Dm, with h(x) being relatively “stable” with varying rain types/intensities. There is considerable latitude in the choice of reference moments {Mi,Mj} in the double-moment scheme depending on the application with N0 expressed as Mi(j+1)/(j-i)Mj(i+1)/(i-j) and Dm as a ratio of (Mj/Mi)1/(j-i). Further, any moment Mk can be expressed as power laws of Mi, Mj, and the kth moment of h(x). RBa showed that the amount of variance in individual DSD moments captured by the normalization scheme depends on the choice of reference moments.

RBb first suggested the use of {M3,M6} as the two reference moments suitable for polarimetric radar retrievals of the DSD. They proposed retrieval of M6 from radar measurements of Zh, while for M3 the retrieval was based on {Zdr,Kdp}. While h(x) can be of any functional form, the G-G model h(x;μ,c), with two positive shape parameters μ and c, from L04 was chosen by RBb. The key to accurate retrievals of Mk depends not only on the retrieval accuracy of the reference moments, but also on h(x;μ,c), which has to be representative of the rain climatology. The estimation of {μ,c} requires a large database of DSD measurements but, more importantly, the small drop end (the fit to which is controlled mostly by μ) of the distributions needs to be measured accurately as discussed in Sect. 2, because otherwise the lower-order moments M0 through M2 will be in error.

In our retrieval the h(x) from 1594 3 min DSDs (with rain rate >0.1 mm h−1) collected by the MPS and 2DVD in Greeley, CO (Easton), during the months from April to October 2015 formed the (spring–summer–fall) “climatological” database. The h(x) for each measured 3 min N(D) was calculated by normalizing using N0 and scaled using Dm. The median values of h(x) in each bin of width δx=0.05 were obtained and fitted to the G-G model through a weighted least-squares minimization leading to optimized values of μ=-0.24 and c=6.03 (see Raupach et al.2019 for details of the fitting procedure). Figure 7 shows empirical h(x) values as a frequency of occurrence plot on which the optimized G-G h(x) is overlaid.

Note that Thurai and Bringi (2018) and Raupach et al. (2019) allowed for μ to be negative in the G-G model, primarily to achieve a better fit for the small drop end of the DSD. The optimized values of c and μ fall within the range of values fitted to 1 min DSDs reported by Raupach et al. (2019). As a result, the analytical Eq. (42; L04) where Mk is derivable exactly in terms of [i,j;μ,c;k] cannot be used. Instead, Eq. (43) of L04, reproduced in (1) below, is employed. The radar estimates of the moments (Mk, k=0, 7) are obtained from the retrieved M3 and M6 and by numerical integration of the following function:

(1) h ( i , j , μ , c ) ( x ) = c Γ i ( j + c μ ) ( i - j ) Γ j ( - i - c μ ) ( i - j ) x c μ - 1 exp - Γ i Γ j c ( i - j ) x c ,

where Γi=Γμ+ic and Γj=Γμ+jc, with i=3 and j=6.

Figure 7The frequency of occurrence plot of h(x) from Greeley, CO, with overlay of the fitted G-G (μ=-0.24, c=6.03). The dashed black line is h(x) based on incomplete spectra using 2DVD data only. Note that the y axis is on a log axis, and therefore many zeros for large values of x are not shown but still affect the per-class median values to which the fits are made.


4.2 Retrieval algorithms

The retrieval algorithms for the reference moments {M3,M6} are based on 2928 3 min averaged complete DSDs from GXY and HSV. The combined DSDs from both locations are used because the frequency of occurrence of significant values of Ah (>1 dB km−1) from GXY alone was not enough to get a good retrieval. The scattering model assumptions are as given earlier in Sect. 3.3. For retrieval of M6 (in mm6 m−3), the obvious choice is Zh, and a power-law fit was derived for three ranges of ZH:


The above ranges of ZH were based on trial and error to minimize the parameterization errors (Fig. 8a). The slight decrease in the exponent from about 1 to about 0.8 as ZH increases is because of the effect of Mie scattering at X-band. The Zh in the above fits is in units of mm6 m−3, and so is M6.

The retrieval of M3 is based on a new multi-step procedure. First, the parameter Dm=(M6/M3)1/3 is retrieved from ZDR, which is reasonable because ZDR is weighted by the axis ratio of the large drops in the distribution and Dm and Dm are related to the drop size. A smoothing spline fit is used as shown in Fig. 8b. Again, the intent was to reduce parameterization errors as much as possible. The spline yields a visibly excellent fit with the Dm0.35 mm as ZDR→0 dB. Next, the Dm=M4/M3 is retrieved from Dm from a DSD-derived linear fit as Dm=0.08+0.8Dm.

Figure 8(a) Retrieval of M6 as a power law of Zh as per Eq. (2); each data point is based on a 3 min averaged complete spectrum from either the Greeley, CO, or Huntsville, AL, site. The simulations of X-band Zh, ZDR, and Ah are based on assumptions in Sect. 3.3. (b) Retrieval of Dm from ZDR along with smoothed cubic spline fit. (c) Same but for retrieval of AhW from Dm, where W is the rain water content. (d) The retrieved M3 versus “true” M3. (e) Box plots of RB for retrieval of M3 and M6 (which is a measure of the deviation of the fitted values from the “true” values because of DSD variability). The inter-quartile range is given by the “height” of the dark blue box, while the orange horizontal line inside the blue box is the median. The outliers are shown as orange circles; <10 % are estimated to be outliers. The orange horizontal lines on top and bottom indicate the upper and lower outer fences, respectively.


The next step is to retrieve AhW from Dm (adapted from Jameson1993, who used a third-order polynomial fit) for which we employ a smoothed spline fit. Here Ah is in dB km−1 and W is the rain water content in gm−3. We restricted the range of AhW at X-band to between 0.02 and 2 to avoid outliers. The smoothing spline fit is shown in Fig. 8c, which again provides a visibly good fit and is robust if the Dm falls outside the specified range. The retrieval of M3 follows from

(3) M 3 = 6000 π W = 6000 π A h f ( D m ) ,

where f(Dm) is the spline fit shown in Fig. 8c. The scatter plot of retrieved M3 versus “true” M3 is shown in Fig. 8d. This multi-step procedure was devised to minimize the parameterization (or algorithm) errors, but we note it is by no means the only way to achieve this.

It is known that the absorption cross section (specifically for X-band used here) depends on the temperature T via the 𝕀m{εr}, where εr is the dielectric constant of water. For a given W, the integral of the extinction cross section weighted by N(D) or Ah increases with colder water temperature but also depends on Dm (Jameson1993). Scattering simulations were performed for 8 and 20C and the spline fits of AhW versus Dm were compared. For low values of Dm, the maximum difference was 35 %, occurring at Dm=0.75 mm (with AhW larger at 8 relative to 20C as expected), but a crossover occurs near Dm=1.8 mm and the deviations increase in the opposite direction, with a maximum deviation of −15 % at Dm=3 mm (AhW at 20C larger than at 8C due to scattering loss). Recall that the National Weather Service (NWS) sounding at Denver about 65 km away showed a surface T of 12C. A lower temperature of 8C was used in the scattering calculations to approximately account for cooling of the atmosphere near Easton due to rainfall. The other temperature dependence is the coefficient α in the relation Ah=αKdp used in the iterative ZPHI method. This method involves finding an optimized α for each beam and is assumed to account for temperature changes. Since the actual drop temperature is not known and the surface T of 12C is close to the assumed T of 8C, the spline fit shown in Fig. 8c is considered to be sufficiently accurate for the retrieval of M3. We note that Diederich et al. (2015) found that the R(Ah) relation at X-band had a relatively “weak” dependence on temperature. Their fitted power law was 45.5Ah0.83 at 10C to 43.5Ah0.79 at 20C. At R=10 mm h−1, the Ah at 10C is larger than at 20C by 6.8 %, while, at 100 mm h−1, the Ah at 10C is lower than at 20C by −10.8 %; this crossover is consistent with our calculations above.

The evaluation of the algorithm error is done by defining the absolute bias of retrieved M, where M=M3 or M6, as Δ=(M(retrieved)M(“true”)) and the relative bias RB=100×Δ/M(“true”) as a percentage. To show the range of the RB and the distribution features (such as the median and 25th and 75th percentiles) in compact form, box plots for M3 and M6 are shown in Fig. 8e. The {25th,median,75th} percentile values for M3 and M6 are {-3.8,1.7,7.2} and {-3.8,0.63,6.2}, respectively. Note that the median relative bias is close to 0 and lies at the center of the box showing very low skewness. The interquartile range (IQR) is nearly the same for both M3 and M6. The extremities of the blue boxes are called hinges which span the IQR or the first (lower hinge) and third (upper hinge) quartiles. The orange line within the blue boxes indicates the median. The outliers (orange circles) lie beyond the first and third quartiles by at least 1.5 times the IQR. In particular, 1.5 and 3 times the IQR limit above (below) the upper (lower) hinge of a box are called the upper (lower) inner fence and upper (lower) outer fence, respectively (Theus and Urbanek2008). A point beyond an inner fence on either side is considered a mild outlier, while a point beyond an outer fence is an extreme outlier. The largest value below the upper inner fence and the smallest value above the lower inner fence are indicated by shorter grey horizontal lines called whiskers, within which lie extreme values that are not considered outliers. If there are no points beyond a whisker, corresponding inner and outer fences are not plotted. Similarly, if there are no samples between the inner and outer fences, only the inner fence is shown in the plot. Otherwise, the inner fence is generally omitted, and only the outer fence is depicted. For example, Fig. 8e shows only outer fence lines on top and bottom. While plotting multiple box plots in the same figure, only a common fence line that is closest and outside of all boxes is shown. The number of outliers for M3 is only 6.5 % of the total number of samples, while for M6 it is 9.5 %. It is interesting to note that the mean RBs for M3 and M6 are, respectively, −1.1 % and 0 %, which are close to the median values, meaning that there is low skewness in the relative error distributions.

Histograms of Δ∕〈M showed Gaussian-like shapes (not shown here). The variances of Δ normalized by M2 were 0.106 and 0.606 for M=M3 and M=M6, respectively; the corresponding fractional standard errors (FSEs) were, respectively, 0.32 and 0.778. These variances are referred to as variances due to parameterization or retrieval algorithm errors which can be added to the variances of the corresponding radar measurement errors to arrive at the total error variances. It is demonstrated (in the Appendix) that the retrieval algorithm error for M6 dominates the total error variance, whereas for M3 the retrieval algorithm and radar measurement errors are comparable.

5 Validation of radar-retrieved moments

The validation procedure essentially follows methods already developed for comparing radar-retrieved rain rates with disdrometers or gages (e.g., Bringi et al.2011b). The mean Lee-filtered Zh, Zdr, and Ah time-series data (see Fig. 6) were used to retrieve time series of {M3,M6}. Using Eq. (1) with the climatological hGG(x;μ,c) shown in Fig. 7, the other moments (0 through 2, 4 through 5, and 7) were computed. Retrievals and performance statistics were calculated for 20:00–21:30 UTC. The period after 21:30 UTC was omitted from this analysis because of the heavy rain observed on radome during the last half-hour of the event.

Figure 9Time series of radar-derived moments and from complete DSDs over the disdrometer site. The radar estimates are mean ±1σ error bars. (a) Moment M0; (b)(c) same but for M3 and M6.


Figure 9 shows the time series of radar-retrieved M0, M3, and M6 with those calculated from the 3 min complete DSDs. The radar retrievals in Fig. 9 show the mean with ±1σ error bars, where σ is the standard deviation. The mean value at each time step is obtained from the Lee-filtered values of {Zh,Zdr,Ah} which are used to retrieve the {M3,M6}. Then, using Eq. (1), the radar retrievals of M0 and other moments are obtained. The error bars or the variances consist of the sum of two terms, namely the parameterization error variances (described above) and the radar measurement errors which are uncorrelated. The Appendix describes the procedure to estimate the total error variances for moments such as M0 and the other non-reference moments in terms of the total error variances of M3 and M6. The last column in Table A1 gives the normalized total error variances for each moment. In Fig. 9, the standard deviation is obtained at each time step by taking the square root of the normalized variances (or the FSE) from the last column of Table A1 for M0, M3, and M6, with respective FSE values being {0.385,0.535,0.805}. The σ at each time step in Fig. 9 is calculated by multiplying the radar-retrieved M0, M3, and M6 at each time step by the corresponding FSE.

Figure 9a illustrates the intercomparison of M0 which is the most difficult to estimate using moments {M3,M6} (see Morrison et al.2019, RBa, and Raupach et al.2019). The error bars on the radar estimates are the total errors with FSE =0.385; see Appendix. The agreement with “ground truth” is visually quite remarkable considering that other error sources such as attenuation correction or point-to-area variance have been neglected (Ciach and Krajewski1999; Bringi et al.2011b). The total concentration (M0) in this event ranges from 100 per m3 to 100000 per m3 at the time of peak rainfall over Easton at 20:45 UTC. Figures 9b and c show time series of M3 and M6, respectively. The M3 retrievals are in excellent agreement with “ground truth” (total FSE =0.535 with 63 % of the variance due to measurement error and 37 % due to algorithm error). Figure 9c compares M6, and now the agreement degrades slightly. But the error bars have also increased substantially (total FSE =0.805 with nearly 93 % of the variance due to algorithm errors). The M6 varies from 10 to 105 mm6 m−3 (or equivalently 10 to 50 dBZ), the peak value at 20:45 UTC being in excess of 50 dBZ. Note that M0 and M3 are the moments prognosed by “bulk” double-moment numerical schemes (actually M0 and mass mixing ratio). So, radar retrievals could potentially play a role in evaluating the microphysical parameterizations in such models (e.g., Meyers et al.1997).

Figure 10Scatter plot of radar-derived moments versus “true” moments from the complete DSD data on log 10 scales. (a) M0, (b) M3, and (c) M6.


Table 2Statistics of the RB (median and IQR range). The correlations are between radar-retrieved moments and directly computed moments from the complete DSD measured by disdrometers.

Download Print Version | Download XLSX

The scatter plots of M0, M3, and M6 are shown in Fig. 10. The high correlation is quite striking and substantiated by both Pearson's and Spearman's rank correlation values in Table 2. Note that Spearman's measures only the monotonic relationship between the correlated variables, while Pearson's provides a measure of both monotonicity and linearity. The RB (%) was defined in Sect. 4.2; for each moment (M0 through M7) the corresponding box plots are shown in Fig. 11. We note that there are very few or no outliers for most of the moments. Table 2 gives the median (%) and the IQR range. The IQR range is smallest for M3. This is expected because it is one of the reference moments. The median of the RB is “best” for M5, with a symmetric IQR range indicating very low skewness. The median RBs for moments M0 through M4 are around −15 %, but the skewness is significant for M0 and progressively less for M1 through M5. It might be unexpected that the retrieval of M6 being one of the reference moments is less accurate than M5. One possible reason is that, at X-band, the larger drops are resonant-sized and the ZH does not vary as M6 but is rather closer to M5 depending on the drop sizes. Figure 7a, in fact, shows that the fit for M6 has a smaller slope for ZH>37 dBZ because of resonant scattering. The median RB for M6 and M7 are <17 %, but the IQR indicates positive skewness (i.e., radar estimates are larger than “truth”).

The difficulty in retrieving M0 from higher-order moments {M3,M6} is clear from the box plot but nevertheless viable with relatively low median values and high correlation coefficients. However, the accuracy of all moment-order retrievals, and especially the lower order, strongly depends on the climatological shape of h(x) for x<0.75 reflecting the shape of the small drop end (concave up for negative μ). This is irrespective of well-constrained measurement and parameterization errors in the retrieval of the reference moments {M3,M6}.

To illustrate this further, the “incomplete” spectra from 2DVD data alone, which are known to underestimate the numbers of small drops, are used to establish the “climatological” h(x) for which the fitted G-G shape parameters are μ=0.54 and c=3.07 (see Fig. 7). The radar moment retrieval steps are the same as before except for the now different h(x). The “true” moments are the same as before, being based solely on the complete DSD spectra. The new box plots of RB are shown in Fig. 12. Note that now the lower-order moments (M0 through M2) are severely underestimated, with median RB slightly less than −100 %, but the IQR is highly compressed, reflecting a distribution of RB which is concentrated as a delta function. The median RB for moments M6 and M7 is very large, >400 % with large IQR. Moments M3 through M5 show more “normal” RB distributions with median values of RB in the range −60 % to 90 %, the minimum occurring for M4. However, Pearson's correlation coefficients, as shown in Table 3, are very low for all moments, implying (practically) no linear variation between the moments. Spearman's rank correlation is higher for moments M4 through M7, implying that a nonlinear monotonic relationship between moments probably exists. The results in Table 3 vis-à-vis Table 2 demonstrate the importance of determining the “climatological” h(x) using the complete DSDs for accurate radar-based retrievals of the lower-order moments. In Table 3, the floor for RB is >-1.

Table 3As in Table 2 except h(x) is from incomplete 2DVD DSDs only.

Download Print Version | Download XLSX

Finally, the radar retrievals are examined from the perspective of identifying coherent “time tracks” as the main echo traversed the Easton site. To this end, Fig. 13 shows tracks in (a) the DmM0 plane, (b) the DmW plane, and (c) the DmM6 plane from 20:12 to 21:12 UTC. Note that we use Dm in panel (c) because it is more closely related to Zdr. For example, panel (a) shows the initial rapid rise in Dm from 1.5 to 2.2 mm (data point numbers 3 to 8 or approximately 2030–2045) with a corresponding increase in M0 (total number concentration) from 1000 to nearly 100 000 per m3. Together with similar behavior in panel (b) where W increases from ≤2 to 6 g m−3 and in (c) where M6 increases from 30 to >50 dB, this suggests the strong echo aloft descending to the surface over Easton. This inference was based on examining successive volume scans from KFTG (WSR-88D in Denver, CO, located about 60 km away) and noting the descent of the echo aloft to the surface at 20:45 UTC. After the peak, the track (data points 8 to 11 or 20:45 to 20:57 UTC) reflects a rapid decrease in M0 and Dm and from panel (b) a rapid decrease in W with a corresponding rapid decrease in M6, reflecting advection of the rain shaft to the north of Easton. Towards the end (last five data points from 20:57 to 21:12 UTC) the Dm decrease is slowed down (from 1.5 to 1 mm), while the M0 increases modestly from 10 000 to 15 000 per m3. This compensatory effect results in the rain water content being more or less steady (at 0.5 g m−3), while M6 decreases from 37 to 33 dB (panels b and c, last five data points). The echo structure during this latter time period was transitioning from earlier descent of the strong echo over Easton to more of a steady rain event.

Figure 11Box plots of relative bias (RB) as in Fig. 8e except between radar-derived moments and completed DSD moments (“truth”). The orange horizontal line on top indicates the upper inner fence.


Figure 12As in Fig. 11 except that the incomplete 2DVD DSDs are used to determine h(x) (see dashed line in Fig. 7). Inset shows magnified box plots for M0 through M2. The orange horizontal line on top indicates the upper inner fence.


6 Discussion

The polarimetric radar retrieval and validation of the lower-order moments of the DSD have not received much attention in the past, except for RBb. However, substantial literature exists in using either the un-normalized gamma model of Ulbrich (1983) or the normalized gamma model of Testud et al. (2001) to estimate the three parameters {N0,μ,Λ} or {Nw,μ,Dm}. As shown by L04, the Testud et al. (2001) formulation falls into double-moment normalization with reference moments M3 and M4, with h(x) being a special case of the G-G with c=1; hence, there is only one shape parameter μ. Note that μ=μULB+1, where μULB is the shape parameter defined in Ulbrich (1983).

Many studies have attempted to retrieve the three parameters {Nw,μ,Dm} using polarimetric measurements {Zh,Zdr,Kdp} at S-, C-, and X-bands, but they are too numerous to discuss herein (e.g., Bringi et al.2003; Brandes et al.2003; Park et al.2005b; Gorgucci et al.2008; Anagnostou et al.2013, to mention a few). Anagnostou et al. (2013) compared three different methods of retrieving Nw but found that validation was very difficult, commenting that “ the estimation of Nw by all algorithms is significantly affected by noise or other factors like radar volume versus point (disdrometer) measurement-scale mismatch and spatial separation.”

However, neither N0 nor Nw is the same as M0, which is simply the total number concentration that scales the gamma pdf. The estimation of either N0 (or Nw) depends on the shape parameter (or the slope parameter Λ). Typically, the μ is assumed to be fixed or empirically derived as f(Λ) or another function of Dm (Schinagl et al.2019). Of course, any moment of the gamma pdf can be derived as functions of the three parameters in the gamma model, but very few validations of, for example, M0 have been conducted. Brandes et al. (2003) used an empirically derived μΛ relation based on 2DVD data. Using Zh and Zdr radar measurements at S-band, they retrieved N0 and Λ and obtained the M0 as N0Λ-(μ+1). They analyzed one convective event (with ZH varying from 10 to 55 dBZ) and showed the mean of log 10(M0) from radar was 2.74 compared with 2.84 from 2DVD measurements (or 550 and 690 per m3, which are much smaller than the values obtained here; see Fig. 10a). The key point is that the mean μ was in the range 3–4, which is caused by truncation at the small drop end of the DSD.

Figure 13Time tracks of radar-derived variables in the (a) Dm versus M0 plane showing the trajectories as a function of time (color coded) over approximately an hour. Each data point reflects the radar-retrieved moments (M0) or ratio of moments (M4M3) (see Fig. 9). (b)(c) Same but time tracks in the Dm versus W plane and Dm versus M6 planes, respectively. Again, all quantities are from radar-retrieved moments.


Wen et al. (2018) describe a different method of estimating the DSD parameters of the gamma distribution and the lower-order moments based on an inverse model where the input is {Zdr,Kdp/Zh} and the output is {μ,Dmax}, where Dmax is the maximum diameter of the retrieved gamma DSD. Their approach follows the well-known k-nearest-neighbor (k-NN) classification from the pattern recognition literature (Shakhnarovich et al.2006). This algorithm stores all input–output associations from the available data as a “training” set. When a new {Zdr,Kdp/Zh} input is presented, the algorithm assigns it the {μ,Dmax} output class that is the most common amongst the k nearest (training set) neighbors of the new input. The k-NN is particularly suitable when large training data are available. Wen et al. (2018) used Euclidean distance to define the closeness of neighbors, although other distance functions are also employed in k-NN algorithms. They applied an empirical μΛ relation based on 2DVD data, while N0 is obtained a posteriori using Zh, μ, Λ, and Dmax. Their training set comprising Zdr and KdpZh was generated using a polynomial function whose inputs μ and Dmax are drawn from 10-year disdrometer data with constraints μ[3,20], Dmax[1.7mm,8mm], and Dmax>Dm. The test stage used S-band radar data from a WSR-88D unit (KTLX) located in Oklahoma City, OK. A large database was analyzed, and the moments M0, M2, M4, and M6 were computed from {N0,μ,Λ,Dmax}. The validation results in terms of what they define as relative absolute error (RAE) ranged from 0.986 (or 98.6 %) for M0 to 0.455 (or 45.5 %) for M6, while Pearson's correlation coefficient between radar-based M0 and 2DVD M0 “truth” was 0.651 (the maximum correlation coefficient for other moments was <0.7). The predictive performance of k-NN was quantified through root relative squared error (RRSE), which computes the difference between the k-NN-predicted values with the actual ones relative to when a simple predictor is used. More than characterizing the accuracy of computation of moments, both RAE and RRSE give an indication of the efficacy of k-NN-based prediction over the most basic mean-value prediction method. Wen et al. (2018) reported low RRSE for M2, M4, and M6, whereas it was large (>1) for M0. They commented that “the inverse model [] produced DSD retrievals with large uncertainties due to the measurement errors, noise, and sampling problems of the instruments.”

The RBb article used X-band radar Zh to retrieve M6 and {Zdr,Kdp} to retrieve M3. Our approach is similar except that we use Ah instead of Kdp. There are several advantages to using Ah (similar to its use in estimating R at X-band; Diederich et al.2015). For instance, Ah is always positive, is highly correlated with Zh variations, thus preserving the spatial resolution, and, at X-band, has decent dynamic range. The article used several networks of the Parsivel disdrometer from three locations to derive h(x) and the G-G fit, but their fit was more similar to that shown in Fig. 7 (black dashed line); in fact, they obtained a larger μ (2.22) and smaller c (1.69). The higher value of μ gives a more convex down shape at small x (relative to black dashed line in Fig. 7), while smaller c results in slower fall of the tail of the distribution. The other issue they had to deal with was the “noisy” Zdr and Kdp measurements when ZH<37 dBZ. They “restored” the noisy Zdr by estimating it using a power law with Zh, while noisy Kdp was restored using power laws of Zh and Zdr. They also commented that the majority of radar-measured ZH was <37 dBZ. So, noise correction dominated the statistics of their moment retrievals. Their radar retrievals of moments were based on a large dataset from three regions (two in Europe and one in Iowa, US). While they obtained median values of RB in the range 4 % to −46 %, their r2 (squared Pearson's correlation) coefficient between radar moments and ground “truth” was very low (0.05 to 0.33), similar to what we obtained in Table 3. They ascribed their poor correlation to spatial representativeness errors, height of the radar pixels above the Parsivel network at longer ranges, and other factors, similarly to Anagnostou et al. (2013).

7 Summary

We demonstrated a proof-of-concept of the viability of radar retrieval of lower-order moments of the DSD using specific attenuation Ah in addition to Zh and Zdr at X-band (an extension of RBb) via a case study approach. The use of specific attenuation (from the iterative ZPHI method) is consistent with its many advantages for rain rate estimation. The multi-cell convective complex which occurred in the area near Greeley, CO, on 23 May 2015 was a target of opportunity as the CSU-CHILL radar system was available to scan the echo complex with a single elevation angle PPI every 90 s over a period of around 90 min. The instrumented site at Easton located 13 km to the south of the radar had an MPS and 2DVD sited inside a DFIR wind shield which made it possible to acquire the “complete” drop spectra with high resolution (50 µm) for the small drop end and good resolution (about 170 µm) for drops ≥0.7 mm. The moment retrieval was based on the double-moment scaling/normalization framework of Lee et al. (2004). Two reference moments {M3,M6} along with a “climatological” estimate of the underlying shape of the scaled/normalized DSDs fitted to the G-G distribution formed the basis of the method. The {M3,M6} retrieval algorithms were trained using scattering simulations of Zh, Zdr, and Ah using 2928 3 min averaged DSDs from Greeley, CO, and Huntsville, AL. The parameterization (or algorithm) errors due to DSD variability about the smoothed spline fits were computed.

Polarimetric X-band radar data (acquired with an exceptionally narrow 3 dB beamwidth of 0.33) were extracted from a small polar box surrounding the instrumented site, and moments M0 through M7 were estimated and validated against ground “truth” from the moments of the complete spectra using MPS and 2DVD. Using a variety of validation measures such as box plots of RB, time-series comparisons, scatter plots, and correlation coefficients, it was determined that good accuracy was obtained for the radar-based moments well beyond what has been possible hitherto. For moments M0 through M2, the RB was <15 % in magnitude, with Pearson's correlation coefficients between radar-derived moments and DSD-based moments exceeding 0.9. A detailed analysis of radar fluctuations or measurement errors propagating to the variance of the moment estimates was performed; in addition, the total variance due to both parameterization and measurement errors was tabulated.

The coherency of “time track” plots of radar-retrieved quantities in the Dm versus M0, Dm versus W, and Dm versus M6 planes as the main 55 dBZ echo passed over the site (as well as 20 min prior to and 20 min after this passage) demonstrated the potential use for precipitation evolution studies for this DSD-retrieval technique. One caveat is that a much larger database is needed before concrete conclusions are drawn. In particular, the possibility of the very narrow beam of 0.33 and the close range (13 km) to the instrumented site contributing to very good validation statistics, found herein relative to other studies, requires investigations with more data.

Appendix A

The error model (Bringi and Chandrasekar2001) we adopt here is an additive one, X^=X+εm+εp, where X^ is the estimated (or retrieved) quantity, X is the “true” value, and εm and εp are, respectively, the radar measurement and parameterization (or algorithm) errors. The εm and εp are zero mean, uncorrelated errors so that E{X^}=X. Thus, it follows that Var(X^-X)=Varεm+Varεp. For different rain rate estimators such as R(Zh), R(Zh,Zdr), and R(Kdp), the Var(R^)/R2 is expressed in terms of the standard deviations of Zh, Zdr, and Kdp which are 1 dBZ, 0.3 dB, and 0.3 km−1, respectively (Thurai et al.2017b).

Errors due to attenuation correction are not considered because most of the attenuation occurred at ranges beyond the instrumented site (see Fig. 5). We refer to Thurai et al. (2017b) for evaluation of such errors.

A1 Radar measurement errors

We consider error variances of the retrieval of M6 and M3 first and then the other moments. Since M6 is retrieved as a power law of Zh0.8 (the exponent is approximate),

(A1) Var ( M 6 ) M 6 2 = 0.8 2 Var ( Z h ) Z h 2 ,

where Zh has units of mm6m−3. Assuming the standard deviation of the radar measurement error is typically 1 dB, we get Var(Zh)/Zh2=0.067. This implies that

(A2) Var M 6 / M 6 2 = 0.043 .

The variance of M3 is more complicated because it is a multi-step procedure as described in Sect. 4 involving smoothed spline fits. We use approximate power-law fits for estimating Var(M3) as follows. We have

(A3) D m 1.18 ( Z dr ) 1.5 ,


(A4) A h W 0.09 D m 2 , for D m > 1 mm .


(A5) Var ( M 3 ) M 3 2 = Var ( A h ) A h 2 + 4 Var ( D m ) D m 2 .

Since Dm is linear with Dm, and assuming the standard deviation of the radar measurement of ZDR is 0.3 dB,

(A6) Var ( M 3 ) M 3 2 = Var ( A h ) A h 2 + 9 Var ( Z dr ) Z dr 2 ,

where Zdr is a ratio and Var(Zdr)/Zdr2=0.0051:

(A7) Var ( M 3 ) M 3 2 = Var ( A h ) A h 2 + 0.0046 .

From Thurai et al. (2017b, Appendix, A5),

(A8) Var A h A h 2 = 0.8 2 Var ( Z h ) Z h 2 + Var K dp K dp 2 ,

assuming that Ah varies as Zh0.8 used in the ZPHI method. The standard deviation of the Kdp measurement is typically 0.3 km−1 and that of ZH is 1 dB. Further, the mean Kdp for our data ≈1 km−1 (Fig. 6c) but is variable in time. In any case,

(A9) Var A h A h 2 ( 0.64 ) ( 0.067 ) + 0.09 = 0.133 .

Substituting above in Eq. (A6) yields

(A10) Var ( M 3 ) M 3 2 = 0.133 + 0.046 = 0.18 .

A2 Variances of the other moments

From Lee et al. (2004), the other moments Mk can be expressed as power laws of the reference moments M3 and M6. They are of the form Mk=CkM3pkM6-qk, where pk and qk are rational numbers and Ck is some constant. The variance of Mk needs more elaboration as XM3 and YM6 are correlated. This correlation arises because Ah is a power law with Zh, i.e., Zh0.8 in the ZPHI method. Together with M6 being a power law with Zh, the M3 and M6 are correlated with a correlation coefficient of 0.93 obtained from radar-derived M3 and M6. For the kth moment Mk, pk=6-k3 and qk=3-k3 for k=0,1,,7, k3,6. The objective is to derive Var(Mk)Mk2 in terms of Var(M3)M32, Var(M6)M62, and Cov(M3,M6).

In the sequel, for notational simplicity, we drop the subscripts k. Then, for certain rational numbers p and q, any moment M is a function of these two random variables as

(A11) M f M 3 , M 6 = C M 3 p M 6 q = C X p Y q .

Consider the parameter vector θ=(X,Y). Then, second-order Taylor series approximation of f(X,Y) around θ produces

(A12) M f ( θ ) + f x ( θ ) ( X - X ) + f y ( θ ) ( Y - Y ) + 1 2 f x x ′′′ ( θ ) ( X - X ) 2 + 2 f x y ′′ ( θ ) ( X - X ) ( Y - Y ) + f y y ′′ ( θ ) ( Y - Y ) 2 ,

where the notations fx(θ) and fxy′′(θ) represent, respectively, the first- and second-order derivatives of the function f with respect to the variables in the subscript and evaluated at θ.

The second-order approximation of the mean M=E{M} is

(A13) M E f ( θ ) + f x ( θ ) ( X - X ) + f y ( θ ) ( Y - Y ) + 1 2 f x x ′′′ ( θ ) ( X - X ) 2 + 2 f x y ′′ ( θ ) ( X - X ) ( Y - Y ) + f y y ′′ ( θ ) ( Y - Y ) 2 = E f ( θ ) + f x ( θ ) E ( X - X ) + f y ( θ ) E ( Y - Y ) + 1 2 f x x ′′′ ( θ ) E ( X - X ) 2 + 2 f x y ′′ ( θ ) E ( X - X ) ( Y - Y ) + f y y ′′ ( θ ) E ( Y - Y ) 2 = f ( θ ) + 1 2 f x x ′′′ ( θ ) Var ( X ) + 2 f x y ′′ ( θ ) Cov ( X , Y ) + f y y ′′ ( θ ) Var ( Y ) ,

where the second equality results because E{(X-X)}=E{(Y-Y)}=0. Note that, from Eq. (A13), the first-order approximation of the mean is simply Mf(θ). Evaluating the function f(⋅) and its derivatives at θ returns


Substituting Eqs. (A14)–(A17) in Eq. (A13) yields the following approximation of the mean:

(A18) M C X p Y q + p ( p - 1 ) 2 X p - 2 Y q Var ( X ) - p q X p - 1 Y q + 1 Cov ( X , Y ) + q ( q + 1 ) 2 X p Y q + 2 Var ( Y ) = C X p Y q 1 + p ( p - 1 ) 2 Var ( X ) X 2 - p q Cov ( X , Y ) X Y + q ( q + 1 ) 2 Var ( Y ) Y 2 .

In order to compute the expression of Var(M), we note that, by definition,

(A19) Var ( M ) = E ( M - M ) 2 E ( M - f ( θ ) ) 2 ,

where we have used the first-order approximation of the mean M. Then, ignoring all the terms above second order, the Taylor series expansion of M around θ gives

(A20) Var ( M ) E f ( θ ) + f x ( θ ) ( X - X ) + f y ( θ ) ( Y - Y ) - f ( θ ) 2 = E f x ( θ ) ( X - X ) + f y ( θ ) ( Y - Y ) 2 = E f x 2 ( θ ) ( X - X ) 2 + 2 f x ( θ ) f y ( θ ) ( X - X ) ( Y - Y ) + f y 2 ( θ ) ( Y - Y ) 2 = f x 2 ( θ ) Var ( X ) + 2 f x ( θ ) f y ( θ ) Cov ( X , Y ) + f y ( θ ) Var ( Y ) .

Again, evaluating at θ, we obtain


Substituting the above in Eq. (A20) leads to the first-order approximation

(A24) Var ( M ) C 2 p 2 X 2 ( p - 1 ) Y 2 q Var ( X ) - 2 p q X 2 p - 1 Y 2 q + 1 Cov ( X , Y ) + q 2 X 2 p Y 2 ( q + 1 ) Var ( Y ) = C 2 X 2 p Y 2 q p 2 Var ( X ) X 2 - 2 p q Cov ( X , Y ) X Y + q 2 Var ( Y ) Y 2 .

From Eqs. (A18) and (A24), we obtain the desired ratio as

(A25) Var ( M ) M 2 p 2 Var ( X ) X 2 - 2 p q Cov ( X , Y ) X Y + q 2 Var ( Y ) Y 2 1 + p ( p - 1 ) 2 Var ( X ) X 2 - p q Cov X , Y X Y + q ( q + 1 ) 2 Var ( Y ) Y 2 2 .

If the correlation coefficient ρXY between X and Y is known, then we replace Cov(X,Y) in Eq. (A25) to obtain

(A26) Var ( M ) M 2 p 2 Var X X 2 - 2 p q ρ X Y Var ( X ) X 2 Var ( Y ) Y 2 + q 2 Var ( Y ) Y 2 1 + p ( p - 1 ) 2 Var ( X ) X 2 - p q ρ X Y Var ( X ) X 2 Var ( Y ) Y 2 + q ( q + 1 ) 2 Var ( Y ) Y 2 2 .

Using Eq. (A26), Table A1 gives the Var(Mk)Mk2 for radar measurement errors Var(M3)M32=0.18 and Var(M6)M62=0.043. From Sect. 4.2, the parameterization errors are Var(M3)M32=0.106 and Var(M6)M62=0.606. The total variances are obtained by adding these parameterization errors to the radar measurement errors, respectively, 0.18 and 0.043 to get 0.286 and 0.649. Thus, for M^6 the parameterization error dominates with 93 % of the total variance, whereas for M^3 the measurement error dominates with 63 % of the total. Using Eq. (A26) and total variances for M^3 and M^6 gives in Table A1 the total variances for the other moments. For M0 through M2, the total variance is smaller than the measurement variance because the covariance term is negative for those moments.

Table A1Variance estimates for the radar-based retrievals of moments M0 through M7. The p and q are the exponents of Mk=M3pM6-q. The normalized variances due to measurement errors are given in the 4th column, while the totals due to measurement and parameterization errors are given in the 5th column.

Download Print Version | Download XLSX

Code availability

The IDL, MATLAB, Fortran, and R codes used in this article are available upon request from the corresponding author.

Data availability

The CSU-CHILL radar, MPS, and 2DVD processed data are available upon request to the corresponding author.

Video supplement

Animation of the radar PPI sweeps for the entire duration as well as the composite DSD spectra can be found in the same web article at (last access: 27 August 2020; Thurai et al., 2015).

Author contributions

Conceptualization was by VB, MT, and THR. Methodology, investigation, and formal analysis were done by VB, KVM, and MT. Data curation was by PCK. Radar analysis was by PCK, MT, and VB. Writing and original draft preparation was by VB and KVM. Writing, review, and editing were done by VB, KVM, PCK, and THR. Supervision was done by VB.

Competing interests

The authors declare that they have no conflict of interest. The funders had no role in the design of this study; in the collection, analyses, or interpretation of its data; in the writing of this paper; or in the decision to publish these results.


Merhala Thurai and Viswanathan Bringi received funding to conduct this research from the National Science Foundation under grant AGS-1901585. The CSU-CHILL radar was made available via a short 20 h project approved by the scientific director Steven Rutledge. The DFIR wind shield which housed the MPS and 2DVD was built at the Easton site near the CSU-CHILL radar under a prior NSF grant (P.I. B. Notaroš). The 2DVD and Pluvio gage were graciously loaned to Colorado State University by Walter Petersen of the NASA/Wallops Precipitation Research Facility. The MPS and 2DVD at Huntsville, AL, continue to be maintained by Patrick Gatlin and Matt Wingo of NASA.

Financial support

This research has been supported by the National Science Foundation (grant no. AGS-1901585).

Review statement

This paper was edited by Gianfranco Vulpiani and reviewed by two anonymous referees.


Abel, S. J. and Boutle, I. A.: An improved representation of the raindrop size distribution for single-moment microphysics schemes, Q. J. Roy. Meteorol. Soc., 138, 2151–2162, 2012. a

Anagnostou, E. N., Krajewski, W. F., and Smith, J.: Uncertainty quantification of mean-areal radar-rainfall estimates, J. Atmos. Ocean. Technol., 16, 206–215, 1999. a

Anagnostou, E. N., Anagnostou, M. N., Krajewski, W. F., Kruger, A., and Miriovsky, B. J.: High-resolution rainfall estimation from X-band polarimetric radar measurements, J. Hydrometeorol., 5, 110–128, 2004. a

Anagnostou, M. N., Kalogiros, J., Marzano, F. S., Anagnostou, E. N., Montopoli, M., and Piccioti, E.: Performance evaluation of a new dual-polarization microphysical algorithm based on long-term X-Band radar and disdrometer observations, J. Hydrometeorol., 14, 560–576, 2013. a, b, c

Barber, P. and Yeh, C.: Scattering of electromagnetic waves by arbitrarily shaped dielectric bodies, Appl. Opt., 14, 2864–2872, 1975. a

Baumgardner, D., Abel, S., Axisa, D., Cotton, R., Crosier, J., Field, P., Gurganus, C., Heymsfield, A., Korolev, A., Krämer, M., Lawson, P., McFarquhar, G., Ulanowski, Z., and Um, J.: Cloud ice properties: In situ measurement challenges, in: Ice formation and evolution in clouds and precipitation: Measurement and modeling challenges, edited by: Baumgardner, D. D., McFarquhar, G. M., and Heymsfield, A. J., Meteorological Monographs (Book 58), p. 320, American Meteorological Society, 2017. a

Belda, M., Holtanová, E., Halenka, T., and Kalvová, J.: Climate classification revisited: From Köppen to Trewartha, Clim. Res., 59, 1–13, 2014. a

Bernauer, F., Hürkamp, K., Rühm, W., and Tschiersch, J.: On the consistency of 2-D video disdrometers in measuring microphysical parameters of solid precipitation, Atmos. Meas. Tech., 8, 3251–3261,, 2015. a

Brandes, E. A., Zhang, G., and Vivekanandan, J.: An evaluation of a drop distribution–based polarimetric radar rainfall estimator, J. Appl. Meteorol., 42, 652–660, 2003. a, b

Bringi, V., Chandrasekar, V., Balakrishnan, N., and Zrnić, D.: An examination of propagation effects in rainfall on radar measurements at microwave frequencies, J. Atmos. Ocean. Technol., 7, 829–840, 1990. a

Bringi, V., Chandrasekar, V., Hubbert, J., Gorgucci, E., Randeu, W., and Schoenhuber, M.: Raindrop size distribution in different climatic regimes from disdrometer and dual-polarized radar analysis, J. Atmos. Sci., 60, 354–365, 2003. a

Bringi, V., Tolstoy, L., Thurai, M., and Petersen, W.: Estimation of spatial correlation of drop size distribution parameters and rain rate using NASA's S-band polarimetric radar and 2D video disdrometer network: Two case studies from MC3E, J. Hydrometeorol., 16, 1207–1221, 2015. a, b

Bringi, V., Thurai, M., and Baumgardner, D.: Raindrop fall velocities from an optical array probe and 2-D video disdrometer, Atmos. Meas. Tech., 11, 1377–1384,, 2018. a

Bringi, V. N. and Chandrasekar, V.: Polarimetric Doppler weather radar: Principles and applications, Cambridge University Press, 2001. a, b, c, d

Bringi, V. N. and Seliga, T. A.: Scattering from non-spherical hydrometeors, Annales Télécommun., 32, 392–397, 1977. a

Bringi, V. N., Keenan, T. D., and Chandrasekar, V.: Correcting C-band radar reflectivity and differential reflectivity data for rain attenuation: A self-consistent method with constraints, IEEE Trans. Geosc. Remote Sens., 39, 1906–1915, 2001. a, b

Bringi, V. N., Hoferer, R., Brunkow, D. A., Schwerdtfeger, R., Chandrasekar, V., Rutledge, S. A., George, J., and Kennedy, P. C.: Design and performance characteristics of the new 8.5-m dual-offset Gregorian antenna for the CSU–CHILL radar, J. Atmos. Ocean. Technol., 28, 907–920, 2011a. a, b

Bringi, V. N., Rico-Ramirez, M. A., and Thurai, M.: Rainfall estimation with an operational polarimetric C-Band radar in the United Kingdom: Comparison with a gauge network and error analysis, J. Hydrometeorol., 12, 935–954, 2011b. a, b

Brunkow, D., Bringi, V. N., Kennedy, P. C., Rutledge, S. A., Chandrasekar, V., Mueller, E. A., and Bowie, R. K.: A description of the CSU–CHILL National Radar Facility, J. Atmos. Ocean. Technol., 17, 1596–1608, 2000. a

Ciach, G. J. and Krajewski, W. F.: On the estimation of radar rainfall error variance, Adv. Water Resour., 22, 585–595, 1999. a

Cocks, S. B., Tang, L., Zhang, P., Ryzhkov, A., Kaney, B., Elmore, K. L., Wang, Y., Zhang, J., and Howard, K.: A prototype quantitative precipitation estimation algorithm for operational S-band polarimetric radar utilizing specific attenuation and specific differential phase. Part II: Performance verification and case study analysis, J. Hydrometeorol., 20, 999–1014, 2019. a

Delanoë, J. M. E., Heymsfield, A. J., Protat, A., Bansemer, A., and Hogan, R. J.: Normalized particle size distribution for remote sensing application, J. Geophys. Res.-Atmos., 119, 4204–4227, 2014. a

Diederich, M., Ryzhkov, A., Simmer, C., Zhang, P., and Trömel, S.: Use of specific attenuation for rainfall measurement at X-band radar wavelengths. Part I: Radar calibration and partial beam blockage estimation, J. Hydrometeorol., 16, 487–502, 2015. a, b, c

Gorgucci, E. and Chandrasekar, V.: Evaluation of attenuation correction methodology for dual-polarization radars: Application to X-band systems, J. Atmos. Ocean. Technol., 22, 1195–1206, 2005. a

Gorgucci, E., Chandrasekar, V., and Baldini, L.: Microphysical retrievals from dual-polarization radar measurements at X band, J. Atmos. Ocean. Technol., 25, 729–741, 2008. a, b

Huang, G.-J., Bringi, V., and Thurai, M.: Orientation angle distributions of drops after an 80-m fall using a 2D video disdrometer, J. Atmos. Ocean. Technol., 25, 1717–1723, 2008. a, b

Hubbert, J. and Bringi, V. N.: An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements, J. Atmos. Ocean. Technol., 12, 643–648, 1995. a

Illingworth, A. J. and Blackman, T. M.: The need to represent raindrop size spectra as normalized gamma distributions for the interpretation of polarization radar observations, J. Appl. Meteorol., 41, 286–297, 2002. a

Jameson, A.: Microphysical interpretation of multi-parameter radar measurements in rain. Part I: Interpretation of polarization measurements and estimation of raindrop shapes, J. Atmos. Sci., 40, 1792–1802, 1983. a

Jameson, A. R.: The meteorological parameterization of specific attenuation and polarization differential phase shift in rain, J. Appl. Meteorol., 32, 1741–1750, 1993. a, b, c

Junyent, F., Chandrasekar, V., Bringi, V. N., Rutledge, S. A., Kennedy, P. C., Brunkow, D., George, J., and Bowie, R.: Transformation of the CSU–CHILL Radar Facility to a dual-frequency, dual-polarization Doppler system, B. Am. Meteorol. Soc., 96, 975–996, 2015. a

Kalogiros, J., Anagnostou, M. N., Anagnostou, E. N., Montopoli, M., Picciotti, E., and Marzano, F. S.: Optimum estimation of rain microphysical parameters from X-band dual-polarization radar observables, IEEE Trans. Geosci. Remote Sens., 51, 3063–3076, 2012. a

Kennedy, P., Thurai, M., Praz, C., Bringi, V. N., Berne, A., and Notaroš, B. M.: Variations in snow crystal riming and ZDR: A case analysis, J. Appl. Meteorol. Climatol., 57, 695–707, 2018. a

Knollenberg, R.: Three new instruments for cloud physics measurements: The 2–D spectrometer probe, the forward scattering spectrometer probe and the active scattering aerosol spectrometer, in: American Meteorological Society International Conference on Cloud Physics, pp. 554–561, 26–30 July 1976, Boulder, Colorado, 1976. a

Knollenberg, R.: Techniques for probing cloud microstructure, in: Clouds, their formation, optical properties and effects, edited by: Hobbs, P. V. and Deepak, A., pp. 15–92, Academic Press, 1981. a

Knollenberg, R. G.: The optical array: An alternative to scattering or extinction for airborne particle size determination, J. Appl. Meteorol., 9, 86–103, 1970. a

Korolev, A., Kuznetsov, S., Makarov, Y. E., and Novikov, V.: Evaluation of measurements of particle size and sample area from optical array probes, J. Atmos. Ocean. Technol., 8, 514–522, 1991. a

Korolev, A., Strapp, J., and Isaac, G.: Evaluation of the accuracy of PMS optical array probes, J. Atmos. Ocean. Technol., 15, 708–720, 1998. a

Kumjian, M. R., Martinkus, C. P., Prat, O. P., Collis, S., van Lier-Walqui, M., and Morrison, H. C.: A moment-based polarimetric radar forward operator for rain microphysics, J. Appl. Meteorol. Climatol., 58, 113–130, 2019. a, b

Larsen, M. L. and Schönhuber, M.: Identification and characterization of an anomaly in two-dimensional video disdrometer data, Atmosphere, 9, 315, 2018. a

Lee, G. W., Zawadzki, I., Szyrmer, W., Sempere-Torres, D., and Uijlenhoet, R.: A general approach to double-moment normalization of drop size sistributions, J. Appl. Meteorol., 43, 264–281, 2004. a, b, c

Lee, J.: Digital image enhancement and noise filtering by use of local statistics, IEEE T. Pattern Anal., PAMI-2, 165–168, 1980. a, b

Low, T. B. and List, R.: Collision, coalescence and breakup of raindrops. Part II: Parameterization of fragment size distributions, J. Atmos. Sci., 39, 1607–1619, 1982. a

May, P., Keenan, T., Zrnić, D., Carey, L., and Rutledge, S.: Polarimetric radar measurements of tropical rain at a 5-cm wavelength, J. Appl. Meteorol., 38, 750–765, 1999. a

McFarquhar, G. M.: A new representation of collision-induced breakup of raindrops and its implications for the shapes of raindrop size distributions, J. Atmos. Sci., 61, 777–794, 2004. a, b

Meyers, M. P., Walko, R. L., Harrington, J. Y., and Cotton, W. R.: New RAMS cloud microphysics parameterization. Part II: The two-moment scheme, Atmos. Res., 45, 3–39, 1997. a, b

Milbrandt, J. A. and Yau, M. K.: A multimoment bulk microphysics parameterization. Part I: Analysis of the role of the spectral shape parameter, J. Atmos. Sci., 62, 3051–3064, 2005. a

Mishra, K. V., Krajewski, W. F., Goska, R., Ceynar, D., Seo, B.-C., Kruger, A., Niemeier, J. J., Galvez, M. B., Thurai, M., Bringi, V., Tolstoy, L., Kucera, P. A., Petersen, W. A., Grazioli, J., and Pazmany, A. L.: Deployment and performance analyses of high-resolution Iowa XPOL radar system during the NASA IFloodS campaign, J. Hydrometeorol., 17, 455–479, 2016. a

Morrison, H., Kumjian, M. R., Martinkus, C. P., Prat, O. P., and van Lier-Walqui, M.: A general N-moment normalization method for deriving raindrop size distribution scaling relationships, J. Appl. Meteorol. Climatol., 58, 247–267, 2019. a, b

Notaroš, B. M., Bringi, V. N., Kleinkort, C., Kennedy, P., Huang, G.-J., Thurai, M., Newman, A. J., Bang, W., and Lee, G.: Accurate characterization of winter precipitation using multi-angle snowflake camera, visual hull, advanced scattering methods and polarimetric radar, Atmosphere, 7, 81, 2016. a

Park, S.-G., Bringi, V. N., Chandrasekar, V., Maki, M., and Iwanami, K.: Correction of radar reflectivity and differential reflectivity for rain attenuation at X band. Part I: Theoretical and empirical basis, J. Atmos. Ocean. Technol., 22, 1621–1632, 2005a. a, b, c

Park, S.-G., Maki, M., Iwanami, K., Bringi, V. N., and Chandrasekar, V.: Correction of radar reflectivity and differential reflectivity for rain attenuation at X band. Part II: Evaluation and application, J. Atmos. Ocean. Technol., 22, 1633–1655, 2005b. a, b, c, d

Park, S.-G., Kim, H.-L., Ham, Y.-W., and Jung, S.-H.: Comparative evaluation of the OTT PARSIVEL2 using a collocated two-dimensional video disdrometer, J. Atmos. Ocean. Technol., 34, 2059–2082, 2017. a

Rasmussen, R., Baker, B., Kochendorfer, J., Meyers, T., Landolt, S., Fischer, A. P., Black, J., Thériault, J. M., Kucera, P., Gochis, D., Smith, C., Nitu, R., Hall, M. Ikeda, K., and Gutmann, E.: How well are we measuring snow: The NOAA/FAA/NCAR winter precipitation test bed, B. Am. Meteorol. Soc., 93, 811–829, 2012. a

Raupach, T. H. and Berne, A.: Invariance of the double-moment normalized raindrop size distribution through 3D spatial displacement in stratiform rain, J. Appl. Meteorol. Climatol., 56, 1663–1680, 2017a. a, b

Raupach, T. H. and Berne, A.: Retrieval of the raindrop size distribution from polarimetric radar data using double-moment normalisation, Atmos. Meas. Tech., 10, 2573–2594,, 2017b. a

Raupach, T. H., Thurai, M., Bringi, V. N., and Berne, A.: Reconstructing the drizzle mode of the raindrop size distribution using double-moment normalization, J. Appl. Meteorol. Climatol., 58, 145–164, 2019. a, b, c, d, e, f, g, h, i

Ryzhkov, A., Diederich, M., Zhang, P., and Simmer, C.: Potential utilization of specific attenuation for rainfall estimation, mitigation of partial beam blockage, and radar networking, J. Atmos. Ocean. Technol., 31, 599–619, 2014. a, b

Ryzhkov, A. V. and Zrnić, D. S.: Radar polarimetry for weather observations, Springer, 2019. a

Schinagl, K., Friederichs, P., Trömel, S., and Simmer, C.: Gamma drop size distribution assumptions in bulk model parameterizations and radar polarimetry and their impact on polarimetric radar moments, J. Appl. Meteorol. Climatol., 58, 467–478, 2019. a

Schönhuber, M., Lammer, G., and Randeu, W. L.: One decade of imaging precipitation measurement by 2D-video-distrometer, Adv. Geosci., 10, 85–90, 2007. a

Schönhuber, M., Lammar, G., and Randeu, W. L.: The 2D-video-distrometer, in: Precipitation: Advances in measurement, estimation and prediction, edited by: Michaelides, S. C., pp. 3–31, Springer, 2008. a, b

Shakhnarovich, G., Darrell, T., and Indyk, P.: Nearest-neighbor methods in learning and vision, MIT Press, 2006. a

Smyth, T. J. and Illingworth, A. J.: Correction for attenuation of radar reflectivity using polarization data, Q. J. Roy. Meteorol. Soc., 124, 2393–2415, 1998. a, b

Stacy, E. W.: A generalization of the gamma distribution, Ann. Mathe. Stat., 33, 1187–1192, 1962. a

Straub, W., Beheng, K. D., Seifert, A., Schlottke, J., and Weigand, B.: Numerical investigation of collision-induced breakup of raindrops. Part II: Parameterizations of coalescence efficiencies and fragment size distributions, J. Atmos. Sci., 67, 576–588, 2010. a

Szyrmer, W., Laroche, S., and Zawadzki, I.: A microphysical bulk formulation based on scaling normalization of the particle size distribution. Part I: Description, J. Atmos. Sci., 62, 4206–4221, 2005. a, b

Testud, J., Le Bouar, E., Obligis, E., and Ali-Mehenni, M.: The rain profiling algorithm applied to polarimetric weather radar, J. Atmos. Ocean. Technol., 17, 332–356, 2000. a, b, c, d

Testud, J., Oury, S., Black, R. A., Amayenc, P., and Dou, X.: The concept of “normalized” distribution to describe raindrop spectra: A tool for cloud physics and cloud remote sensing, J. Appl. Meteorol., 40, 1118–1140, 2001. a, b, c, d, e

Theus, M. and Urbanek, S.: Interactive Graphics for Data Analysis: Principles and Examples, CRC Press, 2008. a

Thurai, M. and Bringi, V. N.: Application of the generalized gamma model to represent the full rain drop size distribution spectra, J. Appl. Meteorol. Climatol., 57, 1197–1210, 2018. a, b, c

Thurai, M., Huang, G. J., Bringi, V. N., Randeu, W. L., and Schönhuber, M.: Drop shapes, model comparisons, and calculations of polarimetric radar parameters in rain, J. Atmos. Ocean. Technol., 24, 1019–1032, 2007. a, b

Thurai, M., Bringi, V., Szakáll, M., Mitra, S., Beard, K., and Borrmann, S.: Drop shapes and axis ratio distributions: Comparison between 2D video disdrometer and wind-tunnel measurements, J. Atmos. Ocean. Technol., 26, 1427–1432, 2009. a

Thurai, M., Bringi, V. N., Carey, L. D., Gatlin, P., Schultz, E., and Petersen, W. A.: Estimating the accuracy of polarimetric radar-based retrievals of drop-size distribution parameters and rain rate: An application of error variance separation using radar-derived spatial correlations, J. Hydrometeorol., 13, 1066–1079, 2012. a, b, c

Thurai, M., Kennedy, P. C., and Bringi, V.: DPWX/Modeling observed drop size distributions: 23 May 2015, available at:, last access: 27 August 2020. 

Thurai, M., Gatlin, P., Bringi, V., Petersen, W., Kennedy, P., Notaroš, B., and Carey, L.: Toward completing the raindrop size spectrum: Case studies involving 2D-video disdrometer, droplet spectrometer, and polarimetric radar measurements, J. Appl. Meteorol. Climatol., 56, 877–896, 2017a. a, b, c, d

Thurai, M., Mishra, K. V., Bringi, V., and Krajewski, W. F.: Initial results of a new composite-weighted algorithm for dual-polarized X-band rainfall estimation, J. Hydrometeorol., 18, 1081–1100, 2017b. a, b, c

Trewartha, G. T. and Horn, L. H.: An introduction to climate, McGraw Hill, 5th Edn., 1980. a

Ulbrich, C. W.: Natural variations in the analytical form of the raindrop size distribution, J. Clim. Appl. Meteorol., 22, 1764–1775, 1983. a, b, c

Wang, Y., Cocks, S., Tang, L., Ryzhkov, A., Zhang, P., Zhang, J., and Howard, K.: A prototype quantitative precipitation estimation algorithm for operational S-Band polarimetric radar utilizing specific attenuation and specific differential phase. Part I: Algorithm description, J. Hydrometeorol., 20, 985–997, 2019.  a

Wen, G., Chen, H., Zhang, G., and Sun, J.: An inverse model for raindrop size distribution retrieval with polarimetric variables, Remote Sens., 10, 1179, 2018. a, b, c

Short summary
The raindrop size distribution and its moments are fundamental in many areas, such as radar measurement of rainfall using polarimetry and numerical modeling of the microphysical processes of rain formation and evolution. We develop a technique which uses advanced radar measurements and complete drop size distributions using two collocated instruments to retrieve the lower-order moments such as total drop concentration and rain water content. We demonstrate a proof-of-concept using a case study.