Retrieval of the raindrop size distribution from polarimetric radar data using double-moment normalisation

A new technique for estimating the raindrop size distribution (DSD) from polarimetric radar data is proposed. Two statistical moments of the DSD are estimated from polarimetric variables, and the DSD is reconstructed using a double-moment normalisation. The technique takes advantage of the relative invariance of the double-moment normalised DSD. The method was tested using X-band radar data and networks of disdrometers in three different climatic regions. Radar-derived estimates of the DSD compare reasonably well to observations. In the three tested domains, in terms of DSD moments, rain rate, and characteristic drop diameter, the proposed method performs similarly to and often better than a state-of-the-art DSD-retrieval technique. The approach is flexible because no specific DSD model is prescribed. In addition, a method is proposed to treat noisy radar data to improve DSD-retrieval performance with radar measurements.


Introduction
The raindrop size distribution (DSD) describes the microstructure of liquid precipitation and is highly variable (Jameson and Kostinski, 2001;Uijlenhoet et al., 2003;Tapiador et al., 2010;Jaffrain and Berne, 2012).The DSD is measured at the point scale by disdrometers.For applications such as numerical weather prediction (e.g.Baldauf et al., 2011) or radar remote sensing (e.g.Bringi and Chandrasekar, 2001), it is often necessary to know the areal DSD at the pixel scale.In other cases, such as studies of the microphysics of precipitation (Pruppacher and Klett, 2000;Tapiador et al., 2014), it would be useful to be able to remotely infer the DSD aloft or in remote locations.For these reasons, retrieval of the DSD from radar data has been a long-standing goal.In this paper we present a new technique for DSD retrieval from polarimetric radar data, which is based on the doublemoment normalisation technique of Lee et al. (2004).
Polarimetric weather radars are particularly useful for remote retrieval of the DSD, because differences between vertically and horizontally polarised electromagnetic waves reflected off hydrometeors in the atmosphere provide information on the particles' concentration, size, and shape.In rainfall, radar reflectivity in horizontal (Z H , dBZ) or vertical (Z V , dBZ) polarisation primarily relates to drop concentration and size, differential reflectivity (Z DR , dB) reflects drop shape, and specific differential phase shift on propagation (K dp , • km −1 ) relates to both the concentration and shape of the drops (Bringi and Chandrasekar, 2001).Seliga and Bringi (1976) showed that Z DR can be linked to the median volume drop diameter, a microphysical property of rain.Since then, many methods for DSD retrieval from radar variables have been proposed.Zhang et al. (2001) introduced the "constrained gamma" method, in which the shape and slope parameters of a gamma DSD model (Ulbrich, 1983) are assumed dependent.This assumption is subject to debate (e.g.Zhang et al., 2003;Atlas and Ulbrich, 2006;Moisseev and Chandrasekar, 2007;Cao and Zhang, 2009).The technique, modified by Brandes et al. (2003), can provide useful DSD information (Brandes et al., 2004a).In the "beta" method (Gorgucci et al., 2002), the effective slope of the drop axis ratio to diameter relationship is retrieved.The slope is used to find parameter values for the normalised gamma model of Willis (1984), which has advantages for use with polarimetric observations Published by Copernicus Publications on behalf of the European Geosciences Union.
T. H. Raupach and A. Berne: Retrieval of the DSD from polarimetric data (Illingworth and Blackman, 2002).Retrieval of the gamma model shape parameter with the beta method is subject to high uncertainty (Gorgucci et al., 2002;Anagnostou et al., 2008).To deal with noisy Z DR and K dp data at low rain rates, Bringi et al. (2002Bringi et al. ( , 2003) ) used the beta method for heavy rain and disdrometer-based regressions on Z H and Z DR for light rain.Brandes et al. (2004b) found that the constrained gamma method was in better agreement with disdrometer data than the beta method, while Anagnostou et al. (2008) reported similar performance from the two techniques, and both studies noted that the beta method is sensitive to errors in K dp .Vulpiani et al. (2006) developed a neural-network DSD-retrieval technique, and spatial correlations of DSD model parameters have been retrieved from radar data (Thurai et al., 2012;Bringi et al., 2015).
X-band polarimetric weather radars are popular due to their portability, small size, and high resolution and sensitivity, but measurements at X band suffer from attenuation by heavy rain (Anagnostou et al., 2013;Kalogiros et al., 2013) and must be corrected (Matrosov et al., 2005;Park et al., 2005a).Several DSD-retrieval algorithms have been developed for X band (e.g. Park et al., 2005b;Gorgucci et al., 2008;Anagnostou et al., 2013;Kalogiros et al., 2013), including some with integrated attenuation correction (e.g.Testud et al., 2000;Yoshikawa et al., 2014).The selfconsistent with optimal parameterization attenuation correction and microphysics estimation (SCOP-ME) algorithm, developed through studies by Anagnostou et al. (2009Anagnostou et al. ( , 2010) ) and Kalogiros et al. (2013), uses relationships calculated for the Rayleigh limit, corrected for Mie scattering at X band.It performs well compared to contemporary algorithms and disdrometer observations (Anagnostou et al., 2013).In this paper we present a new method for DSD retrieval that uses the double-moment DSD normalisation of Lee et al. (2004), we and compare it to SCOP-ME.
The rest of this manuscript is organised as follows: we briefly describe the double-moment DSD normalisation technique of Lee et al. (2004) in Sect. 2. Bulk rainfall variables that we use are introduced in Sect.3. Data used are presented in Sect. 4. In Sect. 5 we propose a new DSD-retrieval method that uses double-moment normalisation to retrieve the DSD from polarimetric radar data.Its performance is compared to that of SCOP-ME using radar variables simulated from DSD measurements in Sect.6.In Sect.7 we introduce a new method to reduce the effects of noise in radar measurements.Using this method, the DSD-retrieval algorithms are compared using radar data in Sect.8. Conclusions are made in Sect.9.

Double-moment DSD normalisation
The DSD is written as N (D) (mm −1 m −3 ) and is defined as the concentration in air of raindrops with equivolume diameter in the interval from D to D + δD mm.The equivol-ume diameter is used because raindrops become oblate with size (e.g.Thurai et al., 2007); it is simply the diameter of a sphere that contains the same volume of water as a drop.M n (mm n m −3 ), the nth-order moment of the DSD, is (1) The double-moment normalisation method of Lee et al. (2004) allows for the DSD to be expressed as a combination of two of its moments M i and M j of arbitrary orders i and j and a double-moment normalised DSD h(x) (-), where x = DM (-) is the second-normalised diameter (Lee et al., 2004).Using the normalisation, the DSD can be written as ( The method is flexible because the function h(x) is not prescribed.Lee et al. (2004) suggested that a generalised gamma model is an appropriate choice for h(x).Following their recommendation, we use the following double-moment normalised DSD (Lee et al., 2004): where is the gamma function, i = (µ + i/c) and j = (µ + j/c), and c (-) and µ (-) are parameters which must be fitted to the generalised gamma model.Since this formulation allows any DSD to be described using only two of its statistical moments, the task of our DSD-retrieval algorithm is to estimate two DSD moments from polarimetric radar data.
The question of whether the double-moment normalised DSD is invariant has been investigated.Compared to previous single-moment normalisation approaches that vary by rainfall type (Sempere-Torres et al., 2000), the doublemoment approach shows more similarity across such changes (Lee et al., 2004).Raupach and Berne (2017) tested the double-moment normalised DSD across spatial displacement and between different climatic regions.They showed that for practical purposes in stratiform rain and with wellchosen input moments, the double-moment DSD can be considered invariant across space with acceptable resulting performance on reconstruction of the DSD. Lee et al. (2007) showed that h(x) derived from time series measurements at one location had low scatter around the average doublemoment normalised DSD.In the DSD-retrieval method proposed here, we make the assumption that the double-moment normalised DSD function h(x) is invariant in space and time over the typical domain of interest and that variance in the DSD is adequately explained through variance in two moments of the DSD.

Bulk rainfall variables
All bulk rainfall variables can be derived from the DSD (a detailed review is provided by Bringi and Chandrasekar, 2001).The mass-weighted mean drop diameter D m (mm), useful as a characteristic drop size, is M 4 /M 3 .Liquid water content W (g m −3 ) is related to the third moment of the DSD and is written as where ρ w (g cm −3 ) is the density of water.The rain rate R (mm h −1 ) is defined as where v(D) (m s −1 ) is the still-air terminal fall speed of a drop with equivolume diameter D. In this study v(D) was calculated using the method of Beard (1976), for site-specific altitudes and latitudes, and an assumed sea-level temperature of 15 • and relative humidity of 0.95.
Radar variables can also be derived from the DSD.In Rayleigh scattering, when the radar wavelength is much larger than the particles being measured and drops are assumed to be spherical, the radar reflectivity is Z = M 6 (Marshall et al., 1947).In Mie scattering, in which the wavelength is of similar size to the particles, reflectivity in horizontal polarisation Z h (mm 6 m −3 ) is defined as (Bringi and Chandrasekar, 2001) where λ (cm) is the wavelength, |K| 2 (-) is the dielectric factor of water, and σ bh (D) (cm 2 ) is the back-scattering cross section at horizontal polarisation of a raindrop of equivolume diameter D. Reflectivity in vertical polarisation, Z v (mm 6 m −3 ), is obtained by replacing σ bh (D) with the vertically polarised back-scattering cross section σ bv (D) (cm 2 ).It is usual practice to deal with radar reflectivities in dBZ, calculated as Differential reflectivity in linear units, ξ dr (-), defined as Z h /Z v , has been shown to relate to the reflectivity-weighted mean drop axis ratio r z (-) (Jameson, 1983).r z is defined as where r(D) is the vertical to horizontal axis ratio of a drop of equivolume diameter D. The relationship found by Jameson (1983) is which is valid for narrow distributions of the raindrop axis ratio (Bringi and Chandrasekar, 2001).Dual-polarisation radars measure specific differential phase shift (on propagation) K dp ( • km −1 ), which is the difference in phase produced between horizontally and vertically polarised waves that pass through rain.It is defined as (Bringi and Chandrasekar, 2001) where Re represents the real part of a complex number and Re(f hh ) (cm) and Re(f vv ) (cm) are the real parts of the forward scattering amplitudes for horizontal and vertical polarisation respectively.Jameson (1985) showed that K dp can be linked to the product of liquid water content and the deviation from unity of the mass-weighted mean raindrop axis ratio r m (-).r m is defined as K dp can be written as with dimensionless value C ∼ 3.75 (Bringi and Chandrasekar, 2001).Various raindrop axis ratio functions are available (e.g.Pruppacher and Beard, 1970;Andsager et al., 1999;Brandes et al., 2002;Thurai et al., 2007).We return to the question of axis ratios and K dp in Sect. 5.The integrals in this section and Eq. ( 1) are idealised because the range of drop sizes is written from zero to infinity.Using measured data, the integrals were calculated over truncated classes of diameter and second-normalised diameter, with D and x as class centres and dD and dx as class widths.Since truncation potentially affects bulk variables (e.g.Willis, 1984;Ulbrich, 1985;Vivekanandan et al., 2004), we used the same truncation limits for compared quantities.When polarimetric variables were calculated from DSDs, the T-matrix codes of Mishchenko and Travis (1998) were used to calculate raindrop scattering properties, with an assumed temperature of 12.5 • C, a Gaussian distribution of raindrop canting angles with zero mean and a standard deviation of 6 • (stated as reasonable by Bringi and Chandrasekar, 2001), and a radar frequency of 9.4 GHz.
T. H. Raupach and A. Berne: Retrieval of the DSD from polarimetric data 4 Data To train and test the new method, data from three networks of OTT Parsivel (Löffler-Mang and Joss, 2000) disdrometers were used.Each network had a nearby X-band weather radar that scanned above the disdrometers.A full description of the data and their treatment, and the coordinates for all stations, are provided in Raupach and Berne (2017), in which the same disdrometer networks were used.Here we provide a summary of the data used in this study.
The first network provided the HyMeX data set.This network was located in Ardèche, France, in the autumns of 2012 and 2013, for the special observation periods of the Hydrological Cycle in the Mediterranean Experiment (HyMeX 1 , Drobinski et al., 2014).In this study we used data from 11 first-generation Parsivel and 5 Parsivel 2 .disdrometers located in the approximately 13 × 7 km 2 network.Also used were data from a METEK GmbH micro rain radar (MRR Peters et al., 2002Peters et al., , 2005;;Tridon et al., 2011) within the network, which provided vertical profiles of estimated DSDs recorded with 100 m vertical resolution and 10 s integration time.MXPol, a transportable Doppler dual-polarisation weather radar (for instrument details see Schneebeli et al., 2013), was located to the north-east of the disdrometer network.In 2013, MXPol recorded "stacked" plan position indicator (PPI) scans above the Parsivel network at elevations of 4, 5, 6, 8, 10, 12, 14, 16, and 20 • above horizontal, with a return time of about 6 min.Six rainfall events in which the MRR and MXPol both recorded data were selected for 2013.The events were from 1.8 to 7.5 h in length.Temperature data from a weather station at Pradel Grainage were used to estimate freezing level cutoff heights, below which precipitation was assumed to be primarily liquid.These heights ranged from 971 to 2386 m a.s.l., and only those MRR data from below the cutoff level per event were used.More network details and the list of identified events are provided in Raupach and Berne (2017).The HyMeX data set was the only set used in which MRR estimates of the DSD aloft were available.
Two more data sets were used in order to incorporate data from different climatologies.The second instrument network was composed of five first-generation Parsivel disdrometers and MXPol, in Payerne, Switzerland, and took measurements from February to July 2014.We used the MXPol PPI scan at 5 • above horizontal, which had a return time of about 5 min.The scans covered the region over three of the disdrometers.The third data set was from a network of 14 Parsivel 2 disdrometers (Petersen et al., 2014) deployed in Iowa, US, during the National Aeronautics and Space Administration Iowa Flood Studies (IFloodS) Global Precipitation Measurement ground validation campaign.Overlooking this network was the University of Iowa's X-band radar XPOL5 (Mishra et al., 2016).We used PPI data recorded at 3 • above horizontal, with a variable return time of about 2 to 8 min, for 3 days of heavy rainfall: the 25th, 26th, and 27th of May 2013.These scans covered the area over 10 of the disdrometers.The three networks were in regions with different climatologies (as described in Wolfensberger et al., 2016).Table 1 provides a summary of the three networks.Radar scans were matched to instrument times by finding scans for which at least half of the typical scanning period was within the instrument integration time.Scans typically lasted about 30 s in HyMeX, about 44 s in Payerne, and about 60 s in Iowa.Table 2 shows the instruments covered by PPI scans, the distance of each station to the PPI radar volumes used, and the number of radar scans that overlapped with 1 min observations.
Disdrometer data, which had raw integration times of either 30 or 60 s, and MRR data, with 10 s integration time, were resampled to 1 min temporal resolution.HyMeX and Payerne Parsivel data were corrected with reference to 2Dvideo-disdrometer (2DVD) measurements from the HyMeX campaign (Raupach and Berne, 2015a, b), with correction factors trained using reprocessed Parsivel data from the HyMeX 2013 campaign.This procedure removed unrealistically large drops and those too far from expected velocities, adjusted velocity measurements, and adjusted drop concentrations so that DSD moments more closely matched those of the 2DVD.These Parsivel data were quality controlled so Table 2. Instrument stations with corresponding PPI volumes, with the number of scans for that volume (S), the volume centre's height above the ground (H , m a.g.l., to nearest 10 m), height above sea level (H , m a.s.l., to nearest 10 m), and horizontal range from the radar (D, km).MI (mm h −1 ) is the maximum 1 min rain intensity recorded by each instrument at a radar scan time.that only error-free time steps containing liquid precipitation were used.Iowa Parsivel data were used as provided without further quality control.Parsivel data are subject to uncertainty due to differences across individual instruments and instrument generations (e.g.Jaffrain and Berne, 2011;Tokay et al., 2014;Thurai et al., 2011, Raupach andBerne, 2015a), and their limited sampling area introduces a bias, as reported by Tapiador et al. (2017).The Iowa data were provided in diameter class definitions that differed from those of the instrument manufacturer (Petersen et al., 2014).The HyMeX and Payerne data sets used the manufacturer's diameter class definitions, which implies the assumption of a raindrop axis ratio to equivolume diameter relationship (Battaglia et al., 2010).Our tests (not shown) showed limited differences made to DSD bulk variables when different axis ratio functions were used to modify the class definitions.Given the uncertainties involved in using modified diameter classes, we decided to use the manufacturer's class definitions for these two data sets.For each of the three regions, the Parsivel data were randomly sampled so that 60 % of records formed a training data set and the remaining 40 % formed an independent validation data set.Sensitivity of the random sampling was evaluated through repeated tests with different sample realisations and was found to be low.
All available disdrometer and PPI data were used, while MRR data were subset to event times so that likely solid precipitation was not considered.MRR data were attenuationcorrected (METEK, 2010;Peters et al., 2010) and contained DSDs retrieved with vertical wind ignored (Strauch, 1976;Peters et al., 2002).Negative concentrations (METEK, 2010) in MRR DSDs were reset to zero.PPI radar reflectivities were compared to measurements from disdrometers (and the MRR in HyMeX), and bias in Z H was corrected on a percampaign basis.Bias in Z DR was estimated using vertical scans (birdbath scans, similar to Grazioli et al., 2015) and was corrected in each of the three data sets.Two days of T. H. Raupach and A. Berne: Retrieval of the DSD from polarimetric data radar data from Payerne (22 March 2014 and 8 April 2014) exhibited higher radar bias due to hardware problems and were not included in this study.Attenuation correction for the PPI data was performed using the ZPHI algorithm (Testud et al., 2000), and K dp was estimated using the method of Schneebeli et al. (2014).PPI scan data were sampled for instrument locations by taking the mean values of radar volumes that horizontally overlapped the instrument coordinates.To discount noise, PPI records were subset to those for which Z H was greater than or equal to 10 dBZ, and the signal-to-noise ratio in horizontal polarisation was greater than or equal to 5 dB.DSD data were treated as in Raupach and Berne (2017): Parsivel DSDs were truncated to 0.2495 (0.2565) to 7 (7.21)mm for HyMeX and Payerne (Iowa) Parsivel data (Raupach and Berne, 2015a); to avoid including overestimated numbers of small drops (Peters et al., 2005), DSDs estimated by the MRR were truncated to 0.6 to 5.8 mm (Raupach and Berne, 2017) and MRR data were further subset to records with R ≤ 150 mm h −1 (thus removing 0.2 % of records); MRR data for altitudes greater than 2250 m were excluded because not enough points were available at those altitudes, and all DSDs were subset to time steps in which R > 0.1 mm h −1 .In each data set, more than 85 % of the DSDs sampled were classified as stratiform type by Raupach and Berne (2017).
To compare measured vs. estimated or retrieved values in this work, we use the median relative bias, the interquartile range (IQR) of relative bias, and the squared Pearson correlation coefficient (r 2 ) between reference and estimated values.If V R is the reference value and V E is the estimated value, the relative bias expressed as a percentage of the reference value is defined as 100(V E − V R )/V R .

DSD retrieval from polarimetric radar data
Raupach and Berne (2017) showed that with reasonably chosen input moments, the double-moment normalised DSD of Lee et al. (2004) can be assumed invariant across spatial displacement in stratiform rain, with a performance loss that is acceptable for practical applications.Results on limited data for non-stratiform rain types suggested that, while the double-moment normalised DSD varies more in these cases, the assumption of its invariance may still lead to acceptable performance with input moments that are not both of low or both of high order.Using the assumption of an invariant double-moment normalised DSD model, the DSD can be estimated using polarimetric radar data.Given a known doublemoment normalised DSD, the task of DSD reconstruction becomes that of estimating from radar information the values of two DSD moments.In this section we present a new DSDretrieval method that uses this idea.The aim of the proposed DSD-retrieval technique is to retrieve two DSD moments using only polarimetric radar data.
The SCOP-ME method was trained with DSDs simulated using a DSD model and a wide range of DSD parameter values.In contrast, we used empirical DSDs measured by Parsivels to train our method to avoid any assumption about the shape of the DSD.A trade-off in using empirical DSDs is that measured DSDs are necessarily truncated by instrumental limitations.However, previous studies have shown that if the considered range of drop diameters is large enough around the median drop diameter D 0 (mm), the effect of truncation on calculated bulk variables is limited (Willis, 1984;Vivekanandan et al., 2004).Willis (1984) concluded that the effect of maximum considered drop size D max on bulk variables is negligible if D max exceeds 2.5D 0 .Using D 0 calculated from the recorded (truncated) Parsivel DSDs, this criteria was met for 99.6 % of the records.The criteria of Vivekanandan et al. (2004) is that, for there to be less than five percent error on bulk variables, the minimum drop size D min should be less than D 0 /2 and D max should exceed 4D 0 .This constraint was met by 90.5 % of the DSDs (93.5 % met this criteria for the upper drop size limit).Calculated D 0 may also be subject to error because of the truncation, but we consider that these calculations give broad confidence in the bulk variables we used to train the method.Further, the truncation on the Parsivel data affects primarily very small drops since large drops are rare, and therefore its influence on the higherorder moments we use is expected to be negligible.
The training data set was sampled as 60 % of each of the three Parsivel data sets and contained 181 829 measured 1 min DSDs.Z H , K dp , and Z DR were calculated for these DSDs for the MXPol stacked PPI incidence angles; temperatures of five, 10, and 15 • C; and each of four drop axis ratio functions: those of Andsager et al. (1999); Brandes et al. (2002); Thurai et al. (2007), and that of Beard and Chuang (1987) in the form shown in Kalogiros et al. (2013).Unusual records with Z DR or K dp less than or equal to zero (0.16 % of all simulated radar records) were excluded.

Retrieval of DSD moment six
Radar reflectivity in linear units, Z h (mm 6 m −3 ), is the sixth moment of the DSD in the Rayleigh scattering regime for spherical drops (Bringi and Chandrasekar, 2001).At X-band frequencies, larger drops enter into the Mie scattering regime and differences appear between M 6 and Z h .We use the observation that Z h departs from M 6 for heavier rain and assume that this departure occurs when Z H is greater than a threshold value.This threshold was determined through comparison of M 6 and Z h for DSDs, classed by Z H in classes of width 2 dBZ between 10 and 40 dBZ, and was set to 28 dBZ.For both smaller and larger reflectivity values, a power law relationship was found using orthogonal least squares fitting in log-log space.The resulting relationship is  13).The Z H threshold of 28 dBZ is shown with a triangle.
On the training set, median relative bias between M 6 and M 6 was 0.1 %, the IQR of relative bias was 2.5 percentage points, and the r 2 value was 0.98.The fitted relationship is shown on samples of training data in Fig. 1.Temperature made only limited difference to the fitted parameters: the prefactor varied from 2.45 to 2.90 for the larger values of Z H , and the other parameters differed by 0.01 or less from the value found for all temperatures combined.

Retrieval of DSD moment three
Retrieving a second, lower-order DSD moment is more difficult than estimating M 6 , because radar variables are more closely linked to the higher-order moments of the DSD.Using theoretical relationships as much as possible, we present a method to estimate the third moment of the DSD from polarimetric data.As shown in Eq. ( 9), the reflectivity-weighted mean drop axis ratio, r z , is related to a negative power of the differential reflectivity in linear units.In Kalogiros et al. (2013), the reflectivity-weighted and mass-weighted drop axis ratios were assumed to be the same and differences were dealt with through fitting of qualitative relationships between radar variables.A similar approach is taken here.Since r z and the mass-weighted mean drop axis ratio r m are both weighted mean drop axis ratios, we assume that r m is also related to differential reflectivity and estimate r m using a polynomial fit to Z DR , such that With our training data, this polynomial of order five produced low relative bias on retrieval of M 3 .Recall from Eq. ( 5) that M 3 relates to W : substituting Eq. ( 5) into Eq.( 12), and solving for M 3 , we have At X band (9.4 GHz, λ = 3.189 cm), assuming that ρ w = 1 g cm −3 , and replacing r m with its estimate based on Z DR , M 3 is predicted by where C is a single representative value for C. K dp is sensitive to the raindrop axis ratio (e.g.Bringi and Chandrasekar, 2001), so values for c i and C were found per axis ratio function.The coefficients c i in Eq. ( 14) were found using least squares polynomial fitting.In rare cases for large values of Z DR the relationships returned unrealistic values of r m (0 ≤ r m or r m > 1).In these few cases, r m was set to 0.75.Estimated r m values were used to find C for each training DSD, and the mean of these values was used as C. The results and their performance statistics are shown in Table 3. Fitted parameters differed across the three tested temperatures.However, parameters fitted using all training data performed similarly on training data for individual temperatures, with the median relative bias remaining within ±1 % of the all-temperatures value and IQR of relative bias varying by about one percentage point.The values fitted using combined training data were used.

Summary of DSD-retrieval technique
The proposed DSD-retrieval method is summarised as follows: the double-moment normalised DSD ĥ(x) with parameters c and µ is assumed trained from data and known.Then, given K dp , Z DR , and Z h , (1) DSD moment six is estimated using Eq. ( 13) and (2) DSD moment three is estimated using Eqs.( 14) and ( 16) and parameters from Table 3.The DSD is then retrieved using Eqs.( 3) and ( 4) with i = 3 and j = 6.

Comparison to an existing DSD-retrieval method
The new DSD-retrieval method was compared to SCOP-ME (Anagnostou et al., 2009(Anagnostou et al., , 2010;;Kalogiros et al., 2013).We implemented SCOP-ME using its description in Anagnostou et al. (2013).SCOP-ME was developed for X band using simulated DSDs and T-matrix simulations of radar variables, and in Anagnostou et al. (2013) it is shown to outperform the algorithms of Anagnostou et al. (2008) and www.atmos-meas-tech.net/10/2573/2017/2013) provides an explicit expression for rain rate using polarimetric variables, but, since we are interested in the whole DSD, in the following we compare R computed from reconstructed DSDs.The comparison of the two methods is first shown using Parsivel data in which the radar values were simulated using T-matrix codes and were therefore free of radar measurement noise.
Comparisons of the two techniques were made using the Parsivel validation data set composed of 40 % of the records from HyMeX, Payerne, and Iowa.For each 1 min DSD, Z h , K dp , and Z DR were calculated using T-matrix codes, for an elevation angle of 4 • above horizontal, and using each of the four drop axis ratio functions.For the double-moment technique, the generalised gamma model ĥ(x) (Eq.4) for i = 3 and j = 6 was used.ĥ(x) was fitted to non-zero median values of h(x) in classes of x with width 0.2, using weighted least squares fitting in log space, with each class weighted by n 4 C , where n C is the number of observations in the class.This is the same technique shown in Raupach and Berne (2017), except that the class weights increase more quickly with the number of class points.We found this different weighting improved the accuracy of the shape of the average reconstructed DSD, which is sensitive to the form of ĥ(x).A1, for Payerne in Table A2, and for Iowa in Table A3.The metrics used were median relative bias, IQR of relative bias, r 2 , and the slope of the linear regression on measured vs. reconstructed points.Differences between the performance metrics for the two techniques were calculated such that a negative difference indicates that the doublemoment technique performed better than SCOP-ME.These differences are shown visually in Fig. 4, in which red colours show negative differences.
In over half of the tested region, axis ratio function, and variable combinations, the double-moment technique produced a better median relative bias than the SCOP-ME technique, with an overall average difference of −0.53 percentage points.IQR of relative bias was usually slightly higher for the double-moment technique, with an average difference of 3.8 percentage points.ter plot slopes were usually similar for both techniques.The average differences across the three tested regions and four tested raindrop axis ratio functions are shown in Table 4. On average, the double-moment technique produced better median relative bias than SCOP-ME on R, D m , and DSD moments two to seven.IQRs were similar on average, with the exception of moments zero, one, and seven for which SCOP-ME produced notably smaller IQRs.As shown in Tables A1,  A2, and A3, the results differed across the different drop axis ratio functions and regions.It was often the case that SCOP-ME produced a less biased estimate of DSD moment zero, but in most of these cases the double-moment technique produced a better r 2 .The double-moment technique's performance variations relate to the accuracy of the prediction of DSD moment three from K dp and Z DR and to the fit of the generalised gamma function ĥ(x).ĥ(x) was trained on data from all data sets combined, in order to have the most general model possible.Our experiments showed that performance for low-order moments could be increased in any one region by training the gamma model on data from that re-gion only.This aligns with the conclusions of Raupach and Berne (2017), who noted that, while the double-normalised DSD can be assumed invariant for practical purposes, some residual variability remains and results in performance loss that depends on the input moments used.We now move to testing the two techniques on measured radar data, in which noise is a problem that must be dealt with.

Reducing the effects of noise
Radar data are noisy at light rain rates, particularly for K dp and Z DR (e.g.Bringi et al., 2002;Schneebeli et al., 2014).
We propose here a method to deal with this noise for the current application of DSD retrieval.Regressions on Z h and ξ dr are used to determine "expected" values for these variables, which can be used when the measured values are likely to be noisy.We found that Z DR can be reasonably predicted from ), R (mm h −1 ), and D m (mm) for the double-moment method, on the HyMeX data set, using the axis ratio function of Beard and Chuang (1987).One-to-one lines are shown in black.Regression lines for the double-moment method are shown in solid red, and dotted red lines show linear regressions for SCOP-ME, for which the densities are not shown.
and K dp can be predicted from Z h and ξ dr using with parameters α Z , β Z , α K , β K1 , and β K2 .Least squares fitting in log-log space, using the training data set described in Sect.5, was used to find best-fitting parameter values per raindrop axis ratio function.Just as for the retrieval of DSD moment six, assumed air temperature made only a small difference (parameter values fitted to individual temperature data sets differed by less than 6 % from those fitted using combined temperatures), whereas different axis ratios produced more diverse parameter values.Resulting parameter values and performance statistics are shown in Table 5. Threshold values are used to determine when K dp and Z DR may be noisy.A threshold value on Z H selects values of Z H for which K dp and Z DR showed large variation around their expected values in the three radar data sets used here.Threshold values on Z DR and K dp are those of Bringi et al. (2002).Although these were established for S band, their influence is rather limited compared to the threshold on Z H .The threshold on Z H is set to 37 dBZ; it is used in addition to the noise thresholds on Z DR and K dp in order to avoid effects Table 5. Fitted coefficients and the performance of the fits on the training data, for Eqs. ( 17) and ( 18), by raindrop axis ratio function (Ratio).Performance is shown in terms of median relative bias (RB, %) and the IQR (% pts) of relative bias.like ground clutter as well as noise that significantly affect any DSD-retrieval algorithm.A drawback of this method is that it reduces the benefit of observed polarimetric variables when Z H is under 37 dBZ, but the average error is kept low.
To reduce the effects of noise, then, if Z H < 37 dBZ or Z DR < 0.2 dB, measured Z DR is replaced by the expected value ẐDR and ξ dr is replaced by 10 ( ẐDR /10) .Likewise, if Z H < 37 dBZ or K dp < 0.3 • km −1 , K dp is replaced by Kdp (calculated with ξ dr if ξ dr was replaced).This treatment method allows radar data with negative or zero K dp or Z DR to be used.The treatment improved DSD-retrieval performance for both the double-moment and SCOP-ME techniques.For example, on PPI data with positive Z DR and K dp , when retrieved DSDs were matched to measured MRR data (Sect.8.1), the median relative bias was reduced by an average (across variables) of ∼ 8 percentage points for SCOP-ME and by ∼ 16 percentage points for the double-moment technique, while average IQRs were reduced more; for example on the comparison with MRR data the IQRs were reduced by ∼ 78 (80) percentage points for the SCOP-ME (double-moment) method.When retrieved DSDs were compared to Parsivel data (Sect.8.2), the noise in the radar data contributed to errors to such an extent that for both techniques the proposed treatment reduced the IQR and at times the median of relative bias by hundreds of percentage points for some variables.We note that because most values of Z H recorded in the PPIs analysed here were lower than 37 dBZ, the noise correction affected the majority of radar records.

Comparisons using radar data
The DSD-retrieval techniques were applied to PPI radar data from the three locations.The double-moment technique was run on noise-corrected data.SCOP-ME was run on uncorrected PPI data (subset to K dp > 0 and Z DR > 0) and noise-corrected data.We used the elevation angles of the stacked PPIs for HyMeX, 5 • for Payerne, and 3 • for Iowa.Measured radar variables Z H , K dp , and Z DR were recovered for volumes corresponding to instrument locations.DSD retrieval was performed using these values, and the resulting DSDs compared to those that were measured by other instruments.All comparisons using PPI data in-T.H. Raupach and A. Berne: Retrieval of the DSD from polarimetric data volved a difference in measurement volume -a change-ofsupport problem that we expect will introduce error spread (e.g.Raupach and Berne, 2016).There were, at times, significant vertical distances between the radar volume and the ground-based Parsivels used in these comparisons (see Table 2).Further, the disdrometers and MRR made integrated measurements, while the PPI scans made instantaneous measurements that could shift in time compared to the instrument's integration period (at worst, this shift could lead to the scan being up to half the typical scan length before or after the instrument integration time).We used the lowest possible radar elevation angles, but the possibility remains that at large distances the radar could have sampled solid precipitation above the ground instruments.These factors and uncertainty in the noise-correction and bias-correction techniques combine to create greater uncertainty in the comparisons of the two techniques made using real data than in those made using simulated radar variables from disdrometer data.
Because the axis ratio of Thurai et al. (2007) produced good results using the double-moment technique on the Parsivel data, the double-moment technique was used with parameters for this axis ratio function.Note that the assumption of axis ratio function affects only parameters of the doublemoment technique, because the radar data used in this section are measured, not simulated, and the SCOP-ME technique is used as presented in Anagnostou et al. (2013).In the HyMeX campaign, the lowest available PPI elevation angle (4 • ) was used to compare results to Parsivels, but there was also an MRR at Pradel Grainage which retrieved estimates of the DSD aloft.MRR-derived DSDs were compared at eight different altitudes using the MXPol stacked PPIs (except 20 • elevation) above the HyMeX instrument network.We first address the comparisons with MRR for HyMeX and then move to the comparisons with the Parsivel networks in all three regions.

Comparisons to MRR DSD estimates aloft
MXPol volume centre altitudes were projected onto MRR altitude classes for comparison.The double-moment DSDretrieval algorithm was used with generalised gamma model ĥ parameters (Eq.4) for MRR data and i = 3 and j = 6.These parameters were found using the same fitting technique as for Parsivel data (Sect.6) but differ since instrumental differences produce different forms of ĥ(x) (Raupach and Berne, 2017).The parameters were set to c = 1 and µ = 5.28.The use of stronger weighting in the fitting procedure reduced the influence of the large numbers of small drops returned by the DSD-retrieval algorithm used by the MRR, and thus the value of µ was smaller than those reported in Raupach and Berne (2017).The reconstructed DSDs were found for classes of drop diameter with class centres from 0.65 to 5.75 mm and a class width of 0.1 mm, so the reconstructed truncation matched that of the MRR data.PPI values from eight 100 m altitude classes between about 900 and 2100 m a.s.l. were compared to MRR estimates of the DSD aloft.Two output pairings are shown here: the first in which both techniques used noise-corrected data and the second in which the SCOP-ME technique used raw data and the double-moment technique used the same raw data set corrected for noise.This second pairing was made to ensure that the performance of SCOP-ME was not compromised by the noise-correction technique.
Results of comparisons between MRR-and PPI-derived DSDs are shown for three example altitudes in Fig. 5.There was good agreement between the recorded radar reflectivity recorded by both instruments, with a median relative bias of −3 %, an IQR on relative bias of 16 percentage points, and a value of r 2 of 0.68.The improvement in SCOP-ME performance made by the noise correction is clear.When both techniques used noise-corrected input, both overestimated DSD moment orders zero to four and underestimated orders six and seven.Rain rate was recovered with a median relative bias of 8 % (IQR 101 % pts) by the double-moment technique and 17 % (IQR 105 % pts) by SCOP-ME.The doublemoment technique showed lower median relative bias than SCOP-ME on moments one to five, seven, and R and smaller IQRs on moments two to seven and R. Similar to some of the Parsivel results, the double-moment technique overestimated moments zero and one of the DSD.r 2 values were low for both techniques (the maximum was 0.35, by SCOP-ME for D m ), but the double-moment technique had the same or a slightly higher value of r 2 in the majority of cases.High best-fit slopes were observed for both techniques for moments five, six, and seven and show the effect of a few outlier points in these cases.Performance differences between the two techniques using noise-corrected data are shown in Table A4.Overall, the double-moment technique for DSDretrieval out-performed SCOP-ME for the retrieval of DSD moments above order zero and rain rate measured aloft by the MRR.

Comparisons to DSDs measured by Parsivels
DSDs retrieved from polarimetric radar data were also compared to those recorded by ground-based Parsivels in the three regions we studied.Unlike in previous sections where training and validation divisions of the Parsivel data were used, here we compared DSDs derived using independent radar data to all available matching Parsivel records.The DSDs were retrieved in truncated Parsivel drop diameter classes, using the Parsivel generalised gamma model parameters quoted in Sect.6. Particularly in the Payerne data set, the noise-correction routine was required in order to retrieve realistic DSDs; the results shown here are thus for the SCOP-ME and double-moment techniques both run on noise-corrected PPI data.Figure 6 shows distributions of DSD-retrieval relative error for each region.
The double-moment technique produced smaller IQRs of relative bias than SCOP-ME for moments four to seven.For moment orders zero to three, the double-moment technique produced better median relative bias than SCOP-ME in the HyMeX and Iowa data sets, but worse in Payerne.Where the double-moment technique produced better median relative bias, the average improvement was of seven percentage points; in cases where SCOP-ME performed better, the average improvement was five percentage points.Values of r 2 and scatter plot slope were similar between the two techniques, with the majority of cases showing differences of less than 0.05 for both variables.Differences in performance between the two techniques are shown in Fig. 7 and Table A4.
The performance of the double-moment technique is reliant on how accurately two DSD moments can be extracted from radar data and in turn on how accurate the radar data are.Both retrieval techniques appear to be similarly affected by radar inaccuracies, and experiments with different reflectivity bias corrections (not shown here) showed similar patterns of results.In Parsivel comparisons, the proposed DSD-retrieval technique was applied using a single doublemoment normalised DSD model in all three tested regions, without significant performance loss between regions.This supports previous findings (Raupach and Berne, 2017) that for practical use with real radar data in primarily stratiform rain, the double-moment normalised DSD may be considered invariant in regions at similar latitudes.Figure 7. Differences in performance between the double-moment technique and SCOP-ME using noise-corrected radar data, for MRR and for Parsivels by region (differences in Table A4).Variables and performance statistics as for Fig. 4. Red indicates that the double-moment technique outperformed SCOP-ME.Grey indicates an r 2 difference greater than 50 on this scale; these slopes were affected by outliers.

Conclusions
Given the assumption of an invariant normalised DSD, and an estimate of that function, the DSD can be predicted using only two of its moments using the double-moment normalisation method of Lee et al. (2004).Two DSD moments are available from polarimetric radar data.At X band, radar reflectivity can be used to accurately predict the sixth moment of the DSD, and moment three can be retrieved relatively accurately using K dp and Z DR .We showed that by estimating these two DSD moments from radar data, the DSD for a radar volume can be predicted using the double-moment formulation.Tests on disdrometer data from three networks in different climatic regions showed that DSD retrieval using this new technique produced similar or slightly better performance than the SCOP-ME DSD-retrieval technique of Kalogiros et al. (2013).The proposed method is also more flexible, because there is no prescribed functional form for the double-moment normalised DSD, and even a non-parametric ĥ(x) could be used.The shape of the retrieved DSD is sensitive to the fitted form of ĥ(x), and better fitting methods may in future improve the performance of the retrieval method.
The flexibility of the technique extends to there being no prescribed method for DSD moment extraction, which means that the moments used could be tailored to the intended purpose.
A new method for treatment of radar data with possibly noisy values of K dp and Z DR was proposed.The method is based on predicting the expected values of these variables from radar reflectivity, and it considerably improved the performance of both the DSD-retrieval techniques when real radar data were used.DSDs were predicted from polarimetric variables in noise-corrected PPI scans measured by X-band radars in each of the three regions.Comparisons of the retrieved DSDs to MRR data for DSDs aloft in the HyMeX region in France, and of radar-retrieved DSDs to disdrometer data from the three regions, showed reasonable agreement but large error spread for both methods.This study provides a proof of concept for DSD retrieval using noise-corrected radar data, the double-moment normalisation method of Lee et al. (2004), and a generalised gamma model for the normalised DSD.Performance improvements may be possible through future work that should test the approach using different instruments and data sets, including testing whether such extensive noise correction is always required.Future work should also address more precise prediction of loworder DSD moments from polarimetric radar data and investigate different models and fitting methods for the doublemoment normalised DSD.
Data availability.HyMeX MXPol, MRR, and Parsivel data are available through the HyMeX database at http://www.hymex.org,and DOIs for individual data sets are listed in Nord et al. (2017).Weather station data, Payerne radar data, and Payerne Parsivel data are publicly available upon request from EPFL-LTE.IFloodS Xband radar data is publicly available in Krajewski et al. (2016).IFloods Parsivel data is publicly available in Petersen et al. (2014).
Table A1.Comparison of double-moment method to SCOP-ME results on all Parsivel data in the HyMeX data set by axis ratio function (Ratio).RB (%) is median relative bias, IQR (% pts) is interquartile range of relative bias (% points), r 2 is squared correlation coefficient.S is the slope on measured vs. reconstructed regression.Difference is the difference in absolute values for RB and IQR and the difference in deviation from unity for r 2 and slope.A negative difference shows that the double-moment method improved on SCOP-ME's performance.

Figure 1 .
Figure 1.A sample of 20 000 points from the training set, showing the relationship between radar reflectivity and DSD moment six in dB scale.The one-to-one line is shown in black; the red dashed line shows the fitted relationship of Eq. (13).The Z H threshold of 28 dBZ is shown with a triangle.
The parameters found for the combined Parsivel training data were c = 1.69 and µ = 2.22.SCOP-ME and the double-moment method were used to retrieve the DSD concentrations N (D) for D in the class centres of the truncated Parsivel diameter classes.For each technique and axis ratio function, retrieved DSDs were compared to measured DSDs by comparing moments zero to seven, D m , and R. Comparisons of relative error distributions by technique are shown in Fig. 2. Example scatter plot results are shown for the HyMeX data set and the drop axis ratio model of Beard and Chuang (1987) in Fig. 3.The Beard model, which

Figure 2 .
Figure 2. Relative bias distributions for the double-moment and SCOP-ME DSD-retrieval methods, by drop axis ratio function and data set (H stands for HyMeX, P for Payerne, and I for Iowa).Variables are moment order n (mm n m −3 ), D m (mm), and R (mm h −1 ).Bold bars show medians, boxes show IQRs, and whiskers show 10th to 90th percentile ranges.

Figure 3 .
Figure3.Density scatter plots of retrieved vs. measured moments M n (mm n m −3 ), R (mm h −1 ), and D m (mm) for the double-moment method, on the HyMeX data set, using the axis ratio function ofBeard and Chuang (1987).One-to-one lines are shown in black.Regression lines for the double-moment method are shown in solid red, and dotted red lines show linear regressions for SCOP-ME, for which the densities are not shown.

Figure 4 .
Figure4.Differences in performance between the double-moment technique and SCOP-ME, using radar variables simulated from Parsivel data, by region and drop axis ratio function (differences in Tables A1, A2, and A3).Reds indicate negative differences, where the double-moment technique outperformed SCOP-ME.Variables are moment order n (mm n m −3 ), D m (mm), and R (mm h −1 ).Differences are shown for median relative bias (RB, % pts), IQR of relative bias (IQR, % pts), r 2 (difference in deviations from unity, multiplied by 100 for display on this scale), and regression slope (S; difference in deviations from unity, multiplied by 100).

Figure 5 .Figure 6 .
Figure 5. Distributions of relative bias on DSD moments zero to seven, comparing DSDs retrieved using PPI data to those measured by the MRR at Pradel Grainage.The results are classed by altitude for a selection of three altitudes across the compared range.Two comparisons are shown: in comparison A, SCOP-ME used raw PPI data and the double-moment technique used the noise-corrected version of the same data set.In B, both techniques used noise-corrected data sets.Variables and symbols as for Fig. 2. Input data for comparison A excluded records with negative or zero K dp or Z DR , while the noise-correction technique deals with these values and they were included in comparison B.

Table 1 .
Summary of instrument networks used.Coordinates for Parsivel networks are bounding boxes.Altitude is above sea level to nearest 10 m.Hours are provided only for non-instantaneous measurements and show total hours of rain data measured by each network.

Table 3 .
Anagnostou et al. (2013)16) and c i (Eq.14), by drop axis ratio function (Ratio).M 3 estimation performance in the training data is shown in terms of median relative bias (RB, %), IQR of relative bias (IQR, % pts), and r 2 .Max Z DR (dB) shows the maximum value of Z DR each relationship can use.ME equations in which K dp is not used, which are to be employed when K dp is absent or close to noise.In this work we used only the version of SCOP-ME that uses K dp , as presented inAnagnostou et al. (2013), and we use it only with positive and non-zero values of K dp and Z DR (tests on our data sets showed better SCOP-ME performance with this configuration).Kalogiros et al. (

Table 4 .
Average differences between double-moment and SCOP-ME techniques, on Parsivel data, over three regions and four axis ratio functions.Negative values show an improvement by the doublemoment technique over SCOP-ME.
(Kalogiros et al., 2013)well to observations(Thurai et al., 2009), is shown because it provided the equilibrium drop shapes around which the SCOP-ME training set was simulated(Kalogiros et al., 2013).Because we are using empirical DSDs, some extreme values of D m are shown in this figure.Values of D m above 5 mm are extremely rare; less than 0.02 % of DSDs in each data set show these values, and they have a negligible influence on the regression lines.Full performance results are shown for the HyMeX data set in Table Correlation coefficients and scat-

Table A4 .
Differences in performance by variable and region, for DSDs retrieved from noise-corrected PPI data using the double-moment technique and SCOP-ME, compared to the MRR at Pradel Grainage (MRR) and Parsivels by region (HyMeX, Payerne, and Iowa).Metrics and differences are defined as for TableA1.An exception is Z, which refers to Z H measured by the radar, not reconstructed through DSD retrieval (hence it is the same for both techniques).