An Inverse Model for Raindrop Size Distribution Retrieval with Polarimetric Variables

Abstract

This paper proposes an inverse model for raindrop size distribution (DSD) retrieval with polarimetric radar variables. In this method, a forward operator is first developed based on the simulations of monodisperse raindrops using a T-matrix method, and then approximated with a polynomial function to generate a pseudo training dataset by considering the maximum drop diameter in a truncated Gamma model for DSD. With the pseudo training data, a nearest-neighborhood method is optimized in terms of mass-weighted diameter and liquid water content. Finally, the inverse model is evaluated with simulated and real radar data, both of which yield better agreement with disdrometer observations compared to the existing Bayesian approach. In addition, the rainfall rate derived from the DSD by the inverse model is also improved when compared to the methods using the power-law relations.

1. Introduction

The raindrop size distribution (DSD) is one of the most important characteristics of a precipitating process, since the raindrop formation and its size evolution are related to the mechanism of cloud microphysics, kinetics, and thermodynamics. Weather radar is suited to obtain backscattering properties of the raindrops on mesoscale and convective scales when the raindrops form and descend toward the surface. The correct retrieval of DSD by radar is valuable for documenting and understanding the microphysical properties of a cloud system, increasing the accuracy of quantitative precipitation estimation, and improving the parameterization of numerical weather prediction. Potential applications include the flash flood warning, water resource management, reservoir construction, aviation, and weather modification.

The polarimetric weather radar has a long and profound impact on the remote sensing of atmosphere. Seliga and Bringi (1976) [1] first indicated that the ratio of the radar reflectivity at the horizontal and vertical polarizations, i.e., differential reflectivity, provided an additional evidence of distorted raindrops. Later, they [2] showed that the polarimetric measurements are in good agreement with model calculations of raindrops with a T-matrix method at S-band frequency. This technology was also transferred to attenuating frequencies, such as C-band [3] and X-band [4]. As the polarimetric radar transmits and receives the electromagnetic waves at two orthogonal polarizations, it also measures polarimetric variables, including differential reflectivity, differential phase shift, specific differential phase, and cross-correlation coefficient, when compared to the conventional radar [5,6]. By combining with absolute reflectivity, these variables can be used to determine the parameters in the truncated Gamma DSD [7,8].

In the literature, the algorithms for the DSD retrieval can be divided into two categories: One is a physical method and the other is a statistical approach. The physical method explores the direct relationship between the DSD parameters and polarimetric measurements. For example, the constraint method [9,10,11] uses radar reflectivity and differential reflectivity to retrieve the Gamma distribution by assuming an empirical relation between the slope and shape parameters. Meanwhile, the beta method [12] derives the effective slope of the axis ratio with respect to the diameter to make a better use of specific differential phase. This method also utilizes a normalized Gamma distribution [13] instead of the nonnormalized one [7], yielding a more realistic DSD [14]. However, it is evident that the constraint method has better consistency with the disdrometer data when compared to the beta method, since the beta method is sensitive to the noise of specific differential phase [15,16]. To reduce the noise effects, a power-law relation was employed for light rain in the beta method [17,18]. Moreover, the relations between the DSD parameters and polarimetric variables can be derived by using a reflectivity-weighted mean diameter [19,20]. This algorithm also considers the Mie effects at X-band frequency, to correct the rain-path attenuation. In addition, Raupach and Berne [21] retrieved the third and sixth moments to obtain a double-moment model for normalized DSD.

In contrast, the statistical approach approximates the nonlinear relation between the Gamma parameters and polarimetric variables with a mathematical process. Vulpiani et al. [22] proposed a regularized neural network method to give a non-parametric mapping from the polarimetric variables to the DSD parameters. Cao et al. [23] developed a Bayesian regression model by assuming the samples of radar measurements are independent and approaching to a Gaussian function. The Bayesian framework consists of a priori and a probability density function. If the probability density function is transferred to a penalty function, a variational method [24,25] can be used to obtain the associated parameters via an optimization procedure.

This study proposes a statistical approach to estimating the parameters of a truncated Gamma model for DSD with the polarimetric radar variables, called an inverse model. This approach designs a forward operator to calculate the polarimetric variables for a given DSD by assuming the shape and orientation of raindrops, ambient temperature, and radar wavelength. Furthermore, it approximates the forward operator with a polynomial function to generate a pseudo training dataset, which is inversely used to obtain the relation between the DSD and polarimetric variables. Finally, it adopts a non-parametric estimator based on a nearest-neighborhood method to retrieve the parameters within the truncated Gamma model. By comparing to the existing algorithms, there are some advantages of the inverse model. First, it uses pseudo training data generated by a forward operator to retrieve the DSD parameters inversely. These pseudo data are adjustable according to the dynamic range of the DSD parameters. Secondly, it can obtain the maximum drop diameter to yield a realistic Gamma model for DSD, whereas the existing algorithms generally assume the maximum drop diameter as a deterministic value or infinity. Thirdly, the model is optimized with regard to the mass-weighted mean diameter and liquid water content, since they are important for the Gamma representation of the DSD.

The paper is organized as follows: Section 2 presents the mathematical representations of the DSD and the microphysical quantities calculated by the DSD. It also introduces the dataset used in this study. Section 3 gives the forward operator for polarimetric variables and its polynomial approximation, which is used to generate a pseudo training dataset. In Section 4, a non-parametric estimator based on a nearest-neighborhood method for the DSD retrieval is proposed and its error characteristics are described. In Section 5, the performance of the proposed approach is evaluated by comparing to the existing Bayesian approach, while the rainfall rate derived from the DSD is also examined against the power-law relations. Finally, Section 6 provides a summary and some further discussions on the proposed method.

2. Representation of Raindrop Size Distribution

The raindrop size distribution is defined as the number density of raindrops per unit size range and per unit volume [26]. It plays a key role in the microphysics and dynamics of raindrops while falling in the atmosphere. Some instruments (e.g., disdrometer) can provide direct measurements of DSD by detecting individual drops with optical light illumination [27,28], whereas the measurements of remote sensing methods (e.g., radar) are indirect. That means the number density at a particular size bin cannot be obtained directly from radar measurements compared to the disdrometer. In this case, a parametric function has to be selected and transferred to represent the DSD of interest. A variety of theoretical distributions could be used for the DSD representation, including Gamma distribution [7,13], log-normal distribution [29], and exponential distribution [30], among others (e.g., [31]). In the literature, the most commonly used one is Gamma distribution, which is described in the following section.

2.1. Gamma Model

The Gamma function [7] is widely used to model the DSD in the precipitating system. It is formulated as

$$N(D)=N_0{D}^{\mu }\exp (-\Lambda D),$$

(1)

where $N(D)$ is the size spectrum in mm⁻¹ m⁻³, N₀ is an intercept parameter in ${\mathrm{mm}}^{-1-\mathsf{\mu }}$m⁻³, $\mu$ is a shape parameter in unitless, and $\Lambda$ is a slope parameter in mm⁻¹. This definition is convenient since the corresponding moment (M_n in mmⁿ m⁻³) of this distribution yields an explicit function other than an implicit integral form:

$$M_n={\displaystyle \int D^{n}N(D)dD=N_0\frac{\Gamma (\mu +n+1)}{{\Lambda }^{\mu +n+1}}}.$$

(2)

It is worth noting that Equation (1) represents the spectrum with a drop size from zero to infinity, thus the DSD moment in Equation (2) becomes an infinite integral accordingly. However, in nature, the measurements are often conducted within a limited sampling volume and time, and the resulting DSD is truncated in a finite size range $\left(D_{\min },D_{\max }\right)$. For example, the two-dimensional video disdrometer (2DVD) can measure drops with a diameter between 0.1 and 8 mm [32]. To reduce the number of parameters in the DSD model, the lower limit is often set to be zero [7], whereas an exception is an unusual DSD associated with a size-sorting effect due to updrafts and wind shear [33]. The Gamma model is then changed to [7]

$$N(D)=N_0{D}^{\mu }\exp (-\Lambda D)\left(0\le D\le D_{\max }\right),$$

(3)

and the corresponding moment is calculated as [8]

$$M_n=\frac{N_0}{{\Lambda }^{\mu +n+1}}\gamma \left(\mu +n+1,\Lambda D_{\max }\right),$$

(4)

where $\gamma (·)$ is a lower incomplete Gamma function.

In Equation (1), the unit of the intercept parameter N₀ varies according to the shape parameter $\mu$. To obtain an independent intercept term, a normalized Gamma function is introduced by scaling the DSD in terms of liquid water content and mass-weighted mean diameter [13,34]. It is formulated as

$$N(D)=N_{w}f(\mu ){\left(\frac{D}{D_m}\right)}^{\mu }\exp \left[-(4+\mu )\frac{D}{D_m}\right],$$

(5)

where N_w is the normalized intercept term in mm⁻¹ m⁻³, D_m is the mass-weighted mean diameter in mm, and $f(\mu )$ is a normalization factor defined in [5]. The new intercept term N_w can be written as

$$N_w=\frac{4^4}{\pi {\rho }_w}\left(\frac{{10}^{3}W}{D_m^4}\right),$$

(6)

where ${\rho }_w$ is the density of liquid water, and W is the liquid water content. The new term N_w is equivalent to the intercept of an exponential distribution with an identical liquid water content [5]. In addition, the slope parameter $\Lambda$ is substituted by an approximation expression:

$$\Lambda D_m=4+\mu ,$$

(7)

which assumes the maximum diameter D_max is sufficiently larger than the mass-weighted mean diameter D_m (D_max/D_m$\ge$2.5 [7,35]). The corresponding DSD moments can then be expressed as

$$M_n=N_{w}F(\mu )D_m^{n+1},$$

(8)

with

$$F(\mu )=\frac{6}{{(4)}^4}{(4+\mu )}^{3-n}\frac{\Gamma (\mu +n+1)}{\Gamma (\mu +4)}.$$

(9)

Another normalization method, called a double-moment generalized Gamma model [21,36,37], takes into account two random DSD moments:

$$N(D)={\left(M_i\right)}^{(j+1)/(j-i)}{\left(M_j\right)}^{(i+1)/(i-j)}h(x),$$

(10)

with

$$x=D{\left(M_i\right)}^{1/(j-i)}{\left(M_j\right)}^{-1/(j-i)}.$$

(11)

Further, the general moment can be derived from the relationship among the DSD moments:

$$M_n=C_n{\left(M_j\right)}^{(n-i)/(j-i)}{\left(M_i\right)}^{(j-n)/(j-i)},$$

(12)

with

$$C_n={\displaystyle \int h(x)x^n}dx.$$

(13)

2.2. Bulk Properties of DSD

The liquid water content W is calculated by using the third moment of DSD:

$$W=\frac{\pi }{6}\times {10}^{-3}{M}_3(\mathrm{g}{\mathrm{m}}^{-3}),$$

(14)

where the water density is considered as 1 g cm⁻³. The median volume diameter, D₀ in mm, is the upper limit of raindrop size that contributes to half of the total liquid water content [5]:

$${\displaystyle \underset{0}{\overset{D_0}{\int }}{D}^{3}N(D)dD=\frac{1}{2}}{\displaystyle \underset{0}{\overset{\infty }{\int }}{D}^{3}N(D)dD}.$$

(15)

Meanwhile, the mass-weighted mean diameter (D_m) is defined as the ratio of the fourth and third moments:

$$D_m=\frac{M_4}{M_3}(\mathrm{mm}).$$

(16)

In the still air, the rainfall rate is determined by the combined effects of drop fall velocity and DSD, which is expressed as

$$R=6\pi \times {10}^{-4}{\displaystyle \int v(D)D^{3}N(D)dD}(\mathrm{mm}{\mathrm{h}}^{-1}).$$

(17)

The terminal velocity at the surface can be approximated in terms of a monotonic increasing function of drop sizes [38,39,40]. By substituting a polynomial fit [38] into Equation (17), the rain rate is then computed via the DSD moments:

$$R=c_1{M}_3+c_2{M}_4+c_3{M}_5+c_4{M}_6+c_5{M}_7,$$

(18)

where the coefficients are $-1.924\times {10}^{-4}$, $9.296\times {10}^{-3}$, $-1.8\times {10}^{-3}$, $1.496\times {10}^{-4}$, and $4.452\times {10}^{-6}$, respectively.

2.3. Dataset

This study analyzed the polarimetric radar measurements collected by the Next Generation Weather Radar (NEXRAD) system located near Oklahoma City (35.333°N, 97.278°W), United States. The KTLX radar is a dual-polarization radar at a frequency of around 2.8 GHz, with a resolution of 250 m in range by 0.5 degree in azimuth below 2.4 degree in elevation. The minimum and maximum elevation ranges from 0.1 to 19.5 degree. We primarily analyzed the first unblocked scan (typically 0.5 degree) in VCP 212 mode updated in 4–6 min. The radar simultaneously measures three base data, i.e., radar reflectivity factor (Z_h in mm⁶ m⁻³), mean Doppler velocity, and spectrum width, as well as three polarimetric variables, namely, differential reflectivity (Z_dr in unitless), differential phase (${\varphi }_{dp}$ in degree), and cross-correlation coefficient at zero lag. The radar dataset is organized and processed using an open-source software package called the Python ARM Radar Toolkit (Py-ART), which has been widely used in weather radar community [41,42]. Prior to performing the DSD retrieval, the data fields were smoothed by a moving window with five range gates to reduce the noise effect, and the ground clutters are removed with a quality control procedure [43,44,45]. Furthermore, the specific differential phase (K_dp in degree km⁻¹), which is the range derivative of ${\varphi }_{dp}$, was computed using an iterative linear regression process [46] and compared with NEXRAD level 3 products [47,48]. Finally, this study utilized three processed radar variables, Z_h, Z_dr, and K_dp, since they are more sensitive to the variability of the DSD when compared to the others [6].

To validate the DSD retrieval method, a second-generation low-profile 2DVD [27] was deployed in the Kessler Atmospheric and Ecological Field Station (KAEFS) in southwest Oklahoma (34.9846°N, 97.5234°W), about 44.7 km away from the KTLX radar. The radar beam at an elevation of 0.5 degree is about 886 m higher than the disdrometer. Strong wind and turbulence may significantly affect the disdrometer measurements of raindrops [49]. Therefore, the disdrometer data were corrected by a terminal velocity filter following the theoretical relation between the diameter and velocity as in [39]. In addition, the time periods with total drop number less than 10 or rainfall rate less than 0.1 mm h⁻¹ were excluded. The processed disdrometer dataset contained 63,806 1-min drop spectra in 852 rain events collected between 0524 UTC 17 June 2006 and 1658 UTC 18 January 2017.

3. Forward Operator for Polarimetric Variables

The polarimetric variables are sensitive to ensemble properties, such as particle size, shape, and different hydrometeor phases [50,51]. For a volume of raindrops, these variables depend on not only the size spectrum, but also the scattering amplitudes, which are determined by the shape and orientation of raindrops, ambient temperature, and radar wavelength. In this section, we first introduce the polarimetric variables used in the DSD retrieval, and then utilize a forward operator to prescribe the scattering amplitude across the size bin. It is simplified with a polynomial approximation, yielding an explicit form for calculating each variable based on the Gamma DSD assumption. Finally, the forward operator is used to generate a pseudo training dataset to optimize a non-parametric estimator that is detailed in the next section.

3.1. Polarimetric Variables

This study used three polarimetric variables on a linear scale, namely, radar reflectivity at horizontal (Z_h) and vertical (Z_v) polarizations, as well as specific differential phase (K_dp). The radar reflectivity can be calculated via the radar cross-section, dielectric constant, and DSD in case of a linear polarization:

$$Z_{h,v}=\frac{4{\lambda }^4}{{\pi }^4|K{|}^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\left|f_{hh,vv}(\pi ,D)\right|}^{2}N(D)dD}({\mathrm{mm}}^6{\mathrm{m}}^{-3}),$$

(19)

where $\lambda$ is the radar wavelength in millimeter, |K|² is the water dielectric constant factor in unitless, and $f_{hh,vv}(\pi ,D)$ is the backward scattering amplitude in mm, of a raindrop with an equivalent-volume diameter (D; hereafter diameter):

$$D=2{(a^{2}b)}^{1/3}(\mathrm{mm}),$$

(20)

where a is the transverse semi-axis in mm and b is the symmetry semi-axis in mm.

Furthermore, the differential reflectivity (Z_dr) is defined as the ratio of the radar reflectivity at the two polarizations:

$$Z_{dr}=\frac{Z_h}{Z_v}.$$

(21)

The specific differential phase is mathematically formulated as the integral of the real part of the propagation constant difference and DSD:

$$K_{dp}=\frac{0.18\lambda }{\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}\mathrm{Re}\left\{f_{hh}(0,D)-f_{vv}(0,D)\right\}N(D)dD}(\mathrm{deg}{\mathrm{km}}^{-1}),$$

(22)

where $f_{hh,vv}(0,D)$ are the forward scattering amplitudes for particles with a diameter D at the horizontal and vertical polarizations, respectively, and $\mathrm{Re}\{·\}$ represents the real part of the scattering amplitude.

In this paper, the scattering amplitudes are computed by using T-matrix codes for non-spherical particles [52,53]. The oblate raindrops with the symmetry semi-axis b less than the transverse semi-axis a are uniformly distributed in a volume with a known size spectrum. Moreover, axis ratio (b/a) is a key parameter for computations of reflectivity and phase shift at the two orthogonal polarizations. To characterize the relationship between the diameter and axis ratio in an analytic form, a universal equation is applied following the observational work in [38], with a modification that the drops with diameters less than 0.5 mm are spherical [26]. This relation is expressed as

$$b/a=\{\begin{array}{cc}0.9951+0.02510D-0.03644D^2+0.005030D^3-0.0002492D^4& D>0.5\mathrm{mm}\\ 1& D\le 0.5\mathrm{mm}\end{array}.$$

(23)

The dielectric constant $|K{|}^2$ is obtained by calculating the complex indices of refraction following the least square fit of the Debye equations [54], with assumptions that the ambient temperature is 10 °C and the radar wavelength is 10.8 cm. It is noted that one typical temperature is selected for T-matrix simulations for raindrops, and the effects of temperature variation have been investigated in previous studies [5,20]. The canting angle distribution has been studied under experimental conditions [55,56]. In this study, we considered a Gaussian distribution with zero mean and small variance (10°).

Figure 1 shows the characteristics of the three polarimetric variables simulated using the disdrometer observations. In Figure 1a, it is clear that Z_dr is gradually rising with Z_h going up. However, this increasing trend is negligible when Z_h is below 0 dBZ, indicating these DSDs are dominated by small spheroid drops. The values of Z_dr then have a rapid growth when Z_h is over 0 dBZ, while the spreads of Z_dr become wider. This might show that the corresponding DSDs contain a significant amount of oblate median-size and large drops. Similarly, Figure 1b illustrates that the values of K_dp remain relatively stable at 0 deg km⁻¹ for Z_h below 37 dBZ, following a boom in K_dp when Z_h is larger than 37 dBZ. When K_dp is very small, it is difficult to accurately estimate the DSD parameters due to the noise effect [12], thus a self-consistency method may be used to provide an additional constraint for the data with Z_h below 37 dBZ [57]. For the scatterplot of Z_dr and K_dp (Figure 1c), it is interesting to note that, when Z_dr reaches 5 dB, K_dp could stay below 0.5 deg km⁻¹, implying the forward scattering amplitude may significantly differ from the backward one due to Mie scattering.

Table 1. Means, standard deviations (STD), minima (min), and maxima (max) for quantities derived from the disdrometer observations. Z_h is the radar reflectivity on a logarithmic scale, Z_dr is the differential reflectivity on a logarithmic scale, K_dp is the specific differential phase, D_m is the mass-weighted mean diameter, W is the liquid water content, and R is the rainfall rate at the surface. Note the means and standard deviations of Z_h and Z_dr are calculated on a linear scale, and then transferred to a logarithmic scale.

	Zh (dBZ)	Zdr (dB)	Kdp (deg/km)	Dm (mm)	W (g/m3)	R (mm/h)
Mean	36.8	0.67	0.076	1.30	0.314	5.4
STD	7.0	0.60	0.251	0.40	0.986	16.1
Min	3.6	0	0	0.36	0.0004	0.1
Max	58.2	5.87	7.031	5.74	35.031	479.5

shows the means, standard deviations, minima, and maxima of microphysical quantities and polarimetric variables simulated using the disdrometer data. It is noticeable that the means of microphysical quantities (D_m, W, and R) are generally small, since the entire dataset is dominated by light rain and drizzle. However, as the moderate and heavy rain cases are taken into account, the probability distributions of these quantities are significantly widened, producing relatively large standard deviations. Furthermore, the dynamic range of Z_h on a linear scale reaches $6.6\times {10}^5$ mm⁶ m⁻³ corresponding to 479.5 mm h⁻¹ for R, whereas the ones for Z_dr and K_dp on a linear scale are 3.86 and 7.03 deg km⁻¹, respectively, which are much smaller than that of Z_h and R. This implies that small measurement errors of Z_dr and K_dp may result in a significant bias on R comparing to Z_h. In contrast, the dynamic ranges of Z_dr and K_dp are comparable to that of D_m and W, indicating a good estimation can be achieved using these variables.

3.2. Forward Operator

In Equation (19), it can be found that the scattering amplitude is independent from the DSD at each size bin. If the scattering amplitude is considered a forward operator, Equation (19) can be transferred as

$$Z_{h,v}=f_{h,v}(D)\otimes N(D),$$

(24)

where $f_{h,v}(D)$ is the scattering amplitudes at the horizontal and vertical polarizations, respectively, $\otimes$ is the integration operator, and $f_{h,v}(D)\otimes$ is thus the forward operators for the radar reflectivity. Similarly, K_dp is computed in terms of a forward operator:

$$K_{dp}=f_k(D)\otimes N(D).$$

(25)

where $f_k(D)$ is the amplitude of the forward operator of K_dp.

The forward operators of Z_h, Z_v, and K_dp can be obtained by performing simulations of monodisperse raindrops with diameters varying from 0.01 mm to 8 mm. As shown in Figure 2a, the forward operator presents a quasi-linear shape in a logarithmic coordinate, while the amplitude of the forward operator monotonically increases with the diameter climbing up. This quasi-linear shape on a logarithmic scale suggests a polynomial relationship between the diameter and the amplitude on a linear scale, which is consistent with the Rayleigh approximation. Similarly, Figure 2b also illustrates a progressively increasing trend for the forward operator of Z_v. A noticeable difference can be seen when the diameter is larger than 1 mm, where the vertical amplitude of a larger drop grows relatively slower than the horizontal dimension. On the one hand, it is because the contribution of the vertical dimension (b axis) to the diameter D in Equation (20) significantly differs from that of the horizontal dimension (a axis). On the other hand, as the diameter gets larger, the axis ratio in Equation (23) gives a large discrepancy from the unity, which indicates that the large raindrops tend to be more oblate than small drops. Under the Rayleigh assumption, the radar cross section at a particular polarization is proportional to sixth order of the corresponding geometric dimension of a raindrop. Therefore, the amplitude at the vertical polarization is generally smaller than that at the horizontal polarization, particularly for larger raindrops. Moreover, Figure 2c shows that the forward operator of K_dp is also characterized by a monotonically increasing curve, but the amplitude of forward operator is equal to zero when the diameter is less than 0.5 mm, due to unity axis ratio at this diameter range in Equation (23). In addition, in this region, the forward operators of Z_h and Z_v are identical, yielding Z_dr equal to 1 (0 dB).

Figure 2 shows that the amplitudes of the forward operators for Z_h, Z_v, and K_dp may have a quasi-linear relation as a function of the diameter on a logarithmic scale, which is equivalent to a polynomial on a linear scale. Therefore, the forward operator may be approximated by a polynomial function to avoid the implicit integral form in Equations (19) and (22), by using the Gamma model of the DSD. Since the amplitude is several orders higher than the diameter, a basis polynomial needs to be selected, and then the residuals are fitted into a high-order polynomial expansion. It can be formulated as

$$f=D^b\left({\displaystyle \sum _{i=0}^p{a}_i{D}^i}\right)={\displaystyle \sum _{i=0}^p{a}_i{D}^{i+b}}$$

(26)

where b is the order number for the basis polynomial, p is the maximum order number for the polynomial fit, and a_i is the polynomial coefficient for the order i.

To obtain the optimal performance for the fitting, the mean square error (MSE) is used to select the polynomial model. It is defined as

$$\mathrm{MSE}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(p_i-a_i\right)}^2},$$

(27)

where p is predicted values and a are actual values.

Figure 3a illustrates that the MSE of Z_h is generally reduced as the basis order increases from 3 to 6, whereas the basis order of 7 produces a significant error by comparing to the others. Therefore, the basis order for Z_h is set as 6 based on the MSE test. Furthermore, for the curve with the basis order of 6, the MSE gradually decreases as the order p goes down, while a sharp fall can be found between 2 and 3. When p is larger than 3, the decline trend becomes steady. Consequently, we set p as 3 to achieve a balance between the model complexity and approximation accuracy. The fitting results are in good agreement with T-matrix computations, as shown in Figure 2a. Moreover, when the polynomial fit is applied to the real DSD, it agrees well with disdrometer data in Figure 4a, with small bias, relative error, and MSE as well as high correlation coefficient (CC) in

Table 2. Coefficients of the polynomial fits for Z_h, Z_v, and K_dp and metrics for the polynomial fits. MSE, mean square error; RE, relative error; CC, correlation coefficient.

	b	p	a0	a1	a2	a3	a4	Bias	MSE	RE	CC
Zh	6	3	1.004	−0.020	0.021	$-1.7\times {10}^{-3}$		−0.417	28.965	$7.4\times {10}^{-4}$	1
Zv	6	4	0.998	0.028	−0.055	$7.8\times {10}^{-3}$	$-3.5\times {10}^{-4}$	−0.265	13.420	$6.6\times {10}^{-4}$	1
Kdp	4	4	$-3.0\times {10}^{-5}$	$6.3\times {10}^{-5}$	$3.6\times {10}^{-6}$	$-2.1\times {10}^{-6}$	$2.0\times {10}^{-7}$	$-6.9\times {10}^{-5}$	$2.7\times {10}^{-7}$	$5.2\times {10}^{-3}$	1

.

As shown in Figure 3b, the basis orders of 4 and 5 for Z_v can yield even better MSEs than that of 6 when p is relatively large (6–9). However, the higher accuracy of fitting results from overfitting of the data samples by the model with high complexity, which is evident from the comparison between the data and fitting curve of higher orders (not shown). Following a way similar to Z_h, the b and p for Z_v are set as 6 and 4, respectively, yielding a good consistency in terms of bias, relative error, MSE, and CC (Table 2), when applied to the monodisperse raindrops (Figure 2b) and real DSDs (Figure 4b).

The selection of order number of polynomial model for K_dp is more straightforward than Z_h and Z_v (Figure 3c). The curve with the basis order of 4 generally gives the lowest MSE among the ones from 3 to 7. Meanwhile, the polynomial order is set to as 4, since the MSEs of p larger than 4 become stable at around $1\times {10}^{-6}$. Note the fitting range for K_dp is from 0.5 to 8 mm rather than from 0.01 to 8 mm for Z_h and Z_v. Figure 2c shows the comparison between the data and the fitting results, where a relatively large difference can be found for the diameter between 0.5 and 0.7 mm. When the real DSDs are considered, an increase of uncertainty is expected for K_dp below 0.01 deg km⁻¹ (Figure 4c). The rest of the coefficients in Equation (26) are listed in Table 2.

3.3. Pseudo Training Dataset

On substituting Equation (26) into Equations (24) and (25), and assuming the Gamma model for the DSD, the polarimetric variables can be re-written as

$$Z_{h,v}={\displaystyle \sum _{i=0}^{p_{h,v}}{a}_i{M}_{i+b}}={\displaystyle \sum _{i=0}^{p_{h,v}}{a}_i{N}_0\gamma \left(\mu +b+i+1,\Lambda D_{\max }\right){\Lambda }^{-\left(\mu +b+i+1\right)}},$$

(28)

and

$$K_{dp}={\displaystyle \sum _{i=0}^{p_d}{a}_i{M}_{i+b}}={\displaystyle \sum _{i=0}^{p_d}{a}_i{N}_0\left[\gamma \left(\mu +b+i+1,\Lambda D_{\max }\right)-\gamma \left(\mu +b+i+1,0.5\Lambda \right)\right]{\Lambda }^{-\left(\mu +b+i+1\right)}},$$

(29)

where M_i+b is the (i + b)th DSD moments obtained by multiplying N(D) in Equations (24) and (25) with D^i+b in Equation (26). To eliminate the factor N₀, the ratios between Z_h and Z_v, and between Z_h and K_dp are taken into consideration, yielding

$$x_1=Z_{dr}=\frac{{\displaystyle \sum _{i=0}^{p_h}{a}_{i,h}\gamma \left(\mu +b_h+i+1,\Lambda D_{\max }\right){\Lambda }^{-\left(\mu +b_h+i+1\right)}}}{{\displaystyle \sum _{i=0}^{p_v}{a}_{i,v}\gamma \left(\mu +b_v+i+1,\Lambda D_{\max }\right){\Lambda }^{-\left(\mu +b_v+i+1\right)}}}$$

(30)

and

$$x_2=\frac{K_{dp}}{Z_h}=\frac{{\displaystyle \sum _{i=0}^{p_d}{a}_{i,d}\left[\gamma \left(\mu +b_d+i+1,\Lambda D_{\max }\right)-\gamma \left(\mu +b_d+i+1,0.5\Lambda \right)\right]{\Lambda }^{-\left(\mu +b_d+i+1\right)}}}{{\displaystyle \sum _{i=0}^{p_h}{a}_{i,h}\gamma \left(\mu +b_h+i+1,\Lambda D_{\max }\right){\Lambda }^{-\left(\mu +b_h+i+1\right)}}},$$

(31)

where x₁ and x₂ are the training data; the subscripts h, v, and d represent the parameter set of b, p, and a_i for Z_h, Z_v, and K_dp, respectively; the minimum diameter for calculation of K_dp is 0.5 mm; and $\mu$, $\Lambda$, and D_max are the parameters in the Gamma model retrieved by x₁ and x₂. Since three dependent variables need to be estimated by two independent variables, the solution of the DSD retrieval may be unstable, thus an empirical constraining assumption has to be made, e.g., the $\mu$-$\Lambda$ relation used in [11]. Although this relation is believed to be valid for limited types of rain events [58] and subjective to the minimum diameter and data filtering [59], it arises from the actual rain microphysical processes [60,61]. An alternative approach to the two-parameter solution is based on the joint probability distribution of the parameters within the raindrop mass spectrum [62].

In this study, we use the thresholds of rainfall rate larger than 5 mm h⁻¹ and total number counts larger than 1000 drops min⁻¹ to avoid the fluctuations arisen from DSD fitting, which is consistent with the studies in [11,63]. It yields the relation as

$$\mu =-0.0279{\Lambda }^2+1.0619\Lambda -2.8281\left({\mathrm{mm}}^{-1}\right).$$

(32)

Figure 5 illustrates the curve for the $\mu$-$\Lambda$ relation obtained by using the ten-year disdrometer data. The result is generally consistent with the previous study in [63]. However, the MSE is 0.184 for our relation, whereas it is 0.727 for the relation in [63].

Finally, the training dataset, Z_dr and K_dp/Z_h, can be generated via Equations (30)–(32), by drawing $\mu$ and D_max. Based on the observations of dynamic ranges of $\mu$ and D_max in the ten-year-disdrometer data, we assume that the drawing poll ranges from −3 to 20 for $\mu$ and from 1.7 mm to 8 mm for D_max, with an additional constraint that D_m is no larger than D_max. Figure 6 presents the simulation results of the training dataset Z_dr and K_dp/Z_h projected onto a two-dimensional plane by assuming a number of fixed $\mu$ and D_max, respectively. In Figure 6a,b, K_dp/Z_h dramatically drops with Z_dr increasing when Z_dr less than 2. The decreasing trend then becomes steady, showing a concave curve. Furthermore, for identical Z_dr, the values of K_dp/Z_h rise as $\mu$ or D_max climbs up, while for identical K_dp/Z_h, the values of Z_dr also increase with $\mu$ or D_max going up. In addition, when Z_dr is relatively small (less than 0.318 dB), the colored curves in Figure 6 may not be disjoint. It means that two distinct parameter sets of $\mu$ and D_max can generate identical training data of Z_dr and K_dp/Z_h, leading to instability of an inverse model. To tackle this problem, we separate the data with instability by applying a threshold of Z_dr less than 0.318 dB. As shown in Figure 6, by applying this threshold, the training dataset is well-separated and the corresponding curves are mutually disjoint. As a result, the problem of inversion of polarimetric variables to Gamma parameters is stable, and the retrieval uncertainty directly relates to measurement errors exclusively.

4. Inverse Model for Gamma Parameter Retrieval

The previous section shows that the polarimetric variables are calculated by the DSDs via the forward operator. Since the generated data are uniquely associated with the physical parameters, the data may be used to retrieve the DSD parameters inversely. In this section, we first introduce a non-parametric model to solve the inverse problem, and then select the optimal model by using real DSD data. Finally, we show the effect of the measurement error on the DSD retrieval.

4.1. Non-Parametric Estimator for the Inverse Model

There are some approaches to solving the inverse problem if the training dataset is available, such as Neural Network [64], Bayesian method [65], and variational method [66]. In this paper, we adopt a nearest neighborhood method, due to its nature of simplicity and high speed.

The idea of the nearest neighborhood method is to use the observations in the training set closest in the feature space (i.e., a space composed by the polarimetric variables) to the input x to estimate the output y. Let $x=(x_1,x_2)=(Z_{dr},K_{dp}/Z_h)$ and $y=(\mu ,D_{\max })$; the k-nearest neighbor estimate $\widehat{y}$ is then defined as

$$\widehat{y}=\frac{1}{k}{\displaystyle \sum _{i=1}^{N_k(x)}{y}_i},$$

(33)

where N_k(x) is the neighborhood defined by the kth closest points to x in the training dataset, and the metric of closeness is the Euclidean distance. Since the two independent variables have very different dynamic ranges, it is necessary to perform the normalization of the variables, which produces a zero-mean and unity covariance matrix. The procedures based on Cholesky factorization [50,67] is described as:

• Subtract the averaged values from the data matrix.
• Calculate the covariance of the residual matrix.
• Calculate an upper-triangular matrix of the covariance using Cholesky factorization.
• Divisde the residual matrix by the upper-triangular matrix.

Afterwards, the k observations closest to x in the feature space are averaged to give the retrievals of $\mu$ and D_max. Finally, we calculate the estimates of the intercept N₀ by substituting the $\mu$ and D_max estimates into Equations (28) and (29) and take an average of the estimates to obtain the N₀.

4.2. Model Selection

The performance of the nearest-neighborhood estimator relates to its capability of prediction on an independent test set of the real data. It is fundamentally important to assess this performance, since it determines the selection of the model, and provides a metric of the model quality. The statistical methods commonly search for the minimum of the MSE of $\widehat{y}$ in Equation (33), yielding an optimal model for estimation of the target variables y. However, as shown in Figure 7, the DSD observed by disdrometer may suffer from significant model errors with the assumption of the Gamma shape. At some stages of a precipitating process, the DSD occasionally exhibits a bi-modal or tri-modal distribution [68,69,70]. In addition, some errors may also arise from the moment method used to fit the Gamma function to the disdrometer observations [60,71]. Due to these problems, the polarimetric variables may result in a retrieval of Gamma parameters that is significantly different from the true DSD (Figure 7).

It is noticeable that the liquid water content (W) and mass-weighted mean diameter (D_m) in Equations (5) and (6) play an important role in the representation of the DSD. Furthermore, the deviation of a Gamma function from the disdrometer observations depends on these two microphysical quantities [17,18,72]. Therefore, the objective functions of the inverse model may be modified to produce a DSD retrieval similar to the disdrometer-derived Gamma parameters. To obtain a number for k neighbors, we search $\mu$ that yields an optimal estimate of D_m by using Equation (32). Furthermore, we determine the k neighbors of D_max by optimizing W via Equation (14), where N₀ is computed using Equations (28) and (29). Figure 8 illustrates the MSEs of W and D_m associated with the model complexity. The optimal numbers are 456 for $\mu$ and 96 for D_max, but they are also affected by the spacing of data points in the pseudo training dataset.

4.3. Error Characteristics

Practically, it is difficult to obtain the polarimetric radar data free of measurement errors. These errors may have a significant impact on DSD retrievals, and sometimes lead to meaningless results. Therefore, it is important to evaluate the magnitude of the errors.

The problem may be formulated as follows. With the assumption of the perfect accuracy of the measurements of Z_dr and K_dp, the uncertainties of D_m and W would correspond to that in Z_h exclusively. Since the radar records Z_h measurements on a logarithmic scale, the impacts of the relative error (ΔZ_h) in decibel are assessed. This analysis can also be applied to Z_dr and K_dp, but K_dp needs to be considered on a linear scale. Note that we study the characteristics of relative square error (RSE), which is defined as

$$\mathrm{RSE}=\frac{{\displaystyle \sum _{i=1}^n{\left(p_i-a_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(a_i-\overline{a}\right)}^2}},$$

(34)

where

$$\overline{a}=\frac{1}{n}{\displaystyle \sum _j{a}_j}.$$

(35)

Figure 9 illustrates the analysis of error characteristics on Z_h, Z_dr, and K_dp, respectively, in the disdrometer dataset. In Figure 9a, the RSE of D_m generally presents an increasing trend associated with either the positive or negative ΔZ_h. In practice, for ΔZ_h < 0.3 dB, the RSE for positive ΔZ_h is not very large (less than 1), and it becomes large beyond this limit. By contrast, the RSE for negative ΔZ_h rises slowly, but it is more sensitive to ΔZ_h between 0.1 and 0.4 dB, leading to a lower threshold of 0.2 dB for the RSE less than 1. Moreover, the RSE of W for positive ΔZ_h increases rapidly when ΔZ_h < 0.2 dB and remains steady at around 0.2 (Figure 9d). It is interesting to note that the RSE of W for negative ΔZ_h decreases slightly when ΔZ_h is below 0.2 dB. This is because Equations (28) and (29) may produce distinct estimates of N₀, where Z_h often gives the largest one. By adding a negative error, the variance of N₀ estimates is reduced, yielding a slightly improved W.

The dependence of the uncertainty of D_m and W on Z_dr errors (Figure 9b,e) generally preserves the trend similar to the previous case (Figure 9a,d), however, some distinct features are noticeable. In Figure 9b, the RSE of D_m for positive ΔZ_dr shows a dramatic rise between 0.2 dB and 0.4 dB. It may indicate a threshold of 0.2 dB for Z_dr to obtain an accurate estimate of D_m, while the RSE is below 0.6 for negative ΔZ_dr less than 0.1 dB. Furthermore, Figure 9e presents a local minimum for ΔZ_dr at around 0.1 dB. Apart from the reduction of the variance of N₀ estimates, this minimum may be due to an effect of the thresholding of Z_dr prior to performing the nearest-neighborhood estimator. With the positive ΔZ_dr adding to the data, the number of points for the calculation of RSE increases more significantly than the squared errors, resulting in a slight decrease of RSE. Moreover, it is clear that D_m is very sensitive to K_dp errors, since a small value (e.g., 0.01 deg km⁻¹) can lead to a boom in RSE, particularly for negative bias (Figure 9c). To determine a threshold for K_dp error, it can be found that the growth rate of RSE for positive ΔK_dp is accelerated when K_dp is larger than about 0.01 deg km⁻¹ (Figure 9f), while the RSEs of the two types of errors are both small (below 0.2).

The important conclusion of this analysis is that, if the D_m and W are retrieved by using the polarimetric variables, we suggest the thresholds of 0.2 dB, 0.1 dB, and 0.01 deg km⁻¹ be applied to the errors of Z_h, Z_dr, and K_dp, respectively.

5. Results

In this section, the inverse model (IM) for the DSD retrieval is first tested with polarimetric data simulated by the ten-year disdrometer observations, and compared to an existing algorithm designed for the same regime [23], called a Bayesian approach (BA). Furthermore, the model is applied to the five-year joint observations of dual-polarized radar and disdrometer to assess its performance on real radar data. Finally, the rainfall rate calculated by the DSD is compared to the results of the power-law relations to evaluate the potential of an operational application. In addition to mean square error (MSE) and relative square error (RSE), we also use mean absolute error (MAE),

$$\mathrm{MAE}=\frac{{\displaystyle \sum _{i=1}^n\left|p_i-a_i\right|}}{n},$$

(36)

relative absolute error (RAE),

$$\mathrm{RAE}=\frac{{\displaystyle \sum _{i=1}^n\left|p_i-a_i\right|}}{{\displaystyle \sum _{i=1}^n\left|a_i-\overline{a}\right|}},$$

(37)

correlation coefficient (CC),

$$\mathrm{CC}=\frac{S_{PA}}{\sqrt{S_P{S}_A}},$$

(38)

where

$$S_{PA}=\frac{{\displaystyle \sum _i\left(p_i-\overline{p}\right)\left(a_i-\overline{a}\right)}}{n-1},$$

(39)

$$S_p=\frac{{\displaystyle \sum _i{\left(p_i-\overline{p}\right)}^2}}{n-1},$$

(40)

$$S_A=\frac{{\displaystyle \sum _i{\left(a_i-\overline{a}\right)}^2}}{n-1},$$

(41)

root mean square error (RMSE),

$$\mathrm{RMSE}=\sqrt{\mathrm{MSE}},$$

(42)

and root relative square error (RRSE),

$$\mathrm{RRSE}=\sqrt{\mathrm{RSE}}.$$

(43)

5.1. DSD Retrieval with Disdrometer-Simulated Data

At the first stage of the assessment, we apply IM to the simulated dataset (see Section 3.1) and compare to BA. Figure 10 illustrates the comparison between the disdrometer observations and the DSD retrievals. As shown in Figure 10a, the D_m results from IM show good consistency with the disdrometer data for D_m less than 0.8 mm, while a boundary exists at around 0.8 mm due to the Z_dr threshold of 0.318 dB. When D_m is between 0.8 mm and 2.5 mm, the retrievals fairly agree with the disdrometer data, following a large spread for D_m larger than 2.5 mm. By contrast, the uncertainty of the retrievals by BA gradually increases as D_m goes up, particularly for D_m larger than 2 mm (Figure 10b). Since the Bayesian approach constructs the joint distribution of polarimetric variables and DSD parameters from the radar data, it gives a better performance on the lower end comparing to IM. By analyzing the statistics in

Table 3. Mean square error (MSE), mean absolute error, relative square error (RSE), relative absolute error (RAE), and correlation coefficient (CC) for D_m and W obtained by the inverse model (IM) and Bayesian approach (BA).

	MSE	MAE	RSE	RAE	CC
Dm (IM)	0.030	0.124	0.183	0.405	0.917
Dm (BA)	0.091	0.239	0.552	0.782	0.821
W (IM)	0.113	0.062	0.128	0.178	0.963
W (BA)	0.218	0.090	0.246	0.257	0.934

, it is clear that the D_m retrievals by IM outperform that by BA. The MSE and RSE of IM are generally one third of BA, while the MAE and RAE of IM are half than of BA. The difference of the CC is less pronounced, with IM at 0.917 and BA at 0.821.

Figure 10b,d shows the comparison of W obtained by IM and BA, respectively. The retrievals of W by IM well match the disdrometer observations between 0.001 g m⁻³ and 10 g m⁻³. When W is above 10 g m⁻³, the points tend to be below the red curve, indicating an underestimation of W at this range. By comparing to IM, BA gives higher estimates when W below 0.01 g m⁻³, while it underestimates W above 5 g m⁻³ (Figure 10d). In addition, it shows high density at around 0.01 g m⁻³ in BA, due to small priori below 0.01 g m⁻³ [23]. As shown in Table 3, the statistics of IM are generally improved in terms of MSE, MAE, RSE, RAE, and CC. Notably, the retrievals of W by IM and BA both achieve a high CC, yielding good agreement between retrievals and observations.

It can be found that Z_h and Z_v are generally proportional to the sixth-order moment of DSD, while the basis polynomial of K_dp is on the fourth order (Section 3). For the simulation data, IM yields high accuracy for the sixth-order moment, Z, with RSE of 2.4$\times {10}^{-3}$ and CC of 0.9989, but very low accuracy for the zero-order moment (or total number concentration), N_t, with RSE of 3.2 and CC of 0.6183. In contrast, the zeroth and sixth moments estimated by BA are less accurate, while RSE for N_t and Z are 77.5 and 0.050, and CC are 0.2417 and 0.6884, respectively.

5.2. DSD Retrieval with Real Radar Data

When the inverse model is applied to real radar data, some problems may arise due to the sampling difference between the instruments. The radar beam often illuminates the hydrometeors at a higher altitude due to the curvature of the earth surface and the radar elevation angle. In our study, we used the first unblocked radar elevation (approximately at 0.5°), thus the altitude difference between the two instruments is about 886 m by considering the altitude of radar location. During the falling of raindrops, their size distribution may evolve as a result of coalescence, collision, breakup, accretion of droplets, and evaporation. Moreover, the polarimetric variables are calculated from the disdrometer observations under the assumption that the raindrops uniformly distributed in the radar volume. However, the radar resolution volume above the disdrometer location is about 3.91$\times {10}^7$ m³, whereas the sensing volume of the disdrometer is around 2.4 m³. Therefore, the variables may be biased in the real case if the radar volume is not filled completely. These sampling problems lead to an uncertainty when real radar data are used for the DSD retrieval.

To compare the radar measurements with the disdrometer observations, a Cartesian coordinate was established with the origin located at the radar center. The perpendicular coordinate of the disdrometer relative to the radar center was then calculated with the haversine formula [73]. To achieve a good consistency, we selected five radar gates closest to the disdrometer coordinate and shift the disdrometer data within $\pm 5$ min for each precipitating process, corresponding to the radar temporal resolution of 4–6 min. The pairs of the radar data and disdrometer data with the highest correlation were processed as a validation dataset for the DSD retrieval. In addition, the mean bias was deduced from the radar measurements for each rain process.

Figure 11 illustrates the comparison between D_m and W obtained by IM and BA using real radar data. It is noticeable that the D_m retrievals by IM concentrate on the range between 1.2 and 2.8 mm, whereas there is a large variation for the retrievals by BA. The BA retrievals also show a lower bound at around 0.01 g m⁻³; by contrast, the IM ones are less biased in this region. Furthermore, the statistics in

Table 4. Same as Table 3, but for real radar data.

	MSE	MAE	RSE	RAE	CC
Dm (IM)	0.062	0.181	0.747	0.882	0.504
Dm (BA)	0.123	0.270	1.470	1.315	0.463
W (IM)	1.081	0.300	0.515	0.492	0.705
W (BA)	1.167	0.317	0.557	0.522	0.679

indicate that IM achieves better agreement than BA in terms of MSE, MAE, RSE, and RAE, but CCs for the two models are generally low. In addition, the W retrievals by IM and BA have a wide spread when compared to the results of the simulated data, due to the measurement errors, noise, and sampling problems.

It is also very important to evaluate the DSD moments obtained by the inverse model, since the moments play an important role in the representation of the DSD. In this study, we compared the zeroth, second, fourth, and sixth moments, which are commonly used in the DSD estimation [21,37,74]. Since the infinite zeroth-moment sometimes results from $\mu$ less than −1, an approximation method has been used. We first calculated the drop concentration N(D) at each size bin using Equation (3), and then integrated the moment concentration $D^{n}N\left(D\right)$. By using this approximation method, the zeroth-moment can yield valid error statistics, while the second, fourth, and sixth moments remain the values calculated by Equation (4). In addition, the Bayesian method has constructed the prior probability with W and D_m, leading to invalid $\mu$ by solving the quadratic equation in Equation (7). Therefore, the DSD moments calculated by the Bayesian retrievals are not presented here.

Figure 12 illustrates the comparison of the moments calculated by the disdrometer observations and retrieved by the radar data. The red curve (reference) crosses the region of high intensity of the blue dots (data) for the second, fourth, and sixth moments. However, the blue dots for the zeroth moment tend to be biased toward the ordinate, leading to an overestimation.

The statistics in

Table 5. Root mean square error (RMSE), mean absolute error (MAE), root relative square error (RRSE), relative absolute error (RAE), and correlation coefficient (CC) for the zeroth, second, fourth, and sixth moments, respectively.

	RMSE	MAE	RRSE	RAE	CC
Zeroth	14,373.8	4602.5	1.233	0.986	0.651
Second	2453.2	660.7	0.746	0.497	0.696
Fourth	2469.1	734.6	0.721	0.480	0.714
Sixth	16,391.1	3800.8	0.855	0.455	0.659

provide quantitative assessment for the performance on the moment estimation. In general, for RMSE, MAE, RRSE, and CC, the second and fourth moments yield better results when compared to the zero and sixth moments, since the inverse model is optimized by W, which is proportional to the third-order moment. By contrast, RAE is constantly reducing as the order increases. In addition, the RMSEs and MAEs of the zeroth and sixth moment are much larger than the others, due to the significance of its magnitude.

5.3. Rainfall Rate Estimation

The rainfall rate estimation is one of the most important applications of polarimetric weather radar [75]. Since the rainfall rate in Equation (18) is related to DSD, a physical relationship between the polarimetric variables and rainfall rate exists. For rainfalls at the experiment site, the curves of Z_h and R in a logarithmic scale present a linear shape, and therefore, yield an empirical relation as

$$R(Z_h)=0.017Z_h^{0.714}.$$

(44)

By applying Z_dr to the Z-R relation, the performance of the rainfall rate estimation is improved, as demonstrated in the JPOLE field experiment [76]. The dual-polarization algorithm is expressed as

$$R(Z_h,Z_{dr})=0.0142Z_h^{0.770}{Z}_{dr}^{-1.67}$$

(45)

In this study, we compared the rainfall rate of the inverse model (IM) calculated by Equation (18) with that of the existing algorithms R(Z_h) in Equation (44) and R(Z_h, Z_dr) in Equation (45).

Figure 13 illustrates the time series of the radar measurements of a convective case that occurred between 0614 and 1519 UTC on 8 June 2014. The precipitating system arrived at the location of the disdrometer at the mature stage, with a maximum reflectivity of 51.5 dBZ. As it propagated eastward, the system gradually decayed and reached an average of 25 dBZ between 1200 and 1500 UTC. Furthermore, there was a large variation in the time series of the polarimetric variables Z_dr and K_dp between 0800 and 0930 UTC, associated with the DSD variability during the mature stage. The two variables then remained steady at around 0.5 dB and 0 deg km⁻¹, respectively, as the system became stratiform rainfall.

Figure 14 shows the results of rainfall rate estimation by IM, R(Z_h), and R(Z_h, Z_dr). In general, the three retrievals well match the disdrometer observations with the correlation coefficient of 0.931, 0.932, and 0.937, respectively, while the MSE of IM is 95.2, about half that of R(Z_h) and R(Z_h, Z_dr). However, all three algorithms tend to give an underestimation of the rainfall rate at the mature stage between 0800 and 0900 UTC. It is more evident from the one-hour rain accumulation (R_1h), as shown in Figure 15. When R_1h reaches its maximum at the third hour of the precipitating process, the bias (relative error) of IM is 1.63 mm (4.5%), comparing to 6.62 mm (18.2%) for R(Z_h) and 11.69 mm (32.2%) for R(Z_h, Z_dr). It can be concluded that the inverse model gives some advantages over the power-law relations, particularly for the periods with high DSD variability.

6. Summary and Conclusions

In this paper, an inverse model for raindrop size distribution (DSD) retrieval has been designed and validated using observations from polarimetric radar and a two-dimensional video disdrometer.

First, the polarimetric variables, Z_h, Z_dr, and K_dp, were calculated via the T-matrix method to study the combined effects of the DSD and scattering properties of raindrops. Subsequently, a forward operator was obtained for each polarimetric variable by performing simulations for monodisperse raindrops under assumptions of ambient temperature, radar frequency, drop shape, and orientation. To yield an explicit form for the forward operator, the high-order polynomial approximation was utilized and optimized according to the mean square error test. The forward operator was then applied to assumed DSDs to generate a pseudo training dataset used by the inverse model. In addition, an empirical relation between the shape parameter $\mu$ and the slope parameter $\Lambda$ was used as an additional constraint to reduce the effect of instability in the forward operator when the maximum diameter was taken into account.

Secondly, a non-parametric estimator based on the nearest-neighborhood method was developed to map the polarimetric variables to DSD parameters. Since the input variables were on different scales, a normalization procedure was introduced to produce a zero-mean and unity-variance data matrix. The inverse model was then trained with regard to the mass-weighted mean diameter (D_m) and liquid water content (W), since these two microphysical quantities were important for the representation of DSD. Moreover, the characteristics of measurement errors of the polarimetric variables were investigated, leading to the conclusion that the error thresholds for Z_h, Z_dr, and K_dp are 0.2 dB, 0.1 dB, and 0.01 deg km⁻¹, respectively.

Finally, the inverse model was validated by ten-year disdrometer data, and five-year joint observations from polarimetric radar. The results were also compared to an existing algorithm called a Bayesian approach. When the inverse model was applied to the simulated data, it generally yielded better agreement with the disdrometer observations when compared to the Bayesian approach. For real radar data, both the inverse model and Bayesian approach produced DSD retrievals with large uncertainties due to the measurement errors, noise, and sampling problems of the instruments. However, the inverse model still resulted in better DSD retrievals than the Bayesian approach in terms of mean square error, mean absolute error, relative square error, relative absolute error, and correlation coefficient. To make an operational use, the rainfall rate derived from the inverse model was compared to the power-law relations R(Z_h) and R(Z_h, Z_dr), showing that the inverse model can improve rainfall rate estimation, particularly for the rainfall regimes with highly varied DSDs.

The inverse model may be extended to other DSD models, such as normalized Gamma model in Equation (5) and generalized Gamma model in Equation (10). In the normalized Gamma model, the normalized intercept is independent from the shape parameter $\mu$ in Equation (6), thus the retrieval procedure is equivalent to the ones proposed in this paper. In the generalized Gamma model, the moments play a key role in calculations of Z_dr and K_dp/Z_h via Equation (12). An inverse mode is then constructed to retrieve the moments within the generalized Gamma model.

The X-band radar has attracted increasing attentions in recent years, due to its low-cost, fine resolution, and high sensitivity to light rain. The inverse model may be adapted to the X-band frequency by considering three variables as a result of Mie effects, namely, specific attenuation (A_H), differential attenuation (A_DP), and differential phase shifts upon scattering (${\mathsf{\delta }}_{\mathrm{b}}$). Therefore, we can simultaneously obtain the rain-path attenuation correction and raindrop size distribution retrieval by using an iterative approach. Some pioneer work for X-band radars can be found in [19,20,77,78].

In the future, we will also apply the inverse model to the radar observations in the complex terrain [79]. With the intervention of the surface topography, the precipitation near the surface cannot be observed, giving a significant challenge on the DSD retrieval as well as rainfall rate estimation. Furthermore, we will investigate the rain characteristics in a tropical regime [72], with the DSD parameters derived by the inverse model.