Doppler spectra from vertically profiling radars are usually considered to retrieve the raindrop size distribution (DSD). However, to exploit both fall velocity spectrum and polarimetric measurements, Doppler spectra acquired in slant profiling mode should be explored. Rain DSD samples are obtained from simultaneously measured vertical and slant profile Doppler spectra and evaluated. In particular, the effect of the horizontal wind and the averaging time are investigated.
The Doppler spectrum provides a way to retrieve the DSD, the radial wind, and a spectral broadening factor by means of a nonlinear optimization technique. For slant profiling of light rain when the horizontal wind is strong, the DSD results can be affected. Such an effect is demonstrated on a study case of stratiform light rain. Adding a wind profiler mode to the radar simultaneously supplies the horizontal wind and Doppler spectra. Before the retrieval procedure, the Doppler spectra are shifted in velocity to remove the mean horizontal wind contribution. The DSD results are considerably improved.
Generally, averaged Doppler spectra are input into this type of algorithm. Instead, high-resolution, low-averaged Doppler spectra are chosen in order to take into account the small-scale variability of the rainfall. Investigating the linear relations at fixed median volume diameter, measured reflectivity-retrieved rainfall rate, for a slant beam, the consistency of the integrated parameters is established for two averaging periods. Nevertheless, the corresponding DSD parameter distributions reveal differences attributed to the averaging of the Doppler spectra.
The new aspects are to obtain the same retrieval quality as vertically profiling and highly averaged spectra in an automated way.
Rain is generally measured by radars in the 1–10-GHz frequency range. Acquiring the raindrop size distribution from radar data is still a challenge. A two-parameter exponential drop size distribution (DSD) can be retrieved using the reflectivity Z and the differential reflectivity Zdr. The specific differential phase Kdp provides a third radar observable in the case of heavy precipitation to strengthen the retrieval technique. For radar slant or vertically profiling, the Kdp values are usually too small. Furthermore, the Zdr measurand cannot be employed when the elevation is near the vertical, or in the case of drizzle and light rain. The Doppler power spectrum is related to the drop size distribution and consequently provides a way to retrieve the required distribution.
A comprehensive review of the information content of the Doppler power spectra of rain and snow can be found in Atlas et al. (1973) in the case of vertically profiling radars. Hauser and Amayenc (1981) propose a least squares fitting of a theoretical Doppler spectrum, depending on a two-parameter exponential DSD and the vertical wind, to the measured one. Spectral broadening of the Doppler spectrum is, however, not considered. Techniques based on wind profiler measurements take advantage of the possibility to measure a bimodal Doppler power spectrum in the case of precipitation (Wakasugi et al. 1986; Rajopadhyaya et al. 1994; May et al. 2001, 2002; May and Keenan 2005). One peak is related to the clear-air echo (Bragg scattering) and the other one corresponds to the Rayleigh scattering of rain. In that case, the vertical wind and the spreading of the Doppler spectrum can be derived from the clear-air peak. Nonetheless, this is not always possible because the two peaks may overlap—that may occur in convective precipitation and moreover, this leads to a minimum observable particle diameter (May et al. 2001). This behooves Williams (2002), and later Moisseev et al. (2006), to develop techniques to retrieve the raindrop size distribution, the radial wind, and the spectral broadening of the Doppler spectra. These approaches are illustrated on highly averaged Doppler power spectra.
So far, these techniques have not been further developed, exploited, and validated. It is our intention to investigate this methodology using the Doppler-polarimetric Transportable Atmospheric Radar (TARA; Unal et al. 2012) to retrieve raindrop size distribution independent of the precipitation regimes, from drizzle to heavy precipitation, with high spatial and time resolutions. Hence, the least squares approach will also be applied on Doppler spectra measured within 2.9 s to take into account possible variations of the DSD and the dynamics. Drop size distribution samples of statistically inhomogeneous rain are mostly obtained in precipitation media (Jameson and Kostinski 2001). Toward high resolution, an attempt can be made to yield the drop size distribution of rain patches. This is highly dependent on the scales of the radar measurement and the rain patch. Because TARA can profile in three directions, we can estimate three profiles of DSD samples to get further insight into the microphysical and dynamical variability of precipitation. Finally, using this wind profiler mode, we can also mitigate the effect of the radial component of the horizontal wind on the retrievals based on slant profiling measurements.
Compared to vertically profiling, slant profiling of precipitation definitively increases the complexity of the raindrop size distribution retrieval procedure because of the impact of variable horizontal winds. However, this observation setup gives the possibility to extract the information contained in both the Doppler spectrum and the polarimetric measurands when the elevation decreases. In addition, this measurement geometry characterizes the range–height indicator (RHI) mode, which is now part of the scanning strategy of an increasing number of weather radars.
In the least squares approach, where a modeled spectrum is fitted to a measured one, the raindrop size distribution is defined by a particular functional form. To prevent the assumption on the form of the DSD, deconvolution methods can be used when the clear-air spectrum is measured (May et al. 2001, 2002) or retrieved (Babb et al. 2000; Moisseev and Chandrasekar 2007). Babb et al. developed a deconvolution method in the case of water clouds. The large-scale vertical velocity is extracted from each spectrum by measuring the zero offset of the smallest particles. As for Moisseev and Chandrasekar, they estimate the radial wind and spectral broadening from the spectral differential reflectivity (sZdr). Because Zdr or sZdr measurement cannot be always exploited, the least squares approach is preferred in an automated procedure for all types of rain.
The paper is structured as follows. Sections 2 and 3 shortly introduce the model of the rain Doppler power spectra and the retrieval algorithm, respectively. Simulation results are provided at the end of section 3 to yield estimation errors. Nevertheless, emphasis is given in the paper to retrievals on real data in sections 4–6. Section 4 describes the measured spectra, input into the optimization algorithm, and the estimation errors of the retrievals. The impacts of the radial wind and averaging are discussed in sections 5 and 6, respectively.
2. Modeling Doppler power spectra of rain
The modeled Doppler spectrum for the horizontal transmit and receive polarization (HH) setting, , can be expressed as the Doppler power spectrum resulting from the DSD, sZHH()d, convoluted by a Gaussian-shaped kernel of width σ0 modeling spectral broadening (Wakasugi et al. 1986). Spectral broadening has several causes: turbulence, wind shear, and raindrop oscillations in the radar resolution volume. Assuming all these contributions are independent, the total spectral broadening, σ02, is the sum of the individual spectral broadening factors (Doviak and Zrnić 1993):
where υ and are Doppler velocities and
The radar wavelength and the equivolume diameter are denoted by λ and D, respectively. The dielectric factor of water is |Kw|2.
Equations (1) and (2) describe a homogeneous medium throughout the radar resolution volume and during the dwell time. This hypothesis has a higher probability to be fulfilled in the case of high-resolution radars. For the analytical expression of the Doppler power spectrum related to inhomogeneous media in the case of vertically profiling radars, the reader is referred to Fang et al. (2012).
The radar cross section σHH is calculated using Rayleigh–Gans scattering of spheroidal hydrometeors (Bringi and Chandrasekar 2001, chapter 2):
The term qHH is related to the integration on the raindrop symmetry axis orientation angles, the axis ratio (r), and relative permittivity (εr). The term r(D) is a combination of Andsager et al. (1999) and Beard and Chuang (1987) relations (see Fig. 1). The axial distribution (Mardia 1972) describes the orientation angle.
The raindrop size distribution N(D) follows the modified gamma distribution:
where D0 is the median volume diameter, Nw is the intercept parameter, and μ is the shape parameter.
The following radial velocity–equivolume diameter relationship (Atlas et al. 1973) is presently used:
where α is the radar beam elevation angle, is the altitude factor, and υt is the terminal velocity at sea level. This relationship is questioned in Montero-Martinez et al. 2009 at large rainfall rates. Finally, because of the contribution of the radial wind, the Doppler spectrum experiences a shift of length υ0 along the Doppler velocity axis, and the radial velocity becomes
which defines υ0 as
where the triplet (W, V, U) represents the standard 3D wind components, and ϕN is the radar beam azimuth angle related to north. In the case of no spectral broadening (σ0 = 0 m s−1), the measured Doppler velocity is −. The model schematic is given in Fig. 1.
When the parameters of the DSD are known, the reflectivity [Eq. (9)]; the liquid water content [Eq. (10)] with ρw = 10−3 g mm−3; the number of raindrops, termed number concentration [Eq. (11)]; and the rainfall rate [Eq. (12)] can be calculated:
3. Retrieval algorithm
The retrieval algorithm obtains the three parameters of the DSD (D0, Nw, μ) and the dynamic parameters (υ0, σ0) by fitting modeled spectra to measured spectra. An optimization procedure minimizes the difference between the fitted spectrum and the measured spectrum by varying the five input parameters:
May et al. (1989) showed that fitting the spectrum on a logarithmic scale requires sufficient signal-to-noise ratio, but it has the advantage that the spectral reflectivity values have the same variance, which leads to a more robust fit.
The five-parameter minimization problem is simplified by separating the retrieval of the intercept parameter Nw:
This allows the direct derivation of the intercept parameter without nonlinear fitting. For Nw estimation, the values of D0, μ, and σ0 have to be known, but the value of the radial wind υ0 is not relevant and hence equals 0 m s−1.
A second simplification of Eq. (13) is done on the estimation of υ0. Assuming the other four parameters of the model are known, the shift between the modeled and measured spectra can be obtained by determining the lag of the cross correlation of the measured and the modeled spectra, .
The practical implementation of the algorithm differs from Moisseev et al. (2006). The cost function [Eq. (13)] is generally not smooth, which requires special attention to the convergence of the minimization in an automated procedure. Because of multiple minima in the cost function, an iterative cascaded retrieval algorithm is preferred to obtain (D0, μ, σ0). This approach is described in the appendix and visualized in Fig. 2.
b. Quality of the retrieval technique based on simulation
To get insight into the quality of the optimization procedure, the optimization is applied on simulated Doppler power spectra. By comparing the input parameters used to create a simulated spectrum with the retrieved parameters, conclusions can be drawn on their errors. To generate signals with real statistical properties, the procedure described by Chandrasekar et al. (1986) is applied. The values of the parameters are randomly selected from the intervals given in Table 1. Like in Moisseev et al. (2006), 30 realizations of the Doppler power spectra are averaged to obtain an estimate of the true spectrum and the retrieval algorithm is only applied on the resulting spectrum when its corresponding reflectivity is between 10 and 55 dBZ (light to heavy rain). Note that the 30 realizations represent the same DSD. The same exercise is carried out on integral parameters. The root-mean-square deviations (RMSD) in Table 1 are estimation errors. The coefficient of variation (CV), RMSD normalized to the mean, is given for the parameters of which the values are positive.
4. Application of the retrieval algorithm on measured Doppler power spectra
a. Input Doppler spectra
Spectral polarimetric processing is performed to obtain noise- and clutter-free dealiased Doppler spectra (Unal and Moisseev 2004; Unal 2009) for the dual-polarized main beam, and classical spectral processing is carried out for the single-polarized beams. This full processing is implemented in real time for TARA (Unal et al. 2012). Negative velocities are downward and the Doppler velocity resolution is 3.1 cm s−1 to obtain a large number of hydrometeor Doppler bins. Prior to the retrieval algorithm, three supplementary processing steps are carried out. Possible missing data in the spectrum top, due to filtering, are replaced using neighbor data. Then an automatic low clipping level is estimated to avoid the partially filtered tails of the spectra. Finally, a light smoothing is performed to reduce the statistical fluctuations of the spectra.
Rain data with TARA were obtained during the Convective and Orographically Induced Precipitation Study (COPS) in Germany on 1 July 2007. The measurement scheme consists of a cycle of five measurements acquired at 2.9 ms [main beam with HH, vertical transmit and horizontal receive (HV), and vertical transmit and receive (VV) data, and two offset beams]. The retrieval algorithm is applied on 31 104 copolar Doppler spectra obtained with the time resolution of 2.9 s. They represent stratiform light rain during 3.5 min between 200 and 1750 m (height resolution of 14.5 m). To achieve sufficient retrieved dataset size despite the averaging of Doppler spectra, the total time is extended to 7 min in section 6. The spectrum top to clipping level is found to be 15 dB on average.
b. Wind-shifted Doppler spectra
The radial contribution of the mean horizontal wind is removed in the Doppler spectra of the slant beams (velocity shifting). The mean horizontal wind is calculated from the averaged mean Doppler velocities measured by the three beams of TARA. Two offset beams are 15° away from the main beam in two orthogonal directions (Fig. 3; Table 2). With this three-beam capability, TARA acts as a wind profiler, which uses the assumption of identical 3D wind at the same altitude. An average is performed over several minutes to reduce the error on the estimate of the horizontal wind due to the violation of this assumption. In particular, the turbulent contribution (updraft/downdraft) should be suppressed. The time average is chosen to be about 10 min for the COPS campaign.
Because the mean horizontal wind is used for the velocity shifting of the slant Doppler spectra, the mean horizontal wind only is assumed to be the same for the two slant beams at fixed altitude. Consequently, the horizontal wind and the vertical wind at the time resolution of 2.9 s are not assumed to be identical for the three beams for the raindrop size distribution retrieval. The radial wind v0 is still retrieved to account for the vertical wind and the variability of the horizontal wind for the slant beams, but we then expect a smaller value of υ0 in nonconvective situations. When velocity shifting is applied, υ0 becomes
Furthermore, using the vertical beam (α = 90°), the υ0 retrieval provides the vertical wind estimate (W). The retrieval technique is carried out on Doppler spectra both velocity shifted [Eq. (15)] and nonshifted [Eq. (8)] for the two slant beams. The corresponding retrievals will be compared to the vertical beam retrievals.
c. Averaged Doppler spectra
The impact of averaging the Doppler power spectra on the retrievals can be investigated when input averaged spectra are provided. For the number of averages, a trade-off has to be made. On the one hand, the number of averages should be small in order to measure spatially and temporally variable precipitation. On the other hand, the number of averages should be large enough to decrease the impact of statistical fluctuations of the Doppler spectrum on the retrievals. The real-time processing of TARA is based on two-averaged Doppler spectra in time (2.9 s). Additionally, 6 of those Doppler spectra are averaged to obtain 12-averaged Doppler spectra (17.5 s). To study the impact of averaging, only wind-shifted Doppler spectra are considered.
d. Quality of the retrieval technique based on measurements of the same radar resolution volume
Because the main beam probes the same medium with two polarization settings, the retrieval algorithm is applied on the Doppler spectra HH and VV. Because these spectra are highly correlated, we expect the same retrieval results. The retrieval results are slightly height smoothed to reduce their variance. Scatterplots of the retrievals are given in Fig. 4. The retrieval algorithm shows a good consistency of D0, μ, υ0, and LWC, whereas the CV of both the intercept parameter and the number of raindrops is large. Applying the same comparison on the 12-averaged Doppler spectra, the same conclusions can be drawn. However, the scales of σ0 and LWC reduce to 0.35 m s−1 and 0.05 g m−3, respectively, compared to 0.75 m s−1 and 0.13 g m−3, respectively, obtained when the two-averaged Doppler spectra are considered.
These scatterplots provide a verification of the algorithm robustness. This does not mean that the retrievals are correct. If we perform this comparison on wind-shifted and nonwind-shifted Doppler spectra, we obtain the same consistency even though the retrieval parameters are improved by using the wind-shifted Doppler spectra. This improvement is demonstrated next.
5. Impact of the radial wind on the raindrop size distribution retrievals
a. Comparison of DSD retrievals from nonwind-shifted and wind-shifted Doppler spectra: Time series at fixed height
At the start of the radar far field (200 m), the separation in distance of the three radar resolution volumes corresponding to the three beams is small. For example, at 300 m from the radar, the radar resolution volumes are separated by 79 m. Therefore, we expect to obtain retrieval results similar for the three probing beams for the light rain stratiform event. The corresponding reflectivity and horizontal wind are plotted versus time in Fig. 5. The reflectivity varies and the horizontal wind can be considered strong and stable. The mean horizontal wind speed and direction are 15.4 m s−1 and 220°, respectively.
To study the impact of the radial wind, the algorithm is applied on Doppler spectra with and without correction for the mean horizontal wind. The parameters of the three-beam DSD are given in Fig. 6 for noncorrected Doppler spectra (top panel). They reasonably agree in trend and in quantitative values for D0 and μ. Differences of 0.2 mm and 1 for D0 and μ, respectively, can be observed. The Nw values may differ by a factor of 10. Are those microphysical differences real? It looks like the time curves of the DSD parameters are shifted along the y axis. To investigate this, a comparison of the 3.5-min-averaged retrieved radial wind (−) with the radial component of the measured mean horizontal wind (−) is shown in Fig. 7,
If the mean vertical wind () is negligible, then the retrieval (−) is underestimated for the slant beams [main beam (MB) and offset beam 2 (OB2)]. In the case of OB2, it means that the retrieved shift (+) applied to the measured Doppler spectra, which is positive, is not sufficient (1 m s−1 underestimation in modulus), and leads to an overestimation of D0 for this beam. In the case of MB, the underestimated retrieved shift (+) is negative, which leads to an underestimation of D0. Consequently, a positive or negative underestimated value of the retrieved Doppler shift will lead to an overestimation or underestimation of D0, respectively, for the slant beam. This explains the top panel of Fig. 6, where the D0 retrieval of the vertical beam is compared with the ones of the slant beams.
If the contribution of the 10.5-min-averaged measured horizontal wind is subtracted from the Doppler velocities of the Doppler spectrum before the retrieval procedure, then the retrieval of the DSD parameters may be improved. This is demonstrated in Fig. 6 (bottom panel), where the agreement on the raindrop size distribution is much better between the different beams. In particular, the retrieval of the intercept parameter becomes meaningful. Note the time dependency of the retrieved parameter D0 (Fig. 6) is similar to the time dependency of the measured reflectivity (Fig. 5).
In this example the trend of the retrieved radial wind is correct but underestimated, which directly decreases the accuracy of the DSD retrievals. In particular, the estimation of the intercept parameter is adversely affected. For vertical profiling radar, Atlas et al. (1973) show that an updraft error of 1 m s−1 produces a large error in particle number density N(D). The same detrimental effect occurs from a radial wind error.
In the case of slant profiling, adding two supplementary beams can overcome this problem because the mean horizontal wind can be measured and be input for the algorithm, as explained previously. Nevertheless, υ0 is still retrieved and varies here between −1.5 and 1.5 m s−1 (see Fig. 4). We can extrapolate that the DSD retrievals may be less accurate in the case of significant updrafts/downdrafts.
b. Comparison of DSD retrievals from nonwind-shifted and wind-shifted Doppler spectra: Vertical profile
Conversely to the horizontal scale, the horizontal wind significantly varies with height (Fig. 8, bottom-right panel). Figure 8 gives an example of the vertical profile of the raindrop size distribution. A light height smoothing has been performed. The top row provides the DSD results without mean horizontal wind correction. One profile of D0 (dotted line, OB2, obtained from the 69° elevation beam) is quite different. Just below (second row), the same retrievals, with correction, are depicted. The three-beam profiles converge to each other. It is clear that the radial wind error is harmful to the D0 estimate. Because they are not directly depending on the raindrops’ mean Doppler velocity, the other parameters, μ and σ0, characterizing the shape of the modeled Doppler power spectrum, are less affected. To obtain the correct reflectivity, which is related to D07, the Nw value drastically changes following the D0 error. The two bottom rows of the figure concern estimates with mean horizontal wind correction.
From the height 1300 to 600 m, the shape parameter approaches 5. The median volume diameter increases, while the number concentration diminishes indicating growth processes via coalescence. The liquid water content () and the reflectivity () increase with D0 until the number concentration decreases too much. Because of the difference in dependency with D0, LWC first starts to diminish (1100 m) shortly, followed by a decrease in reflectivity (800 m). We note that the retrieved spectrum width broadening strongly decreases from 1300 m to reach a constant value approaching 0 m s−1 between 1000 and 600 m. Above 1300 until 1600 m, the measured horizontal wind speed increases from 16 to 26 m s−1 (mean horizontal wind from 18 to 20 m s−1). The retrievals are noisier in this area. Examples of input Doppler power spectra are displayed in Fig. 9. Larger Doppler velocities (from −6 to −8 m s−1) are measured at lower heights and there is a signal decrease at −2 m s−1, which confirms the hypothesis of coalescence.
c. Comparison of rainfall-rate retrievals from nonwind-shifted and wind-shifted Doppler spectra
Finally, mixing all heights and times, the retrieval of rainfall rates is depicted with histograms in Fig. 10. With mean horizontal wind correction, the histograms are similar for the three beams. Without wind correction, a broad histogram with large values is obtained for the main beam. Compared to the vertical beam, the median volume diameter has the tendency to be underestimated, which is compensated by large values of Nw to get the right reflectivity. Being less sensitive to D0 than the reflectivity, the rainfall rate is consequently overestimated. For the other slant beam, OB2, D0 tends to be overestimated, leading to small values of Nw and rainfall-rate underestimation results from this.
6. Impact of the average on the raindrop size distribution retrievals
Probing the rain medium at different heights, times, and directions gives the possibility to study precipitation variability. However, DSD retrievals should be carefully interpreted. Jameson and Kostinski (2001) indicate that the DSD retrievals from radar Doppler spectra measured with large sampling volumes and times are probably samples from statistically inhomogeneous rain. Range and time resolutions can be increased to attempt to characterize rain patches, and this requires estimating the DSD at high resolution. That is done in sections 4 and 5. Hereby, a few consistency checks are made to strengthen the proposition of high-resolution retrievals. For these verifications, only retrievals of wind-shifted Doppler spectra are considered.
a. Quality of the fit
Similar results to those acquired in Figs. 4 and 6 are obtained when using averaged Doppler spectra (over 17.5 s). Despite the fit on Doppler spectra with less statistical fluctuations, the retrievals do not show improvement in RMSD or CV. Nevertheless, the quality of the fit (Rust 2001), given by the coefficient of determination RsZ varying from 0 to 1 (perfect fit), is better:
where Nυ is the number of Doppler bins considered for the fit and the superscript “opt” labels the results of the optimization. The mean and standard deviation of RsZ are (0.81, 0.11) for the two-averaged Doppler spectra and (0.91, 0.07) for the 12-averaged spectra.
Prior to averaging, the Doppler spectra are shifted to have the same mean Doppler velocity to reduce unwanted broadening due to wind variability (Giangrande et al. 2001). However, the DSD is generally not steady because of the rain patchiness. An averaged Doppler spectrum may represent a mixture of DSD.
b. Consistency check: Z–R relations at fixed D0
For identifying statistically homogeneous rain, Jameson and Kostinski (2001) propose to display the reflectivity versus the rainfall rate. When linearity exists between these two moments of the distribution [Eq. (19)], the DSD may be unique and steady,
For Eq. (19), the expressions of Z [Eq. (9)] and R [Eq. (12)] are used. The ratio Z/R increases with D0 and decreases when μ increases at fixed D0 (see Fig. 11). A simulated scatterplot between the reflectivity and rainfall rate is depicted in Fig. 12. The range of retrieved D0 values, from 0.4 to 1.2 mm, is selected. The shape parameter and the intercept parameter data are comprised in the interval [−1, 5] and [102, 105], respectively, in the left panel. The right panels represent the same scatterplot when both μ and Nw are fixed (top) and when only Nw is constant (bottom) to understand the building of the left panel. In the simulation, the maximum reflectivity is constrained to 75 mm6 m−3 (18.8 dBZ), which is the largest measured value. Consequently, the maximum possible median diameter reaches 0.8 mm when the intercept parameter equals 8000 mm−1 m−3. Lower values of Nw are necessary to extend the range of D0 data.
Next, let us consider real data and retrievals. Figure 6 shows that the retrieved DSD gradually varies versus time. The same occurs versus height (Fig. 8). The liquid water content (Fig. 8) and the rainfall rate (Fig. 10) exhibit small but clear variations. The precipitation may be statistically inhomogeneous at the space and time scales considered (height interval [200, 1750] m and a few minutes, respectively). For the whole raw data measurement (2 × 3.5 min), the measured reflectivity is plotted against the retrieved rainfall rate in Fig. 13, considering different intervals of retrieved D0 from 2-averaged (left panel) and 12-averaged Doppler spectra (right panel). The different Z/R ratios resulting from different D0 intervals are clearly visible and they show an increase with D0 as expected. The spread of the Z/R ratio at fixed D0 is caused by the μ values. In addition, the different Nw values lead to more points on the linear relations toward small and large values of rainfall rates. The rationale for this plot is a consistency check of the retrievals. The left panels display the elementary drop size distributions leading to linearity between Z and R. They are acquired at high time and spatial resolutions. With more averaging the linearity between Z and R is conserved, but the largest data of rainfall rate and reflectivity disappear (see also rainfall rates histograms in Fig. 14). Note how height smoothing of the retrievals can distort the linearity, especially for the 12-averages case. The relations Z–R bend for the large values of Z and R in the two-averages case. This evaluation suggests that no smoothing should be carried out for statistical comparisons although the variances of the retrievals are large.
c. Consistency check: Retrieval distributions
Histograms of the retrieval data, not smoothed, are shown in Fig. 14. Only the averaging number of Doppler power spectra differs. With more averaging, the median volume diameter distribution shifts toward higher values, which leads to lower intercept parameter values. Therefore, the rainfall rate decreases.
d. Concluding remarks
These first results on light rain indicate that high-resolution DSD retrievals with satisfactory performance are possible. Even with careful averaging, as soon as the Doppler spectrum width varies, the resulting averaged spectrum will be broadened, which reduces the number of small D0 values. Consequently, averaging of the Doppler power spectra is not recommended when the fit quality is acceptable. The variances of the retrievals do not decrease anyway. Finally, smoothing of the retrievals may adversely affect the interpretation of the results and its use should be restricted to visualization when necessary.
In the case of slant profiling of precipitation when Zdr is about 0 dB, samples of the raindrop size distribution can still be retrieved from the Doppler spectra using the described technique. This technique is automated for study cases. Simulation and comparison on measured Doppler spectra of the same radar resolution volume yield the estimation errors of this method. The median volume diameter and the shape parameter are well retrieved, but the error on the intercept parameter may be large. A way to evaluate the DSD results is to probe stratiform rain in three different directions and to investigate the time series of retrievals at low heights when the three radar resolution volumes are near each other. We have shown a case where the retrieved radial wind is underestimated, which adversely affects all the DSD parameters and in particular the intercept parameter. Adding two extra beams mitigates this effect because we can input the mean horizontal wind in the retrieval algorithm. Note that homogeneity of only the mean horizontal wind is assumed across the three beams. The estimation of the radial wind is still carried out to account for the vertical wind and the variability of the horizontal wind. The obtained raindrop size distributions are considerably improved at all heights. The same conclusion is obtained from the distributions of rainfall rates for the three beams and from performing the analysis on a vertical profile of DSD where the three radar resolution volumes become distant and the horizontal wind varies with height. We can speculate that an important updraft/downdraft will have a detrimental influence on the raindrop size distribution calculations.
For a consistency check, the measured reflectivity is related to the retrieved rainfall rate at fixed D0 for 2-averaged and 12-averaged Doppler spectra, wind corrected. The relations show the expected linearity in both cases. Noteworthy is the degradation of the linearity when the retrievals are smoothed for the case of 12 averages. The corresponding distributions of retrievals exhibit differences that are attributed to the averaging. Careful averaging of the Doppler spectra can still lead to broadening of the resulting Doppler spectrum when their Doppler widths vary. Consequently, the number of small D0 values decreases, which diminishes the number of large Nw values and leads to a smaller number of large rainfall rates.
Based on this study case of stratiform light rain, this technique gives satisfactory performances for the retrieval of raindrop size distributions comparing vertical and slant profiling. Furthermore, the method can be directly applied on the Doppler spectra without strong averaging, giving the possibility to obtain high-resolution DSD retrieval samples when the rain is not homogeneous.
The author gratefully acknowledges the contributions of the reviewers; Herman Russchenberg’s support; and the COPS team, Yann Dufournet and Christine Brandau.
Iterative Cascaded Retrieval Algorithm
A three-parameter nonlinear least squares optimization is carried out on the median volume diameter D0, the spectral broadening σ0, and the shape factor μ using a cascaded approach. The iterative optimization procedure is divided into five steps (see also Fig. 2).
Iterative selection procedure
The values of D0, σ0, and μ are selected for consecutive optimizations. The values of D0, σ0, and μ are bounded.
Boundaries for the search of the median volume diameter D0 are estimated to avoid exotic fitting when Zdr values are about 0 dB. Scattering computations have been carried out using as input a set of 73 017 DSD that displays the natural variability of rain by varying the drop size distribution parameters as 2 ≤ log[Nw (mm−1 m−3)] ≤ 5, 0.109 mm ≤ D0 ≤ 3.5 mm, and −1 ≤ μ ≤5. The obtained reflectivity values have been classified into intervals of 10 dBZ, which correspond to intervals of D0. For example, the D0 search intervals are [0.3, 0.9] and [0.4, 1.25] mm for reflectivity values ranging from 0 to 10 and from 10 to 20 dBZ, respectively.
Consequently, based on the reflectivity value, the parameter D0 varies from Dmin(Z) to Dmax(Z) (ND values), μ varies from −1 to 5 (Nμ values), and σ0 varies from 0 to 1 m s−1 (Nσ values).
Estimation of the intercept parameter Nw uses Eq. (14). The retrieval of Nw depends on D0, σ0, and μ. Next, the ambient wind velocity υ0 is calculated as the lag of the cross correlation of the measured and the modeled spectrum . Hence, the retrieval of υ0 depends on D0, σ0, μ, and Nw.
Output of the model
A set of five parameters is obtained with the corresponding modeled spectrum.
Optimization procedure from cost function study
Optimization on μ
The cost function for the optimization of the shape factor is given by
In practice, steps 2 and 3 are repeated Nμ times for each value of μ. The terms D0 and σ0 are fixed. We obtain thus Nμ values of Nw and υ0, which depend on μ. The minimization of L(μ) is carried out using the obtained Nμ values of Nw and υ0. The result is μopt (σ0, D0), Nw (μopt, σ0, D0), and υ0 (Nw, μopt, σ0, D0) for fixed D0 and σ0.
Optimization on σ0
Steps 2 and 3 and step 4a are repeated using the entire range of values of σ0. Here D0 is fixed. We obtain thus a set of Nσ values of μopt (σ0, D0), Nw (μopt, σ0, D0), and υ0 (Nw, μopt, σ0, D0). The cost function
is minimized using the set of Nσ values. The result is σ0opt (D0), μopt (σ0opt, D0), Nw(μopt, σ0opt, D0) and υ0 (Nw, μopt, σ0opt, D0) for fixed D0.
Optimization on D0
Steps 2 and 3 and steps 4a and 4b are repeated using the entire range of values of D0. We obtain thus a set of ND values of σ0opt (D0), μopt (σ0opt, D0), Nw (μopt, σ0opt, D0), and υ0 (Nw, μopt, σ0opt, D0). The cost function
is minimized using the set of ND values. The solution of the cascaded retrieval algorithm is D0opt, σ0opt (D0opt), μopt (σ0opt, D0opt), Nw (μopt, σ0opt, D0opt), and υ0 (Nw, μopt, σ0opt, D0opt).
The solution consisting of the optimized five parameters, the output of step 4c, and the final fit between the measured and modeled radar observables is obtained.