The Quality Control and Rain Rate Estimation for the X-Band Dual-Polarization Radar: A Study of Propagation of Uncertainty

: In this study, the speciﬁc differential phase ( K dp ) is applied to attenuation correction for radar reﬂectivity Z H and differential reﬂectivity Z DR , and then the corrected Z H , Z DR , and K dp are studied in the rain rate ( R ) estimation at the X-band. The statistical uncertainties of Z H , Z DR , and R are propagated from the uncertainty of K dp , leading to variability in their error characteristics. For the attenuation correction, a differential phase shift ( Φ dp )-based method is adopted, while the statistical uncertainties σ ( Z H ) and σ ( Z DR ) are related to σ ( K dp ) via the relations of K dp -speciﬁc attenuation ( A H ) and K dp -speciﬁc differential attenuation ( A DP ), respectively. For the rain rate estimation, the rain rates are retrieved by the power-law relations of R ( K dp ) , R ( Z h ) , R ( Z h , Z dr ) , and R ( Z h , Z dr , K dp ) . The statistical uncertainty σ ( R ) is propagated from Z H , Z DR , and K dp via the Taylor series expansion of the power-law relations. A composite method is then proposed to reduce the σ ( R ) effect. When compared to the existing algorithms, the composite method yields the best performance in terms of root mean square error and Pearson correlation coefﬁcient, but shows slightly worse normalized bias than R ( K dp ) and R ( Z h , Z dr , K dp ) . The attenuation correction and rain rate estimation are evaluated by analyzing a squall line event and a prolonged rain event. It is clear that Z H , Z DR , and K dp show the storm structure consistent with the theoretical model, while the statistical uncertainties σ ( Z H ) , σ ( Z DR ) and σ ( K dp ) are increased in the transition region. The scatterplots of Z H , Z DR , and K dp agree with the self-consistency relations at X-band, indicating a fairly good performance. The rain rate estimation algorithms are also evaluated by the time-series of the prolonged rain event, yielding strong correlations between the composite method and rain gauge data. In addition, the statistical uncertainty σ ( R ) is propagated from Z H , Z DR , and K dp , showing higher uncertainty when the large gradient presents.

The major limitation of the X-band frequency is strong rain attenuation effects on the conventional Z H , but the X-band radars have some advantages over the longer wavelengths, including finer resolution with a smaller antenna, easier mobility, and lower cost. The radars with dual-polarization capacity can also measure the backscattering properties of hydrometeors at two orthogonal

Experimental Site and Data
This study analyzes the rainfall measurements collected by the MZZU radar and rain gauges (RGs) during the Missouri Experimental Project to Stimulate Competitive Research (EPSCoR). In this section, we first discuss the significance of weather radar surveillance in central Missouri, and then describe the radar system and data.

Significance of Weather Radar Surveillance in Central Missouri
Central Missouri has a typical continental climate characterized by strong seasonality in precipitation due to the inland location. In summer, moist and warm air masses from the Gulf produce abundant rain by frontal and convective processes. The total amount of summer precipitation is twice than that of winter precipitation on average. In winter, dry and cold air masses invade from northern plains, producing rain, snow and mixed-phase precipitation. Spring and fall are transitional seasons when the precipitation abruptly changes due to fast-moving cold fronts. Additionally, hail occurs throughout the entire year, but it is more likely to be observed in spring. In all, central Missouri is a good experimental location to study a variety of precipitating events using the radar observations.
In Missouri, the amount of precipitation water averages 1083 mm per year, a large portion of which is evaporated back to the atmosphere or transpired by plants. In the summer months when the loss of the water is high, a drought may occur due to the lack of rainfall. According to [38], if a month has a precipitation amount of less than 40% of the climatological average, the average probability of drought is about 15%. Furthermore, for the months of April through October, the probability of drought is 13% when three consecutive months have a precipitation amount of less than 60% of the historical average. The drought may occur as a result of precipitation deficiency, indicating the importance of precipitation surveillance in central Missouri.
Central Missouri is drained by tributaries of the Missouri River, which flows eastward into the Mississippi River.
Tributary flooding is expected once or twice per year, while thunderstorm-producing flash flooding occurs more frequently in the minor streams in spring. In contrast, the main steam flooding caused by prolonged periods of heavy rain is less frequent, but it may occur in the summer. On the other hand, there are about 30 tornadoes in Missouri each year, including 8 strong or violent tornadoes. Records show that about 70% tornadoes occur between March and June, during which 25% occur in March. Notably, tornadoes have been reported in Missouri every year. Because of the high spatiotemporal resolution, the radar can provide timely warnings to reduce the fatalities and economic losses of flooding and tornadoes.

System and Dataset Description
As part of the EPSCoR program, a new X-band dual-polarization Doppler radar (MZZU) was deployed in South Farm Research Center (38.906 • N, 92.269 • W) in central Missouri in mid-2015. The data have been available in real-time since September 2015. During the study period between April 2016 and June 2018, the MZZU radar operated in a volumetric scanning mode at 9 elevation angles to obtain instantaneous measurements of precipitation every 3-5 min. The radar range is equally sampled every 260 m between 1.3 and 94.64 km, yielding 360 gates for each ray. Moreover, the radar uses a parabolic dish antenna with a beam width of 1.27 • , leading to an angular resolution of about 1 • . For a distance of 60 km, the radar beam gives a width of approximately 1 km. Overall, the radar can obtain the data satisfying the criteria of a minimum temporal resolution of 1-5 min and a minimum spatial resolution of 1 km, as recommended by [39,40].
From 6 April 2016 to 2 June 2018, the radar collected 19,698 volumetric scans. To analyze the performance of radar quality control, we selected two representative events occurred on 7 March 2017 and 2-4 July 2016. Table 1 presents the storm characteristics, including movement, size, shape, duration and type. On one hand, the convective case records the phenomena of rain and hail accompanied with tornadoes. When the storm is propagating southeastward, it is manifested as a straight-line shape, lasting 5-6 h. On the other hand, the stratiform case presents prolonged rain on a relatively large scale. , respectively. Due to radar side lobes, the radar measurements at Bradford and Sanborn suffer from clutter contamination, giving a reduction of correlation coefficient of the MZZU-RG comparison data. Therefore, we consider the second elevation at 2 • at Bradford and Sanborn, with the beam heights of 410 m and 465 m ASL, respectively. The rain gauges are carefully maintained and calibrated at least once a year for hydrological applications, and a simple quality control is performed to remove zeros and anomalously high values (250 mm in an hour).
It is noted that the radar collects instantaneous measurements every 4-5 min, whereas RGs obtain the precipitation accumulations over 60 min. Therefore, it is necessary to average 12-15 consecutive radar scans to derive the hourly rain amounts. Moreover, the radar and RG data may both include some uncertainties, such as spatial variability of rainfall, measurement errors, calibration errors, resolution errors, K dp -rain rate conversion errors and errors due to winds and local disturbances (e.g., [41,42]). The uncertainties related to the spatial variability and measurement errors become even larger during the events of flash flooding [2]. To have a better agreement between the two types of instruments, the radar estimates are averaged across an area of 2 • ×0.5 km, i.e., 1 azimuthal angle and 1 range gate on either side of an RG. The selected smoothing region is consistent with the previous studies at other frequencies (e.g., [43]).

Quality Control
The quality control of the weather radar data is crucial for the radar products and model analysis.
In this section, we first remove the clutter using a simple clustering procedure and then estimate the K dp data fields with a Gaussian mixture method (GMM). Finally, we derive the attenuation-corrected Z H and Z DR and their statistical uncertainties with the GMM K dp .

Clutter Removal
The noise and residual clutter in Z H , Z DR , and Φ dp need to be removed prior to the rain rate estimation. As this step is independent of the K dp estimation and attenuation correction, the existing techniques may be adopted with respect to the radar specifications and research purposes [44][45][46]. In this study, a simplified Gaussian mixture method is applied to the MZZU radar data [37]. The range profile of Φ dp is first separated into multiple clusters, each of which corresponds to a data segment. Secondly, a segment is classified into weather and clutter by comparing the weight of each data segment with the pre-set thresholds. Thirdly, the clusters are evaluated based on the types of the segments. Fourthly, azimuthally isolated data points are removed to yield the clean Φ dp data. Finally, the Z H and Z DR data are processed referring to the clean Φ dp .

K dp Estimation
The prevalent method for K dp estimation is a linear regression model (LRM), which fits a first-order polynomial function to the range profile of Φ dp [12]. It is formulated as where Φ dp is denoted as y, the radar range is denoted as x, β 0 is the intercept, β 1 is the slope, and is the error. The error is individually and independently distributed, and drawn from a zero-mean Gaussian distribution with variance σ 2 . The specific differential phase K dp is then estimated by the slope β 1 , yielding: where n is the number of gates andx is the mean value in the n gates. As is the only stochastic variable, i.e., σ 2 (y) = σ 2 , the variance of K dp is estimated as For linear polarization, σ 2 is related to the variance in phase for each of the two polarizations [47], leading to σ(Φ dp ) = 2.61 • using 32 pulses for the MZZU radar. By substituting to Equation (3), it yields σ(K dp ) = 0.1 deg km −1 if we apply 31 range gates with a gate spacing of 260 m. Moreover, as the variance of K dp is a constant inherited from σ 2 (Φ dp ), lower K dp has higher coefficient of variation σ(K dp )/K dp compared to higher K dp . For example, if the radar collects a datum of K dp = 0.1 deg km −1 , its coefficient of variation is about 100%, whereas the coefficient of variation is only 10% for K dp = 1.0 deg km −1 .
Nevertheless, it is clear that the assumption of σ 2 (ŷ) = σ 2 may not always hold when the randomness of y is considered in the data generation model. If the variable y follows a probability distribution with meanȳ and variance σ 2 (y), the estimateŷ iŝ From Equations (4) and (5), we can see that the data generation model leads to a better representation of the random error of Φ dp . Therefore, [37] proposed a probabilistic method based on a Gaussian mixture model for K dp estimation using the MZZU radar data. The Gaussian mixture method (GMM) can not only estimate the expected values of K dp by differentiating the conditional expectation of Φ dp , but also calculate the variance of K dp by regarding the errors in the calculation of the first derivative of Φ dp .
In this method, the range gate and the Φ dp data are fitted into a Gaussian mixture, i.e., where x is the range gate, y is the Φ dp , m is the number of components in the mixture, w i is a weight with ∑ m i=1 w i = 1, and N (x, y; µ i , Σ i ) represents a multi-variant Gaussian distribution with mean . The parameters are iteratively computed by an Expectation-Maximization algorithm. Subsequently, the probability of y conditioned on x is also a Gaussian mixture, yielding where w y|x i , µ y|x i , and Σ y|x i are the conditional weight, mean, and covariance, respectively. The mean differential phase shift Φ dp is then retrieved by the mathematical expectation of Equation (7).
, the expected value of y conditioned on x is expressed as Meanwhile, its variance is given as Substituting Equation (14) to Equation (5), the variance of Φ dp is finally obtained by For simplicity, we assume that σ is 2.61 deg km −1 for the MZZU radar. By taking the derivative of Equation (11) with respect to the radar range x, K dp is retrieved by Moreover, the variance of K dp is approximated by a Taylor series expansion, i.e., GMM σ 2 (K dp ) = σ 2 (ŷ |x) ≈ [E (y|x)] 2 σ 2 (φ dp ) = (K dp ) 2 σ 2 (φ dp ), where σ 2 (Φ dp ) is given in Equation (15). The derivative of K dp is expressed as where It can be found that K dp is estimated by the derivative of the mathematical expectation of Φ dp conditioned on the range x, while the variance of K dp is related to the range derivative of K dp and the variance of Φ dp .
It is known that the X-band radar data are affected by ambiguous Φ dp and backscattering differential phase shift δ co ([37], Figure 5). To correct these errors, we compare the means of two consecutive clusters along the range. If the latter cluster has a mean larger than the former one by 80 • , the latter one is labeled as ambiguous Φ dp and unfolded by adding 180 • . On the other hand, the latter cluster is removed due to δ co if its mean is 85 • smaller than the former one or if its weight is less than 0.0501.
Moreover, a finite impulse response (FIR) filter is applied to reduce the K dp variances [48]. In [37], the filtering for each ray is optimized in terms of relative square error, leading to a window length between 29 and 33 gates for the MZZU radar. As discussed in [49], the spatial resolution gives a strong correlation with the temporal resolution, giving a relation of r = 5t 0.3 where t in min and r in km. If the radar data are updated every 4.5 min, the spatial resolution requires a minimum of 7.85 km, equivalent to 31 range gates.

Attenuation Correction
It is important to do the attenuation correction before the rain rate estimation, as the Mie scattering significantly affects the power measurements at X-band frequency. The attenuation effects at X-band are significant due to the scattering and absorption of raindrops along the wave propagation path [12].
The T-matrix simulations show that A H and A DP at X-band are about 10 times and 5 times larger than those at S-and C-band, respectively [13,50]. Nevertheless, K dp is insensitive to the attenuation effects, and therefore useful for determining A H and A DP as a function of range. In this section, we apply GMM K dp to the attenuation correction for the MZZU radar.
The specific attenuation is related to the integral of the raindrop size distribution (DSD) and the imaginary part of forwarding scattering amplitude, which is defined as where D is the raindrop size in millimeter, N(D) is the drop spectrum in m −3 mm −1 , {·} is the imaginary part and f hh,vv (0, D) are the forward scattering amplitudes at the horizontal and vertical polarizations, respectively. The specific differential attenuation is then given as As K dp is also related to the forward scattering amplitude, A H and A DP may be calculated by an exponential function of K dp , respectively, i.e., where α, β, γ 1 , and γ 2 are constant coefficients. The variables Z H and Z DR are then corrected by where Z att H and Z att DR are corrected Z H and Z DR , respectively, and Z raw H and Z raw DR are raw Z H and Z DR , respectively. For simplicity, we assume the coefficients γ 1 and γ 2 are unity, leading to a φ dp -based attenuation correction method [13], where φ dp (0) is the initial φ dp , and φ dp (r) is the re-constructed φ dp at the rth gate. The coefficients α and β depend on the radar frequency, DSD and drop shape [23,51]. With the assumption of the drop shape of [13,52] calculate α and β using the T-matrix simulation based on observed DSDs. Following this work, we derive α and β as 0.25 and 0.05, respectively, using the MZZU-RG comparison dataset. The coefficients are similar to the ones presented in the previous studies [14,53]. From Equations (29) and (30), the variance of corrected Z H and Z DR are related to the variance of φ dp and the variance of raw Z H and Z DR , respectively. By assuming that Z H , Z DR , and φ dp are independent, we have where σ Z raw H and σ Z raw DR are 1.36 dB and 0.436 dB for the MZZU radar, respectively. Equations (31) and (32) indicate that the attenuation correction may increase the uncertainties for Z H and Z DR , and consequently produce higher errors in the Z H and/or Z DR -based rain rates. Finally, we use the Z H and Z DR after the attenuation correction in the rain rate estimation, while the variables in linear units are calculated by where Z h is in mm 6 m −3 and Z dr is unitless.

Rain Rate Estimation
In this section, the rain rate estimation algorithms are first derived from the MZZU-RG comparison data, and the statistical uncertainties σ(R) are then calculated by considering the propagation of the uncertainties of Z h , Z dr , and K dp in the power-law relations. The algorithms of R(Z h ), R(Z h , Z dr ), R(K dp ), and R(Z h , Z dr , K dp ) are optimized in terms of the root mean square error (RMSE). It is defined as where N is the sample size, R i is the radar hourly rain amount, and G i is the hourly gauge data. In addition to RMSE, the normalized bias (NB) and Pearson correlation coefficient (CORR) are also calculated to analyze the accuracy of the algorithms. They are given as whereR andḠ are the sample means for radar and gauge, respectively.

Retrieval of R
The linear radar reflectivity Z h is the sixth order moments of DSD [12,54], while the rain rate R is about 3.67th order moment. Consequently, the rain rate may be retrieved via a power-law relation of Z h and R [55], where Z h is in units of mm 6 m −3 , and the coefficients of 0.238 and 0.411 are optimized for the MZZU-RG data. To reduce the sensitivity to the variability of DSD, Nevertheless, the R(Z h ) and R(Z h , Z dr ) algorithms suffer from the absolute radar calibration errors and the Mie effects at X-band [56]. As K dp is a phase measurement, the R(K dp ) algorithm is often used to improve the accuracy of the rain rate estimation at X-band, yielding R(K dp ) = 17.33K 0.92 dp .
It is noted that the R(K dp ) algorithm may have a lower range resolution than R(Z h ) and R(Z h , Z dr ) due to the window smoothing. Therefore, Z h , Z dr , and K dp are combined to obtain a tradeoff between the accuracy and the resolution. The R(Z h , Z dr , K dp ) algorithm is given as Figure 1 gives a comparison between the radar rain rates and the RG data obtained by Equations (38)- (41). The scatterplots allow an easy examination of the algorithms, especially when the one-to-one line (red line) is included. The rain rates that perfectly match the RG data would fall on this line.
In Figure 1, the radar data of 1151 h contain 1940 points available for the four RGs. The majority of the points are associated with small R (light rain), followed by a sizeable number of points with medium R (moderate rain). The points with high R (heavy rain) only take a small portion but make a significant contribution to the total rain amounts. Looking at the scatterplots in more detail, we can see that R(Z h ) ( Figure 1a) gives a clear underestimation with a considerable amount of the points below the red line. As a result of the effects of the Mie scattering, negative outliers are clearly identified in the region of RG larger than 20 mm and radar R less than 10 mm. In addition, there are a small number of points above the red curve, leading to large RMSE of 2.97 and small CORR of 0.63. Moreover, the R(Z h , Z dr ) algorithm ( Figure 1b) has a similar trend to R(Z h ) due to the Mie effect, giving RMSE of 2.87, NB of 0.05 and CORR of 0.67. The R(K dp ) and R(Z h , Z dr , K dp ) algorithms improve the estimation of R when compared to R(Z h ) and R(Z h , Z dr ). The scatterplots in Figure 1c,d show a better concentration on the red line with noticeable changes for higher R. According to Figure 1c, R(K dp ) has improved the R in the bottom-right region, yielding remarkable RMSE and ρ RG at 2.22 and 0.81, respectively. The discrepancy still exist in the region of heavy rain, particularly for R larger than 20 mm. Furthermore, R(Z h , Z dr , K dp ) in Figure 1d are similar to R(K dp ) in Figure 1c, while the points with higher R concentrate on the red line. Consequently, the RMSE and CORR of R(Z h , Z dr , K dp ) are better than those of R(K dp ). In contrast, the NB of R(Z h , Z dr , K dp ) is a little worse, leading to a value of −0.05.

Retrieval of σ(R)
By using the perturbation analysis [12], the statistical uncertainty σ(R) in Equation (38) is given as σ(Z H ) = 10 log 10 where Z H = 10 log 10 (Z h ). Similarly, the variance of R(Z h , Z dr ) is expressed as σ(Z DR ) = 10 log 10 where Z DR = 10 log 10 (Z dr ). It is noticeable that σ(R) in Equations (42) and (44) are independent of the absolute values of Z h and Z dr . In contrast, K dp is in a linear scale, and thus the statistical uncertainty σ(R) is related to K dp and σ(K dp ), leading to the expression of Finally, the variance of R in Equation (41) is derived as It is known that R(K dp ), R(Z h , Z dr ), and R(Z h , Z dr , K dp ) are less sensitive to the DSD variability than the conventional R(Z h ). However, the polarimetric algorithms yield higher statistical uncertainty of R due to the uncertainties of Z dr and K dp [31]. This problem is obvious when R is small. As shown in Figure 1, R(K dp ) and R(Z h , Z dr , K dp ) give a significant improvement over R(Z h ) and R(Z h , Z dr ) for R exceeding 10 mm h −1 , whereas R(Z h ) and R(Z h , Z dr ) have a better estimation for R less than 5 mm h −1 . The explicit solution is to combine all the available algorithms by considering their statistical errors at a particular location in a certain period of time. For example, the rain rate estimation algorithms can be selected based on the reflectivity thresholds to reduce the effects of the statistical uncertainty for lower rain rates [57].
Nevertheless, it is difficult to determine the hard thresholds for various rain rate estimation algorithms, and the thresholds may be biased due to the statistical uncertainties of Z h , Z dr , and K dp , particularly for X-band radars. The rain rate algorithms can then be combined softly by adding a weight for each rain rate algorithm, leading to a composite method [58,59]. In this method, the weight is calculated by multiplying the inverse of the statistical uncertainty σ(R) by a normalized constant a i , i.e., , R(K dp ), and R(Z h , Z dr , K dp ). It can be seen that this method uses the statistical errors as a function of range gates to yield the variability of R. The previous research often analyzes the errors of Z h , Z dr , and K dp from a view of signal processing [47], leading to a constant σ(R)/R except for R(K dp ). In this study, we also consider the propagation of the statistical uncertainties of Z h , Z dr , and K dp in the rain rate estimation. To determine the coefficient a i in Equation (48), it is straightforward to set a i as the value of σ(R i ) if the corresponding algorithm has the smallest σ(R i ) weighted by the absolute values A i of the polarimetric values. Therefore, the composite method selects the rain rate with the smallest σ(R) as derived from the four algorithms: where A i is equal to Z h , Z dr , and K dp , respectively. Figure 2 compares the rain rates estimated by radar algorithms in Equations (38)- (41) and Equation (49) to the RG data. The accompanying colors show the radar algorithms of R(K dp ) (blue), R(Z h ) (cyan), R(Z h , Z dr ) (yellow), R(Z h , Z dr , K dp ) (green), and the composite method (red). Consistent with the RG data, the entire body of rain rates by the composite method is predominantly made up of light rain with only a small quantity of moderate rain and heavy rain making up the remainder. Among the five radar algorithms, the composite method yields the best performance for light rain, showing a narrow spread around the reference line. Furthermore, the composite method has a slightly higher variation for moderate rain when compared to R(K dp ), as the red dots tend to less close to the reference line than the blue dots. Nevertheless, the rain rates by the composite method well concentrate on the reference line for heavy rain. Comparison between hourly rain amounts of RG and radar rain rates. The blue, cyan, yellow, green, and red dots represent the rain rates by R(K dp ), R(Z h ), R(Z h , Z dr ), R(Z h , Z dr , K dp ), and the composite method, respectively.
When we study the associated RMSE, NB, and CORR, it is apparent that R(K dp ), R(Z h , Z dr , K dp ), and the composite method all have their advantages. The RMSE and CORR for the composite method are the best among the five algorithms, yielding improved results of 2.08 and 0.84, respectively. In particular, the RMSE for the composite has a significant improvement over R(K dp ) and R(Z h , Z dr , K dp ). However, R(K dp ) and R(Z h , Z dr , K dp ) still outperform the composite method in terms of NB.

Case Studies
The detailed analyses of individual cases are important for the accuracy assessment for the quality control and rain rate estimation for the MZZU radar. We have selected two rainfall events, including a squall line event on 7 March 2017 and a prolonged rain event on 2-4 July 2016.
corrected Z DR , (f) σ(Z DR ), (g) LRM K dp , (h) GMM K dp , and (i) σ(K dp ). The data were collected at 0411 UTC on 7 March 2017. The axes are in kilometers from the radar. Figure 3g and h both show strong convective cells (K dp > 3 deg km −1 ) with horizontal widths of a few kilometers surrounded by moderate echoes (1 < K dp < 3 deg km −1 ). Notably, GMM K dp (Figure 3h) presents a series of intense cells in the center of the line, with a maximum of greater than 8 deg km −1 . The edge of the line center appears to be jagged with forward extending protrusions. By taking a closer look at the leading edge, we can see a bow echo with a horizontal scale of 10 km, where corrected Z H is more than 60 dBZ (Figure 3b), corrected Z DR is 4-6 dB ( Figure 3e) and GMM K dp is about 2 deg km −1 . These radar signatures signify the occurrence of hail, which is consistent with the surface hail report at Columbia, Boone County (38.95 • N, 92.33 • W) at the time of 9 min later than the radar PPIs. At either side of the center, a number of weaker convective cells are aligned along the line, showing a symmetric shape of the storm. In contrast, LRM K dp (Figure 3g) tends to yield an extended region of the convective cores, but the maximum K dp values are considerably lower. When compared to GMM K dp , LRM K dp obscures the manifestations of the serrated leading edge and elongated rain cells, increasing difficulties for the analysis of the storm structure.
Moreover, the radar echoes for the trailing stratiform cannot be observed due to strong attenuation at X-band caused by the leading edge. Nevertheless, GMM K dp shows a narrow sector of the stratiform rain behind a notch-like concavity at the leading edge. The concave structure is also associated with the bow echo, which may be linked to the mid-level inflow that inserts into the back of the convective line. The trailing stratiform rain is characterized by low K dp and relatively small K dp gradients.
Likewise, LRM K dp depicts a uniformly distributed signature for the stratiform rain. The blank regions within LRM K dp are produced by missing Φ dp data, since the signal-to-noise ratio is below the noise threshold. In contrast, the Gaussian mixture method can provide meaningful K dp via the joint probability density function of Φ dp at these regions, though the statistical uncertainty σ(K dp ) tends to be increased (Figure 3i).
The statistical uncertainty σ(K dp ) gives a further insight into the radar signatures of GMM K dp . Figure 3i shows that σ(K dp ) is increased in the convective regimes due to significant K dp gradients. However, the largest σ(K dp ) is not well associated with the highest K dp value, but it occurs at a location between the leading edge and the trailing stratiform, with a large K dp gradient. It may imply that the significant σ(K dp ) is an indicator of the transition regimes [60], where the DSDs are unique. In addition, the trailing stratiform regimes have smaller σ(K dp ) compared to the leading convection due to relatively stable K dp values.
As illustrated in Figure 3a,d, raw Z H and Z DR are affected by the rain attenuation. In Figure 3a, the rear part of the leading edge shows raw Z H of 20-30 dBZ, weaker than the typical signatures of convective rain. Moreover, the attenuation effects lead to a Z H gap between the leading edge and the rain cell at the far range. On the other hand, raw Z DR (Figure 3d) is unrealistically lower than −1 dB after passing the front part of the leading edge. As K dp is immune to the attenuation, it is applicable for the attenuation correction via Equations (29) and (30). Figure 3b,e depict that corrected Z H and Z DR in the leading edge have been intensified by about 10-20 dBZ and 1-3 dB, respectively, while the gap in the trailing part is transformed to a region of stratiform rain. Moreover, it is notable that the strong attenuation produces a "V"-shape echo in the hail region at X: −20~−15 km and Y: 15~20 km. The high K dp values have given significant compensation to corrected Z H and Z DR , leading to intense echoes in this region. However, the correction using GMM K dp has the radial effects, particularly for the regions with missing values in raw Z H and Z DR . It may be due to the K dp smoothing and the constant α and β in Equations (29) and (30), respectively.
The statistical uncertainties σ(Z H ) and σ(Z DR ) inherit σ(K dp ) via Equations (31) and (32), and therefore σ(Z H ) and σ(Z DR ) for corrected Z H and Z DR are larger than those for raw data. According to Figure 3c,f, σ(Z H ) and σ(Z DR ) are significant at the regions of high K dp gradients, as σ(Z H ) and σ(Z DR ) are positively related to σ(K dp ). In Figure 3c, it can be seen that σ(Z H ) is negligible relative to the magnitude of corrected Z H . When σ(Z H ) reaches the maximum, it is less than 1% of the corresponding Z H magnitude. In contrast, σ(Z DR ) is on the same order as the magnitude of corrected Z DR . At the point with the largest σ(Z DR ), the ratio of σ(Z DR ) and Z DR can reach more than 20%, which is much higher than the upper limit for Z DR error. Overall, the attenuation correction using GMM K dp introduces relatively large errors for Z DR , leading to higher uncertainty when corrected Z DR is used to retrieve the rain rate.
To give a further evaluation of the attenuation correction, the corrected Z H and Z DR are compared to the self-consistency relations of Z H , Z DR , and K dp [12,21,23,53,61]. Self-consistency relations are related to the drop size distribution [62]. In this study, we adopt the Pruppacher-Beard linear fit [63], yielding the following relations for the X-band radar: K dp = 1.37 × 10 −3 × 10 0.068Z H × 10 −0.042Z DR , where Z H = 10 log Z h and Z h is in mm 6 m −3 . Figure 4 compares the scatterplots of corrected Z H , corrected Z DR and GMM K dp to the self-consistency relations in Equations (50)- (52). As depicted in Figure 4a, K dp increases greatly when raw Z H is above 30 dBZ, showing a large difference from the self-consistency relations. In contrast, in Figure 4b, there is a low and steady increase in K dp when Z H is between 10 and 40 dBZ, following a significant K dp jump when Z H is larger than 40 dBZ. The most noteworthy increase is recorded for Z H greater than 50 dBZ, where K dp rises dramatically from 4 to 8 deg km −1 . Overall, the Z H -K dp points after the correction have good agreement with the self-consistency relations. Figure 4. Scatterplots of (a) raw Z H and K dp , (b) corrected Z H and K dp , (c) raw Z DR and K dp , (d) corrected Z DR and K dp , (e) raw Z DR and the ratio of K dp and raw Z h , and (f) corrected Z DR and the ratio of K dp and corrected Z h for the data presented in Figure 3, where Z H = 10 log(Z h ). The color scale is the number of points, and the black lines show the theoretical self-consistency relations.
In Figure 4c,e, the points of Z H -Z DR and Z DR -K dp /Z h can be separated into two data clusters with high populations, which remain after the attenuation correction in Figure 4d,f. As shown in Figure 4d, the cluster concentrated on 2 dB has a steady increasing trend, which agrees with the self-consistency relations. However, the higher slope of the cluster indicates that the coefficient β in Equation (30) is not optimal for the MZZU radar. On the other hand, the cluster centered at Z H of 50 dBZ and Z DR of 4 dBZ are likely produced by a mixture of rain and hail, as melting hail or raindrops with ice cores make a significant contribution to Z DR . In addition to the Z H -K dp and Z H -Z DR relations, it is evident that the ratio of K dp and Z h also shows a dependency on Z DR . As Z DR rises, K dp /Z h tends to drop steadily, following a similar trend with the scattering simulation [54].
Through Equations (38)-(49), we can retrieve the rain rates and their uncertainties. As might be expected, R(Z h ) (Figure 5a) and R(Z h , Z dr ) (Figure 5b) give an underestimation at the center of the leading edge, with the maximum rain rates of no more than 40 mm h −1 . In contrast, R(K dp ) (Figure 5c) and R(Z h , Z dr , K dp ) (Figure 5d) produce the values as high as 100 mm h −1 , indicating intensive rain cells. What remains truly remarkable is that K dp -based algorithms yield the meaningful rain rate estimation in the hail-contaminated regions at X-band. Somewhat surprisingly however, according to Figure 5b-d, R(Z h , Z dr ), R(K dp ), and R(Z h , Z dr , K dp ) present some missing data in the trailing stratiform due to negative values in Z DR and K dp .  According to Figure 5e, the composite method shows the rain rate results similar to R(K dp ) in the leading edge. Likewise, the statistical uncertainty (Figure 5f) for the composite method has given higher values in the transition regime, which is consistent with the uncertainty of Z H , Z DR , and K dp . However, it is clear the trailing stratiform in the composite method is characterized by lower R, giving the improvement of rain rate estimation in this region.
It is found that the elongated shapes in the convective cores may be resulted from the radial effect of R(K dp ) as the K dp is estimated independently for each ray. Furthermore, the K dp is applied to the attenuation correction along the radial, and thus R(Z h , Z dr , K dp ) gives the discontinuity in the azimuthal direction. However, the strong vertical winds have a large contribution to the raindrop size distribution as well as fall speed, leading to uneven rain in space. The updrafts may lift the smaller raindrops and enhance the coalescence process during the falling. Therefore, the rain rates increase when the updrafts exist in the area. On the other hand, the radar ray may spread after a number of range gates, while the radar variables and the resulting rain rates are averaged over the volume. If the small convective cells present in a few adjacent radar gates, the rain rates within these gates may be higher than the surrounding gates, giving the discontinuity in space.
In conclusion, the quality control is effective for Z H , Z DR , and K dp in the rain regions of the squall line event. The data of GMM K dp , corrected Z H and Z DR can provide an insight into the storm structure, and the scatterplots are consistent with the self-consistency relations at X-band. Furthermore, it is evident that R(K dp ) and R(Z h , Z dr , K dp ) have better results than R(Z h ) and R(Z h , Z dr ) at X-band, and the composite method gives some improvements in the trailing stratiform regime.

Prolonged Rain Event: 2-4 July 2016
From 2 July 2016, a weak lee-side trough formed at the southern Rocky Mountains and moved toward Missouri. It persisted over central Missouri and produced large areas of precipitation. Meanwhile, a cold front trailing from the low advanced southeastward across the radar surveillance area, causing convective rain and strong winds. Figure 6 illustrates the radar PPIs at 0339 UTC on 3 July 2016, when the convection embedded in the stratiform is approaching. By briefly glancing at the PPIs for this event, the radar intensity is considerably weaker than that for the squall line event in Figure 2. According to the data in Figure 6e,h, corrected Z DR and GMM K dp are characterized by low values in the widespread stratiform regime, indicating small and less-oblate raindrops. Moreover, the evaporation below the melting layer may further decrease K dp as a result of mass flux changes [64]. In the west and northwest of Figure 6h, GMM K dp shows some pockets of convective cells with moderately large values (>2 deg km −1 ). As the updrafts in the convection are narrow, the convective cells are manifested as a small number of adjacent range gates of intensified K dp . Notably, a small number of range gates are discarded at the edge due to a low signal-to-noise ratio. Nevertheless, GMM K dp has reduced the effects of the missing data when compared to LRM K dp (Figure 6g), leading to continuous regions of precipitation. Furthermore, GMM K dp tends to produce larger but fewer pockets of small and negative values, which gives a slight improvement over LRM K dp . The ring echo of negative K dp near the radar center results from the switching of short and long pulse repetition frequencies (PRFs). In addition, stable K dp yields negligible σ(K dp ) (Figure 6i) in the stratiform regime, while negative K dp (<−1 deg km −1 ) increases σ(K dp ) at northeast and south of the PPI. Figure 6b and a compare corrected Z H to raw Z H for the stratiform event. The major changes for Z H occur in the west and northwest of the PPI, where Z H has been increased by 10-25 dBZ. In the raw data (Figure 6a), the convective cells appear to be a number of isolated pockets with maximum Z H of no more than 45 dBZ. After the correction (Figure 6b), there is an extended region of Z H larger than 45 dBZ, while the cores of the cells reach about 60 dBZ. It is noteworthy that Z H at X: −80~−70 km and Y: 0~20 km have been significantly intensified as a result of high K dp . Moreover, the uncertainty Figure 6c is very similar to the constant error of 1.36 dB. In the south, σ(Z H ) has reached about 1.8 dB due to increased σ(K dp ).
Likewise, Figure 6d and e show Z DR before and after the correction, respectively. The main difference between the two data is that raw Z DR (Figure 6d) has small and negative values in the sector from west to north, whereas corrected Z DR (Figure 6e) shows the values larger than 1 dB in this region. Looking at the data in more detail, the attenuation correction has also enhanced Z DR in the stratiform region, increasing from −1 to 0.5 dB. Figure 6f shows that σ(Z DR ) is very small, except for the echoes at the edge of the southern part.
Moreover, Figure 7 shows the scatterplots of Z H , Z DR , and K dp . According to the data in Figure 7a,b, we can discern a clear trend in the relationship between Z H and K dp . When Z H is below 35 dBZ, K dp has a slow and steady increase, ranging from near zero to a few tenths' deg km −1 . The specific differential phase then gives a dramatic rise when Z H is between 35 and 45 dBZ. The high probability falls in the region of Z H between 20 and 38 dBZ and K dp between 0 and 0.5 deg km −1 for the stratiform rain. Furthermore, it can be seen that the attenuation correction has intensified Z H , giving better agreement with the reference curve when K dp is above 1.7 deg km −1 . Although K dp can reach as low as −0.5 deg km −1 due to the PRF switching and phase fluctuation, the negative K dp has a relatively low number concentration.
The scatterplots in Figure 7c,d give a positive correlation between Z H and Z DR , which is consistent with the reference curve. In contrast with the squall line event, the points in the stratiform event concentrate on one rain cluster with Z DR between 0 and 2 dB. By comparing corrected Z DR (Figure 7d) to raw Z DR (Figure 7c), it can be found that the attenuation correction has reduced the number of points with negative Z DR and also widened the distribution of corrected Z DR . These effects are apparent in the scatterplots of Z DR and K dp /Z h as illustrated in Figure 7e,f. Moreover, corrected Z DR has been shifted toward the positive direction, indicating more oblate raindrops. However, there is a small mismatch between the high probability regions of K dp /Z h in Figure 7f and the reference curve due to low coefficient α. In addition to Z H , Z DR , and K dp , Figure 8 depicts the PPIs of the rain rates retrieved by R(Z h ), R(Z h , Z dr ), R(K dp ), and R(Z h , Z dr , K dp ) and the composite method for the stratiform rain event. It is clear that R(K dp ) (Figure 8c) has a similar distribution to K dp (Figure 6h), as R(K dp ) is almost linearly proportional to K dp . The principal improvement of the composite method (Figure 8e) is that it has filled gaps within the regions where R(K dp ) presents the missing values, and therefore the rain accumulations are increased in these regions. By taking a closer look at the embedded convection, we can see that R(K dp ) (Figure 8c), R(Z h , Z dr , K dp ) (Figure 8d) and the composite method (Figure 8e) produce values as high as 60 mm h −1 , whereas R(Z h ) (Figure 8a) and R(Z h , Z dr ) (Figure 8b) are generally below 50 mm h −1 . On the other hand, the composite method can present uniformly distributed rain in the stratiform region, which is characterized by σ(R) between 0.5 and 5 mm h −1 as shown in Figure 8f. When compared to the squall line event, σ(R) is generally small for this stratiform event. In addition, the negative K dp is not used for attenuation correction or rain rate estimation in this paper. However, the negative values are frequently shown in the stratiform rain due to the fluctuation of phase measurements. To reduce this effect, we have set the rain rates with negative K dp to zero, leading to a continuous region of rain.
Indeed, the discontinuity of the attenuation correction and rain rates by K dp is inevitable. In Figure 8, we present the instantaneous data to give the best comparison between the composite method and the power-law methods. When comparing with the rain gauge data, we take an average over the area of 3 gates × 3 azimuths. There are a number of sophisticated solutions for rain rate estimation, such as the areal rain rate estimation using negative K dp [65] and X-band radar network [66].
To validate the algorithms for the rain rate estimation, we also present the time series of the prolonged rain event that occurred on 2-4 July 2016. Figure 9 depicts an example of hourly accumulated rain amounts using Equations (38)- (41) and (49). This is a portion of the time series at Bradford (Figure 9a  . Time series of hourly rain amounts of RG (black), R(K dp ) (blue), R(Z h ) (cyan), R(Z h , Z dr ) (yellow), R(Z h , Z dr , K dp ) (green) and the composite method (red) at (a) Bradford, (b) Sanborn, (c) Auxvasse and (d) Williamsburg. The error bars represent the statistical uncertainties associated with the composite method. The data were collected between 0000 UTC and 2300 UTC on 3 July 2016.
At Bradford (Figure 9a), the rain rates remain steady during the first five hours, followed by a sudden rise at 06 UTC before a slight decline at 07 UTC. The time series of the rain rates arrive at a peak at 08 UTC, and then fall sharply over the following five hours. There is a slight increase up to a few millimeters when the rain rates drop steadily and finally reach the lowest point at 18 UTC. Overall, the radar algorithms for the rain rate estimation have given an underestimation when compared to the RG data. As the magnitudes of the rain rates go up, the gap between the radar and RG becomes wider, reaching about 10 mm at the peak. By observing the individual PPIs, we can find that the radar data may be affected by strong signal attenuation. However, the attenuation correction cannot sufficiently compensate the radar signal at some time periods, producing smaller total amounts when the rain rates are accumulated over an hour. Moreover, the composite method (red) yields the best agreement with the RG data among the five algorithms, but its statistical uncertainties (error bars) are considerably large. In contrast, R(Z h ) and R(Z h , Z dr ) present very small values at the peaks, and relatively larger values for the rest of the time interval. In addition, R(Z h ) and R(Z h , Z dr ) have very low correlations with RG, indicating a poor performance. Figure 9b,c show the time series of the rain rates at Sanborn and Auxvasse, respectively. According to the graph, the trends for the rain rates at both sites are very similar, and strong correlations with RG can be found. It can be seen that the rain rates over the first three hours remain stable before rising steadily between 03 and 05 UTC. Remarkably, however, the rain rates shoot up to 40 mm at Sanborn and 25 mm at Auxvasse at 07 UTC. As then the rain rates decrease dramatically, and remain at a few millimeters until 19 UTC. Moreover, the time series at Sanborn and Auxvasse show clear gaps between the radar and RG, but the gaps have been reduced when compared to that at Bradford. Among the five algorithms, the composite method still gives the best performance in terms of RMSE, NB, and CORR. Notably, the rain rates produced by the composite method at the peaks of Sanborn and Auxvasse are higher than R(K dp ) and R(Z h , Z dr , K dp ), yielding better agreement with the RG. In contract, R(Z h ) and R(Z h , Z dr ) show a moderate underestimation at this point. Overall, R(K dp ), R(Z h , Z dr , K dp ), and the composite method is consistent with the RG for higher rain rates, as K dp is more sensitive to heavy rain. On the other hand, R(Z h ), R(Z h , Z dr ), and the composite method have good agreement for lower rain rates, but R(Z h ) and R(Z h , Z dr ) produce significant bias at the plateaus, leading to considerable RMSE, NB, and CORR.
As illustrated in Figure 9d, the time series of the rain rates at Williamsburg have weak correlations with the RG. It is clear that the time series exhibit three local peaks at 06, 10, and 18 UTC, respectively. However, the first peak of the rain rates is found at a time one hour later than the RG though the magnitude at this time is consistent with the RG. At the second and third peaks of RG, the rain rates generally produce much smaller values than the RG with relatively stable trends. These differences may be due to the radar beam spread at the far range. As Williamsburg is 46.2 km away from the radar center, the width of the range gate reaches about 700 m, covering an area of 700 × 260 m 2 . In contrast, the measurement area of a tipping bucket gauge is only a few hundred cm 2 . The difference between the areal rainfall and the point rainfall leads to significant bias and uncertainty, which can reduce the accuracy of the radar rain rate estimation [67]. Furthermore, referring to the timing of the peaks at Sanborn and Bradford, it can be found that the rain cell passed the radar center between 06 and 08 UTC, which may have coincided with heavy rain at Williamsburg. However, the signal attenuation at X-band reduces the magnitude of the measurements at Williamsburg at this time, although K dp may not very sensitive to the attenuation. In addition, the peak may be shifted to the subsequent hour if the heaviest rain is right at 07 UTC since the rain gauges collect the samples every hour.
To give a further assessment, Figure 9 illustrates the statistical uncertainties σ(R) (error bars) associated with the composite method. By briefly glancing at the time series, it is apparent that there are considerable variations in σ(R)s at the four sites between 02 and 19 UTC. At Bradford (Figure 9a), σ(R) remains large in most of the times, and then gives a moderate increase at 06 UTC, reaching about 25 mm. Standing a similar trend, σ(R) at Sanborn (Figure 9b) peaks at the same time, but it has even larger magnitudes when compared to Bradford. Furthermore, the average σ(R) at Sanborn is generally higher than the others. In contrast, σ(R) at Auxvasse exhibit narrow plateaus at 06 UTC due to the steady trend at this time. In addition, there is a steady trend for σ(R)s at Williamsburg (Figure 9d), with an average of a few tenths of a millimeter.
If we look more closely at Figure 9, we can discern that the largest σ(R) often occurs with the high variations of the rain rates. During the periods of 05 to 09 UTC, all the four sites present one or more peaks, leading to larger σ(R) correspondingly. However, the increased σ(R) is not always related to the peaks of the rain rates. For example, σ(R) at Sanborn at 06 UTC is about three times larger than the peak at 07 UTC.
In conclusion, the quality control can effectively retrieve K dp and correct the attenuations in Z H and Z DR for the prolonged rain event on 2-4 July 2016. The comparison between the scatterplots and the self-consistency relations confirms that the attenuation correction yields a good performance. Furthermore, the radar algorithms for the rain rate estimation have been derived from the MZZU-RG dataset, yielding a fairly good agreement with the RG for the prolonged rain event. The major improvement of the composite method is that it can fill the missing data in the stratiform regions, and provide the lower statistical uncertainties for the rain rates. From the time series data, it is clear that σ(R) is often associated with the variation of the rain rates rather than its magnitude.

Conclusions
In this study, the quality control and rain rate estimation were investigated using the MZZU X-band dual-polarization radar in central Missouri. For quality control, a simple data masking technique was adopted to remove the clutter data in the radar reflectivity (Z H ), differential reflectivity (Z DR ) and differential phase shifts (Φ dp ). The Gaussian mixture method (GMM) was then applied to the specific differential phase (K dp ) estimation, which can not only retrieve the expected value of K dp but also calculate the statistical uncertainty σ(K dp ).
For the attenuation correction, the data of K dp and σ(K dp ) were used to correct Z H and Z DR . The Φ dp -based method was adopted with the assumption that the specific attenuation (A H ) and specific differential attenuation (A DP ) were linearly proportional to K dp . Meanwhile, the statistical uncertainties σ(Z H ) and σ(Z DR ) were propagated from σ(K dp ), leading to the variability in their error characteristics. The corrected Z H , Z DR , and K dp were finally investigated in the rain rate (R) estimation.
For the rain rate estimation, the radar rain rates were retrieved by Z H , Z DR , and K dp via the power-law relations, while the statistical uncertainty σ(R) was propagated from σ(Z h ), σ(Z dr ), and σ(K dp ). The coefficients in R(K dp ), R(Z h ), R(Z h , Z dr ), and R(Z h , Z dr , K dp ) were optimized with the comparison data of the MZZU radar and rain gauge (RG). It was clear that R(Z h ) and R(Z h , Z dr ) gave an underestimation for higher rain rates, but R(K dp ) and R(Z h , Z dr , K dp ) improved the estimation of higher rain rates, giving remarkable root mean square error (RMSE), normalized bias (NB) and Pearson correlation coefficient (CORR). Furthermore, the statistical uncertainty σ(R) was approximated by a first order Taylor series expansion of the power-law relations. To reduce the effects of σ(R), the composite method was proposed by combining the existing rain rate algorithms based on the inverse of their σ(R). It was found that the composite method yielded the best performance in terms of RMSE and CORR, but showed slightly worse NB when compared to R(K dp ) and R(Z h , Z dr , K dp ).
Two rainfall cases were evaluated, including a squall line event on 7 March 2017 and a prolonged rain event on 2-4 July 2016. In the squall line event, the serrated leading edge, the elongated convective cells and the trailing stratiform were manifested in the K dp , giving a clear picture of the storm structure. It was notable that σ(K dp ) was increased with the K dp gradient, while the significant σ(K dp ) indicated the transition region between the regimes of the leading convection and the trailing stratiform. The radar reflectivity and differential reflectivity were intensified as a result of the attenuation corrections, leading to realistic radar signatures. It was found that σ(Z H ) was very small relative to corrected Z H , whereas σ(Z DR ) was significant relative to corrected Z DR . Furthermore, the scatterplots of corrected Z H , corrected Z DR and K dp had fairly good agreement with the self-consistency relations at X-band, indicating the effectiveness of the attenuation corrections. The rain rate retrieved by the composite method was compared to R(K dp ), R(Z h ), R(Z h , Z dr ), and R(Z h , Z dr , K dp ). It was clear that the composite method produced the meaningful rain rates in the regions with missing data, giving the best performance among the five rain rain algorithms. In the composite method, the statistical uncertainty σ(R) was fairly high in the transition regime but relatively small in the leading edge and the trailing stratiform.
The prolonged rain event was dominated by the widespread rain, showing weaker radar echoes than the previous case. In this event, the K dp was characterized by low values in the stratiform regime and by moderately large values in the convective regime, due to differences in the mean vertical velocity and particle growth mechanisms. It was notable that the GMM K dp reduced the effects of the missing values and also slightly improved negative K dp . The statistical uncertainty σ(K dp ) was negligible, except for the negative K dp regions at the northeastern and southern parts. By applying the K dp to the attenuation corrections, Z H and Z DR were given a considerable rise in the convective regime, while Z DR was also enhanced in the stratiform regime, yielding realistically positive values. The statistical uncertainties σ(Z H ) and σ(Z DR ) were very close to the constant errors, but they were increased in the negative K dp regions. To further validate the attenuation corrections, the scatterplots of Z H , Z DR , and K dp were compared to the reference curves obtained by the self-consistency relations at X-band. The high probability regions were consistent with the reference curves, indicating reliable correction results for the prolonged rain event. Moreover, Z H , Z DR , and K dp were applied to the retrievals of the rain rates via R(K dp ), R(Z h ), R(Z h , Z dr ), and R(Z h , Z dr , K dp ). It was clear that the composite method filled the missing data as presented in R(Z h , Z dr ), R(K dp ), and R(Z h , Z dr , K dp ), and also provided the smallest σ(R). In the case study of a prolonged rain event, clear gaps were shown in the time series of the four RG sites within the radar coverage. Nevertheless, strong correlations were found between the radar rain rates and RG, indicating a fairly good performance. Among the five algorithms, the composite method had the best agreement with RG. Furthermore, the statistical uncertainty σ(R) was considerable in the time series, and σ(R) was related to the variation of the rain rate rather than its magnitude.
It was apparent that the relations of A H − K dp and A DP − K dp were related to the raindrop size distribution, shape and temperature, leading to moderate bias in the attenuation corrections if constant coefficients α and β were used. In the future, the self-consistent ZPHI method will be applied to the attenuation corrections for the MZZU radar, following the work of [21,22]. The statistical uncertainties σ(A H ) and σ(A DP ) will also be derived from σ(Φ dp ) via the Gaussian mixture model. For rain rate estimation, the R(A H ) and R(A V ) algorithms will be derived to retrieve the rain rate [68], while the statistical uncertainties σ(R AH ) and σ(R AV ) will also be calculated by σ(A H ) and σ(A V ), respectively. Some discussions have been provided based on the Z-PHI method in Appendix A.
Author Contributions: N.I.F. designed the experiment, provided the radar data and reviewed the paper, G.W. developed the methods and prepared the manuscript, and P.S.M. reviewed the paper. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors would express the appreciations for constructive comments from two anonymous reviewers, under the effective and efficient management of the editor Joan Bech and associate editor Mladen Radosavljević. The MZZU radar data can be made available upon request to the authors. The rain gauge data are available online: http://agebb.missouri.edu/weather/stations/index.php.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Attenuation-Based Rain Rate Estimation
From the examples in Section 5, it can be seen that the attenuation corrections of Z H and Z DR using GMM K dp introduce moderate bias when compared to the self-consistency relations in Equations (50)- (52). This is because the relations of A H -K dp and A DP -K dp are subject to drop shape, DSD, and temperature. First, the coefficients α and β in Equations (25) and (26) are assumed to be constants, but α and β are closely related to the axis ratio [12]. The change of axis ratio may lead to the variations of α and β at each range gate. Secondly, K dp is less sensitive to raindrops with spherical shapes (D < 0.5 mm), whereas A H and A DP are related to the liquid water content contributed by all sizes of drops [69]. Furthermore, for the median diameter D 0 larger than 2.5 mm, A H and A DP are proportional to D 0 × K dp , and therefore α and β increase as D 0 rises. If raindrops with D 0 larger than 2.5 mm exist along the path, the attenuation corrections based on constant α and β may not be reliable. Thirdly, the scattering simulations which derive α and β depend on the ambient temperature [70]. As given in [17], α and β under various temperatures fall into a band rather than a straight-line. For accurate attenuation correction, α and β cannot be assumed as constants, but derived from a range of values corresponding to the band.
To tackle these problems, [21] extended the methods of [19,20], which use the increase of Φ dp to calculate the range integrals of A H and A DP . They also optimize α(r) by minimizing the difference between the measured Φ dp and the Φ dp reconstructed from A H , and then obtain β(r) following the self-consistency relation of Z H and Z DR . This method is later adapted to the X-band [22]. In addition, the specific attenuation A H at each range gate can be used to retrieve the rain rate R via a power-law relation between A H and R.