Triple‐Frequency Doppler Retrieval of Characteristic Raindrop Size

A retrieval for characteristic raindrop size and width of the drop size distribution (DSD) based on triple‐frequency vertical Doppler radar measurements is developed. The algorithm exploits a statistical relation that maps measurements of the differential Doppler velocities at X and Ka and at Ka and W bands into the two aforementioned DSD moments. The statistical mapping has been founded on 7,900 hr of disdrometer‐observed DSDs and their simulated Doppler velocities. Additionally, a retrieval of Dm based only on DDVX−W measurements is also presented, and its performance is compared to the analogous algorithm exploiting DDVKa−W data. The retrievals are tested using triple‐frequency radar data collected during a recent field campaign held at the Juelich Observatory for Cloud Evolution (JOYCE, Germany) where in situ measurements of the DSD were carried out only few meters away from the vertically pointing radars. The triple‐frequency retrieval is able to obtain Dm with an uncertainty below 25% for Dm ranging from 0.7 to 2.4 mm. Compared to previously published dual‐frequency retrievals, the third frequency does not improve the retrieval for small Dm ( <1.4 mm). However, it significantly surpasses the DDVKa−W algorithm for larger Dm (20% versus 50% bias at 2.25 mm). Also compared to DDVX−W method, the triple‐frequency retrieval is found to provide an improvement of 15% in terms of bias for Dm=2.25 mm. The triple‐frequency retrieval of σm performs with an uncertainty of 20–50% for 0.2


Introduction
Knowledge of the raindrop size distribution (DSD) and its vertical evolution is paramount for a wide variety of applications ranging from improving remote sensing quantitative precipitation estimates (e.g., Bringi & Chandrasekar, 2001, and references therein) to better understand rain processes and their description in atmospheric models (Kumjian & Prat, 2014). In particular, the characterization of DSDs for large characteristic size is of specific relevance for the parameterization of breakup which impacts upon storm properties like the propagation speed, the cold pool strength, and the total rain accumulation (Morrison et al., 2012).
The last decade has seen a wealth of DSD studies involving multifrequency Williams et al., 2016) and polarimetric Doppler radars (Kumjian, 2018). While vertically pointing multifrequency radar observations provide a detailed description of DSD variability at very fine spatial and temporal scales, scanning polarimetric radars are unrivaled for mapping the DSD properties over wide areas but at much coarser resolutions . An increasing number of research sites are routinely running vertically pointing multifrequency Doppler radars with frequencies ranging from X to W band. Dual-frequency radars (Ka/W) have been operated within the Atmospheric Radiation Measurement program (Mather & Voyles, 2013) at permanent supersites such as Southern Great Plains or North Slope of Alaska for several years. Several European sites, such as the Jülich Observatory for Cloud Evolution Core Facility (JOYCE-CF, X/Ka/W; Dias Neto et al., 2019;Löhnert et al., 2015), the Chilbolton Observatory (S/Ka/W; Mason et al., 2019), or the Hyytiälä Forestry Field Station (C/Ka/W; Falconi et al., 2018;Kneifel et al., 2016) started to routinely record triple-frequency radar data. Similarly, suites of multifrequency (three or more) radar observations are now available from the NASA DC8 and ER-2 airborne platforms (and references therein Battaglia et al., 2016;Kulie et al., 2014;. Multifrequency approaches with Doppler capabilities are currently pursued for spaceborne configurations as well (Durden et al., 2016;Leinonen et al., 2015).

10.1029/2019EA000789
Multifrequency radar rain microphysics retrievals exploit measurements from radars with wavelengths ranging from decimeters to millimeters by capitalizing on the diversity of scattering properties of raindrops (characterized by maximum sizes ≤9 mm). While raindrops can be considered Rayleigh targets at S, C, and X bands, the backscattering cross sections significantly deviate away from the ∝ D 6 behavior at Ka and W bands due to non-Rayleigh effects (Kollias et al., 2007). In addition to backscattering effects, attenuation becomes increasingly important with higher frequencies (Lhermitte, 1990). Both features have been widely used in reflectivity-based retrievals, with the major obstacle remaining the separation of attenuation from non-Rayleigh effects (Tridon et al., 2013). Such separation can be more easily achieved when Doppler information is available because the raindrop terminal velocity is an increasing function of the raindrop size (Atlas et al., 1973), which allows the return from different sizes to be separated in the spectral domain. Sophisticated single-and dual-frequency retrievals of the full rain microphysics and the entangled dynamics (turbulence and vertical wind) have been developed over the last two decades from Doppler spectral moments (Williams et al., 2016) and full Doppler spectra (Giangrande et al., 2012;Tridon & Battaglia, 2015;. A specific subset of multifrequency Doppler radar DSD retrievals is represented by those based on differential Doppler velocity (DDV) measurements (Liao et al., 2008;Matrosov, 2017;Tian et al., 2007). Such retrievals are aimed at estimating a characteristic raindrop size parameter (e.g., the mean mass-weighted equivolume diameter, D m ) only. Compared to reflectivity-based methodologies, DDV techniques have the clear advantage of being (1) much simpler in their implementation, (2) immune to radar miscalibration and attenuation, and (3) unaffected by the vertical wind speed. However, careful alignment of the different antennas must be performed for a proper Doppler cross calibration between different frequencies. The idea underpinning DDV-based retrievals is generally based on a relationship between DDV and the characteristic raindrop size. Such relationships have been thoroughly discussed for different pairs of frequencies: Tian et al. (2007) derived results for the X-W band pair based upon the assumption of a three-parameter Γ function DSD (see their Figure 1), whereas Matrosov (2017) exploited disdrometer observations to build a statistical relation for the Ka-W pair (see their Figure 1). In both cases, precision of the retrievals is hindered by the unknown DSD width and double solutions of D m that appears as DDV becomes large (typically when D m exceeds 1.7-1.8 mm). Such ambiguity can be mitigated by using the absolute reflectivity or the Doppler velocity measurement at one frequency. For the DDV Ka−W technique, Matrosov (2017) concluded that D m can be retrieved with a precision of approximately 20% but only for values in the range between 0.5 and 2 mm and no extensive validation with in situ data has yet been performed. Similarly, no assessment of the DDV X−W technique has been established so far. Liao et al. (2008) have shown that retrieval of characteristic sizes from DDV X−W techniques is consistent with the dual-frequency reflectivity ratio retrievals. Therefore, this work has three main objectives: 1. to refine DDV techniques by using an extensive data set of DSDs acquired during the GPM validation field program and exploiting the climatological relationship between D m and width of the DSDs; 2. to assess retrieval performances (precision and accuracy) of D m in DDV X−W and in DDV Ka−W by expanding upon the results in Matrosov (2017); 3. to evaluate the potential improvements of a triple Doppler velocity (TDV) approach in retrieving both D m and in the width of the DSD.
The paper is structured as follows: first, the theory of underpinning DDV retrieval is shortly revisited (section 2). In section 3 the application of the TDV methodology to a case study measurements and a validation of different retrievals is shown. Conclusions and future work are presented in section 4.

Methodology
The DSD fully describes rain microphysics by providing the number of drops with a given equivolume diameter per unit volume of air. The changes of the DSD in time and space give insight into processes that control the evolution of precipitating systems, such as evaporation, condensation, collision-coalescence, and drop breakup. DSDs are measured either via in situ or remote sensing techniques. Most in situ DSD observations are obtained via optical disdrometers, where drops which pass through a sampling cross section are counted over some time period and then converted to volumetric drop concentration indicated as N (Battaglia et al., 2010;Raupach & Berne, 2015). Such measurements provide DSD binned according to the necessity of the 10.1029/2019EA000789 user (generally, the 2-D video disdrometer is sampling at 0.2 mm bin widths resolution), but they represent only a small volume and are restricted to the location of the instrument.
Remote sensing techniques offer an alternative approach to measure DSDs. Vertically pointing radars recording Doppler spectra are capable to provide a detailed description of the DSD at high spatiotemporal resolution (e.g., . Conventional scanning radars have the advantage of covering larger areas, but retrieving information about the DSD requires the use of more concise multiparameter models. The three-parameter (N 0 , , and Λ) Gamma function (Γ) is the most widely used multiparameter model (Testud et al., 2001). Within this framework, N 0 or any normalized concentration parameter accounts for the magnitude of the DSD, whereas two additional parameters describe its shape. Equivalently, for a parametrization of the Γ model, more physically meaningful parameters can be used such as the mass-weighted moments (Williams et al., 2014): where m denotes the mass spectrum, that is, m(D) ∝ D 3 N(D). The equivalency of ( , Λ) or ( m , D m ) modeling is given by the following relations: D m = ( + 4)∕Λ and m = √ + 4∕Λ.

DDV Methods
The mean Doppler velocity, v D , is the line-of-sight velocity of targets relative to the radar, v LOS , within the radar volume weighted by their frequency-dependent backscattering cross sections, b ( , D): Throughout this manuscript, positive velocities correspond to motions approaching the radar. The backscattering cross-section term can be calculated with the T-matrix method (Mishchenko & Travis, 1994), whereas in quiet air and for a vertically pointing radar, v LOS equals the raindrop terminal fall speed as determined by Gunn and Kinzer (1949) and revised further by many others (e.g., Atlas et al., 1973;Best, 1950). The notationṽ D refers to the mean Doppler velocity in the absence of vertical wind. Since the total number of drops in the backscattering volume cancels out in equation (3),ṽ D depends only on the shape of the DSD (i.e., on D m and m for the Γ model). The impact of D m and onṽ D is demonstrated in Figure 1a for various Γ DSDs simulated for Ka band (35 GHz). It can be seen that the variability ofṽ D due to D m is much stronger than that of , especially for small drops (D m < 2 mm). This implies that only under ideal condition (e.g., no vertical air motion) single-frequency Doppler velocity can be used to determine the mean drop size in the radar like done for Micro Rain Radar systems (Peters et al., 2005).
However, Doppler measurements have two main drawbacks: first, the Doppler velocity is affected by the vertical wind; second, any deviation from the vertical pointing direction of the radar antenna results in a contamination by horizontal winds. These factors can cause v D to significantly deviate fromṽ D , with a large impact on the retrieval of D m . For instance, vertical motions in the range of ±0.25 m/s are likely to occur even in stratiform rain due to turbulent motions or gravity waves. The impact of those vertical air motions on D m assuming also a variability of is illustrated by the gray-shaded areas in Figure 1a. Due to the saturation of the D m -ṽ D relationship, retrieval errors, already significant for small D m , strongly increase for D m larger than 2-3 mm.
DDV measurements are not affected by vertical air motion because the contribution of the wind is canceled out when taking the difference between the two Doppler velocities. The connection of DDV to the DSD shape is caused by increasing non-Rayleigh effects of larger rain drops for the higher frequency (e.g., compare the yellow curves in Figure 1). This idea has been exploited to develop a D m retrieval using pairs of frequency for the radar bands X and W (e.g., Liao et al., 2008;Tian et al., 2007) and Ka and W (Matrosov, 2017). For the latter, Figure 1b illustrates the relationship between D m and DDV Ka−W for Γ DSDs with different . Uncertainties in D m are clearly driven by uncertainties in the width of the DSD. Moreover, the decrease of DDV Ka−W for D m > 1.5 mm means that two possible solutions for D m exist, and hence, the retrieval becomes ambiguous regardless of the choice of . In order to address the uncertainties due to variable , Matrosov (2017) derived a statistical relation (plotted as a black line in Figure 1b) between D m and DDV Ka−W based on a disdrometer data set of DSDs. The problem of the double solution was bypassed by limiting the analysis to the DSDs characterized by Ka Band Doppler velocities smaller than 6.9 m/s (D m ≈ 2 mm, assuming the Γ model). For D m in the range between 0.6 and 1.7 mm, such a methodology is expected to perform with an accuracy of 20%. Note that DDV techniques are not immune to the problem of superterminal and subterminal drops in the radar volume. Deviations from the expected falling velocity, due to hydrodynamical changes lagging behind the microphysical processes, introduce unpredictable anomalies on the mean Doppler velocity. The effect of the turbulence can introduce a positive or negative biases depending on the shape of DSD. Furthermore, the magnitude of this deviation, as the backscattering cross-section area, depends on the radar frequency. Due to the complexity of this problem, we do not attempt to compensate for these effects and we treat superterminal and subterminal velocities as a source of random errors in the retrieval. Having said that, the radars used in this study are characterized by very narrow beam widths (see Table 1) which effectively limits aforementioned turbulence effects corresponding to subterminal/superterminal drops.
An alternative approach to effectively retrieve D m from DDV data is to constrain the parameters of the DSD model. This can be achieved either by fixing one parameter (e.g., or m as done in Tian et al. (2007)) or Lowest clutter-free range (m) 300 400 300

Radome No No Yes
Note. Note that the W band radar is a FMCW system for which chirp repetition frequency, number of spectral average, Nyquist velocity, and range resolution are changing for different range intervals (see also Dias Neto et al. (2019)); values provided here are for the lowest range gate region from 220 to 1,480 m.
by imposing a climatological relation between DSD parameters (Williams et al., 2014;Zhang et al., 2001). For instance, the black dashed line in Figure 1b corresponds to the D mm relation of Williams et al. (2014). Unfortunately, this approach is still affected by the double solution problem and it is restricted to DDV values where the inversion procedure is applicable (e.g., for the Ka-W pair DDVs have to be smaller than 1.8 m/s). Moreover, the derived relationship is based on the Γ model assumption (Thurai & Bringi, 2018) which ignores any natural variability of m for a given D m .

TDV Method
Including a third frequency has the potential of improving the retrieval of the characteristic size and the shape of the DSD without the need for a constraint between the model parameters (e.g., and D m ). This is demonstrated in Figure 2 for the X, Ka, and W band frequencies: Γ-DSDs with different D m and different tend to be well separated in the DDV Ka−W -DDV X−Ka space, at least for D m exceeding 1 mm. In theory, measurements of DDV Ka−W and DDV X−Ka could be used to retrieve and D m . Such an approach has the clear limitation of being based on the Γ parameterization; in this work, instead, the retrieval capitalizes on a comprehensive DSD database and on a radar Doppler forward model to generate inversion LUTs (lookup tables) that map TDV measurements into DSD properties (D m and m ). The approach is similar to Matrosov (2017) and based on mapping DDV Ka−W to D m .

Inversion LUTs
In situ measurements of the rain DSD gathered during field campaigns and from permanent sites of the Ground Validation program of the Global Precipitation Measuring Mission (Hou et al., 2014) are exploited here (for more detail see ; Dolan et al., 2018). The analysis is restricted to the data from the two-dimensional video disdrometer (2DVD; Kruger & Krajewski, 2002). Although the integration time of the instrument was originally set to 1 min, five consecutive rainy measurements are combined using a running average to obtain DSD samples that are more representative of large radar volumes. This leads to a data set of almost 190,000 samples of rainy measurements across the globe, thus thoroughly covering natural variability. Any observation contaminated by frozen precipitation is discarded from the analysis. Furthermore, an additional threshold on minimum rainfall rate (R ≥ 0.1 mm/hr) and number of drops (N T ≥ 100) is set, as well as on the mean mass-weighted diameter (D m ≤ 4 mm). These thresholds are set in order to have statistically and physically meaningful DSDs and decreased the data set sample size to approximately 150,000 DSDs for analysis.
For each DSD, Doppler velocities,ṽ D , are simulated at the three frequencies. The backscattering cross sections are modeled with the T-matrix method (Mishchenko & Travis, 1994) assuming an ellipsoidal shape with the aspect ratio from Brandes et al. (2005). The refractive index of water is calculated at the temperature of 280 K following Ellison (2007). The raindrop terminal velocities are calculated according to Gunn The mapping between the pair of observables (DDV Ka−W , DDV X−Ka ) and the pair of DSD parameters (D m , m ) is derived using Bayesian theory, that is, for any hypothetical measurement the weights of the DSD parameters from our database are calculated taking into account the distance to the simulated DDV values and the uncertainty of the measurement: In equation (4) a subscript i denotes ith element of the DSD data set, Σ is a covariance matrix of the measurement error, and DDV m and DDV i are vectors composed of the measured and simulated (DDV Ka−W , DDV X−Ka ) pairs, respectively. The expected value and the uncertainty of the estimate for a given measurement are computed as a weighted mean and a standard deviation of the DSD data, respectively. In this study, the measurement error of 0.1 m/s is assumed for all frequency channels which corresponds to In Figure 3 the black contour lines represent the expected values for D m (a) and m (b), whereas the color of the background corresponds to the relative uncertainty of the estimate. White areas represent the regions of insufficient sampling. The white hatching corresponds to the retrievals where a strong overestimation is expected due to the characteristics of the data set. Sampling deficiency of the 2DVD instrument causes that the training data set is missing DSDs characterized by the D m smaller than 0.3 mm. This reduced variability of D m and m occurs for small DDVs and therefore results in the overestimation of the DSD moments and the underestimation of the uncertainty of the retrieval whenever DDV X−Ka and DDV Ka−W are close to 0 m/s. In this case, the retrieval can be only considered as an upper limit of the real characteristic size. Note that the largest retrievable size is limited to 2.5-3 mm, which is a consequence of two facts: first, DSDs with a characteristic size exceeding 3 mm are very rare, so they are washed out in the statistical retrieval whenever their occurrence domain in the DDV-DDV space overlaps with smaller sizes that occur more often; second, very small variability in the terminal velocity of raindrops exceeding 3 mm (Gunn & Kinzer, 1949) introduces a large ambiguity in the retrieval of the characteristic size even when the full Doppler spectra information is available, all the more if only the mean Doppler velocity is measured. In order to expand the methodology of Matrosov (2017) to the X-W pair of frequencies, a mapping between DDV X−W and D m is derived using the same training data set as for TDV retrieval (plot not shown). Again, Bayesian statistics and the measurement error of 0.1 m/s are used. The best fit curve and the retrieval error estimate are given by the following: where DDV is an abbreviation of DDV X−W for this formula only. By virtue of formulae (5) and (6), the theoretical retrieval error varies between 11% and 22% for DDV X−W from 0 to 4 m/s. Note that this retrieval is also affected by the lack of small characteristic sizes in the training data set; thus, for DDV X−W < 0.2 m/s the provided size estimate should be considered only as the upper limit.

TRIPEx-pol Field Campaign
The "TRIple-frequency and Polarimetric radar Experiment for improving process observation of winter precipitation" ( Although the Parsivel distrometer is utilized as a reference, it must be noted that it is not an error-free instrument and it has its own disadvantages. Biases in the measurements can be a result of drop splashing, partial raindrop recordings, multiple drops appearing at the same time, and contamination with insects or spiderwebs. The second generation of the instrument operates at the JOYCE site, which has been shown to perform better than its predecessor (Tokay et al., 2014). According to the manufacturer reports, this instrument measures rainfall rate with the accuracy of ±5% which has been confirmed by Tokay et al. (2014) where the absolute bias between the disdrometer and the gauge was 6%. In comparison to the other available instruments, Parsivel2 exhibits good agreement with Joss Waldvogel disdrometer (Tokay et al., 2014).
On the other hand, some underestimation of raindrops between 1.38 and 3.25 mm in diameter has been reported by Raupach and Berne (2015) compared to a 2DVD. The same author pointed out that the numbers of drops larger than 3.5 mm is overestimated by both generations of the Parsivel unit. These differences in the shape of the PSDs measured by different instrument result in difference in the integrated distribution moments that has been quantified to be as large as 39% and 47% for third and fourth moments, respectively.
The three radars are installed on the same roof platform within a maximum distance of 10 m in order to maximize radar volume matching (main technical specifications of the radars are provided in Table 1). The absolute pointing accuracy of the scanning Ka band radar has been estimated to be smaller than ±0.1 • in elevation and azimuth using a Sun-tracking method similar to Muth et al. (2012). In addition, "bird bath scans" (i.e., azimuth scans with antenna pointing zenith) have been performed in order to check the proper zenith pointing of the Ka band. Even in situations with strong horizontal winds in high-altitude ice clouds (up to 50 m/s), no signal of a misalignment of the antenna could be detected. The X/W band radar has been compared to the Ka band system for several cases of different horizontal wind velocities and directions similar to the method described in Kneifel et al. (2016). The relative misalignment of the three radars in elevation was found to be in the range of 0.1 • .
The beam width of the W and Ka band instruments is matched reasonably well, but the lowest frequency radar has its sampling volume twice as large as the others. On the other hand, for the sake of higher sensitivity the W band radar samples the atmosphere at the lower rate with longer integration time than the other instruments. Both of these factors introduce disparity in measurements at different frequencies. In order to minimize these effects, the data are averaged over 1 min for all the radars. This procedure removes differences in the integration time and also minimizes differences in the sampling volumes and additionally reduces random fluctuations due to random noise. The overall uncertainty of the measurements is quantified by the standard deviation of DDVs for Rayleigh particles. Threshold values of −10 dBZ and 2 m/s on the Ka band reflectivity and the velocity are used for downselecting the matching domain. Furthermore, only the data with a high signal-to-noise ratio are considered; that is, all the reflectivity channels are required to be at least 9 dB above the sensitivity threshold. The standard deviation of DDVs is 0.25 and 0.1 m/s for Ka-W and X-Ka pairs, respectively. Similar standard deviations have been observed for rain. Moreover, for all days included in this analysis, the mean DDV Ka−W and DDV X−Ka for Rayleigh targets within ice clouds are −0.067 and −0.004 m/s respectively, which confirms the very tight alignment of the antennas. Figure 4 presents time-height plots of the radar measurements for a 4 hr rain event that occurred on 8 December 2018. The top panel depicts the radar reflectivity factor at X band. The melting layer can be well identified at approximately 1.5 km in the X band reflectivity. This is consistent with the Doppler measurements which clearly show an acceleration due to melting at the same height. Between 18:00 and 22:00 UTC, three periods of heavy rainfall are noticeable in the reflectively data: two showers around 18:30 and 21:00 UTC and one longer-lasting event from 19:45 until 20:30 UTC. This is also reflected in the Doppler measurements where the X band velocity reaches values close to saturation, which suggests characteristic sizes exceeding 1.5 mm. As expected from Figure 1, strongly enhanced DDV is found for the rainfall events with large Doppler velocities. Surprisingly, the DDV X−Ka are found during moderate rainfall (Z X < 35 dBZ, e.g., 20:10-20:30) to be slightly negative. Rather than being a measurement artifact, this effect can be well explained by a super-Rayleigh scattering property of raindrops at the Ka band: for D m < 1 mm the backscattering cross sections exceed the Rayleigh values before dropping to their first minimum, which causes negative DDV Ka−X (see Figure 2).

Case Study
The retrieval of D m and m for the DDV X−W and TDV methods is presented in Figure 5. All Doppler measurements are corrected for changes in the density of air prior to the retrieval. The gray area in panels (a)-(c) indicates the region where the retrieval is invalid; that is, for panels (a) and (c) this corresponds to the hatched area in Figure 3; for panel (b) it is where DDV X−W < 0.2 m/s. Both retrieval methods provide a very similar pattern of characteristic diameters, but there are some subtle differences between them. For example, for two intensification periods at 18:30 and 19:45-20:00 the triple-frequency technique provides a clearly larger characteristic size estimate than the DDV X−W technique. On the other hand, an opposite is true for the shower after 21:00. Although it is expected that the TDV retrieval provides a more reliable representation of the observed processes (because it is based on more measurements), an independent validation is necessary.
The validation is performed by comparing the retrieval for the lowest available range gate with the measurement of the Parsivel optical disdrometer which operated in close proximity to the radars. In order to minimize the difference in sampling volume, 3 min averaging is applied. For this comparison, the DDV Ka−W retrieval of Matrosov (2017) is also included. All considered algorithms provide almost exactly the same D m estimate whenever D m is smaller than 1.2 mm. This is not surprising because for such small particles the DDV X−Ka oscillates around 0 m/s, and it is mainly the DDV Ka−W (≈ DDV X−Ka ) that drives the retrieval. For larger sizes, differences between algorithms become more noticeable. For the whole duration of the event, the triple-frequency retrieval outperforms the other algorithms, which is particularly evident for the periods with most intense rainfall and associated large D m . The DDV X−W retrieval performs slightly worse for the largest D m . This can be caused by the previously mentioned problem of a double solution in the DDV − D m space. Because the same DDV is observed for a range of diameters, the retrieval provides the most probable estimate corresponding to smaller sizes that are more frequent in the training data set. The same issue affects the DDV Ka−W retrieval, but in this case, the saturation of the signal takes place for even lower sizes and thus the underestimation of D m is more severe.
In situ measurements and the retrieved values of the m are presented in Figure 5d as dashed lines. The estimate based on DDV X−W assumes a fixed D mm relation of Williams et al. (2014); that is, the retrieved width of the DSD is directly computed from the characteristic size. As it is for D m , the retrievals of m differ only for the periods with broadest DSDs, where DDV X−Ka provides additional information to the algorithm. During these time slots, the DDV − DDV retrieval follows the measured m values more closely, though a slight underestimation of m is still visible.

Validation With Disdrometer
The proximity of the Parsivel to the three radars at JOYCE-CF provides a very valuable data set to validate the Doppler-based retrieval of D m and m for several rainfall cases. For this study, Parsivel observations from 13 rainy days in late autumn, early winter, and late spring are analyzed providing ≈ 55 hr of precipitation in total. A complete list of the dates used for validation is given in Table 2.
The Doppler data are corrected for changes in the density of air prior to the retrieval. For this, a pressure adjustment factor of ( ∕ 0 ) 0.4 is used as suggested by Foote and Du Toit (1969), where 0 (=1.2434 kg/m 3 ) is the air density used for the drop velocity calculations and is the density at the level of observation. In order to match the temporal resolution of the instruments, the radar measurements are averaged over 1 min which corresponds to the sampling interval of the disdrometer. Although the radar data acquired at the lowest range gate are used for comparison, the distance between the in situ and remote measurements is approximately 300 m. In order to match them more precisely, the X band reflectivities are simulated based on the disdrometer PSD. Then, for each 15 min time window, the optimal time lag that maximizes the correlation between the simulated and the measured reflectivity values is derived. The retrieval results and the PSD derived estimates of the D m and m are compared using this optimal matching. Changes in the DSD along the 300 m path from radar volume to disdrometer due to evaporation, condensation, and raindrop interaction are difficult to estimate and are therefore treated as a random error. Because the algorithms show different behaviors for different characteristic sizes, all retrieved diameters are clustered in logarithmically spaced bins according to their disdrometer counterpart and then the performance is assessed in each diameter class. One of the metrics used for the assessment of the retrieval performance is the mean normalized bias (MNB) defined as follows: where D r m is the mass-weighted diameters in mm retrieved from the radar data, D d m denotes the disdrometer class, and the angle bracket is an average over the class. The uncertainty of the retrieval is quantified in terms of the normalized root-mean-square error (NRMSE): which is a measure of the scatter in the data around the "ground truth." Figure 6 presents the MNB and NRMSE for the two DDVs and our new TDV retrieval over several D m classes.
In general, all the algorithms perform nearly identical for D m < 1.2 mm. This similarity of all three retrievals was already observed for the case study. For small D m , all rain drops are Rayleigh scatterers at X and Ka bands (DDV X−Ka = 0). Consequently, the main information on the DSD comes from the W band transitioning into non-Rayleigh scattering. The NRMSE of all retrievals is largest for the smallest D m where all retrievals overestimate the D m measured by the Parsivel. It should be noted that the Parsivel sorts the drops into predefined 32 size bins (Angulo-Martínez et al., 2018, Table 1). Especially for very narrow DSD, the fixed Parsivel size bins might introduce biases and should be interpreted with care. On the other hand, theoretically derived uncertainties of the D m and m in Figures 3a and 3b predict an enhanced retrieval ambiguity for D m smaller than 0.8 mm appearing as a noticeable peak around the origin of the DDV-DDV space. As it was discussed in section 2.2.1, the predicted uncertainty is clearly underestimated compared to the value derived during 10.1029/2019EA000789 Figure 7. The same as Figure 6 but for m . the validation. The difference is a consequence of a limited number of DSDs characterized by a small mean size in the training data set due to sampling limitations of disdrometers. This deficiency of the training data set leads to an underestimate of the theoretical variability of D m and m for DDVs close to 0.
For D m between 0.8 and 1.6 mm, the MNB of all retrievals does not exceed 20% in the absolute value which is close to the predicted retrieval uncertainty. For D m exceeding 1.4 mm, differences in the three retrievals become noticeable. The DDV Ka−W retrieval performs worse than the others; both the precision and the accuracy are decreasing with increasing D m reaching a maximum underestimation of −50% and a NRMSE of 50% for the largest considered size of 2.5 mm. DSDs characterized by the mean mass-weighted diameter exceeding the threshold value of 2.5 mm are not considered in this study. It should be noted that Matrosov (2017) suggests to filter out Doppler velocities at Ka band which are larger than 6.9 m/s. For this comparison, it was decided not to apply this filter in order to quantify potential biases in the DDV Ka−W approach when large drops are observed. From Figure 6 it is clear that the DDV Ka−W retrieval becomes increasingly inaccurate and imprecise when moving to D m above 1.7 mm, so that precision of 20% for characteristic sizes reaching 2 mm cannot be achieved as previously suggested in Matrosov (2017). The retrieval based on DDV X−W measurements performs well for the D m between 1 and 1.9 mm, where the MNB changes only from 10% to −20%, whereas the uncertainty of the retrieval does not exceed 25% (Figure 6b). For larger sizes, as it is for the DDV Ka−W retrieval, there is a progressive reduction in the reliability of the algorithm reaching the NRMSE of 45% for D m = 2.5 mm. The TDV algorithm performs better than DDV algorithms with a reduced bias for large drops and a lower NRMSE for the entire D m range. The normalized bias does not exceed 20% up to 2.25 mm, whereas the NRMSE is lower than 30% for the whole range of considered sizes, with a majority of the domain being below 20% uncertainty threshold. However, the disdrometer data should be interpreted with care. A comparison of the Parsivel derived and the actual X band reflectivity reveals a systematic overestimate of approximately 2.5 dB of the disdrometer simulated radar measurements for Z X reaching 40 dBZ (for stronger echos the overestimate is even larger). This discrepancy suggests that the ground "truth" might be subject to a positive bias that can be caused by the shadowing effect, that is, when two drops appear as one that is much larger.
The MNB and NRMSE of the DSD width are also assessed in Figure 7. The DSD width parameter m is retrieved directly in the TDV approach; for the dual-frequency retrievals the climatological relation from Williams et al. (2014) is applied to the retrieved D m values. Note that using the "retrieval" term for DDV techniques is a bit of a stretch because this methodology is based on the assumption that there is one-to-one correspondence between D m and m . Nevertheless, the assessment of the performance of TDV and DDV techniques presented below is essential to quantify the added value of the third frequency channel on top of the information given by two radars and the DSD climatology. The performance of the retrievals for m is similar to that for D m ; there is large positive bias up to 20% for low m that continuously decreases reaching up to −65% for the worst performing technique (DDV W−Ka ). The uncertainty of the DDV Ka−W retrieval varies between 25% and 55%, and it generally increases with m . The TDV technique is superior to the DDV retrievals for m > 0.5 mm. The ability of the TDV algorithm to account for the variability in the triple-frequency DDV space rather than using a single DDV results in approximately 20% (10%) reduction of the MNB and NRMSE for m larger than 0.8 mm compared to DDV Ka−W (DDV X−W ). The inability of the DDV Ka−W algorithm to derive large widths of the DSD is a consequence of the underestimation of large characteristic sizes with the Matrosov (2017) approach; the underestimation of D m is then translated to a negative bias in the width of the DSD via the m − D m formula. The MNB for the TDV retrieval of m changes nearly linearly from a positive bias of 20% for small widths to approximately 50% underestimation for m = 1.3 mm. The uncertainty of the estimate varies between 25% and 50% with a clear upward trend as the width increases. The mean normalized bias is found to be within the theoretically derived uncertainty values for m below 0.7 mm (see supporting information Figure S2); for larger widths the error increases linearly and is approximately 2.5 times larger than expected when m reaches 1.3 mm. This is not very surprising considering, for example, the imperfect collocation of in situ and radar observation and the aforementioned issue of the shadowing effect in the disdrometer data. Figure 8 presents a histogram of characteristic sizes and DSD widths as measured by the disdrometer at the ground for the whole validation data set. For reference, the climatological relation between D m and m derived by Williams et al. (2014) is overlaid as a magenta line. Clearly, data collected in Cologne differ from those observed in Huntsville, Alabama, where the study of Williams et al. (2014) was performed. On average, for the same characteristic size larger DSD width is observed at the German site. This discrepancy may be a result of two factors: first, different DSD climatology and second, differences between the measurement instruments. The state of Alabama is characterized by a humid subtropical climate, Huntsville experience, ∼1,390 mm of rainfall over a year; the mean annual temperature is 16.6 • C. Cologne has a temperate oceanic climate with a much lower average annual temperature of 10.1 • C. In a year, the average rainfall is 774 mm, 10.1029/2019EA000789 nearly a half of Huntsville amount. Therefore, differences in DSD statistics would not be a surprise. In his study, Williams et al. (2014) used 2DVD for PSD measurements, whereas the Parsivel disdrometer operates at JOYCE site. Differences between these instruments have been described and quantified by, for example, Raupach and Berne (2015) and Tokay et al. (2014) where a potential overestimation of a number of large drops in Parsivel measurements has been pointed out. This could, at least partially, explain differences in the distributions of D m and m for different locations. The TDV retrieval results (black dots in Figure 8) follow much better the Parsivel data than the climatological formula of Williams et al. (2014), especially for PSDs where D m is larger than 1.1 mm. This improvement indicates an added value of the third radar frequency for retrieving the DSD width when non-Rayleigh scattering is experienced at two radar bands. Nevertheless, a clear underestimation of extreme D m and m is apparent.

Conclusions and Discussion
Collocated measurements of the difference in the mean Doppler velocities at two or more frequency bands can be effectively used to estimate mean mass-weighted diameter (D m ) and mass-weighted standard deviation ( m ) of the DSD in the radar volume. The methodology relies on the presence of non-Rayleigh scattering raindrops which cause the DDV to deviate from zero. This implies that these methods are performing best for DSDs whose characteristic size exceeds a certain threshold. This threshold depends on the frequency combination adopted for the retrieval. In this study, combinations of X, Ka, and W band radar observations are investigated. For these frequency combinations, retrievals are reliable for D m exceeding 0.7-0.8 mm. For smaller sizes, an implicit assumption on the vertical winds must be made to perform the inversion directly from the Doppler velocity measurements. On the other hand, the methodology described in this paper is limited to DSDs whose characteristic size is smaller than 2.5 mm, which results from the saturation of the raindrop terminal velocity exceeding 2.5-3 mm.
For the first time, a combination of three frequencies was included in the retrieval framework. This TDV retrieval of D m outperforms the dual-frequency DDV approaches for a wide range of considered D m . The TDV algorithm shows the greatest advantage compared to the DDV Ka−W methodology, which is affected by large uncertainty already for D m exceeding 1.7 mm. For these characteristic sizes, the combination of X and Ka bands with W band reduces the uncertainty of the retrieval by 20-30 percentage points, compared to DDV Ka−W technique. The dual-frequency algorithm based on DDV X−W is comparable to the TDV retrieval for D m < 2 mm; for larger sizes the TDV technique has a clear advantage of a reduced bias by 10 percentage points and smaller uncertainty of the retrieval (up to 15 percentage points gain). Worse performance of the dual-frequency approaches results from ambiguous inversion region for large D m (see Figure 1) which can be circumvented using a combination of three frequencies.
An additional advantage of the TDV technique is its ability to derive also the width m in addition to D m . For the dual-frequency retrievals, this parameter can only be derived using climatological relations between D m and m . With the TDV method, m can be derived with smaller bias and uncertainty as compared to the dual-frequency methods.
The principles of the TDV technique presented in this study can be also used for other combinations of the radar bands. In particular, if the retrieval of drizzle-like precipitation is of interest, radar frequencies higher than 94 GHz must be used, as suggested in Battaglia et al. (2014). Inclusion of the G band (1.5 mm wavelength) data should extend the range of plausible retrievals to about 0.4 mm, because non-Rayleigh effects at this frequency occur at much smaller sizes than at W band. Moreover, simulations with Gamma function DSDs indicate that the DDV X−G is very sensitive to the shape of the DSD; that is, for D m larger than 1 mm, DDV is strongly dependent on but it is less variable in terms of D m (see Figure S3 in the supporting information). However, the possibility of using this characteristic may be significantly hampered by strong attenuation at 200 GHz channel. An alternative modification of our technique would be a replacement of the X by the Ku band and of the Ka by the K band (which is adopted by the relatively low cost Micro Rain Radars). These new sets of frequency pairs should offer performances similar to the ones discussed in this work. Because the majority of raindrops scatter according to the Rayleigh approximation for frequencies smaller or equal to X band, wind profilers and vertically pointing S band radars, used, for example, at Atmospheric Radiation Measurement sites, can effectively substitute X band instrument. Nevertheless, it must be noted that formula (5) can introduce a positive bias when applied to DDV S−W or DDV C−W measurements because the super-Rayleigh behavior of large raindrops can cause DDV X−W to larger than DDV Ra leigh−W by up to 0.4 10.1029/2019EA000789 m/s according to our simulations. Therefore, corresponding formulas for DDV L−W , DDV S−W , and DDV C−W are presented in the supporting information.
Retrievals using the information of the full Doppler spectrum (e.g., 2017) remain superior to dual-wavelength and triple-wavelength DDV retrievals in terms of uncertainties. However, they are complex and computationally expensive and they rely on high-quality Doppler spectra which are not always available. With the clear advantage of being unaffected by the attenuation, radome attenuation uncertainties, and receiver saturation, TDV retrievals can be synergistic with reflectivity-based multifrequency radar retrievals. The integration of both Doppler and reflectivity-based rain microphysics retrieval in a unified scheme is the topic of our future research. In addition, an assessment of the potential of the full Doppler spectra retrieval incorporating X and Ka band measurements in extreme rainfall conditions is planned for future studies.