the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Drop size distribution retrieval using dualpolarization radar at Cband and Sband
Daniel Durbin
Yadong Wang
PaoLiang Chang
Having knowledge of the drop size distribution (DSD) is of particular interest to researchers as it is widely applied to quantitative precipitation estimation (QPE) methods. Polarimetric radar measurements have previously been utilized to derive DSD curve characteristics frequently modeled as a gamma distribution. Likewise, approaches using dualfrequency measurements have shown positive results. Both cases have relied on the need to constrain the relationship between the DSD parameters based on location or assumed weather conditions. This paper presents a methodology for retrieving the DSD parameters using the dualfrequency and polarimetric nature of measurements from a unique data set taken at colocated Sband and Cband dualpolarization radars. Using the reflectivity and differentialphase measurements from each radar, an optimization routine employing particle swarm optimization (PSO) and Tmatrix computation of radar parameters is able to accurately retrieve the gamma distribution parameters without the constraints required in previous methods. Retrieved results are compared to known truth data collected using a network of OTT Parsivel disdrometers in Taiwan in order to assess the success of this procedure.
 Article
(5388 KB)  Fulltext XML
 BibTeX
 EndNote
Knowledge of the drop size distribution (DSD) for a given location is an incredibly valuable piece of information. The most precise rainfall estimates can be achieved with an accurate assessment of the DSD. The applicability of deducing localized DSDs to quantitative precipitation estimation (QPE) is evident. QPE provides society with many benefits. Being able to precisely identify or even forecast heavy precipitation can help emergency and water management services better deploy their resources and can alert others of dangerous flash floods, which can ultimately save lives.
The DSD provides critical information about the composition of a given volume of atmosphere in a region of interest. It is expressed as the number of particles for each drop size given as an equivalent sphere with a reference diameter, typically measured in millimeters. By knowing the DSD, important measurements such as reflectivity (Z), rainfall rate, and total water content can be derived. These parameters are only as accurate as the DSD representation, which highlights the value of flexibility in the model.
The Marshall–Palmer distribution is a wellknown model for describing the DSD of rain (Marshall and Palmer, 1948). It is based on an exponential distribution and assumes that the DSD is spatially and temporally homogeneous. The intercept parameter of the Marshall–Palmer distribution, which is set to 8000, determines the overall scale of the distribution.
However, research has shown that the Marshall–Palmer distribution is limited in its applicability to different types of precipitation and atmospheric conditions. Specifically, it has been found that the Marshall–Palmer distribution is only accurate in the stratiform precipitative region, which is a region of moderate, steady rainfall. As the rainfall rate increases and diverges from this mode of precipitation, the intercept factor exhibits a large degree of variability (Sauvageot and Lacaux, 1995). While it can be a useful approximation, it has limitations that should be taken into account when using it in applications such as QPE.
A larger number of DSDs can be expressed generically in the form of a gamma function given by Eq. (1), where N(D) (m^{−3} mm^{−1}) is the concentration for each diameter D (mm) and μ (unitless), Λ (mm^{−1}), and N_{0} are the shaping parameter, the exponential slope parameter, and the intercept parameter of the distribution, respectively (Ulbrich, 1983).
While the gamma DSD provides more flexibility than the exponential or Marshall–Palmer distribution, its complexity can make it difficult to use for retrieval purposes, as three unknown parameters must be solved, and the model is still limited in its ability to account for the smallest drop sizes. In the past, researchers have attempted to simplify the process by assuming a relationship between these parameters.
Extensive research has been conducted to estimate DSDs, with many studies utilizing measurements taken at two frequencies. This is because using data from a single frequency is typically insufficient for accurately estimating the DSD and can only provide information about a single parameter at a specific location (Meneghini et al., 1997). Gorgucci and Baldini (2016) demonstrated this approach using global precipitation measurement (GPM) dualwavelength radar measurements. The surface reference technique is used to calculate pathintegrated attenuation for each frequency so that a set of reflectivity and attenuation integral equations is formed. The system can then be solved for the DSD parameters that best match the measurements. The disadvantage of the approach is that it relies on assuming μ is a function of the median volume diameter (Gorgucci and Baldini, 2016).
Williams et al. (2014) also worked to show dualwavelength approaches could be performed with GPM data. The underconstrainment issue was overcome by creating a distribution of standard deviations of the mean drop diameter using a large number of surface disdrometer measurements. Superior results were achieved relative to approaches that assume μ is a constant in order to overcome the illconstrained nature of the retrieval (Williams et al., 2014). Many others have shown dualfrequency measurements are exploitable for discerning these properties (e.g., Mardiana et al., 2004; Chandrasekar et al., 2005; Eccles, 1979; Kozu et al., 1991; Kummerow et al., 1989; Marzoug and Amayenc, 1994; Meneghini et al., 1992).
In addition to leveraging the scattering variances across different frequencies, there have also been efforts to retrieve DSD parameters using dualpolarization measurements. These efforts capitalize on the characteristic oblateness of raindrops, which would obviously not be available in the case of the previously mentioned vertical profiles measured by GPM. Brandes et al. (2002) used Sband measurements of reflectivity and differential reflectivity to seek the gamma parameters. An empirical relation is assumed between μ and Λ. While this relationship was applicable to both convective and stratiform rains, it was localized to the Florida region of the United States that was used for the study. Similar approaches using horizontal and differential reflectivity have also been successful in experiments in central Oklahoma but still rely on constrained relationships of the parameters (Cao et al., 2010).
The aim of this present work is to solve for the three variables of the gamma drop size DSD. To achieve this, measurements taken at two frequencies with both vertical and horizontal polarizations are incorporated into optimization routines. This approach allows for the determination of DSD solutions that accurately represent the input reflectivity and specific differentialphase information.
This paper is organized in the following manner. Section 2 conveys the methodology of this work. Highlevel design characteristics of the instruments used for data collection are provided, including the disdrometers and the two ground radars. Key elements of the data preprocessing routines are provided. Prior to concluding thoughts in Sect. 4, Sect. 3 demonstrates the results through multiple examples, comparing the disdrometercollected DSDs with the retrieved results from the proposed algorithm and an evaluation of the error when applied to QPE. A notable increase in accuracy is found when estimating the rainfall rate using the proposed algorithm when compared to the traditional Z–Rderived rainfall rate. Furthermore, we discuss and demonstrate a supplementary validation approach, utilizing Cband reflectivity data previously excluded from the proposed optimization routines. This approach serves as an external validation mechanism for the proposed algorithm and provides a comparison of the results obtained from the algorithm with those obtained from the Cband reflectivity data. The results demonstrate the effectiveness of the proposed algorithm in retrieving accurate DSDs and its potential usefulness in radar systems that utilize more than one frequency.
2.1 Instrumentation
In the current work, the measurements of two colocated polarimetric radars, RCWF (Sband) and RCMD (Cband), are used in the algorithm development and validation. RCWF is an Sband dualpolarization radar that is part of Taiwan's operational MultiRadar MultiSensor QPE system (Chang et al., 2021). To achieve a QPE accuracy within 10 %, the reflectivity bias must be kept within 1 dBZ. This is accomplished by regularly calibrating the RCWF's reflectivity using a selfconsistency algorithm (Le Loh et al., 2022). In contrast, RCMD is a Cband radar used experimentally and for research purposes rather than operationally.
A photo of these two radars is shown in Fig. 1, where RCWF and RCMD are located on the right and left, respectively. The Central Weather Bureau of Taiwan operates these two radars that provide realtime observations used for severe weather surveillance, quantitative precipitation estimation, hydrometeor classification, and precipitation microphysics studies. These radars play a critical role in monitoring and forecasting precipitation in Taiwan and provide key data that are useful in the research of precipitative processes. The manufacturer models of RCWF and RCMD are the Weather Surveillance Radar 1988 Doppler (WSR88D) and Gematronik, respectively. More technical specifications about these two radars can be found in Table 1. As shown in Fig. 1, these two radars are adjacent to each other, and the distance between them is much less than the range resolution (250 m). Since these two radars have the same polarization and location, the divergence of the observed differential phase (ϕ_{DP}) and reflectivity (Z) can therefore be attributable to only the radar frequency difference. These measurement differences are the key variables utilized in the retrieval method.
To validate the performance of the proposed algorithm, DSD data collected by a laserbased optical system, the OTT Particle Size Velocity (PARSIVEL or Parsivel) disdrometer, are used as the ground truth (Table 3). The Parsivel disdrometer derives the DSD through the measured particle's size and velocity. There are 32 diameter and velocity bins available to measure particles with sizes between 0.062 and 24.5 mm and with velocities between 0.05 and 20.8 m s^{−1}. The accuracy of the Parsivel disdrometer has been studied by Sheppard (2007) and Jaffrain and Berne (2011).
Sampling with the OTT Parsivel disdrometer may exhibit some degree of statistical variance, which is a common characteristic of all measurements. Microphysical events are inherently stochastic in nature, and the physical sampling effects and noise also play roles in the deviations (Sheppard, 2007). The Parsivel data used in this data set have been compared to more accurate 2D disdrometer data in order to gauge the variance. The findings of Jaffrain and Berne (2011) demonstrated that sampling uncertainty is minimal for small to moderate drop sizes but starts to escalate for larger classes (greater than 2.0 mm). Notably, the data set used in this study lacks a large representation of drop sizes exceeding 2 mm, which is advantageous as the smaller class measurements can be assumed to be more accurate (Jaffrain and Berne, 2011). Possible error sources of OTT Parsivel were investigated by AnguloMartínez et al. (2018). It was found that uneven power distribution over the beamwidth or any time variation can adversely affect the accuracy. Other factors, such as the angle of the drop trajectory, coincidentally observed particles, and particles that intersect with only the edge will also lead to biases.
The locations of the two dualpolarization radars and OTT Parsivel disdrometers are depicted in Fig. 2 with black and red markers, respectively. The complex terrain of Taiwan presents a challenge due to radar beam blockage. The Central Mountain Range (CMR) is visible in Fig. 2, running from the north to the south of the island, with the highest peak exceeding 3800 m. To mitigate the effects of vertical variability in radar measurements, only data from the two lowest unblocked tilts are used in the DSD retrieval. Additionally, two rings with ranges of 25 and 70 km are shown in Fig. 2. The performance of the proposed approach was validated using only the disdrometers located within these two ranges. More details related to the validation data selection are provided in the following section.
The final sources of error we attempted to mitigate relate to the radar measurement itself. Observation bias is a primary concern, stemming from the inherent differences between radar and disdrometer measurements – where radars capture data over a large volume and disdrometers provide pointspecific observations, leading to potential discrepancies in data interpretation. Additionally, vertical variations in DSDs pose a significant challenge due to the radar's observation volume being several hundred meters above the ground. This discrepancy is particularly problematic for all radarbased DSD retrieval methods. The various equipment heights and locations are shown in Table 2. To address this issue, we constrained our analysis to data from the two lowestelevation tilts of the radar scans. Although it does not fully eliminate the issue, this approach helps to reduce the impact of this error source. These aspects underscore the complexities involved in accurately retrieving DSDs and highlight the necessity of carefully considering these factors when analyzing radar and disdrometer data.
2.2 Radar and disdrometer data
2.2.1 Preprocessing of radar and disdrometer data
The Z and ϕ_{DP} fields from both Cband and Sband radars are the proposed parameters for DSD retrieval, and the qualities of both fields are examined and processed through a set of quality control procedures. The quality control process pertaining to the reflectivity field includes identifying and removing nonprecipitation radar echoes and smoothing along the radial direction. Any gates associated with correlation coefficient (ρ_{HV}) less than 0.98 are considered to be possibly produced by nonprecipitative clutter, and reflectivity is excluded from the subsequent average operation. The obtained reflectivity is smoothed with a 4 km smoothing window along each radial direction. The raw ϕ_{DP} field of RCMD is processed with the new Selex–Gematronik family of digital receiver and signal processor (GDRX) (Bringi et al., 2005). The GDRX processes the raw field using the field unwrapping, “good data” mask application, and finite impulse response (FIR) filtering. The details of the procedure can be found in Bringi et al. (2005). A similar procedure is also applied on the raw ϕ_{DP} field from RCWF. Examples of processed reflectivity and differentialphase fields are shown in Fig. 3, where panels Fig. 3a (c) and Fig. 3b (d) show the reflectivity (raw differential phase) measured by Sband and Cband, respectively. The yellow arrow in Fig. 3a indicates the radial path from the radars to the disdrometer location of interest. It should be noted that the reflectivity fields from both frequencies are consistent, and the difference is mainly caused by the attenuation. The differentialphase fields, on the other hand, show significant difference, which indicates that the differential phase is more sensitive to the radar frequency. A detailed analysis of the impact from frequency on the differential phase and specific differential phase is presented in Sect. 3.1.
Figure 3 serves as a visual representation of the data used in the DSD retrieval algorithm and demonstrates the differences in the data between the Sband and Cband radars. Figure 4 shows the preprocessed (dashed) and postprocessed Z (solid), without any attenuation correction, and ϕ_{DP} fields along the yellow arrow of Fig. 3, contrasting the difference between the postprocessed differentialphase profiles (solid) and the profiles prior to unfolding and smoothing (dashed). The measurements from Sband and Cband are shown in blue and red, respectively.
The DSD parameters are derived from the Parsivel disdrometer observations through the approach proposed by Raupach and Berne (Raupach and Berne, 2015). In this approach, data concerning individual raindrops, including their diameters and fall velocities, and the effective sampling areas of the instrument are recorded. The drops are binned into diameter classes, and the concentration is then calculated using Eq. (2), where N_{i} is the drop concentration for the ith equivolume diameter class, S is the effective sampling area, V is the particle velocity, and ΔD_{i} and Δt are the class width and sampling period, respectively.
2.2.2 Data selection
The algorithm development and validation were based on 10 d of data collected by RCMD and RCWF. These days included 1 June and 11–17 June 2017 and 7 January and 7 May 2018. The majority of these days correspond to the Meiyu season in Taiwan and represent events of light to moderate precipitation which were observed in the region covered by the two radars. These days provided data representative of multiple precipitation intensities.
One challenge in using both RCWF and RCMD is that they operate under different volume coverage patterns (VCPs). As a result, there can be a slight time lag between scans from these two radars when observing the same location. To minimize the retrieval biases caused by the DSD variation during the time lag, the time stamp differences between scans from the two radars were limited to within 1 min. This time limitation helped ensure that the radar data adhere to an acceptable degree of synchronization.
In this work, the differences in both reflectivity and differentialphase fields obtained from the Sband and Cband radars play critical roles in the DSD retrieval. Sufficient differences are expected from the two frequencies for differential phase. If the differences are too small, the algorithm may result in a biased or inaccurate retrieval of the DSD. Thus, it is crucial to carefully select the observation range such that enough differentialphase difference has been experienced by the radar return.
Another important factor to consider is that biases in the retrieved DSD can accumulate along the range. This means that the farther the distance is between the radar and the target area, the larger the potential error in the retrieved DSD will be. The presence of atmospheric attenuation is also a predictable issue in radar data processing, and it also underscores the importance of carefully selecting the range of interest when estimating the DSD.
To achieve reasonable results, the following criteria for candidate data are therefore used for the terminal gate of the retrieval that contains the disdrometer:

25 km < range < 70 km,

Z^{S} > 25 dBZ.
This set of criteria is designed to strike a balance between creating sufficient deviation between C and Sband differential phases for an accurate DSD retrieval, while also preventing excessive error accumulation and minimizing the vertical separation between the radar observation volume and ground location. Additionally, a reflectivity threshold of 25 dBZ is imposed to ensure that there is enough observable precipitation in the terminal gate where the disdrometer is located. These criteria are primarily chosen to increase the quality of the data for development and validation; however, similar standards would need imposed if the algorithm were to be operationally applied to address the error accumulation and elevation differences.
2.3 Drop size distribution with an artificial intelligence method
The flowchart of the proposed DSD retrieval algorithm is presented in Fig. 5. Three variables, Sband reflectivity (Z^{S}(r,θ)) and differential phase from both Sband and Cband (${\mathit{\varphi}}_{\mathrm{DP}}^{\mathrm{S}}(r,\mathit{\theta})$, ${\mathit{\varphi}}_{\mathrm{DP}}^{\mathrm{S}}(r,\mathit{\theta})$), are implemented as inputs, where r and θ are the coordinate of range and azimuthal angle, respectively. The analysis exclusively employs Sband reflectivity as it is considered a more dependable variable than Z^{C} in the data set since RCWF is utilized operationally and is well calibrated. Although Sband reflectivity does experience atmospheric attenuation, Z^{C} is much more vulnerable to this effect and is therefore not used as an input to the algorithm. Since Z^{C} displays high correlation with Z^{S}, it is unlikely that the inclusion of Cband reflectivity would furnish any additional useful insights into the retrieval. Moreover, the exclusion of Z^{C} from the process serves as an additional validation parameter, as discussed further in Sect. 3.
The radar variables from a given gate are first preprocessed with the routine described in Sect. 2.2.1, and the processed data are then used to retrieve three parameters as described in Sect. 2.3.2.
2.3.1 Tmatrix computation of radar parameters
In the retrieval procedure, a set of DSD parameters (Eq. 1) is first initialized within commonly observed ranges of the following parameters (Zhang et al., 2001):

${\mathrm{10}}^{\mathrm{2}}<{N}_{\mathrm{0}}<{\mathrm{10}}^{\mathrm{10}}$,

$\mathrm{0}<\mathit{\mu}<\mathrm{10}$,

$\mathrm{0}<\mathrm{\Lambda}<\mathrm{15}$.
With the initial parameters, the DSD (N(D)) is calculated, and radar variables (Z^{S}, ${K}_{\mathrm{DP}}^{\mathrm{S}}$, ${K}_{\mathrm{DP}}^{\mathrm{C}}$, and specific attenuation, A) are then calculated with the following equations integrated from 0 to a maximum diameter (D_{max}) of 8 mm (Ryzhkov et al., 2013):
In these equations, vertical polarization scattering amplitudes, ${f}_{a}^{\mathrm{0}}\left(D\right)$ (${f}_{a}^{\mathit{\pi}}\left(D\right)$), and horizontal polarization scattering amplitudes, ${f}_{b}^{\mathrm{0}}\left(D\right)$ (${f}_{b}^{\mathit{\pi}}\left(D\right)$), are calculated with the Tmatrix method (Waterman, 1965), where 0 and π indicate forward and backward directions, respectively. These formulations make use of a simplification with zero canting angle, which will introduce only a small error due to the limitation of elevation tilts. The dielectric constant of water is referenced at 10 °C, which is a reasonable temperature for the radar volume situated between the melting layer and the groundlevel temperature of Taiwan during these dates.
With the obtained K_{DP} field, the ϕ_{DP} fields from both the Sband and Cband radars were then calculated through Eq. (6), where ΔR represents the range difference between the ith gate and the previous gate. ${\mathit{\varphi}}_{\mathrm{DP}}^{\mathrm{sys}}$ is the system ϕ_{DP} found in Table 1 for both radars.
2.3.2 Optimization routine
The estimation of radar variables can be achieved by adjusting a parameterized DSD in order to reduce the difference between the estimated and observed values. As this difference decreases, the optimization problem gradually approaches a minimum value, ultimately resulting in the retrieval of the desired information. Essentially, the retrieval problem can be considered an optimization problem where the goal is to find the optimal set of parameters that minimize the difference between the estimated and observed values.
Multiple methods exist for minimizing the error between the simulated radar variables and the measured variables. A relatively simple approach that was first used is the Gauss–Newton method. While it can quickly converge to a solution, the technique will often only find local minima rather than the global minimum of the solution space. Another early attempt used the genetic algorithm (GA), which very reliably found better solutions. The GA, however, was very computationally intensive and relied on finetuning of the crossover and mutation factors to efficiently solve for the DSD.
Particle swarm optimization (PSO) was ultimately used for this work since it is comparatively more efficient at seeking the most representative DSD. Figure 6 shows the organization of the PSO application. For a gate with preprocessed reflectivity and differentialphase measurements, particles with random N_{0}, μ, and Λ are initialized. Ranges for the particle positions are chosen according to commonly observed intervals (Zhang et al., 2001). In this manner, a coordinate space is created in which each DSD parameter is a dimension.
The three coordinates of each particle position collectively define a DSD that is used to calculate Z^{S} and ${\mathit{\varphi}}_{\mathrm{DP}}^{\mathrm{S},\mathrm{C}}$. Various cost functions that measure the distance from the truth data were tried. The fitness function given by Eq. (7) allows for relative weightings to be applied to each variable. Experimentally determined values (α = 5, β = γ = 1) led to reliable retrievals. It should also be noted that Z^{S} is expressed as a logarithmic value rather than in linear units.
For a gate's retrieval, the cost of every particle is calculated, and the iteration's current best solution and the global best solution of all iterations are recorded. The particle positions are then updated according to Eq. (8), where ϵ is the local acceleration factor, ζ is the global acceleration factor, r_{1} and r_{2} are random numbers between zero and one that are generated for each particle, and i denotes the current iteration. This allows the particles to move towards the current and global best solutions and possibly find better solutions along the path. The convergence speed must be weighed against the possibility of “overshooting” viable candidate solutions. ϵ and ζ were determined experimentally, and values of 0.15 and 0.0015 were used for processing the overall data set.
The swarm consists of 5000 particles which are allowed 400 iterations for each retrieval. It should be noted that the setting of 400 iterations is intended solely for prototype algorithm development, and computational efficiency has not yet been addressed. A simple iteration control algorithm could be implemented to terminate the computation once the rootmeansquare error reaches the predefined threshold. Any embedded solution should take advantage of speed improvements gained from parameter tuning or even substituting a more efficient technique in place of PSO.
2.3.3 Retrieval along radial
The goal of the retrieval process is to obtain a DSD similar to the disdrometer measurement. To achieve this, the PSO retrieval algorithm is applied iteratively starting from the gate nearest the radar. Initially, the first gate is assumed to be unattenuated, and its measured reflectivity and differential phase are directly input into the PSO algorithm. After completing the retrieval for a gate, attenuation is calculated using Eq. (5), and the reflectivity for the next farthest gate is adjusted. In this manner, each gate's input reflectivity accounts for attenuation experienced between the radar and that gate. This process, illustrated in Fig. 7, continues at 4 km intervals until reaching the disdrometer location.
2.4 Simulation
The role of the additional frequency in providing extra information is one of the key questions addressed in this study. We demonstrate the viability of employing dualpolarization measurements at Cband and Sband for rainfall estimation through simulations. The algorithm was evaluated under ideal conditions by extracting disdrometerrecorded distributions and fitting each to a DSD with gamma parameters as discussed in Sect. 2.3.1. For each distribution, Z and K_{DP} values were calculated using Eqs. (3) and (4). A total of 700 distributions were generated.
These calculated Z and K_{DP} values were treated as true observations, devoid of confounding factors such as elevation differences, noise, or unaccounted attenuation discussed in Sect. 2.1. Rainfall rates were calculated using the fitted DSDs according to Eq. (9) with v(D) given by Eq. (10), where D is given in mm and v(D) is in mm s^{−1} (Ulbrich, 1983). These rates are considered true if perfectly observed.
The cost function, defined in Eq. (7), was modified to operate directly on specific differential phase rather than differential phase, according to Eq. (11) with simplified weights: $\mathit{\alpha}=\mathit{\beta}=\mathit{\gamma}=\mathrm{1}$.
The retrieval process was conducted using this study's dualfrequency approach on a single gate. For comparison, retrievals were also performed using a single frequency by setting γ to zero, thus discarding the Cband contribution. Additionally, retrievals incorporated a μ−Λ constraint (Eq. 12, which Seela et al., 2018 demonstrated to be accurate for Taiwan rain systems during summer months). Rain rates for each retrieval method were calculated and compared to true values using relative absolute error (RAE), as defined in Eq. (13), where R_{d} represents the truth rain rate and R denotes the rainfall rates from the methods being compared.
As a final comparative benchmark, rain rates were also calculated using varying Z–R relationships. Equation (14) represents perhaps the most common meteorological radar relationship (Ulbrich and Lee, 1999), Eq. (15) is specific to Taiwan as proposed by Chang et al. (2021), and Eq. (16) is derived from fitting calculated Z and R values found in the initial steps of the simulation.
Figure 8 illustrates the cumulative distribution of errors for each method. In this plot, the cumulative portion of errors is shown on the y axis, while the sorted error values of each method are displayed on the x axis. As an example interpretation of the plot, 90 % of current study's errors (blue) are less than 0.2 in terms of RAE, while 90 % of the Taiwan Z–Rbased errors (dashed black) are not contained until an error level of 0.78 RAE. Interpreting the plot from the constant RAE perspective, 65 % of the current study's errors are below 0.1, while 60 % of μ–Λ (purple) method's errors and 30 % of the Taiwan Z–R (black) method's errors are below 0.1 RAE. The conclusion of the simulations can be drawn from the plots: the dualfrequency approach provides a modest improvement over singlefrequency retrievals, even with a relevant μ–Λ constraint, and a significant improvement compared to all Z–Rbased rates. Under these ideal simulated conditions, the median RAE of the dualfrequency approach was 0.0623, while the μ–Λ singlefrequency approach was 0.0725. The median RAE of the bestperforming (at the 50th percentile) Z–R relationship (Z=207R^{1.45}) was 0.1861.
At the parameter level, the benefit of an additional frequency can also be investigated while still using the gamma model assumption. Distributions with varying μ and Λ values, while keeping N_{0} constant, can be used to calculate specific differentialphase values at two wavelengths (10.0 cm for Sband and 5.0 cm for Cband). These values are depicted in Fig. 9. This simulation highlights that specific differentialphase values exhibit distinct separations at the two wavelengths. This simulation underscores the potential of using dualpolarization measurements at C and S radar bands for DSD retrieval.
The evaluation of the proposed algorithm was carried out through both qualitative and quantitative methods. Initially, the retrieved DSDs were subjectively compared to the DSDs measured by the disdrometers at the times nearest to the radar scans. This qualitative comparison provides a preliminary measure of the algorithm's overall accuracy as adjustments were made to its parameters. A subsequent discussion of the role of using two frequencies rather than a singlefrequency approach highlights the value of using the Cband reflectivity (unused in the retrieval) for further validation of the alongradial results. Following the finetuning of the cost function and PSO parameters, a more rigorous quantitative analysis was performed. This involved calculating the rainfall rates based on the disdrometerrecorded DSDs, which were considered the benchmark, and comparing these rates to those obtained from the retrieved DSDs. Additionally, the calculated rainfall rates were contrasted with those derived using the traditional Z–R relationship methodology, offering a comprehensive evaluation of the algorithm's performance.
For the performance validation of this study, data spanning 10 d from 2017 and 2018 were utilized. A potential case is defined as one plan position indicator (PPI) scan of each radar. The stringent criteria for time synchronization narrowed the data set to 167 cases that not only met the synchronization requirements of the radar scans but also had the requisite data quality and fell within the disdrometer range requirement.
3.1 Qualitative assessment
Figure 10 showcases a representative selection of the retrieved DSDs the algorithm generated as compared to the disdrometer data. Cases were systematically sampled across the data set's collection period to accurately reflect the algorithm's retrieval performance. These cases strictly adhere to the time synchronization criteria established by our data quality standards and exhibit reflectivity values at the terminal gate that exceed the minimum threshold. The radar scan time tags are specified with precision down to the second, whereas disdrometer measurements are recorded on a perminute basis. Consequently, both the preceding and subsequent disdrometer timestamps are presented (red and green, respectively), in addition to the closest time (blue). This approach acknowledges instances where these adjacent timestamps may more accurately represent the radar data.
Six of these eight cases show a high degree of agreement across the spectrum of drop sizes, with particularly strong correlation observed for the measurements of moderate drop sizes. It can be found from Fig. 10 that the retrieved DSD fits the disdrometer observations very well for larger drops (D > 0.5 mm). Deviations can be found in cases such as 20170613 – 08:04:49 UTC and 20180107 – 09:26:32 UTC for smaller drops (D < 0.5 mm). This is a predictable result given the retrieval input parameters are heavily dominated by contributions of larger diameters. Not all disdrometer sizes need to be equally prioritized for accurate fitting. Examination of Eqs. (3) and (4) reveals that the contribution to radar parameters from each diameter increases with the size of the drop. Furthermore, previous research evaluating the accuracy of disdrometers across various drop sizes indicated that drops of 0.6 mm and larger are the first reliably measured sizes by optical disdrometers (Tokay et al., 2001). This finding supports the notion that inaccuracies in measuring smaller drop sizes do not significantly impact the calculations of reflectivity or DSDderived metrics such as rain rate and attenuation. While some researchers have proposed that the goal of a retrieval should be to fit the predominance of the size spectrum (Adirosi et al., 2013), it is clearly more important to represent medium to large diameters when the integrated parameters are the focus.
To ensure the accuracy of the retrieval algorithm, one can utilize the Cband parameters obtained from the retrieved DSDs to validate the Sband parameters. To accomplish this, the attenuation that is calculated from the retrieved DSDs can be applied to correct the Cband reflectivities affected by the large attenuation experienced over the path. Figure 11 suggests that the Sband reflectivity is less affected by attenuation as expected. The dashed blue lines in the plots represent the raw reflectivity values measured by the Sband radar, which demonstrates the relatively low attenuation at this wavelength. The solid blue line in the plot represents the attenuationcorrected reflectivity values measured by the Sband radar. The corrected Cband reflectivities are shown as the solid red line in each plot and are seen to match well with the Sband values. This indicates that the attenuation factor derived from the retrieved DSD is effective in converting the Cband reflectivities to equivalent Sband values. By using this correction factor, we can confirm that the retrieved DSD accurately predicts the atmospheric effects at both radar bands, and thus the corrected Cband reflectivity serves as a “sanity check” for the alongradial retrieval. The input (dashed) and retrieved (solid) differentialphase profiles are also provided for each example.
3.2 Quantitative assessment
A comprehensive quantitative analysis was conducted on the 167 cases, utilizing the disdrometers in Table 2, with a particular focus on the algorithm's accuracy in QPE, where the importance of accuracy is most prominently showcased.
At the parameter level, Fig. 12 demonstrates the retrieval routine produces DSDs corresponding to radar parameters which effectively match the inputs. The x axis represents the inputs to the retrieval based on measurements, taking into account attenuation in the case of Z. The y axis represents the parameters output by the algorithm retrievals for the final gate containing the disdrometer. Metrics for the retrieval pairs are indicated for each parameter: the mean bias (MB), rootmeansquare error (RMSE), and the correlation coefficient (CC). These are defined as follows: the mean bias MB $=\frac{\langle {X}_{\mathrm{Retrieved}}\rangle}{\langle {X}_{\mathrm{Observed}}\rangle}$; the rootmeansquare error RMSE $={\u2329{\left({X}_{\mathrm{Retrieved}}{X}_{\mathrm{Observed}}\right)}^{\mathrm{2}}\u232a}^{\mathrm{1}/\mathrm{2}}$; and the correlation coefficient CC $=\frac{\langle \left({X}_{\mathrm{Retrieved}}\langle {X}_{\mathrm{Retrieved}}\rangle \right)\left({X}_{\mathrm{Observed}}\langle {X}_{\mathrm{Observed}}\rangle \right)\rangle}{{\mathit{\sigma}}_{{X}_{\mathrm{Retrieved}}}{\mathit{\sigma}}_{{X}_{\mathrm{Observed}}}}$, where angle brackets indicate the mean, X_{Retrieved} is the algorithm output for the given parameter (Z^{S}, ${K}_{\mathrm{DP}}^{\mathrm{S}}$, ${K}_{\mathrm{DP}}^{\mathrm{C}}$), X_{Observed} is the observed input value for the parameter, and σ is the standard deviation of the parameter.
In the context of QPE, rainfall rates were calculated using the retrieved DSDs according to Eq. (9) with v(D) given by Eq. (10) as in Sect. 2.4. For comparison, the Sband reflectivity data were applied using Eq. (15), which is the regionspecific relationship also referenced in Sect. 2.4 (Chang et al., 2021). The Parsivel data were used to calculate the ground truth via Eq. (17), where D_{j} is the equivolume diameter of the jth recorded drop, S is the effective sampling area, and Δt is the sampling period (Raupach and Berne, 2015).
The evaluation involved comparing the rainfall rates derived from the Z–R relationship, that of a μ−Λconstrained (Eq. 12) Sbandonly retrieval, and that of this study's retrieval algorithm against the ground truth values from the Parsivel disdrometer, whose biases were discussed in Sect. 2. This comparison was based on RAE following the same approach as discussed in Sect. 2.4.
Figure 13 features a plot of the cumulative distribution of RAE for all three methods. The retrieval algorithm's accuracy is depicted by the blue line with circles, while the μ–Λ constrained retrieval is shown in purple, and the Z–R relationship benchmark is indicated with the dashed black line. This visualization highlights that the proposed method of estimating rainfall rates using retrieved DSD parameters significantly enhances accuracy over a regionspecific Z–R relationship (Chang et al., 2021). Specifically, the median RAE for the Z–R method is 0.76, while the retrieval results correspond to a median of 0.53, marking a notable improvement of 30.3 % in this study's context. The dualfrequency approach is comparable to the singlefrequency constrained approach, and the only claimed relative benefits are not needing to ascertain a regional μ–Λ relationship or to use cases in which the precipitation deviates from the assumed relationship.
It is important to acknowledge a notable limitation observed in the performance of the retrieval algorithm, as reflected in the resulting cases with large RAE. Specifically, the algorithm has the potential to significantly overestimate the rainfall rate. Closer examination reveals that such discrepancies arise from illposed conditions for retrieval. Although incorporating multiple frequencies enhances the capability to retrieve the DSD, the specific differentialphase and reflectivity values must still exhibit a consistent correlation. Deviations from this correlation can mislead the algorithm, yielding suboptimal outcomes. Therefore, incorporating selfconsistency relations between the input variables stands out as a promising direction for enhancing algorithmic accuracy in future research endeavors (Park et al., 2005; Giangrande et al., 2013; Gorgucci et al., 1999; ReinosoRondinel et al., 2018).
A novel approach to retrieving the DSD using PSO has been discussed. The outlined dualfrequency dualpolarization method can provide significantly improved estimates for rainfall compared to Z–R relationship estimates and offers modest improvement to retrievals utilizing relevant μ–Λ constraints. While the retrieval is unconstrained regarding the gamma distribution parameters, the price is the additional data needed. The authors predict that radar systems utilizing more than one frequency will continue to become more commonplace and that producing the required synchronized measurements will be more feasible with the adoption of phased array systems. The value of dual polarization in radars is already universally accepted. Algorithms such as the one prototyped in this work will become more valuable as radar systems produce data with this type of increased diversity.
There are several limitations in this research that have been briefly mentioned and that future work should address. A larger data set or alternative preprocessing criteria should be assessed for longterm evaluation of the approach. While 10 d of data were screened, the selection process for useful cases was highly discriminatory. Only measurements from the two radars that were synchronized within 1 min became candidates. Of these, only data with corresponding disdrometer truth data and minimum terminal gate reflectivity were included in the processing pool. Preprocessing the measurements led to the greatest reduction in potential data. The majority of these cases were excluded because of data quality issues that were primarily due to unreliable differentialphase profiles. Future work should exhaustively assess the approach by processing a larger amount of highquality data.
One focus of this project was to show that N_{0}, μ, and Λ are independently solvable when incorporating both multifrequency and multipolarization information. Restricted relationships between the variables are still highly useful. If the intended application of any retrieval algorithm can utilize a constraint with high confidence, the DSD retrieval process becomes more efficient. The tradeoff is that shaping characteristics may not be captured, which is why study of more flexible retrieval processes should continue.
The optimization approach involves various adjustable parameters, including swarm size, number of iterations, and acceleration coefficients. Finetuning these parameters could result in faster and more optimal results. Moreover, adapting the algorithm to function as an embedded application for field testing is also a promising area for further development.
The data sets and source code used in this study are available from the corresponding author upon request.
YW and DD created the initial concept of utilizing dualfrequency, dualpolarization measurements for DSD retrieval and developed preliminary algorithms to demonstrate its feasibility. DD further applied the PSO technique, processed the data set, and prepared the drafts of the paper. YW and DD collaborated on the development of the final algorithm. PLC provided and processed radar data from the Central Weather Bureau (CWB) and was further involved in algorithm discussion and article authoring.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The authors thank the radar engineers at CWB for helping collect and process the radar data used in the study. The CWB provided the unique data used in the project.
This paper was edited by Jorge Luis Chau and reviewed by Ricardo ReinosoRondinel and three anonymous referees.
Adirosi, E., Baldini, L., Lombardo, F., Russo, F., and Napolitano, F.: Comparison of different fittings of experimental DSD, in: AIP Conference Proceedings, AIP Conf. Proc., 1558, 1669–1672, https://doi.org/10.1063/1.4825850, 2013. a
AnguloMartínez, M., Beguería, S., Latorre, B., and FernándezRaga, M.: Comparison of precipitation measurements by OTT Parsivel^{2} and Thies LPM optical disdrometers, Hydrol. Earth Syst. Sci., 22, 2811–2837, https://doi.org/10.5194/hess2228112018, 2018. a
Brandes, E. A., Zhang, G., and Vivekanandan, J.: Experiments in Rainfall Estimation with a Polarimetric Radar in a Subtropical Environment, J. Appl. Meteorol., 41, 674–685, https://doi.org/10.1175/15200450(2002)041<0674:EIREWA>2.0.co;2, 2002. a
Bringi, V. N., Thurai, M., and Hannesen, R.: DualPolarization weather radar handbook, 2nd edn., Neuss, 163 pp., 2005. a, b
Cao, Q., Zhang, G., Brandes, E. A., and Schuur, T. J.: Polarimetric Radar Rain Estimation through Retrieval of Drop Size Distribution Using a Bayesian Approach, J. Appl. Meteorol. Clim., 49, 973–990, https://doi.org/10.1175/2009JAMC2227.1, 2010. a
Chandrasekar, V., Li, W., and Zafar, B.: Estimation of raindrop size distribution from spaceborne Radar observations, IEEE T. Geosci. Remote, 43, 1078–1086, https://doi.org/10.1109/TGRS.2005.846130, 2005. a
Chang, P.L., Zhang, J., Tang, Y.S., Tang, L., Lin, P.F., Langston, C., Kaney, B., Chen, C.R., and Howard, K.: An Operational MultiRadar MultiSensor QPE System in Taiwan, B. Am. Meteorol. Soc., 102, E555–E577, 2021. a, b, c, d
Eccles, P. J.: Comparison of Remote Measurements by Single and DualWavelength Meteorological Radars, IEEE T. Geosci. Elect., 17, 205–218, https://doi.org/10.1109/TGE.1979.294650, 1979. a
Giangrande, S. E., McGraw, R., and Lei, L.: An Application of Linear Programming to Polarimetric Radar Differential Phase Processing, J. Atmos. Ocean. Tech., 30, 1716–1729, https://doi.org/10.1175/JTECHD1200147.1, 2013. a
Gorgucci, E. and Baldini, L.: A SelfConsistent Numerical Method for Microphysical Retrieval in Rain Using GPM DualWavelength Radar, J. Atmos. Ocean. Tech., 33, 2205–2223, https://doi.org/10.1175/jtechd160020.1, 2016. a, b
Gorgucci, E., Scarchilli, G., and Chandrasekar, V.: Specific Differential Phase Estimation in the Presence of Nonuniform Rainfall Medium along the Path, J. Atmos. Ocean. Tech., 16, 1690–1697, https://doi.org/10.1175/15200426(1999)016<1690:SDPEIT>2.0.CO;2, 1999. a
Jaffrain, J. and Berne, A.: Experimental Quantification of the Sampling Uncertainty Associated with Measurements from PARSIVEL Disdrometers, J. Hydrometeorol., 12, 352–370, https://doi.org/10.1175/2010jhm1244.1, 2011. a, b, c
Kozu, T., Nakamura, K., Meneghini, R., and Boncyk, W.: Dualparameter radar rainfall measurement from space: a test result from an aircraft experiment, IEEE T. Geosci. Remote, 29, 690–703, https://doi.org/10.1109/36.83983, 1991. a
Kummerow, C., Mack, R. A., and Hakkarinen, I. M.: A SelfConsistency Approach to Improve Microwave Rainfall Rate Estimation from Space, J. Appl. Meteorol., 28, 869–884, https://doi.org/10.1175/15200450(1989)028<0869:ASCATI>2.0.CO;2, 1989. a
Le Loh, J., Chang, W.Y., Hsu, H.W., Lin, P.F., Chang, P.L., Teng, Y.L., and Liou, Y.C.: LongTerm Assessment of the Reflectivity Biases and WetRadome Effect Using Collocated Operational S and CBand DualPolarization Radars, IEEE T. Geosci. Remote, 60, 1–17, https://doi.org/10.1109/TGRS.2022.3170609, 2022. a
Mardiana, R., Iguchi, T., and Takahashi, N.: A dualfrequency rain profiling method without the use of a surface reference technique, IEEE T. Geosci. Remote, 42, 2214–2225, https://doi.org/10.1109/TGRS.2004.834647, 2004. a
Marshall, J. S. and Palmer, W. M. K.: THE DISTRIBUTION OF RAINDROPS WITH SIZE, Journal of Meteorology, 5, 165–166, https://doi.org/10.1175/15200469(1948)005<0165:TDORWS>2.0.CO;2, 1948. a
Marzoug, M. and Amayenc, P.: A Class of Single and DualFrequency Algorithms for RainRate Profiling from a Spaceborne Radar. Pad I: Principle and Tests from Numerical Simulations, J. Atmos. Ocean. Tech., 11, 1480–1506, https://doi.org/10.1175/15200426(1994)011<1480:ACOSAD>2.0.CO;2, 1994. a
Meneghini, R., Kozu, T., Kumagai, H., and Boncyk, W. C.: A Study of Rain Estimation Methods from Space Using DualWavelength Radar Measurements at NearNadir Incidence over Ocean, J. Atmos. Ocean. Tech., 9, 364–382, https://doi.org/10.1175/15200426(1992)009<0364:asorem>2.0.co;2, 1992. a
Meneghini, R., Kumagai, H., Wang, J., Iguchi, T., and Kozu, T.: Microphysical retrievals over stratiform rain using measurements from an airborne dualwavelength radarradiometer, IEEE T. Geosci. Remote, 35, 487–506, https://doi.org/10.1109/36.581956, 1997. a
Park, S.G., Bringi, V. N., Chandrasekar, V., Maki, M., and Iwanami, K.: Correction of Radar Reflectivity and Differential Reflectivity for Rain Attenuation at X Band. Part I: Theoretical and Empirical Basis, J. Atmos. Ocean. Tech., 22, 1621–1632, https://doi.org/10.1175/jtech1803.1, 2005. a
Raupach, T. H. and Berne, A.: Correction of raindrop size distributions measured by Parsivel disdrometers, using a twodimensional video disdrometer as a reference, Atmos. Meas. Tech., 8, 343–365, https://doi.org/10.5194/amt83432015, 2015. a, b
ReinosoRondinel, R., Unal, C., and Russchenberg, H.: Adaptive and HighResolution Estimation of Specific Differential Phase for Polarimetric XBand Weather Radars, J. Atmos. Ocean. Tech., 35, 555–573, https://doi.org/10.1175/JTECHD170105.1, 2018. a
Ryzhkov, A. V., Kumjian, M. R., Ganson, S. M., and Khain, A. P.: Polarimetric Radar Characteristics of Melting Hail. Part I: Theoretical Simulations Using Spectral Microphysical Modeling, J. Appl. Meteorol. Clim., 52, 2849–2870, https://doi.org/10.1175/jamcd13073.1, 2013. a
Sauvageot, H. and Lacaux, J.P.: The Shape of Averaged Drop Size Distributions, J. Atmos. Sci., 52, 1070–1083, https://doi.org/10.1175/15200469(1995)052<1070:TSOADS>2.0.CO;2, 1995. a
Seela, B. K., Janapati, J., Lin, P.L., Wang, P. K., and Lee, M.T.: Raindrop size distribution characteristics of summer and winter season rainfall over north Taiwan, J. Geophys. Res., 123, 11602–11624, 2018. a
Sheppard, B. E.: Sampling Errors in the Measurement of Rainfall Parameters Using the Precipitation Occurrence Sensor System (POSS), J. Atmos. Ocean. Tech., 24, 125–140, https://doi.org/10.1175/JTECH1956.1, 2007. a, b
Tokay, A., Kruger, A., and Krajewski, W. F.: Comparison of Drop Size Distribution Measurements by Impact and Optical Disdrometers, J. Appl. Meteorol., 40, 2083–2097, https://doi.org/10.1175/15200450(2001)040<2083:CODSDM>2.0.CO;2, 2001. a
Ulbrich, C. W.: Natural Variations in the Analytical Form of the Raindrop Size Distribution, J. Clim. Appl. Meteorol., 22, 1764–1775, https://doi.org/10.1175/15200450(1983)022<1764:NVITAF>2.0.CO;2, 1983. a, b
Ulbrich, C. W. and Lee, L. G.: Rainfall measurement error by WSR88D radars due to variations in Z–R law parameters and the radar constant, J. Atmos. Ocean. Tech., 16, 1017–1024, 1999. a
Waterman, P.: Matrix formulation of electromagnetic scattering, P. IEEE, 53, 805–812, https://doi.org/10.1109/PROC.1965.4058, 1965. a
Williams, C., Bringi, V., Carey, L., Chandrasekar, V., Gatlin, P., Haddad, Z., Meneghini, R., Munchak, s., Nesbitt, S., Petersen, W., Tanelli, S., Tokay, A., Wilson, A., and Wolff, D.: Describing the Shape of Raindrop Size Distributions Using Uncorrelated Raindrop Mass Spectrum Parameters, J. Appl. Meteorol. Clim., 53, 1282–1296, https://doi.org/10.1175/JAMCD13076.1, 2014. a, b
Zhang, G., Vivekanandan, J., and Brandes, E.: A method for estimating rain rate and drop size distribution from polarimetric radar measurements, IEEE T. Geosci. Remote, 39, 830–841, https://doi.org/10.1109/36.917906, 2001. a, b