This paper presents a method to retrieve raindrop size distributions (DSD) from slant profile dual-polarization Doppler spectra observations. It is shown that using radar measurements taken at a high elevation angle raindrop size distributions can be retrieved without making an assumption on the form of a DSD. In this paper it is shown that drop size distributions can be retrieved from Doppler power spectra by compensating for the effect of spectrum broadening and mean velocity shift. To accomplish that, spectrum deconvolution is used where the spectral broadening kernel width and wind velocity are estimated from spectral differential reflectivity measurements. Since convolution kernel is estimated from dual-polarization Doppler spectra observations and does not require observation of a clear-air signal, this method can be used by most radars capable of dual-polarization spectra measurements.
To validate the technique, sensitivity of this method to the underlying assumptions and calibration errors is evaluated on realistic simulations of radar observations. Furthermore, performance of the method is illustrated on Colorado State University–University of Chicago–Illinois State Water Survey radar (CSU–CHILL) measurements of stratiform precipitation.
Rainfall radar observations depend on underlying raindrop size distributions. Since rain-rate retrieval from radar measurements is a problem of current interest there is a great importance in independent measurements of naturally occurring drop size distributions. Often in situ measurements of raindrop size distributions are performed using disdrometers (Bringi et al. 2003). Despite wide acceptance of these measurements as a ground truth, the small sampling volume limits accuracy of these measurements (Ulbrich and Atlas 1998). Jameson and Kostinski (2001) have argued that raindrop size distributions (DSDs) are sampling volume dependent and that observed raindrop size distributions should be compared only if observed at comparable scales.
The relation between raindrop fall velocity and equivolumetric diameter (Atlas et al. 1973) brought a possibility of converting vertical Doppler radar measurements to raindrop size distributions. This approach allows for increasing an observation volume to radar resolution volume scales. Since Doppler power spectra of rainfall are influenced by turbulence and wind, the DSD retrieval procedure is not straightforward. Atlas et al. (1973) have shown that if effects of turbulence and wind are not accounted for, the retrieved DSD would be largely erroneous.
Parameterization of raindrop size distributions (Ulbrich 1983; Willis 1984) allowed for using optimization approach to find parameters of DSD. Hauser and Amayenc (1981) have shown that under certain assumptions intercept parameter and median volume diameter of exponential DSD can be retrieved from Doppler spectra observations. The proposed method, however, did not include effects of turbulence and wind. Williams (2002) have extended this approach by including spectrum broadening and wind effects and assuming that raindrop size distributions follow a gamma functional form that allowed DSD parameters retrieval from profiler measurements. A similar procedure was also implemented by Moisseev et al. (2006) for slant profile radar measurements. Russchenberg (1993) used spectral moments to retrieve parameters of a drop size distribution. The main limitation of these methods, though, is dependence on one or another functional form of the DSD.
The sensitivity of VHF wind profilers to both clear-air and precipitation signals has brought another possibility for estimating drop size distributions. Gossard (1988) and Rajopadhyaya et al. (1993) have proposed to use deconvolution techniques to remove the effect of turbulence and wind from Doppler spectra measurements. In this procedure the broadening kernel width and wind velocity were obtained from clear-air measurements. The drawback of this approach is that clear-air spectrum has a finite length and often contaminates measurements of the precipitation spectrum. To mitigate this problem Currier et al. (1992) has proposed the use of dual-frequency profiler retrieval, by combining clear-air measurements from a VHF profiler with precipitation observations taken by a UHF profiler. May et al. (2001, 2002 and May and Keenan (2005) have shown that the dual-frequency approach can successfully be applied in variety of measurement conditions.
Introduction of dual-polarization radar measurements (Seliga and Bringi 1976, 1978) has brought a new possibility to radar rainfall observations. Dual-polarization rain-rate retrieval techniques have been demonstrated to be more accurate than reflectivity-based ones. The dual-polarization methods rely on changes of raindrop shapes with diameter. Generally, raindrop shape is approximated by a spheroid with axis ratio being a function of raindrop diameter. There are a number of theoretical and experimental studies on raindrop shapes (Pruppacher and Beard 1970; Beard and Chuang 1987; Andsager et al. 1999; Thurai et al. 2006).
Kezys et al. (1993) and Unal et al. (1998) have demonstrated measurements of the spectral differential reflectivity in precipitation. Recently Moisseev et al. (2006), Spek et al. (2005), and Unal et al. (2001) have combined dual-polarization and Doppler observations to retrieve microphysical properties of precipitation. Wilson et al. (1997) have also shown that a combination of dual-polarization and Doppler observations can be used to constrain the shape parameter of a gamma DSD. Based on simulations, Moisseev et al. (2006) have shown that an optimum combination of Doppler and dual-polarization measurements can be achieved if observations are taken at elevation angles ranging between 30° and 60°. Furthermore, it was shown that spectral differential reflectivity measurements can be related to raindrop shapes and air motion. In this work we propose to use spectral differential reflectivity observations to estimate a spectrum broadening kernel width and wind velocity. Then by applying this information one can remove influence of spectral broadening and wind from a Doppler power spectrum and therefore retrieve a raindrop size distribution.
This paper is structured as follows. In section 2, dual-polarization spectral observations are introduced. Furthermore, in this section the connection between spectral differential reflectivity and precipitation microphysical properties is described. The DSD retrieval method is outlined in section 3. In section 4 a detailed error analysis of the proposed method is given. Implementation of the method to the Colorado State University–University of Chicago–Illinois State Water Survey radar (CSU–CHILL) data is shown in section 5. Finally discussion and conclusions are given in sections 6 and 7.
2. Differential reflectivity spectrum
Doppler power spectra of precipitation echoes can be written as follows (Wakasugi et al. 1986):
where Spr(υ) is the power spectrum due to scattering from hydrometeors, Sb(υ) is the spectrum broadening kernel, the asterisk (*) denotes circular convolution, and υ0 is the wind radial velocity component. The spectrum broadening is caused by turbulence, crosswind, and spectral window. The spectrum broadening kernel can be approximated as following the Gaussian shape (Doviak and Zrnic 1993; Wakasugi et al. 1986; Rajopadhyaya et al. 1993):
where σb is the broadening kernel width given in m s−1 and includes contributions from all causes of spectral broadening (Doviak and Zrnic 1993, section 5.3). The precipitation spectrum, Spr (υ), can be written in terms of scattering cross section, σ(D), and drop size distribution, N(D), as follows:
where |dD/dυ| is the Jacobean of diameter to velocity transformation.
For observations taken at some elevation angle, θ, the radial velocity of a raindrop can be written as (Atlas et al. 1973)
where ρ0 and ρ are the air densities at the surface and the measurement altitude, respectively.
If radar measurements are taken at an elevation angle ranging between 30° and 60°, one can combine dual-polarization observations with Doppler spectra observations to get more insight into precipitation microphysics (Moisseev et al. 2006). In cases of no spectrum broadening and no wind, the spectral differential reflectivity can be defined as
One can see that in this particular case the spectral differential reflectivity depends only on two copolarized-scattering cross sections. This ratio is fully defined by raindrop shapes and therefore can be related to an existing raindrop size—shape relation, as can be seen in Fig. 1a. It is a continuing debate about a naturally occurring raindrop size–shape relation. As a result there is a number of known relations that are either based on experimental observations (Andsager et al. 1999; Pruppacher and Beard 1970; Chandrasekar et al. 1988; Bringi et al. 1998; Thurai et al. 2006), theoretical considerations (Beard and Chuang 1987), or a combination of both (Brandes et al. 2002). In Fig. 2 we have depicted some of the known relations and the corresponding spectral differential reflectivities are shown in Fig. 1a). From these two figures one can observe that differences in the shapes for diameters between 0.5 and 4 mm are the most important for the spectral differential observations. The relative unimportance of raindrop shapes for diameters larger than 4 mm is caused by nonlinearity of the expression (4). The effect of an assumed axis ration relation on the retrieval will be discussed later.
In reality, Zdr(υ) observations are influenced by spectral broadening and the resulting observation can be expressed as
The resulting Zdr(υ) values depend not only on parameters of the broadening kernel and raindrop shape–size relation but also on a DSD. Assuming an exponential relation (Bringi and Chandrasekar 2001) with Nw = 8000 mm−1 m−3 and D0 = 1.2 mm we have plotted Zdr (υ) in Fig. 1b for different values of the spectrum broadening kernel width. One can observe a strong dependence of Zdr(υ) on the spectrum broadening width. It should be noted that a nonzero wind radial velocity would translate the spectral differential reflectivity, defined by (6), along the velocity axis.
3. DSD retrieval
To retrieve drop size distributions from Doppler spectra measurements, one needs to compensate for the effect of spectrum broadening and wind. Gossard (1988) and Rajopadhyaya et al. (1993) have shown that VHF profiler measurements of clear-air and precipitation returns can be used to estimate raindrop size distributions. The wind-profiling DSD retrieval procedure is based on deconvolution of an observed precipitation spectrum (1) by using a clear-air signal spectrum as the broadening kernel (Gossard 1988; Rajopadhyaya et al. 1993). This method gives a direct way of retrieving raindrop size distributions from VHF Doppler radar observations. For shorter wavelength radars, scattering from hydrometers is much larger than Bragg scattering, and therefore this method is not always applicable.
Let us consider dual-polarization slant profile spectral measurements taken at an elevation angle ranging between 30° and 60°. In the previous section it was shown that in the absence of spectral broadening and wind the spectral differential reflectivity is completely defined by a drop size–shape relation. Assuming that this relation is known, one can use expression (5) as the reference differential reflectivity spectrum that corresponds to Doppler power spectra due to scattering from hydrometeors, as given by expression (3). The expression (5) is independent of a DSD. Therefore, the estimation procedure of the spectrum broadening kernel width and wind radial velocity component can be formulated as a least squares problem of finding parameters σb and υ0 that minimize the sum square residuals, SSR:
Here Zdecdr(υ) is obtained by applying the deconvolution to the observed hh and vv power spectra with a convolution kernel Sb(υ − υ0, σb) defined by the expression (2), and υmax is the value of the radial component of the largest raindrop terminal velocity (here it is assumed that raindrops with diameters of 0 to 8 mm are present in the radar volume). The coherency spectrum Whv(υ) is defined as follows:
where Shv(υ) is the cross spectrum between hh and vv signals, and Shh(υ) and Svv(υ) are observed copolar power spectra. From expression (7) it can be seen that the coherency spectrum is used to eliminate parts of the spectral reflectivity that are affected by noise and clutter, as discussed in Moisseev et al. (2000).
b. Practical implementation and simulation results
To illustrate feasibility of the procedure we have simulated dual-polarization radar observations using the method described by Chandrasekar et al. (1986). For these simulations the shape of spectra was assumed to be determined by expressions (1) and (3), where the drop size distribution follows gamma functional form (Bringi and Chandrasekar 2001, section 7.1.4). The input parameters to this simulation are given in the Table 1. For each simulation run Doppler spectra were estimated by averaging over 30 simulated spectra.
The DSD retrieval procedure consists of three steps. At the first step parameters of the convolution kernel are found. As described above this step is carried out by solving the following minimization problem:
To deconvolve Doppler spectra we have used the iterative technique proposed by Lucy (1974). This deconvolution technique is based on Bayes’ theorem on conditional probabilities and maximizes likelihood of an observed sample. An example of the cost function for one of the simulation runs is shown in Fig. 3. It was observed that in most of the cases the cost function (7) is well behaved. However, in some cases the function can have several local minima, as can be seen in Fig. 3, which can be associated with instabilities in the deconvolution process. To obtain a spectral broadening width and wind velocity component that correspond to a global minimum, the optimization procedure is initiated using several seed values. The seed values are selected randomly from a range of values in which the spectrum broadening width and wind velocity are expected to lie. One can observe that the solution of the problem (9) gives the wind velocity component and spectrum broadening width values that are close to the input ones. For this simulation the spectrum broadening width was equal to 0.5 m s−1 and ambient wind velocity 0 m s−1. It should be noted that the difference between input and retrieved kernel width values can be attributed to the effect of spectrum window. The Chebyshev window with the sidelobe level of −50 dB is used for spectral estimation.
At the second step the convolution kernel is used to recover Spr(υ). In Fig. 4 simulated and recovered power spectra, spectral differential reflectivity, and coherency spectra are shown. It can be seen that the recovered spectrum, solid line with crosses, closely follows the modeled spectrum where σb = 0. It is interesting to note that the deconvolution procedure did not affect much the coherency spectrum values in the precipitation region. That comes from the fact that the values of the coherency spectrum in the precipitation region of the spectrum are not influenced by the spectrum broadening. Similar to the copolar correlation coefficient the copolar coherency spectrum is related to the statistical correlation between hh and vv radar signals.
In the final step the drop size distribution is estimated using the following expression:
To assess the accuracy of this method 100 simulation runs were carried out. The resulting bias and standard deviation of the retrieved log N(D) is shown in Fig. 6. One can observe that under current assumptions bias in the log N(D) is less than 0.2 and the standard deviation in log N(D) is less than 1 for diameters less than 7 mm. Also, estimated D0, rain-rate values, and corresponding standard deviations are given in Table 2.
The sensitivity of the proposed DSD retrieval method on the input DSD parameters was also tested. For this test input D0 values were varying from 0.5 to 2.5 mm with a step of 0.5 mm, and μ were varying from 0 to 5 with a step of 1. The remaining input parameters are given in the Table 1. Comparing results of the simulation, given in the Table 2, with an error analysis of the dual-frequency technique shown in Schafer et al. (2002), one may conclude that these two techniques have comparable performances.
Based on these simulations it was observed that a minimum drop diameter that can successfully be resolved is 0.3–0.4 mm and the maximum resolvable drop diameter is 7–8 mm. Due to increased errors in N(D) for small equivolumetric diameters, retrieved D0 values have larger errors in cases where input D0 values were smaller than 1 mm, as shown in Fig. 7. This effect is caused by the edge effects of the deconvolution procedure, where edges of deconvolved spectra have larger errors, and not by dependence of Zdr on D0. As a matter of fact, the reference spectral differential reflectivity is independent of DSD parameters. In Fig. 7 the bias and standard deviation of the retrieved D0 are shown. It was observed that the bias and the standard deviation do not exhibit dependence on μ; as a result the presented results are averaged over all simulations with different input μ values.
4. Evaluation of the methodology
a. Effect of drop size–shape relation
In the previous section we have shown that the proposed method is able to retrieve a DSD under the assumption that a raindrop size–shape relation is known. Currently there are a number of known relations that would result in different Zdr(υ) values especially for the diameters ranging between 0.5 and 4 mm, as shown in the Fig. 1. In this diameter range relations by Pruppacher and Beard (1970) and by Andsager et al. (1999) differ the most. Therefore, to study dependence of the proposed method on the raindrop shape assumption we have simulated 100 spectra using the Andsager et al. (1999) size–shape relation; the remaining input parameters are given in the Table 1. Then the spectrum broadening kernel width and wind velocity radial component were retrieved by solving (8), where the ratio [σhh(D)/σvv(D)] was calculated assuming the (Pruppacher and Beard 1970) relation. The results of this study are given in the Fig. 8. It can be seen that as a result of the uncertainty in the drop size–shape relation, one might expect an increase in the log N(D) bias of the estimate for diameters less than 5 mm. In Table 2 retrieved values and standard deviation of D0 and rain-rate estimates are shown. One can observe that an error in the assumption has a strong result on the retrieved median volume diameter and rain rate.
b. Effect of Zdr calibration
In the proposed method the retrieval of the spectrum broadening kernel width and wind velocity depends on Zdr(υ) observations and hence can be expected to be influenced by the calibration errors. To investigate this effect we have simulated three datasets of 100 spectra each, with Zdr offsets of 0.1, 0.2, and 0.5 dB, using parameters given in the Table 1. Then the DSD retrieval method was applied to each dataset. The results of this study are summarized in Fig. 9. It can be observed that Zdr offset has a strong influence on the bias in log N(D). It ranges between 0 and 0.5 for the offset of 0.1 dB, −0.2 and 1 for the offset of 0.2 dB, and between −1 and 1.7 for the offset of 0.5 dB.
In the Table 2 the retrieved values of D0 and rain rate are shown. One can see that 0.1-dB Zdr offset results in about 20% standard error in the retrieved D0 and 40% in the retrieved rain rate. A 0.2-dB offset would result in 50% standard error in D0 and 120% error in the rain rate. And a 0.5-dB Zdr offset would result in underestimation of D0 by almost 4 times and overestimation of rain rate by almost 6 times.
Another way of approaching this problem would be to consider Zdr offset as an additional unknown parameter and to solve (8) for three parameters, σb, υ0, and the offset. This approach was applied to the simulated radar observations with Zdr offset of 0.2 dB. The resulting DSD, bias, and standard deviation are shown in Fig. 9 as black solid lines. It can be seen that this approach reduces bias in the DSD estimate, which as a result ranges between −0.25 and 0.5, but it also increases the standard deviation of the estimate. And the retrieved values of D0 and rain rate have standard errors of 10% and 15%, respectively.
5. CSU–CHILL data
On 23 July 2004 time series data measurements were collected during a stratiform rain event by the CSU–CHILL radar. The radar measurement settings are given in Table 3. The observed reflectivities for this event were around 35 dBZ. To investigate a performance of the DSD retrieval method, the observed spectra for each range gate were averaged over complete dataset. Prior to averaging, all spectra were shifted to the same mean velocity (Moisseev et al. 2006). Then the method was applied to the averaged spectra. For the estimation of the convolution kernel (Beard and Chuang 1987) raindrop shapes were assumed. The retrieved DSD were further averaged over five neighboring range gates, resulting in the effective range resolution of 250 m. The resulting raindrop size distributions for selected range gates are shown in Fig. 10. Moisseev et al. (2006) have used the same dataset to estimate precipitation DSD parameters by fitting the parametric model to the observed spectra. It was concluded that this measurement can be characterized by a gamma DSD with log Nw = 3.54, D0 = 1.2 mm, and μ = −0.4. Moisseev et al. (2006) have also shown that the retrieved DSD parameters yield Zdr values comparable to the measured ones. The resulting log N(D) plot is shown in Fig. 10. One can observe that DSD retrievals from these two methods are in good agreement. The proposed nonparametric DSD retrieval method, nonetheless, shows more detail in the observed DSD.
6. Discussion and conclusions
In this paper we have proposed a new nonparametric method for the retrieval of drop size distributions. This method does not use any assumption on the form of a DSD. It was shown that this method allows for an accurate DSD retrieval over large range of diameters. Since this is a radar-based retrieval method the accuracy of a retrieved DSD does not suffer from a small observation volume. This method is based on slant profile dual-polarization radar observations and does not require radar sensitivity to the Bragg scatter. As a result the proposed retrieval technique can be used by most dual-polarization radars.
The research was supported by the National Science Foundation (ATM-0313881).
Corresponding author address: Dmitri N. Moisseev, Department of Electrical and Computer Engineering, Colorado State University, Fort Collins, CO 80523. Email: email@example.com