Retrieving Rain Drop Size Distribution Moments from GPM Dual-Frequency Precipitation Radar

A novel method for retrieving the moments of rain drop size distribution (DSD) from the dual-frequency precipitation radar (DPR) onboard the global precipitation mission satellite (GPM) is presented. The method involves the estimation of two chosen reference moments from two specific DPR products, namely the attenuation-corrected Ku-band radar reflectivity and (if made available) the specific attenuation at Ka-band. The reference moments are then combined with a function representing the underlying shape of the DSD based on the generalized gamma model. Simulations are performed to quantify the algorithm errors. The performance of methodology is assessed with two GPM-DPR overpass cases over disdrometer sites, one in Huntsville, Alabama and one in Delmarva peninsula, Virginia, both in the US. Results are promising and indicate that it is feasible to estimate DSD moments directly from DPR-based quantities.


Introduction
The global precipitation measurement (GPM) core satellite hosts the dual-frequency precipitation radar (DPR), which is the first spaceborne radar operating at Ku-and Kabands for precipitation mapping [1]. The dual-frequency radar measurements [2] provide a more complete depiction of precipitation globally (±65 degree latitude; [3,4]). The DPR also facilitates better estimates of the two main parameters governing the drop size distributions (DSDs), namely the mass-weighted mean diameter (D m ) and the normalized intercept parameter (N W ) [5,6]. The DPR precipitation retrieval algorithm [7] employs the normalized representation of the drop size distribution based on the gamma distribution with shape parameter fixed (µ = 3) and uses a rather complex adjustment parameter to estimate the DSD and precipitation rate [8].
A number of recent studies have analyzed the statistics of [D m , N W ] from the DPR and compared these with similar statistics from a diverse network of ground-based polarimetric radars [9] and surface disdrometers on a seasonal or global scale [10,11]. On the global scale, the higher rain rate (>8-10 mm/h) DSD parameters [D m , N W ] are consistent with disdrometer data from continental convection, deep oceanic convection, and oceanic shallow rain. However, on a seasonal basis, the N W estimate is not consistent with disdrometer data [10]. Such statistical comparisons have to encompass many years of data to overcome the "snapshot" view from DPR over a footprint of 5 km × 5 km versus nearly continuous [M 3 , M 6 ] two-moment predictions. Thus, it is not clear if [M 3 , M 6 ] is the "best" choice for the prognostic moments when considering in general all warm rain processes.
The choice of reference moments will be limited to higher orders for both polarimetric and dual-frequency radars. For dual-polarization radar algorithms, the physical basis relies on the correlation between raindrop shape and size, while dual-wavelength weather radar algorithms rely primarily on non-Rayleigh scattering at the shorter wavelength. The equations for estimating parameters of the DSD are nearly identical in the presence of attenuation [29].
Recently, retrieval of DSD moments from ground-based polarimetric radar measurements has been evaluated using X-band data, with special emphasis on the prediction of lower-order moments. Two examples are [30,31]. In both cases, the procedure initially entailed the retrieval of two "reference" moments, in both cases [M 3 , M 6 ], followed by reconstructing the DSDs and calculating other moments assuming a specific function to represent the underlying shape (denoted h(x), as explained in the next section). In [30], the two reference moments are estimated using the co-polar reflectivity (Z h for horizontal polarization), differential reflectivity (Z dr ), and the specific differential propagation phase (K dp ), whereas in [31], they are estimated using Z h , Z dr , and specific attenuation, A h .
In this paper, we examine a conceptually similar approach for retrieving the DSD moments from the DPR products. As with the polarimetric radar retrievals, we first estimate two (chosen) reference moments [M 3 , M 6 ]. Then, using the "most probable" underlying shape function h(x), we determine other DSD moments. The methodology and concepts are given in Section 2. Section 3 describes the DSD datasets and simulation results to demonstrate which moments are the best reference moments and which DPR-based quantities should be used to estimate those. Parameterization errors are quantified based on "true" versus retrieved moments from our simulation results. In Section 4, we present two GPM-DPR overpass cases and assess the accuracy of the retrievals, but noting an important limitation, namely that one of the quantities needed for our proposed method is not available as an official DPR product. Error sources are considered in Section 5, along with some discussion of results, caveats, and conclusions.

Methodology and Concepts
The scaling normalization framework of Lee et al. [18] expresses any moment of the DSD (inclusive of the lower order moments, M 0 -M 2 ) in terms of the product of the power laws of the chosen reference moments and the underlying function h(x) describing the intrinsic shape of the DSD. The invariance of h(x) has been demonstrated in a number of articles ( [32,33] and Figure 2 in [20]). In compact notation, the DSD is expressed as N(D) = N 0 h(x), where (a) the intrinsic shape of the DSD, h(x) is a function of x, the scaled diameter D/D m , and (b) N 0 is the normalized "intercept" parameter (see later in Equations (4) and (5)).
The nth moment, M n , of the DSD (units are mm n /m 3 ) is given by: where D is the equi-volume drop diameter, and N(D) is the drop concentration per unit volume in the diameter interval D to D + δD. When n equals zero, i.e., the zeroth moment, M 0 will represent the total number of drops per unit volume. In remote sensing applications at microwave frequencies, M 3 will be related to specific attenuation, and M 6 will be related to the radar reflectivity, but in the non-Rayleigh scattering region, the relationships will show more scatter depending on the frequency-band.

Exponential and Gamma Distributions
Characterization of the DSD dates back to the 1940s, with one of the most-quoted references being that of Marshall and Palmer [34] who used manual DSD measurements in stratiform rain and used an exponential distribution (with a rainfall intensity-dependent slope parameter) to represent their data. Decades later, a three-parameter, un-normalized gamma distribution was introduced [35] to better represent DSD measurements with smaller time intervals. This has been widely used in numerous studies on DSD variability, e.g., [36][37][38][39]. Retrievals of DSD parameters, e.g., from polarimetric radar data and from spaceborne radar data, have also assumed the gamma model [40][41][42][43]. The model will be referred to as the "standard gamma (SG) model".

Double-Moment Normalization and Generalized Gamma Model
A more novel and sound approach for DSD representation was introduced by Sempere-Torres et al. [44] and by Testud et al. [5] after several studies had shown that the underlying shape of the DSD can be revealed if N(D) is normalized by the parameter N W , the normalized intercept parameter, together with the normalization D by the parameter D m , the mass-weighted mean diameter [5,45]. Often, the underlying shape is denoted as h(x) where x is the so-called "scaled diameter", [46], given by: where D m is defined later in Equation (5).
Further to the scaling law of Sempere-Torres [44] and the scaling normalization of Testud [5], a unified approach was proposed in [18] based on double-moment normalization with any two reference moments, [M i , M j ], resulting in: The choice of the two reference moments will depend on the application. Lee et al. [18] also considered the use of the so-called 'generalized gamma' (GG) model, which had previously been shown to be applicable and useful for DSD studies in, e.g., [47]. They showed that the GG was far more suitable to represent h(x), since it uses two shape parameters, µ GG and c, and as a result, yielded more flexibility. The term h(x) will then be given by: Combining Equations (3), (4), and (6) and rearranging these, we obtain: Equation (7) forms the basis of our analyses presented here. As seen, N(D) is dependent on the two chosen reference moments and the two shape parameters µ GG and c. Note also that, if we choose the third and the fourth moments as our reference moments, i.e., set i = 3 and j = 4, and further, c is set to 1, then Equation (6) will simplify to the standard gamma model, with the shape parameter (denoted as µ SG in this paper). Under these conditions, µ SG = µ GG − 1, and further, when µ GG is set to 1, i.e., µ SG = 0, an exponential distribution is obtained.

Choice of Reference Moments
The choice of the two reference moments in Equation (6), as mentioned earlier, depends on the application. For rain microphysical studies, M 3 and M 4 are often used, whereas for remote sensing applications, M 3 and M 6 are more commonly used.
Since the moments of rainfall rate, R and liquid water content, W, are close (i.e., M 3 vs. M 3.67 , e.g., [40]), we can expect a near linear relation between W and k Ka . For Rayleigh scattering reflectivity, Z is the sixth moment of the DSD, which is approximately the case for Ku-band. Thus, we have selected k Ka for retrieval of M 3 and Z Ku for M 6 , i.e., [M 3 , M 6 ] have been chosen as the reference moments. Although this quantity is currently not available as an official DPR product, it can be 'reproduced' by the DSD parameters that are currently accessible. We now examine these variabilities via scattering simulations using DSD measurements.

DSD Data for Simulations
The DSD data used herein were obtained from three observation campaigns in three different locations, namely (i) Greeley, Colorado, (ii) Huntsville, Alabama, and (iii) Wallops Island, Virginia, three climatically very different locations in the US. They represent semiarid continental climate, a sub-tropical continental climate, and a mid-latitude coastal region, respectively. In all three cases, the DSDs were obtained from (a) Meteorological Particle Spectrometer (MPS; [48]) and (b) 2D video disdrometer (2DVD; [49,50]). While the MPS, with its 50-micron resolution, provided accurate measurements of concentrations of small drops (especially for D < 1.2 mm), the 2DVD provided better measurements for larger drop sizes. In all three locations, the MPS and the 2DVD had been installed within a two-third scaled double wind fence (DFIR; [51]). The overlap of the DSD measurements between the two instruments has been investigated extensively in [52] and found to be reasonably close to within acceptable accuracy in the 0.7 to 1.2 mm drop diameter range. The composite DSD, which we refer to here as the "full DSD spectra" [53], ranged from D = 0.15 mm up to large and very large drops, the maximum recorded in our datasets being D = 8 mm in the outer rain-bands of a hurricane event (Dorian) over Wallops site [54]. MPS measurements were used for D ≤ 0.8 mm and 2DVD measurements for D > 0.8 mm.

Scattering Simulations
A total of around 3000 3 min DSDs from the three locations were used as input to T-matrix scattering calculations at Ku and Ka bands. A 20 deg C temperature was assumed. This computational method uses the spherical vector wave functions for the farfield scattering matrix with unknown expansion coefficients, and the total field inside the dielectric is also expanded the same way. The incident wave is also expanded in the same way but with known coefficients. The extended boundary condition is the main principle whereby the surface electric and magnetic currents on the dielectric radiate the negative of the incident field throughout the particle and the far field scattered external to the minimum sphere containing the particle. Thus, the far field expansion coefficients and the incident field are cast into a matrix called the transition matrix (T-matrix), which depends on the dielectric shape (must have an axis of rotational symmetry) and the dielectric constant. The method is very fast compared to the method of moments, which use patch elements (a large number compared to the entire basis functions used in the T-matrix). A comprehensive review of this method can be found in [55].
The computed radar reflectivity (Z, dBZ) and the specific attenuation (k), together with the DSD moments (from M 0 to M 7 ), showed that the best estimates (i.e., with the lowest uncertainties) for [M 6 and M 3 ] were Z Ku and k Ka , respectively. The results are shown in Figure 1a,b as color-contoured frequency of occurrence (2D) plots. Panel (a) shows very little scatter, indicating that the uncertainties associated with the estimation of M 6 will be very little, whereas panel (b) shows somewhat more significant scatter, indicating the Remote Sens. 2021, 13, 4690 6 of 22 uncertainties in M 3 estimation will be somewhat larger. Note, however, that the color scale of number N in Figure 1 is on a log 10 scale.
Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 23 very little, whereas panel (b) shows somewhat more significant scatter, indicating the uncertainties in M3 estimation will be somewhat larger. Note, however, that the color scale of number N in Figure 1 is on a log10 scale. The dashed lines in Figure 1a,b represent the fitted equations based on second-order polynomial equations: where Zku is in dBZ units in Equation (8a). The fitted coefficients and their standard errors are given in Table 1.

Stability of the Underlying DSD Shape
Next, we consider the uncertainties associated with the underlying shape of the DSDs using Equations (2)-(4). Figure 2a,b show the double-moment normalization (i.e., h(x) versus x) applied to the Greeley (GXY) and the Huntsville (HSV) datasets. In each case, the 25th, the 50th (median), and the 75th percentile curves are superimposed. The two sets of curves appear to be similar to each other.
Each of the 3 min DSDs from the GXY and HSV campaigns had been fitted to Equation (6) with i = 3 and j = 6 [53]. The fitted μGG versus c are shown in panel (c) as a colorintensity plot. μGG shows a narrow range, but the fitted c values range a wider span. Note, once again, that the color scale of N is on a log10 scale. The maximum value of N is reached for [μGG, c] = [−0.25, 3.67]. We will use these as the 'most-probable' combination and will substitute these values into Equation (6) for i = 3 and j = 6. To test this most-probable h(x), we show in panel (d) the color-intensity plot of h(x) versus x derived from the 3 min DSDs recorded during category-1 Hurricane Dorian (outer bands only) at the Wallops site [56] as well as the most-probable h(x) from the GXY-HSV analyses. The fit appears to be very representative, capturing the 'high-intensity' trend from the DSD data very well.  (7) and (8). The dashed lines in Figure 1a,b represent the fitted equations based on second-order polynomial equations: where Z ku is in dBZ units in Equation (8a). The fitted coefficients and their standard errors are given in Table 1.  (7) and (8)

Stability of the Underlying DSD Shape
Next, we consider the uncertainties associated with the underlying shape of the DSDs using Equations (2)-(4). Figure 2a,b show the double-moment normalization (i.e., h(x) versus x) applied to the Greeley (GXY) and the Huntsville (HSV) datasets. In each case, the 25th, the 50th (median), and the 75th percentile curves are superimposed. The two sets of curves appear to be similar to each other.
Each of the 3 min DSDs from the GXY and HSV campaigns had been fitted to Equation (6) with i = 3 and j = 6 [53]. The fitted µ GG versus c are shown in panel (c) as a color-intensity plot. µ GG shows a narrow range, but the fitted c values range a wider span. Note, once again, that the color scale of N is on a log 10 scale. The maximum value of N is reached for [µ GG , c] = [−0.25, 3.67]. We will use these as the 'most-probable' combination and will substitute these values into Equation (6) for i = 3 and j = 6. To test this most-probable h(x), we show in panel (d) the color-intensity plot of h(x) versus x derived from the 3 min DSDs recorded during category-1 Hurricane Dorian (outer bands only) at the Wallops site [56] as well as the most-probable h(x) from the GXY-HSV analyses. The fit appears to be very representative, capturing the 'high-intensity' trend from the DSD data very well. Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 23

Simulations and Algorithm Errors
As seen in Sections 3.3 and 3.4, our proposed retrieval method is expected to have uncertainties associated with each of the retrieved moments. To quantify these "algorithm errors", forward simulations and retrievals were performed using the same set of 3 min DSDs from the three locations. Figure 3 shows the block-diagram for this procedure. The steps can be summarized as follows: i.
Use the 3 min DSDs for scattering calculations at Ku and Ka bands; ii. Use the Zku and kka outputs from (i) to estimate M6 and M3 using Equations (7)  These steps are numbered in red in Figure 3. For step (iii) (with green-box outline), the retrieval of other moments can be achieved analytically using equation (42) of Lee et al. [18]; however, that expression is not valid when μGG is negative, which is the case for our h(x). An alternate solution is to estimate the other moments numerically by first constructing the full DSD spectra from each of the retrieved M3 and M6 in step (ii) and our h(x) with i = 3 and j = 4 for μGG = −0.25 and c = 3.67 in Equation (7). This is followed by the calculation of other moments using Equation (1).
Results from step (v) are shown in Figure 4 as scatterplots of the retrieved moments versus the 'true' moments, i.e., those derived directly from the 3 min DSDs. Only the cases with kKa > 1 dB/km were chosen. The 1:1 line is included in each plot. The higher-order moments, viz. M3 to M7 show very good correlation and lie close to the 1:1 line. The lower-

Simulations and Algorithm Errors
As seen in Sections 3.3 and 3.4, our proposed retrieval method is expected to have uncertainties associated with each of the retrieved moments. To quantify these "algorithm errors", forward simulations and retrievals were performed using the same set of 3 min DSDs from the three locations. Figure 3 shows the block-diagram for this procedure. The steps can be summarized as follows: i.
Use the 3 min DSDs for scattering calculations at Ku and Ka bands; ii. Use the Z ku and k ka outputs from (i) to estimate M 6 and M 3 using Equations (7) and (8) Compare (iii) and (iv).
These steps are numbered in red in Figure 3. For step (iii) (with green-box outline), the retrieval of other moments can be achieved analytically using equation (42) of Lee et al. [18]; however, that expression is not valid when µ GG is negative, which is the case for our h(x). An alternate solution is to estimate the other moments numerically by first constructing the full DSD spectra from each of the retrieved M 3 and M 6 in step (ii) and our h(x) with i = 3 and j = 4 for µ GG = −0.25 and c = 3.67 in Equation (7). This is followed by the calculation of other moments using Equation (1).
Results from step (v) are shown in Figure 4 as scatterplots of the retrieved moments versus the 'true' moments, i.e., those derived directly from the 3 min DSDs. Only the cases with k Ka > 1 dB/km were chosen.  Table 2. They range from nearly 11% for M 0 to just over 3% for M 7 , gradually decreasing. These values represent the overall algorithm errors of the retrieval method.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 23 order moments, especially M0, show more scatter. The errors are quantified in terms of the fractional standard error (FSE) in Table 2. They range from nearly 11% for M0 to just over 3% for M7, gradually decreasing. These values represent the overall algorithm errors of the retrieval method.    Table 2. They range from nearly 11% for M0 to just over 3% for M7, gradually decreasing. These values represent the overall algorithm errors of the retrieval method.    Algorithm errors are sometimes referred to as parameterization errors. Bringi et al. [31] considered an additional source of uncertainty, namely those due to radar measurement errors. They found that, for lower-order moments, the parameterization error dominates the overall uncertainties, whereas the radar-measurement errors dominate those for higher-order moments. Although their analyses were related to retrieval of moments from polarimetric radar data, the same trend is likely to apply for the DPR-based retrievals also. Next, we preform basic tests for our retrieval method for two GPM-DPR overpass cases.

The Huntsville Event
The first of the two events considered here occurred over the Huntsville site, where the MPS and two 2DVDs had been installed ( [53]; see also Section 3.2 earlier). It consisted of precipitation associated with a semi-organized line-convection on 11 April 2016, which moved across northern Alabama from~17 to 24 h UTC. The line convection did not directly pass over the disdrometer site; instead, it was approximately 10 km to the west. The GPM-DPR overpass over the site occurred at 23:31 UTC.
The near-surface DPR products include measured and attenuation-corrected radar reflectivities both at Ka and Ku bands as well as N w and D m for each of the 4 km by 4 km radar pixels. The two quantities needed for our retrievals, as mentioned earlier, are the Ku-band radar reflectivity and the specific attenuation at Ka-band. The former is a 'direct' product, but of course, Z Ku should be the reflectivity values after attenuation-correction is used. The latter, however, is not readily available, but on the other hand, it can be recreated by (i) utilizing N w and D m to derive the DSDs for each pixel, assuming the DPR assumption of standard gamma model with µ = 3, followed by (ii) T-matrix calculations to derive k Ka for each of the DSDs. These are shown in panels (a) and (b) in Figure 5. In each panel, the color scale is included. The diamonds represent the location of the disdrometers in terms of latitude and longitude. Z ku (attenuation-corrected) over the site was 34.6 dBZ, and k Ka was 0.60 dB/km. The estimated M 6 and M 3 are shown in panels (c) and (d), respectively. The color scales (included within each panel) show the log 10 of the moments. Using these two moments, together with the h(x) in Equation (6), the DSDs per pixel were derived, and the other moments, M 0 , M 1 , M 2 , M 4 , M 5 , and M 7 , were estimated. Four are shown in Figure 6.
Validation of the retrieved moments is not straightforward, but in Table 3, we show the mean and standard deviation of all moments M 0 to M 7 derived from the DPR data in (22) pixels over and surrounding the disdrometer site, which are compared with those derived directly from the (20) 3 min DSD within ±30 min time interval around the DPR overpass.      From Table 3, we observe the following: i. Overall good agreement between the two sets of mean values, especially for the zeroth moment. ii. However, the DSD-based moments show somewhat lower mean values that the DPR-retrieved moments; although, when the standard deviations are included, the overlaps are considerable. iii. In both cases, M 2 shows the lowest mean values and M 7 the highest; the standard deviations also show a similar trend.
Furthermore, values of D m , given by the ratio of the fourth moment to the third ranged from 0.9 to 1.8 mm for the DPR-based retrievals and 1.1 to 1.8 mm for the DSD-based estimates. The mean was 1.5 mm for both cases. Hence, it is likely that our retrievals of moments represent the 'true' values over the 4 km by 4 km pixel areas. Further comparisons are given in Appendix A.

Remnants of Storm Sally
The second event we consider in this paper is the remnants of Hurricane Sally over the Wallops area. This storm made landfall in the southern part of the US as category-1 hurricane but weakened considerably as it headed north-east. Remnants of this storm passed over the disdrometer location, and the GPM-DPR overpass occurred north of the site on 18 September 2020 at 04:40 UTC. The closest approach for the nadir beam was around 170 km NE of the disdrometer site, but the closest off-nadir beam was around 65 km north-east. Panel (a) in Figure 7 shows the (near-surface) attenuation-corrected Z Ku , and panel (b) shows k Ka . The retrieved M 6  Direct validation of the retrieved moments is, once again, not possible, especially because the DPR did not directly go over the disdrometer site, but in Figure 8, we compare the histograms of [M 0 . . . M 7 ] from DPR and the DSDs. For the latter, ±4 h around the DPR overpass time was used, whereas for the former, only the data shown in Figure 7 were used (a total of 3269 pixels). There seems considerable overlap in the two sets of histograms in all panels, but the DPR-based moments generally appear to be somewhat higher. On the other hand, if we restrict the DSD-based moment calculations to within 2 h of the DPR overpass, improvements in the agreement were observed. The green circle in each panel of Figure 8 represents the averaged value from the DSD data within the 2 h interval. They correspond very well to the mode of the DPR-based histograms.

Evaluation of DSD Parameters
The two main DSD parameters often used in the radar retrievals are N W and D m . There is also another parameter, namely the standard deviation of mass-spectrum, denoted as σ M [37], considered useful for rain microphysics studies. While the first two parameters can be defined in terms of the third and the fourth moments, the third parameter additionally requires the fifth moment. The following equations can be derived: Note, however, that Equation (11) is only valid for µ GG > 0. For negative values, it is more appropriate to derive σ M from the DSDs constructed using the two reference moments and the most-probable h(x). Direct validation of the retrieved moments is, once again, not possible, especially because the DPR did not directly go over the disdrometer site, but in Figure 8, we compare the histograms of [M0 … M7] from DPR and the DSDs. For the latter, ±4 h around the DPR overpass time was used, whereas for the former, only the data shown in Figure 7 were used (a total of 3269 pixels). There seems considerable overlap in the two sets of histograms in all panels, but the DPR-based moments generally appear to be somewhat higher. On the other hand, if we restrict the DSD-based moment calculations to within 2 h of the DPR overpass, improvements in the agreement were observed. The green circle in each panel of Figure 8 represents the averaged value from the DSD data within the 2 h interval. They correspond very well to the mode of the DPR-based histograms.

Evaluation of DSD Parameters
The two main DSD parameters often used in the radar retrievals are NW and Dm. There is also another parameter, namely the standard deviation of mass-spectrum, denoted as σM [37], considered useful for rain microphysics studies. While the first two parameters can be defined in terms of the third and the fourth moments, the third parameter additionally requires the fifth moment. The following equations can be derived:  Remnants of storm Sally shows narrower D m histograms (both DPR-retrieved and DSD-based) having lower D m values with a maximum value of only 1.6 mm whereas the whereas the HSV event shows D m 's ranging up to 2.1 mm. This is to be expected, because storm Sally originated as a Hurricane, and it is well known that such storms contain an abundance of small drops (higher concentration) compared with other rain regimes [57,58]. ii. The HSV event shows two peaks in the D m , and these bimodal peaks are evident in both the DPR-retrieved and the DSD-based histograms. It is very plausible that the two peaks arise from the (semi-organized) line convection being embedded within a larger widespread (probably stratiform) rain region. The bimodal peaks are also noticeable in the σ M histograms. Storm Sally, on the other hand, shows only one peak in both D m and σ M . iii. For the HSV event, there is considerable overlap between the DPR-retrieved and the DSD-based histograms, whereas for storm Sally, the DSD-based histogram has a higher number of cases with lower D m values (i.e., <0.6 mm). This may well be due to the DPR sensitivity, which has a lower limit of approximately 12 dBZ for the radar reflectivity at Ku-band, which, in turn, indicates that light rainfall cases will not be detected often enough. On the other hand, the disdrometer-based DSDs include the MPS measurements with good accuracy for the concentration of small drops. iv. Another feature that is different between the two panels concerns the proportion of stratiform to convective rain. From the two clear peaks in D m and σ M histograms observed in the HSV event (as noted earlier in point (ii)), we can infer that the proportion of the two rain types are somewhat comparable. For storm Sally, by contrast, we have ascertained from NPOL radar quasi-vertical profiles (QVP; [59]) that it was largely stratiform rain. The single peak D m histograms from both DPR and DSD data in panel (b) of Figure 9 support this.
the whereas the HSV event shows Dm's ranging up to 2.1 mm. This is to be expected, because storm Sally originated as a Hurricane, and it is well known that such storms contain an abundance of small drops (higher concentration) compared with other rain regimes [57,58]. ii. The HSV event shows two peaks in the Dm, and these bimodal peaks are evident in both the DPR-retrieved and the DSD-based histograms. It is very plausible that the two peaks arise from the (semi-organized) line convection being embedded within a larger widespread (probably stratiform) rain region. The bimodal peaks are also noticeable in the σM histograms. Storm Sally, on the other hand, shows only one peak in both Dm and σM. iii. For the HSV event, there is considerable overlap between the DPR-retrieved and the DSD-based histograms, whereas for storm Sally, the DSD-based histogram has a higher number of cases with lower Dm values (i.e., <0.6 mm). This may well be due to the DPR sensitivity, which has a lower limit of approximately 12 dBZ for the radar reflectivity at Ku-band, which, in turn, indicates that light rainfall cases will not be detected often enough. On the other hand, the disdrometer-based DSDs include the MPS measurements with good accuracy for the concentration of small drops. iv. Another feature that is different between the two panels concerns the proportion of stratiform to convective rain. From the two clear peaks in Dm and σM histograms observed in the HSV event (as noted earlier in point (ii)), we can infer that the proportion of the two rain types are somewhat comparable. For storm Sally, by contrast, we have ascertained from NPOL radar quasi-vertical profiles (QVP; [59]) that it was largely stratiform rain. The single peak Dm histograms from both DPR and DSD data in panel (b) of Figure 9 support this.  Further improvement in the agreement between DPR and DSD histograms in Figure 9 may be obtained if we restrict the off-nadir angles for the DPR to be less than 10 degrees and/or by utilizing the 'type_Precip' flag product from DPR, in particular, the second digit of the flag. This will be investigated in the near future.

Evaluation of (Stratiform and Convective) Rain Types
Regarding stratiform and convective rain types, it is also of interest to examine the applicability of the classification based on N w versus D m separation technique. This observational-based empirical method has been used for disdrometer data in the past, and our own recent work has shown that the two rain types can be separated in the N W -D m space [60][61][62]. The same method had also been applied to the S-band NPOL radar data and compared with a texture-based method described in [63]. Considerable agreement of over 86% was obtained [64]. Similar comparisons [65] were made for the C-band CPOL radar based in Darwin, Australia [66].
Disdrometer data represent "point-measurements" (albeit over a 1 or 3 min time interval), and ground-based radar data represent measurements over pulse volumes (but instantaneous). Despite these differences, the stratiform-convective rain separation technique seemed applicable to both. Here, we perform similar tests for the DPR-retrieved N W -D m for the two cases, noting that the DPR footprint is much larger (4 km by 4 km) than the ground-based radar data.
The separation technique is based on whether a given pair of [N W , D m ] lies above or below or specified line. A simple index i given by the following equation is used to quantify the rain-type likelihood [64,65]: where 'sep' corresponds to the separation line, and 'est' represents the estimated value.
Values of c 1 and c 2 may be location dependent, but to be consistent with our previous studies, we have used −1.682 and 6.541, respectively. Broadly speaking, when i is negative, stratiform rain is indicated, and when i is positive, convective rain is indicated. The pixel-by-pixel DPR-based index value derived for the HSV event is shown in panel (a) of Figure 10. The red color represents i > 0 (hence, convective rain regions), and the lighter colors indicate i < 0 (hence, stratiform rain regions). Panel (b) of Figure 10  pared with a texture-based method described in [63]. Considerable agreement of over 86% was obtained [64]. Similar comparisons [65] were made for the C-band CPOL radar based in Darwin, Australia [66].
Disdrometer data represent "point-measurements" (albeit over a 1 or 3 min time interval), and ground-based radar data represent measurements over pulse volumes (but instantaneous). Despite these differences, the stratiform-convective rain separation technique seemed applicable to both. Here, we perform similar tests for the DPR-retrieved NW-Dm for the two cases, noting that the DPR footprint is much larger (4 km by 4 km) than the ground-based radar data.
The separation technique is based on whether a given pair of [NW, Dm] lies above or below or specified line. A simple index i given by the following equation is used to quantify the rain-type likelihood [64,65]: where 'sep' corresponds to the separation line, and 'est' represents the estimated value.
Values of c1 and c2 may be location dependent, but to be consistent with our previous studies, we have used −1.682 and 6.541, respectively. Broadly speaking, when i is negative, stratiform rain is indicated, and when i is positive, convective rain is indicated. The pixelby-pixel DPR-based index value derived for the HSV event is shown in panel (a) of Figure  10. The red color represents i > 0 (hence, convective rain regions), and the lighter colors indicate i < 0 (hence, stratiform rain regions). Panel (b) of Figure 10    Storm Sally, by contrast (as inferred earlier), was mostly stratiform rain when it passed over Wallops. Most of the DPR-retrieved index values (from Figure 7) were negative. Panel (a) of Figure 11 shows the variation of the index values with the PIA from the normal scan (NS; [7]) scan. Only data within a specified range of (off) nadir angles were used. Interestingly, one can see an increase in the i values with PIA. In other words, the N W -D m points move closer to the separation line as the PIA increases, although they lie mostly on the stratiform rain side of the separation line. Panel (b) of Figure 11 shows the QVP plot (a conical scan at 20 deg elevation, height up to 10 km, and range of around 30 km) of the radar reflectivity from the S-band NPOL volume scans on 18 September 2020. The clear radar bright-band seen at 4-4.2 km height throughout the event confirms that the event was dominated by stratiform rain. This validates the predominantly negative i values from panel (c), which were derived from DPR-retrieved [N W , D m ]. The increase in i values with PIA in panel (a) is also consistent with another case event analyses using disdrometer data and vertically pointing X-band radar during a cold-rain event in Ontario, Canada (see Section 2 in Thurai et al. 2016), which showed an increase in the (negative values of) i with increasing bright-band thickness. Both the Ontario event analysis and the HSV event in Figure 11a show a very similar trend. the event was dominated by stratiform rain. This validates the predominantly negative i values from panel (c), which were derived from DPR-retrieved [NW, Dm]. The increase in i values with PIA in panel (a) is also consistent with another case event analyses using disdrometer data and vertically pointing X-band radar during a cold-rain event in Ontario, Canada (see Section 2 in Thurai et al. 2016), which showed an increase in the (negative values of) i with increasing bright-band thickness. Both the Ontario event analysis and the HSV event in Figure 11a show a very similar trend.

Summary and Discussion
Our results have shown that it is feasible to retrieve DSD moments directly from DPR-based quantities. The retrieved moments can act as an initial step, prior to deriving DSD parameters such as NW and Dm. Additionally, the DSD moments are very helpful in understanding the microphysical processes involved in various types of rain regimes. The basic premise of this article is that the lower order moments (M0, M2, M3) of the DSD, which are involved in the formulation of the warm rain process rates (such as evaporation, sedimentation, and binary collisions), can be expressed as products of power laws of the a priori chosen reference moments of higher order (M3, M6) along with moments of h(x), where h(x) is the intrinsic shape of the distribution. The form of h(x) can be any functional shape, provided its moments are finite and positive. The latter formulation also enables estimation of the standard deviation of the mass spectrum σM [67], as well as NW and Dm.
In our study, the retrieval of moments used only two DPR products, namely the attenuation-corrected radar reflectivity at Ku-band and the specific attenuation at Ka-band. They were used to retrieve the sixth (M6) and the third moments (M3), respectively. Using these two quantities as our chosen reference moments, and assuming an underlying shape function, h(x), for the DSDs, it is possible to retrieve other moments, ranging from M0 to M7. Note that h(x) is dependent on the pair of chosen moments, and we have assumed the generalized gamma model to represent h(x). Using 3000 3 min full DSD spectra from two

Summary and Discussion
Our results have shown that it is feasible to retrieve DSD moments directly from DPR-based quantities. The retrieved moments can act as an initial step, prior to deriving DSD parameters such as N W and D m . Additionally, the DSD moments are very helpful in understanding the microphysical processes involved in various types of rain regimes. The basic premise of this article is that the lower order moments (M 0 , M 2 , M 3 ) of the DSD, which are involved in the formulation of the warm rain process rates (such as evaporation, sedimentation, and binary collisions), can be expressed as products of power laws of the a priori chosen reference moments of higher order (M 3 , M 6 ) along with moments of h(x), where h(x) is the intrinsic shape of the distribution. The form of h(x) can be any functional shape, provided its moments are finite and positive. The latter formulation also enables estimation of the standard deviation of the mass spectrum σ M [67], as well as N W and D m .
In our study, the retrieval of moments used only two DPR products, namely the attenuation-corrected radar reflectivity at Ku-band and the specific attenuation at Ka-band. They were used to retrieve the sixth (M 6 ) and the third moments (M 3 ), respectively. Using these two quantities as our chosen reference moments, and assuming an underlying shape function, h(x), for the DSDs, it is possible to retrieve other moments, ranging from M 0 to M 7 . Note that h(x) is dependent on the pair of chosen moments, and we have assumed the generalized gamma model to represent h(x). Using 3000 3 min full DSD spectra from two (climatically) different locations, we found that the underlying function corresponding to [M 3 , M 6 ] pair seemed reasonably stable.
Simulations using the 3000 3 min DSDs showed that the algorithm errors for the retrievals would be low; however, for lower-order moments, especially for M 0 , the errors were found to be more significant. When applying the technique to DPR data, other error sources also need to be considered. They include variance of radar measurement errors and 'point-to-area' variance, which arises from the fact that disdrometer data represent point measurements (although over a 3 min period), whereas DPR data are over a much larger footprint (although instantaneous), and factors such as non-uniform beam filling also need to be considered. We also need to bear in mind that the DPR has a lower threshold of around 12 dBZ for Ku-band measurements, which will put a lower limit of the estimated M 6 used in this method.
Two GPM-DPR overpass cases were used to examine the validity of our retrieval technique, one over the Huntsville site on 11 April 2016 at 23:31 UTC and one over the Wallops site on 18 September 2020 at 04:42 UTC during the passage of remnants of storm Sally. We compared the mean and standard deviation of all moments (M 0 to M 7 ) derived from the DPR data in pixels in the 0.05 degree of latitude and 0.05 degree longitude around the disdrometer site and compared them with those derived directly from the DSD within a 30 min time interval around the DPR overpass. The agreement for all moments (in terms of log 10 ) was well within 10%. Furthermore, values of D m , given by the ratio of the fourth moment to the third, ranged from 0.9 to 1.8 mm for the DPR-based retrievals around the disdrometer site and 1.1 to 1.8 mm for the DSD-based estimates over a ±30 min period. For the second case, histograms of the moments were derived from DPR and DSDs, and when the DPR-based moments were restricted to areas over and surrounding the disdrometer location, their modes showed very good agreement with the DSD-based mean values, which determined when the disdrometer data were restricted to within 2 h of the DPR overpass. D m and σ M histograms were also compared for both events. Of the two, remnants of hurricane Sally showed narrower D m histograms, both DPR-retrieved and DSD-based, having lower D m values with a maximum value of only 1.6 mm, whereas the HSV event shows D m values ranging up to 2.1 mm. The latter, in fact, showed two peaks in the D m , and these bimodal peaks are evident in both the DPR-retrieved and the DSD-based histograms. They were attributed to the semi-organized convective line being embedded within a larger widespread (and seemingly stratiform) rain region. The bimodal peaks were also noticeable in the σ M histograms. Sally (around the Wallops disdrometer area), on the other hand, shows only one peak in both D m and σ M . The tropical environment may also be the reason that no larger drop mode was observed by the disdrometers or DPR.
Finally, a previously published method for separating stratiform and convective rain types was tested for the two events. The method is based on N w versus D m and has been tested using disdrometer data in several locations, e.g., [60][61][62], and against the texturebased method used in [63]. For the two events considered here, the DPR-based N W versus D m (again note considerably larger footprint) were used to identify regions of the two rain types. For the HSV event, the convective rain regions corresponded well with regions of high PIA from the matched scan. On the other hand, DPR data for Sally showed mostly stratiform rain during the overpass. This was confirmed by QVP constructed from NPOL radar for the event.
The two initial test cases show encouraging results for more informative and accurate DSD retrievals from DPR. Our retrieval method for DSD moments followed a procedure similar to previous work by Bringi et al. [31], using X-band ground-based polarimetric radar data. The third and the sixth moments had also been used for that study. It is likely that these are also applicable for the DPR based retrievals, but further investigations are needed. Another aspect to consider is the stability of h(x), specifically, whether having two sets of h(x), one for stratiform rain and one for convective rain would improve the overall retrievals. This will also be examined in the future.
The novelty of our method lies in the conceptual fact that the rain DSDs are better characterized in terms of moments rather than the DSD parameters governing the distributions. As mentioned earlier in the Introduction, DSD moments are far more informative about the microphysical processes involved in the observation regions. Given that the DPR is capable of providing vertical profiles of attenuation corrected radar reflectivities as well as specific attenuations at Ku and Ka bands, the height variations of the DSD moments will lead to better understanding of the dominant processes. One drawback to our method is that the Ka-band specific attenuation, as mentioned several times earlier, has not been made available as one of the 'official DPR products'. It is hoped that this will be included in the near future. The specific attenuation can also be used as one of the rain-rate estimators, as has been the base for S-band ground-based radars [68] as well as X-band ground-based radars [69].    Figure A1 and those from the ground measurements in maroon. Figure A1. Histograms of moments (units are mm n /m 3 for the nth moment) shown in green derived from DPR data over and around the disdrometer site for the 11 April 2016 event in Huntsville ( Figure 5). The maroon arrows represent those derived from the ground-based DSD measurements within ±30 min period around the overpass time.
Remote Sens. 2021, 13, x FOR PEER REVIEW 20 of 23 Figure A1. Histograms of moments (units are mm n /m 3 for the n th moment) shown in green derived from DPR data over and around the disdrometer site for the 11 April 2016 event in Huntsville (Figure 5). The maroon arrows represent those derived from the ground-based DSD measurements within ±30 min period around the overpass time. Figure A2. DSDs shown in green reconstructed from the moments in Figure A1 and those from the ground measurements in maroon. Figure A2. DSDs shown in green reconstructed from the moments in Figure A1 and those from the ground measurements in maroon.