Retrieval of Lower-Order Moments of the Drop Size Distribution using CSU-CHILL X-band Polarimetric Radar: A Case Study

. The lower-order moments of the drop size distribution (DSD) have generally been considered as difﬁcult to retrieve accurately from polarimetric radar data because these data are related to higher-order moments. For example, the 4 . 6 th moment is associated with speciﬁc differential phase, 6 th moment with reﬂectivity and ratio of high-order moments with differential reﬂectivity. Thus, conventionally, the emphasis has been to estimate rain rate ( 3 . 67 th moment) or parameters of the exponential or gamma distribution for the DSD. Many double-moment “bulk” microphysical schemes predict the total number concen- 5 tration (the 0 th moment of the DSD, or M 0 ) and the mixing ratio (or equivalently, the 3 rd moment M 3 ). Thus, it is difﬁcult 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 ﬁtted by the generalized gamma (G-G) distribution. The two reference moments are M 3 and M 6 which are shown to be retrievable using the X-band radar reﬂectivity, differential reﬂectivity and speciﬁc attenuation (from the iterative 10 ZPHI method). Along with the climatological shape parameters of the G-G ﬁt 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 15 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 M 0 through M 2 was < 15 % in magnitude, with Pearson’s correlation coefﬁcient > 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 M 0 resulted in coherent “time tracks” that can potentially lead to studies of precipitation evolution that have not been possible so far. SHV 3-dB three separate feed orthogonal mode transducers (a) an S-band feed width far-ﬁeld 1°) whose performance is described in Bringi et al. (2011a), an 180 band our retrievals of the reference moments { M 3 ,M 6 } are based on the X-band polarimetric data { Z h ,Z dr , Φ 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 speciﬁc 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 speciﬁc 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 the 171 . 25 ◦ azimuth). Details of the terrain variation between the radar and the Easton site are in Kennedy et al. ﬁlter versus ﬁltering along a ﬁxed showed that they The Lee-ﬁltered values of the radar data show the time evolution of the main echo over Easton site.


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 25 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 (A h ) is much larger than at S-band but not linear with R. The development of R(A h ) algorithms rests on a large body of work since the early 1990s and is related to attenuation-correction using differential propagation phase as a constraint (Bringi et al. 1990;Smyth and Illingworth 1998;Testud et al. 2000;and 30 references therein).
The retrieval of drop size distribution (DSD) parameters has also been a strong impetus for radar polarimetry. In this context, 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 Berne 2017a;Thurai and Bringi 2018). 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 60 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 {M 2 , M 3 , M 6 } was one possible choice mentioning lidar backscatter (M 2 ), microwave attenuation (M 3 ) and radar reflectivity (M 6 ). While the combination of M 3 and M 6 was not optimal for estimating the lower-order moments (in particular, M 0 ), it was better than using M 6 alone.
Recently, using the double-moment approach of L04, Raupach and Berne (2017a;RBa) showed that measured DSDs in 65 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 statiform 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. 70 The polarimetric (X-band) radar-retrieval of moments using reference moments M 3 , M 6 and an "invariant" h(x) of the G-G form were first described in Raupach and Berne (2017b;RBb). Their retrieval of M 6 was based on Z h while M 3 was retrieved in a two-step procedure using Z dr and K dp . 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" Z dr and K dp radar 75 data (for Z h < 37 dBZ), their validation of the moments (M 0 through M 7 ) using radar measurements was not conclusive but sufficient to demonstrate that their approach gave results similar to other methods based on normalized gamma model using "more classical" radar-retrievals of {N w , µ, D m } (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 M 3 and M 6 , there has been a reverse moment-based polarimetric forward operator . This reverse approach 80 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 Z h , Z dr , and K dp 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 {M 6 , M 9 } was optimal in terms of lowest variability in {Z h , Z dr , K dp } but that the pair {M 3 , M 6 } was suboptimal but still "useable". Thus, RBb's 85 choice of {M 3 , M 6 } as the two reference moments was validated by Kumjian et al. (2019). The rationale through which M 9 (whose sampling error would be very large using available disdrometers) entered the moment pair is not entirely clear. It could be because of correlating Z dr with absolute moments (M 0 through M 9 ) as opposed to the more physically-based ratio of moments such as D m = (M 6 /M 3 ) 1/3 in RBb or M 7 /M 6 as in Jameson (1983).
This work further develops on RBb but, instead of K dp , we employ specific attenuation (A h ) given its operational use in estimating R. The iterative ZPHI algorithm, which is a variant of Testud et al. (2000), is used here to estimate A h Park et al. 2005a,b). The reference moment M 3 , which is proportional to rain water content (W ) is retrieved using a modification of Jameson (1993) by fitting A h /W as a smoothed cubic spline with D m . The prior step is the retrieval of D m from Z dr and then retrieving D m from D m . This multi-step procedure was found to minimize the parameterization errors (also referred as algorithm errors) in the estimation of M 3 . As in RBb, the reference moment M 6 is derived as power law fit to Z h .

95
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 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 100 the instrumented site were additional factors that differed from RBb. We also show coherent "time tracks" in the D m versus M 0 , D m versus W , and D m versus M 6 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 Section 3. The retrieval of the two reference moments {M 3 , M 6 } follows in Section 4.

105
Several different ways of validating the moment retrievals are presented in Section 5. We follow this by a discussion in Section 6 and summarize the case study in Section 7. We cannot draw firm conclusions from just one case study even though the analysis is quite detailed. Rather this work may be considered as 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 110 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 {M 3 , M 6 } and then the variances of the other moments. The estimates of the variances of the ratio of correlated variables of the form X p Y −q are derived using a Taylor expansion to 2 nd order. The parameterization error variances are estimated for {M 3 , M 6 } and summed with the radar measurement error variances to yield the total error variance 115 for each moment retrieved.
Throughout this paper, we use "H" as a subscript for reflectivity Z H to denote units of dBZ at horizontal polarization. The lower case "h" in Z h means units of mm 6 m −3 . The same applies to Z DR (in dB) or Z dr (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 E{·} is used for the statistical expectation; · for the average 120 of its argument; Γ(·) for the gamma function; and Im {·} for the imaginary part of its complex argument.

Surface Instrumentation
The principal surface-based instruments used in this study are the MPS (or Meteorological Particle Spectrometer, manufactured by Droplet Measurement Technologies) and 3 rd 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 125 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 Section 4) of the reference moments M 3 and M 6 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 130 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 (MSL) as compared with 1.4 km MSL for Greeley. According to the Köppen-Trewartha climate classification system (Trewartha and Horn 1980), 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 (1970Knollenberg ( , 1976Knollenberg ( , 1981 and mea-135 sures drop diameter in the range from 0.05−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 140 et al. 1991(Korolev 140 et al. , 1998Baumgardner 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).

145
The 3 rd 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;Schönhuber et al. 2008;Thurai et al. 2007Thurai et al. , 2009Huang 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 150 is < 5% for fall velocity < 10 m s −1 (Schönhuber et al. 2008). A comparison of fall speeds between the MPS and 2DVD have 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 ms −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önhuber 2018). Thus, caution is necessary when using the 2DVD fall speeds for sizes < 0.6 mm (about 3 − 4 pixels).
In our application, we utilize the MPS for measurement of small drops with 0.1 ≤ D < 0.75 mm and the 2DVD for larger sized drops (see Raupach et al. 2019 for the rationale). This is termed here as the "complete" size spectrum and 2, 928 3-min 160 averaged spectra were available from the two sites with minimum rain rate of 0.1 mm h −1 and 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 Boutle 2012). The "shoulder" is the diameter range where the N (D) either remains steady or falls off more "slowly" (generally found under equilibrium conditions (McFarquhar 2004)).

165
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 shows 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 170 at the instrumented site during the 23 May 2015 event. The shape is equilibrium-like but with a single "drizzle" mode, a welldefined shoulder and faster (than exponential) fall off in the tail (Low and List 1982;McFarquhar 2004;Straub et al. 2010).
These features cannot be captured by either standard gamma or exponential fits, as shown in Fig. 1c.

CSU-CHILL Radar
The CSU-CHILL (Colorado State University-University of Chicago/Illinois State Water Survey) radar is described in Brunkow  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 gain of 55 dB. There are three separate feed or orthogonal mode transducers (OMTs) available: (a) an S-band feed (beam width in far-field is 1°) whose performance is described in Bringi et al. (2011a), (b) an 180 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 {M 3 , M 6 } are based on the X-band polarimetric data {Z h , Z dr , Φ 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 185 multiple cells with strong azimuthal and elevation gradients across the "echo cores" which generally precludes accurate dualwavelength 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    The general echo movement near Easton was estimated at 10 ms −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 210 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

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 Z h and Z dr 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). For correcting the measured Z h , we apply an iterative version of the ZPHI method, which uses a Φ dp constraint (Testud et al. 2000;) that was originally 225 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 (A h ) and specific differential phase (K dp ) is determined by minimizing a cost function based on least squares [we refer to Bringi and Chandrasekar 2001 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 A h = b 2 Z b1 h is assumed where b 1 = 0.78 and b 2 are constants (Park et al. 2005a). The method gives the estimate of A h at each resolution volume 230 in the selected range interval (here 0-40 km). The upshot of using A h instead of K dp is that the former closely follows the variations in Z h without the smoothing needed for estimating the latter, but with all the advantages of K dp such as immunity to calibration offsets and partial beam blockage Ryzhkov et al. 2014. Note that using A h for retrieval of W is restricted to precipitation comprising pure rain. In contrast, using K dp (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 Z H . In our multi- of 2 • ). The Easton instrumented site is located at range of 13 km, (b) measured and attenuation-corrected ZDR at X-band, (c) measured and filtered Φ dp , (d) specific attenuation (A h ). step retrieval procedure, it is reasonable to not 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 backscatter differential phase from the propagation phase. In essence, the estimate of A h may be considered as a by-product of attenuation correction of the measured Z h using the differential propagation phase over the selected path interval as a constraint.

240
The correction of the measured Z dr for differential attenuation is based on an extension of the method proposed by Smyth and Illingworth (1998) for C-band, which is described in  as a "combined Φ dp -Z dr " 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 A h determined by the Φ dp constraint is scaled by a factor ν and the measured Z dr is corrected for differential attenuation (A dp = νA h ) such that a desired value is reached at the end of the beam. The desired value is the 245 intrinsic or "true" Z dr at the end of the beam, which is estimated from the corrected Z h using a mean Z h -Z dr 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 Z dr at the end of the beam (generally Z dr ≈ 0 dB because of light rain at the end of the beam or because of ice particles above the 0 • C 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 Z H and Z DR , the measured and filtered Φ dp (which is used as constraint from 0-40 250 km), and A h are shown in the four panels of Fig. 4 at 20:27 UTC along the radial to the instrumented site located at Easton.
The Z H profiles show that very minor attenuation-correction is needed at this time, while the Z DR 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 A h peak is 1.5 dB/km coinciding with the Z H peak at 25 km. At the Easton location (13 km range) the A h is negligible. Figure 5a shows the PPI (at elevation angle of 1.5 • ) of the measured X-band Z H at 20:43 UTC at which time the peak 55 255 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 Figure 5b where the attenuation has increased dramatically with Z H corrected by 35 dB and Z DR 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 A h values have increased with peaks of 3 dB/km. At the Easton site the A h ≈ 1.5 dB/km. From a moment-retrieval 260 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).

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 Z h , Z dr , K dp , and A h 265 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 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 270 elevation angle of 1.5 • at 13 km range is 340 m AGL. The radar data from each pixel is 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 mins. The radar time series were shifted backward in time by 60 s as is common to match the peak in Z h (e.g., May et al. 1999). A more general 275 analysis of the error characterization of radar-gauge 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 Easton site. The composite 3-min averaged DSDs (an example was shown in Fig. 1b) were used to simulate radar observables as a 280 time series using the T-matrix scattering code (Barber and Yeh 1975;Bringi and Seliga 1977). The time resolution of 3mins corresponds to spatial scale of 1.8 km (using echo movement speed of 10 m s −1 ) which is less than the echo cell sizes respectively, for corrected ZDR, K dp and A h . which studied stratiform rain with embedded weak convection using 4 s samples; the 1/e-folding time was around 200 s, where e is the Napier's constant. For a highly convective case of our present case study, the decorrelation time would be substantially 285 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 radar sampling time of 40 s was 2.5 mins. In our case, the echo speed of 10 ms −1 together with the "effective" Lee filtered radar resolution of 2.5 mins 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).

290
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 Z H is about 5 dB lower than that measured by radar for time period 20:00-20:30 UTC. The measured Z H was 18-25 dBZ, implying very low rain rates (~0.5 295 mmh −1 ) and 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 AGL to SSE at surface. Given the unsteady conditions in this complex of echoes, it is not surprising that the disdrometer-based Z H calculation is biased low by around 5 dB relative to low radar Z H values of 18-25 dBZ. These problems are mitigated when heavier rain rates traverse the site about 15 mins later.

300
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 wavelength of 3 cm and temperature of 8 • C (Ray 1972) were 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 Z dr relative to radar measurements by 305 around 0.5 dB at~20:45 UTC. The discrepancy at 21:40 UTC noted in Figs. 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 show good calibration of Z h and Z dr . The radar-retrieved specific attenuation closely follows the Z h due to the A h -Z h power law assumption in the ZPHI method (only the fixed exponent is relevant) whereas the K dp does not, as expected.

310
As mentioned in Section 1, the methodology we used here follows RBb except for the use of specific attenuation rather than K dp in the retrieval of M 3 (however, both methods use Z dr in a multi-step retrieval described below). There are several advantages to this approach. First, "noisy" A h is strictly positive, as opposed to K dp 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 K dp from DSD measurements, meaning that all simulated K dp values used to train retrieval algorithms are 315 positive. Second, the smoothing in range necessary for K dp is not needed for A h which closely follows the spatial variability in Z h . 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 summarise it in the next subsection.

Overview
Moment retrievals from polarimetric radar data are a relatively recent application of the scaling/normalization of the DSD.

320
There are several aspects in this scaling as described by L04 namely, 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 N 0 and D m with h(x) being relatively "stable" with varying rain types/intensities. There is considerable latitude in the choice of reference moments RBb first suggested the use of {M 3 , M 6 } as the two reference moments suitable for polarimetric radar retrievals of the DSD.
They proposed retrieval of M 6 from radar measurements of Z h while for M 3 the retrieval was based on {Z dr , K dp }. 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 M k not only depends 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 Section 2, because otherwise the lower-order moments M 0 through M 2 will be in error.

350
The retrieval algorithms for the reference moments {M 3 , M 6 } 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 A h (> 1 dB km −1 ) from GXY alone was not enough to get a good retrieval. The scattering model assumptions are as given earlier in Section 3.3. For retrieval of M 6 (in mm 6 m −3 ), the obvious choice is Z h and a power law fit was derived for three ranges of The above ranges of Z H were based on trial and error to minimize the parameterization errors (Fig. 8a).  The next step is to retrieve A h /W from D m (adapted from Jameson 1993 who used a 3 rd order polynomial fit) for which we employ a smoothed spline fit. Here A h is in dB km −1 and W is the rain water content in gm −3 . We restricted the range of A h /W 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 D m falls outside the specified range. The retrieval of M 3 follows from where f (D m ) is the spline fit shown in Fig. 8c. The scatter plot of retrieved M 3 versus "true" M 3 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 375 Im {ε 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 A h increases with colder water temperature, but also depends on D m (Jameson 1993). Scattering simulations were performed for 8 • C and 20 • C and the spline fits of A h /W versus D m were compared. For low values of D m , the maximum difference was 35% occurring at D m = 0.75 mm (with A h /W larger at 8 • relative to 20 • C as expected) but a cross-over occurs near D m = 1.8 mm and the deviations increase in the opposite direction, with maximum deviation of -15% at D m = 3 380 mm (A h /W at 20 • C larger than at 8 • C due to scattering loss). Recall that the National Weather Service (NWS) sounding at Denver about 65 km away showed surface T of 12 • C. A lower temperature of 8 • C 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 A h = αK dp 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 385 surface T of 12 • C is close to the assumed T of 8 • C, the spline fit shown in Fig. 8c is considered to be sufficiently accurate for the retrieval of M 3 . We note that Diederich et al. (2015) found that the R(A h ) relation at X-band had a relatively "weak" dependence on temperature. Their fitted power law was 45.5A 0.83 h at 10 • C to 43.5A 0.79 h at 20 • C. At R = 10 mm h −1 , the A h at 10 • C is larger than at 20 • C by 6.8% while, at 100 mm h −1 , the A h at 10 • C is lower than at 20 • C by -10.8%; this cross-over is consistent with our calculations above. it is 9.5%. It is interesting to note that the mean relative bias for M 3 and M 6 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 M 2 were 0.106 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 M 6 dominates the total error variance, whereas for M 3 the retrieval algorithm and radar measurement errors are comparable.

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 Z h , Z dr , and A h time series data (see Fig. 6) were used to retrieve time series of {M 3 , M 6 }. Using Eq. (1) with the climatological h GG (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.

420
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.  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 430 normalized variances (or, the FSE) from the last column of Table A1   with "ground truth" (total FSE= 0.535 with 63% of the variance due to measurement error and 37% due to algorithm error).
440 Figure 9c compares M 6 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 M 6 varies from 10 to 10 5 mm 6 m −3 (or equivalently 10 to 50 dBZ); the peak value at 20:45 UTC being in excess of 50 dBZ. Note that M 0 and M 3 are the moments prognosed by "bulk" double-moment numerical schemes (actually M 0 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). was defined in Section 4.2; for each moment (M 0 through M 7 ) 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 450 range is the smallest for M 3 . This is expected because it is one of the reference moments. The median of the relative bias is "best" for M 5 with symmetric IQR range indicating very low skewness. The median RB for the moments M 0 through M 4 are around -15% but the skewness is significant for M 0 and progressively less for M 1 through M 5 . It might be unexpected that the retrieval of M 6 being one of the reference moments is less accurate than M 5 . One possible reason is that, at X-band, the larger drops are resonant-sized and the Z H does not vary as M 6 but rather closer to M 5 depending on the drop sizes. Fig. 7a, in fact,  shows that the fit for M 6 has a smaller slope for Z H > 37 dBZ because of resonant scattering. The median RB for M 6 and M 7 are < 17% but IQR indicates positive skewness (i.e., radar estimates are larger than "truth").
The difficulty in retrieving M 0 from higher-order moments {M 3 , M 6 } 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 460 drop end (concave up for negative µ). This is irrespective of well-constrained measurement and parameterization errors in the retrieval of the reference moments {M 3 , M 6 }.
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 Figure 11. Box plots of relative bias (RB) as in Fig. 8(e) except between radar-derived moments and completed DSD moments ("truth").
The orange horizontal line on top indicates the upper inner fence. 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" 465 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 (M 0 through M 2 ) 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 M 6 and M 7 is very large > 400% with large IQR. The moments M 3 through M 5 show more "normal" RB distributions with median values of RB in the range -60 to 90%, the minimum occurring for M 4 . However, 470 the Pearson's correlation coefficients, as shown in Table 3, are very low for all moments implying (practically) no linear variation between the moments. The Spearman's rank correlation is higher for moments M 4 through M 7 implying a non-linear 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 relative bias is > −1.  (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 M 0 and D m and from panel (b) a rapid decrease in W with corresponding rapid decrease in M 6 , reflecting advection of the rainshaft to the north of Easton.

485
Towards the end (last 5 data points from 20:57 to 21:12 UTC) the D m decrease is slowed down (from 1.5 to 1 mm) while the M 0 increases modestly from 10, 000 to 15, 000 per m 3 . This compensatory effect results in the rain water content being more or less steady (at 0.5 g m −3 ) while M 6 decreases from 37 to 33 dB (panel b and c, last 5 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.

490
The polarimetric radar-retrieval and validation of the lower-order moments of the DSD has 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)  h(x) being a special case of the G-G with c = 1, hence there is only one shape parameter µ. Note that this µ = µ ULB + 1, where 495 µ ULB is the shape parameter defined in Ulbrich (1983).
Many studies have attempted to retrieve the three parameters {N w , µ, D m } using polarimetric measurements {Z h , Z dr , K dp } 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 N w but found that validation was very difficult commenting that, ". . . the estimation of N w by all 500 algorithms is significantly affected by noise or other factors like radar volume versus point (disdrometer) measurement-scale mismatch and spatial separation." However, neither N 0 nor N w are the same as M 0 which is simply the total number concentration that scales the gamma pdf.
The estimation of either N 0 (or N w ) depends on the shape parameter (or the slope parameter Λ). Typically, the µ is assumed to be fixed or empirically derived as f (Λ) or other function of D m (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, M 0 have been conducted. Brandes et al. (2003) used an empirically derived µ-Λ relation based on 2DVD data. Using Z h and Z dr radar measurements at S-band they retrieved N 0 and Λ and obtained the M 0 as N 0 Λ −(µ+1) . They analyzed one convective event (with Z H varying from 10 to 55 dBZ) and showed the mean of log 10 (M 0 ) from radar was 2.74 compared with 2.84 from 2DVD-measurements (or, 550 and 690 per m 3 which are much smaller than the values obtained here, see Fig.10a). The key 510 point is that the mean µ was in the range 3-4 which is caused by truncation at the small drop end of the DSD. Wen et al. (2018) describe a different method of estimating the DSD parameters of the gamma distribution and the lowerorder moments based on an inverse model where the input is {Z dr , K dp /Z h } and the output is {µ, D max } where D max is the maximum diameter of the retrieved gamma DSD. Their approach follows the well-known k-nearest neighbour (k-NN) classification from pattern recognition literature (Shakhnarovich et al. 2006). This algorithm stores all input-output associations 515 from the available data as a "training" set. When a new {Z dr , K dp /Z h } input is presented, the algorithm assigns it the {µ, D max } output class that is the most common amongst the k nearest (training set) neighbours 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 neighbours although other distance functions are also employed in k-NN algorithms. They applied an empirical µ-Λ relation based on 2DVD data while N 0 is obtained a posteriori using Z h , µ, Λ, and D max . Their training set comprising Z dr and K dp /Z h was 520 generated using a polynomial function whose inputs µ and D max are drawn from ten-year disdrometer data with constraints The RBb article used X-band radar Z h to retrieve M 6 and {Z dr , K dp } to retrieve M 3 . Our approach is similar except that we use A h instead of K dp . There are several advantages to using A h (similar to its use in estimating R at X-band, Diederich et al. (2015). For instance, A h is always positive, is highly correlated with Z h variations thus preserving the spatial resolution and, 535 at X-band, has decent dynamic range. The article used several networks of 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 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" Z dr and K dp measurements when Z H < 37 dBZ. They "restored" the noisy Z dr by estimating it using a power law with Z h 540 while noisy K dp was restored using power laws of Z h and Z dr . They also commented that the majority of radar-measured Z H was < 37dBZ. 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 r 2 (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 545 errors, height of the radar pixels above the Parsivel network at longer ranges and other factors, similar to Anagnostou et al. (2013).

Summary
We demonstrated a proof-of-concept of the viability of radar retrieval of lower-order moments of the DSD using specific attenuation A h in addition to Z h and Z dr at X-band (an extension of RBb) via a case study approach. The use of specific attenuation 550 (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 mins.
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 555 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 {M 3 , M 6 } 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 {M 3 , M 6 } retrieval algorithms were trained using scattering simulations of Z h , Z dr , and A h 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 560 computed.
Polarimetric X-band radar data (acquired with an exceptionally narrow 3 dB beamwidth of 0.33 • ) data were extracted from a small polar box surrounding the instrumented site and the moments M 0 through M 7 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 relative bias, time series comparisons, scatter plots and correlation coefficients, it was determined that good 565 accuracy was obtained for the radar-based moments well beyond than possible hitherto. For the moments M 0 through M 2 , the relative bias 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 were tabulated.

570
The coherency of "time track" plots of radar retrieved quantities in the D m versus M 0 , D m versus W , and D m versus M 6 planes as the main 55 dBZ echo passed over the site (as well as 20 mins prior to and 20 mins 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, 575 requires investigations with more data.

Appendix A
The error model  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) = 580 Var (ε m ) + Var (ε p ). For different rain rate estimators such as R(Z h ), R (Z h , Z dr ), and R (K dp ), the Var( R)/R 2 is expressed in terms of the standard deviations of Z h , Z dr , and K dp 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 M 6 and M 3 first and then the other moments. Since M 6 is retrieved as a power law of Z 0.8 h (the exponent is approximate), Var(M 6 ) where Z h has units of mm 6 m −3 . Assuming the standard deviation of the radar measurement error is typically 1 dB, we get 590 Var(Z h )/Z h 2 = 0.067. This implies Var M 6 /M 6 2 = 0.043.
The variance of M 3 is more complicated because it is a multiple step procedure as described in Section 4 involving smoothed spline fits. We use approximate power law fits for estimating Var(M 3 ) as follows. We have Thus, Since D m is linear with D m , and assuming the standard deviation of the radar measurement of Z DR = 0.3 dB, where Z dr is a ratio, and Var(Z dr )/Z dr 2 = 0.0051, From Thurai et al. (2017b; Appendix, (A5)),

605
assuming that A h varies as Z 0.8 h used in the ZPHI method. The standard deviation of the K dp measurement is typically 0.3 • km −1 and that of Z H is 1 dB. Further, the mean K dp for our data ≈ 1 • km −1 (Fig.6c) but variable in time. In any case, , where p k and q k are rational numbers and C k is some constant. The variance of M k needs more elaboration as X ≡ M 3 an Y ≡ M 6 are correlated. This correlation arises because A h is a power law with Z h , i.e., Z 0.8 h in the ZPHI method. Together with M 6 being a power law with Z h , the M 3 and M 6 are correlated with a correlation coefficient of 615 0.93 obtained from radar-derived M 3 and M 6 . For the k-th moment M k , p k = 6−k 3 and q k = 3−k 3 for k = 0, 1, · · · , 7, k = 3, 6. The objective is to derive Var(M k ) 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

Var (A h )
Consider the parameter vector θ = (X, Y ). Then, second-order Taylor series approximation of f (X, Y ) around θ produces where the notations f x (θ) and f xy (θ) represent, respectively, the first and second-order derivatives of the function f with respect to the variables in the subscript and evaluated at θ.
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 θ) Var(X) + 2f x (θ)f y (θ) Cov(X, Y ) + f y (θ) Var(Y ).
Again, evaluating at θ, we obtain Substituting the above in (A20) leads to the first-order approximation From (A18) and (A24), we obtain the desired ratio as If the correlation coefficient ρ XY between X and Y is known, then we replace Cov(X, Y ) in (A25) to obtain Using (A26), Table A1  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 (A26) and total variances for M 3 and M 6 gives in Table A1 the total variances for the other moments. For M 0 through M 2 , the total variance is smaller than the measurement variance because the covariance term in negative for those moments.

665
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 data are available by request from either the corresponding author or by request to http://www. chill.colostate.edu/w/Contacts. The MPS and 2DVD processed data are available upon request to the corresponding author. Table A1. Variance estimates for the radar-based retrievals of moments M0 through M7. The p and q are the exponents of M k = M p 3 M −q 6 . The normalized variances due to measurement errors are given in the 4 th column while the total due to measurement and parameterization errors are given in the 5 th column.

670
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: http://www.chill.colostate.edu/w/DPWX/Modeling_observed_drop_size_distributions:_23_May_2015.