The Precipitation Occurrence Sensor System (POSS) is a small Doppler radar originally designed by the Meteorological Service of Canada (MSC) to report the occurrence, type, and intensity of precipitation in automated observing stations. It is also used for real-time estimation of raindrop size distributions (DSDs). From the DSD, various rainfall parameters can be calculated and relationships established, such as between the radar reflectivity factor (Z) and the rainfall rate (R). Earlier work presented first-order estimates of the sampling errors for some POSS rainfall parameter estimates. This work combines a Monte Carlo simulation and “inverse problem” analysis to better estimate errors due to the specific sampling problems of this disdrometer type. The uncertainties are necessary to determine the statistical significance of differences between DSD estimates by the POSS and other collocated disdrometers, or between POSS measurements in different climatologies. Additionally, confidence limits can be assigned to regression coefficients for rainfall parameter relationships determined from POSS estimates. An example is given of the uncertainties in the coefficients of measured Z–R relationships.
The Precipitation Occurrence Sensor System (POSS) is a bistatic, X-band, continuous wave (CW) radar (Sheppard 1990, hereafter S90) originally designed by the Meteorological Service of Canada (MSC) to report the occurrence, type, and intensity of precipitation for use in automated observing stations (Sheppard and Joe 2000). Researchers have also used this instrument as a disdrometer (e.g., Campos and Zawadzki 2000). The term disdrometer is used here to refer to any instrument that estimates a raindrop size distribution (DSD), that is, the number of raindrops per unit volume of air in a specified diameter range. From the DSD, precipitation parameters such as rainfall rate and the radar reflectivity factor can be calculated and relationships established. This paper determines the uncertainty due to sampling statistics in the POSS estimates of DSDs and derived parameters.
In principle, there are two types of disdrometers: the surface flux and the volumetric types. One example of a volumetric type is a raindrop camera developed by the Illinois State Water Survey (Jones 1992). Most disdrometers are of the surface flux type; that is, they estimate the number of drops crossing a defined surface area in a given time interval. Some count the number of impacts on a surface (e.g., Joss and Waldvogel 1967) and some image drops as they cross an optical surface (e.g., Kruger and Krajewski 2002). In surface flux–type disdrometers it is necessary to assume a relationship between terminal velocity and raindrop size in order to convert the flux measurement to a volumetric number concentration.
Radars, like POSS, are basically a volumetric type because they estimate the DSD from the characteristics of the scattered power received from drops in a defined measurement volume. As discussed in section 3, for POSS, the size and shape of this volume depends on the raindrop’s diameter and velocity vector.
There have been previous studies of the statistics of sampling raindrops using disdrometers (Joss and Waldvogel 1969; Gertzman and Atlas 1977; Smith et al. 1993). The assumption is usually made that the measurement technique of the instrument is without error, and that any variability in the measurement is due to sampling variability and inhomogeneity in the precipitation. In general, when modeling the sampling statistics of a disdrometer, a DSD model is assumed (e.g., Marshall and Palmer 1948, hereafter MP48), to determine the mean number of drops in each sample volume. A statistical distribution, such as Poisson, is also assumed for the variations in the number of drops from sample to sample. More complex statistical models assume that higher-order moments, such as the mean value, are also Poisson distributed (e.g., Kostinski and Jameson 1997).
A distinction between the present work and previous POSS studies is the approach used to represent the sampling error. Previous work used a “forward model” to determine the conditional probability distribution of the sensor estimate for a given population model DSD. Here, an “inverse problem” analysis estimates the conditional probability distribution of a population parameter for a given sensor estimate (measurement). Examples of forward versus inverse models are given in the appendix.
Compared to many disdrometers, the sampling volume of POSS is several orders of magnitude larger. Therefore, the range of raindrop number concentrations, for which sampling errors due to Poisson variability are significant, is much less for POSS than for other disdrometers with smaller measurement volumes. However, POSS has a different specific sampling problem caused by the large variability of the Doppler signal received from a scatterer of a fixed cross section at different locations in the measurement volume. A first-order estimate of this sampling error was given in Figs. 12 and 13 of S90, and Fig. 2 of Sheppard and Joe (1994, hereafter SJ94). Those analyses assumed that each drop size produced power at a single Doppler frequency corresponding to its terminal velocity. In this paper, a Monte Carlo simulation incorporates the effect of sampling on the Doppler spectral distribution of power, which results in a more realistic estimate of the measurement uncertainties.
2. Instrument description
The transmitter and receiver of the bistatic X-band radar are housed separately and mounted on a frame so that the centers of the antenna radomes are about 25 cm apart. The antennas are oriented 20° from the vertical so that the beam axes intersect midway between them, and about 34 cm above the horizontal plane through the center points of their radomes.
The radar characteristics and signal processing parameters are given in Table 1. The sensor’s processing unit outputs a real-time 1-min-average Doppler velocity spectrum. In 2003, the POSS real-time processing algorithms were upgraded so that the number of spectra in the 1-min average increased from 380 to 960.
3. Modeling the Doppler signal
A computer model is used to simulate the Doppler signal produced by any hydrometeor of known size, shape, refractive index, and velocity vector while traversing the measurement volume during the 64 m s−1 sampling window. A Cartesian coordinate system is defined with origin midway between the centers of the radomes of the transmitter and receiver. The x axis passes through the radome centers; the y axis is in the horizontal plane and the z axis in the zenith direction.
The characteristics of the POSS radar signal are different from those of large weather surveillance radars. In the latter, the range gates and antenna beam pattern define the measurement volume boundaries. Because the POSS is a CW radar, there is no ranging. The combined transmitter and receiver antenna beam patterns, and the scattering cross section and velocity vector of the hydrometeor, define the size and shape of the measurement volume. The perimeter of the volume is defined by the location of the scatterer at the start of a 64 m s−1 sampling window if the received power during the sampling window is detectable, that is, the amplitude of the Doppler voltage exceeds the threshold corresponding to the ½ bit resolution of the analog-to-digital converter used by POSS. Therefore, not only the scattering cross section, but also the velocity vector, determines the measurement volume’s size and shape. For liquid hydrometeors, with no wind, Fig. 1 shows the size of the measurement volume for a single 64 m s−1 sampling window as a function of the equivolumetric drop diameter. The largest diameter raindrops are detected at a maximum distance of about 3 m above the sensor when there is no wind. In general the measurement volume includes both the near and far fields of the antennas.
In the far field, the amplitude of the Doppler signal is calculated from the bistatic radar equation for a single scatterer, assuming no atmospheric transmission losses, adapted from Skolnik (1980):
where and are the average transmitted and received powers, respectively; Lw is the combined transmission factor for both radomes; Gt(θi, ϕi) is the transmitter antenna gain in the direction of the scatterer at the incident zenith and azimuth angles (θi,ϕi), respectively; Gr(π − θs, π − ϕs) is the receiver antenna gain in the direction of the scatterer at the scattered zenith and azimuth angles (θs, ϕs), respectively; λ is the transmitted wavelength; σ is the scattering cross section of a single particle; and Rt and Rr are the distances from the scatterer to the transmitter and receiver antenna phase centers, respectively.
A far-field theoretical model of Sletten (1988) for smooth-walled pyramidal horns is used to determine the antenna gains in (1). In the near field, a laboratory calibration procedure (S90) is used. Distilled water drops of known size, and therefore scattering cross section, are released from a hypodermic needle so that they fall through the POSS measurement volume at known locations in the near field. The resultant Doppler signal is used to calculate the combined radar system constants in the bistatic radar equation (1) that are independent of the scattering cross section. In this paper the average of near-field measurements from 11 different POSS instruments is used in the volume with coordinates between −30 cm <= x <= 30 cm, −30 cm <= y <= 30 cm, and z < 115 cm. Beyond this range, the far field model is used in (1). Figure 2 represents the measured received power as a grayscale in dB in the vertical x–z plane. The maximum 0-dB power occurs at coordinates (0, 0, 28) cm.
A second important difference from weather surveillance radars is that a single raindrop falling at constant velocity generates a range of Doppler frequencies and amplitudes. Because of the proximity of the drop to the transmitter and receiver, the Doppler velocity component and amplitude change significantly during the sampling window, which results in spectral spreading not seen in larger Doppler radars.
The Doppler frequency generated at any point in the volume is determined by the rate at which the pathlength from the transmitter antenna phase center, to the scatterer, and back to the receiver antenna phase center, changes by one wavelength. Surfaces of constant scattering pathlength define ellipsoids with foci at the phase centers of the antennas. These are represented in Fig. 2 by the dashed lines.
Figure 2 also shows isopleths of the ratio of the Doppler velocity measured by the POSS, to the velocity of the scatterer when it is falling vertically (i.e., no horizontal wind). If there is also no vertical wind, this is the terminal velocity. The terminal velocity model used here is from Beard (1977) and is a ninth-order polynomial fit between velocity and raindrop diameter. The POSS Doppler velocity measurement always underestimates the terminal velocity when there is no vertical wind. The underestimation increases the farther the trajectory is away from central vertical z axis and also as the scatterer falls along any given vertical trajectory (i.e., for constant x and y, and decreasing z).
The scattering cross section, σ in (1), is calculated using the T-matrix method as implemented in software from Mishchenko (2000). All of the preceding theory is general for any type of scatterer, but from here on, this paper will deal exclusively with sampling errors for liquid scatterers. Falling raindrops distort from spheres due to aerodynamic effects. For this work we have assumed an oblate shape at terminal velocity, with axial ratio given by Brandes et al. (2002):
where r is the ratio of the axis of rotational symmetry to horizontal axis, and D is the equivolumetric diameter in millimeters. The T-matrix scattering theory allows the calculation of the scattering cross section of nonspherical raindrops at different incident and scattering angles in (1).
From the models of this section, a simulated Doppler signal can be generated for any raindrop size traversing any path through the measurement volume (e.g., see in Fig. 2 the signal generated by a 2.5-mm drop falling vertically at x = 0 and y = 0).
In summary, this simulation differs from earlier work (SJ94) in four aspects: the near-field radar system constants are determined from measurements rather than assuming far-field theory; the measurement volume size is determined using an analog-to-digital converter detection threshold of ½ bit instead of 1 bit; the scattering cross section is calculated using the T-matrix theory rather than using Mie theory, and at each scattering angle rather than a fixed “weighted average” scattering angle for the volume; and the drops are assumed to be oblate spheroids rather than spheres.
4. Calibration for DSD estimation
In operation, because the radar is CW, the location of scatterers detected in the volume is unknown. The calibration must be based on the mean of the population of Doppler spectra, generated by a single drop size at all locations in the measurement volume.
The measured Doppler power density spectrum S( f ), at frequency f, from an ensemble of drops of different equivolumetric diameters D (mm), is
where N(D) (m−3 mm−1) is the number of the drops per unit volume of air, per unit diameter size interval; V(D) (m3) is the measurement volume; S( f, D) is the volume-averaged Doppler power density at frequency f from a single drop of diameter D. The minimum detectable raindrop diameter is Dmin = 0.34 mm and the maximum raindrop diameter present in the DSD is Dmax.
The products V(D) S( f, D) are referred to as “weighting functions.” They represent the fundamental calibration values necessary to invert the measured S( f ) in order to estimate N(D) [S90, his Eq. (7)]. The weighting functions are not sensitive to small model errors in the estimation of the exact size of the measurement volume because these are compensated for by an approximately equal reciprocal error in the volume-averaged Doppler spectrum such that the product of these factors remains nearly constant.
5. Sampling statistics analysis
In earlier work, the sampling errors were analyzed parametrically by assuming that a raindrop of a given diameter generated power at only one Doppler frequency component (i.e., that the weighting functions at each Doppler frequency were zero for all diameters but one; see S90, his section 5c and SJ94, their section 3). This approach ignored the effect of the Doppler spectral spreading for each drop size on the sampling error.
In this work the sampling errors are evaluated using Monte Carlo sampling of the population composed of the Doppler spectra associated with each location in the measurement volume and for the diameter of the midpoint of each measurement channel. Each Doppler spectrum is calculated from the simulated signal generated by a drop falling vertically, using the signal processing parameters given in Table 1. These are the same parameters used in the sensor’s real-time processing software. A separate spectrum is found for each 64 m s−1 sampling window starting at different points spanning the measurement volume with a resolution of 5 cm in the horizontal and 1 cm in the vertical. Each trajectory is assumed to have an equal probability of occurring.
The sampling error in the estimation of the mean value of any rainfall parameter x, can be represented by the conditional probability density function (pdf) of the population value of x, given a measured Doppler spectrum S( f ). The measured S( f ) can always be inverted to give an estimate of the DSD (section 4), from which an estimate of x can be calculated. For simplicity of notation, will be written as x̂. The sampling statistics can therefore be represented by the conditional pdf p(x|x̂) instead of p[x|S( f )].
It is relatively straightforward to determine the sampling distribution of POSS estimates x̂ given a population mean x as described in steps 1–7 below. This distribution p(x̂|x) is referred to as the “forward” problem in this paper. However, the distribution of interest is the “inverse” problem p(x|x̂); that is, what is the distribution of population means x that might generate a given POSS estimate x̂? Tarantola (2005) gives a comprehensive study of inverse problem theory. The simple approach used here is based on the multiplicative law of probability:
where p(x̂) = Σx p(x̂, x) and p(x) = Σx̂ p(x̂, x) are the marginal pdfs of the joint pdf p(x̂, x).
To estimate the pdfs in (4), the following simulation steps are performed.
Assume a DSD model and calculate a corresponding mean rainfall parameter xi, where the index i represents the DSD model used. This analysis uses a Marshall–Palmer (MP) DSD model with different precipitation rates (index i). Here xi can be the number concentration or any derived parameter (e.g., rainfall rate, radar reflectivity factor, etc.).
For the diameters associated with the midpoint of each of the 34 POSS channels, calculate the mean number of drops sampled for a given number of measurements in the averaging period. This requires determination of the measurement volume as a function of diameter (Fig. 1). The sensor’s current real-time processing speed is 960 Doppler spectra in 1 min. Figure 3 shows contours of the log of the number of drops in each channel, sampled in 960 measurements, as a function of the midchannel diameter and the log number concentration per unit diameter size interval. The contours are valid for any DSD shape; the MP DSDs for different precipitation rates are superimposed only for reference.
Assume, for a stationary DSD, that the number of drops sampled in the averaging period has a Poisson distribution (Joss and Waldvogel 1969) with mean equal to that found in step 2. Randomly sample this Poisson distribution to determine the total number of drops in each channel in the averaging period.
For each drop in each POSS diameter channel during the averaging period, randomly select a starting location in the measurement volume for the 64 m s−1 time domain signal, and then determine the corresponding Doppler velocity spectrum. The seed for the random number generator was reinitialized before each new MP rate in step 1.
Calculate a composite Doppler spectrum from the individual spectra found in step 4.
Invert the composite Doppler spectrum using S90 [his Eq. (7)], to give the corresponding estimate of the DSD. Figure 4 shows 10 estimates of DSD assuming an MP population model in step 1 with a rate of 1 mm h−1.
Calculate x̂ from the DSD estimate in step 6. Note that all estimations are in the linear domain. Because of the large dynamic range some parameters are then transformed to the log domain. The estimates and population means are then written log x̂ and log x, respectively. Note that all logarithms represented in this paper by log are to the base 10. Any natural logarithms are represented by ln.
Repeat steps 3–7 for a large number of simulations.
Populate a two-way contingency table of the number of occurrences (nij) of x̂ in the jth subrange of x̂ values for a population value xi. Figure 5 uses interpolated contour lines to represent the values in the contingency table for the example of x representing the MP log number concentration of 3.14-mm-diameter drops and the number of spectra in the average is 160.
Return to step 1 and choose a different rate for the MP DSD model with a corresponding population mean value of xi+1. Repeat steps 1–9 until an appropriate range of x values is spanned. A total of 90 rates from 0.04 to 160 mm h−1 is used here.
From the contingency table, estimate the conditional pdf p(x̂|x). The bottom panel in Fig. 5 is an example of this pdf for x = logN = −2 for a 3.14-mm-diameter drop and the number of spectra in the average is 160. The marginal pdfs p(x̂) and p(x), used in (4), are calculated from the contingency table’s columns and rows, respectively. This simulation assumes that the 90 MP rates used in step 1 are equally probable and therefore p(x) = 1/90. In fact, p(x) would depend on the climatology of the site. This is discussed in more detail in section 6.
6. Results and discussion
Sampling statistics are presented for the following rainfall parameters: number concentration, rainfall rate, radar reflectivity factor, differential reflectivity, and specific differential phase in sections 6a–e, respectively. The simulation in section 5 was run for different numbers of Doppler velocity spectra in the averaging period. The speed of the POSS processor prior to 2003 allowed a maximum of 380 spectra in the 1-min average. The current version has 960 spectra in a 1-min average. Because of interest in using POSS to study DSD evolution with higher temporal resolution (e.g., 10s), an averaging number of 160 spectra was also analyzed.
The previous POSS evaluation studied only the forward problem [i.e., the conditional pdf p(x̂|x)], whereas this paper considers the inverse pdf p(x|x̂). Other papers studying the sampling statistics of disdrometers have also only dealt with the forward problem (e.g., Smith et al. 1993). The difference between forward and inverse conditional pdfs can become significant if the forward distribution is strongly skewed. Disdrometers that are discrete raindrop counters with uniform sensitivity to drops at all locations in their sampling volume are subject to Poisson variability only. (In the appendix an example is provided in Fig. A1 comparing the forward versus inverse conditional Poisson probability distributions.)
The statistics used to represent the forward pdf of the estimated DSD parameters in previous analyses were the mode and standard deviation (S90) and the mode and quartiles (SJ94). Here the inverse pdf is characterized by the median and interquartile difference for the log number concentrations, and by the mean and standard deviation for the other parameters derived from the DSD. The difference between the POSS estimate of number concentration and the median of the inverse pdf, or in the case of the derived rainfall parameters, between the estimate and the mean of the inverse pdf, is referred to here as the “bias” in the estimate.
The DSD population was assumed to follow the MP model in step 1 of the simulation in section 5. The MP model was used because its formulation provides a simple way to generate a wide range of raindrop number concentrations in a realistic manner. The choice of model is unimportant in evaluating sampling errors for number concentration provided that it is continuous and moderately “smooth.” This was verified by repeating the simulation using a modified gamma distribution (Ulbrich 1983).
The sampling error of the rainfall parameters that are integrated over the DSD will depend on its shape. The results for an MP DSD are presented in this section. For other shapes of DSD, the sampling error will depend on the relative contribution of each diameter channel to the integral parameter.
The effect of horizontal wind on the sampling errors was also analyzed. The simulation described in section 5 was repeated, but in step 4 the trajectory of the raindrop falling through the volume was determined from the resultant of the terminal and wind velocity vectors for each drop size. Also, in step 6, the Doppler spectra were inverted using weighting functions appropriate to the simulated wind speed. Because the resultant speed of the drops is always greater than without horizontal wind, more drops are sampled in a given averaging period. This results in a reduction of the inverse pdf dispersion for all rainfall parameters. Therefore, the results of this section, which have assumed no wind, represent the maximum limits expected.
As mentioned in step 11 of the simulation in section 5, the assumption is made that all rainfall rates are equally probable [i.e., p(R) is uniform]. In reality p(R) will depend on the climatology of the location. For the MP DSD model, the pdf p(x) of any derived parameter x can be determined if the distribution of rainfall rate p(R) is known. A lognormal pdf has been shown to be a realistic model for some situations (e.g., Cho et al. 2004). A frequency of occurrence histogram of 1-min-average POSS estimates of rainfall rate, measured at Egbert, Ontario, Canada, during 1998, was well represented by a lognormal distribution with a mean and standard deviation of the logarithm (base 10) of the rate of −0.56 and 0.79, respectively. The effect of using a nonuniform versus uniform pdf p(x) on the inverse pdf calculation in (4) depends on the magnitude of the change in p(x) over the nonnegligible range of p(x|x̂). If the dispersion of p(x|x̂) is small enough that variations in p(x) are insignificant, then the results will be similar to the uniform case presented below. The dispersion of p(x|x̂) depends on the parameter x and the number of spectra in the average as will be shown in the following sections.
Last, apart from uncertainties caused by sampling effects, there is another uncertainty introduced by the algorithm used to invert the Doppler spectra to DSD [S90, his Eq. (7)]. Although for most parameters this uncertainty is small, it is worth noting that it is included as part of the sampling errors presented here. The inversion algorithm bias (x̂ − x), where x is any rainfall parameter, was quantified using the same procedure described in section 5, but with the following modifications. First, steps 2 to 4 are omitted. Second, step 5 is replaced by a calculation of the “population mean” Doppler spectrum using (3). This spectrum was then inverted to produce an estimate x̂ of the population x used in step 1.
Further investigation determined that the magnitude of the bias depends on the convergence criterion used to terminate the iterations of the general relaxation inversion algorithm. The convergence criterion is presently based on the fractional change in the estimated precipitation rate from the last iteration. In this paper a value of 0.001 is used. If convergence is not achieved in 1000 iterations the process is terminated. It was beyond the scope of this paper to investigate different inversion algorithms and convergence criteria.
a. Number concentrations
The conditional pdf p(logN|logN̂) (see the example for logN̂ = −2 in the right-hand panel in Fig. 5) is used to determine percentile statistics for each POSS diameter channel. Figures 6a–c present contours of the difference of logN̂ minus the median of p(logN|logN̂), as a function of diameter D and logN̂, for different numbers of spectra in the averaging period. As expected, the underestimation of the population median decreases at a given diameter as both the number concentration and the number of spectra averaged increase. Using Figs. 6a–c, the bias in POSS estimates could be corrected in a statistical sense.
As mentioned above, p(logN|logN̂) will depend on the assumed climatological distribution of logN. The lognormal pdf p(logR) described above was converted to the corresponding pdf p(logN) using the MP48 relationship. The inverse conditional pdf p(logN|logN̂) was recalculated. The result of the lognormal weighting of logN is a decrease in the median of p(logN|logN̂). Figure 6d shows the corresponding reduction of the underestimation bias compared to the uniform pdf used in Fig. 6a. Both analyses use 160 spectra in the average. For estimates near logN̂ = 1 and diameters in the range of 3–4 mm in Fig. 6d, the POSS overestimates the median of the population. The effect is less as the number of spectra in the 1-min-average increases (not shown here).
Figure 7 represents contours of the interquartile range between the third and first quartiles of p(logN|logN̂) as a function of D and logN̂, for 960 spectra in the average. This represents a measure of the dispersion of the population pdf that could generate a given POSS estimate. The larger the contour value, the greater the uncertainty in the population value. The uncertainty increases with decreasing number concentration and decreasing number of spectra in the average (not shown here).
Figure 8 plots the inversion algorithm bias (logN̂ − logN) as a function of D and logN̂ over the range of MP rates from 0.1 to 128 mm h−1. The region indicated by the unpatterned area between the ±0.02 contours represents the part of estimation space where the agreement is better than 5% (approximately antilog 0.02). The agreement degrades at the extremes of the diameter ranges. At the smallest diameter classes, logN is increasingly overestimated with rate. This is the source of the overestimation at small diameters seen in Fig. 6. For the largest diameter class at 5.34 mm, logN is overestimated at low rates and underestimated at high rates.
b. Rainfall rate
The rainfall rate R (mm h−1) is defined by
where vT (D) is the terminal velocity (m s−1) of a drop of equivolumetric diameter D (mm), and N(D) is the number concentration (m−3 mm−1) per unit diameter, and Dmin and Dmax are the minimum and maximum raindrop diameter in the DSD, respectively.
Because of the large dynamic range of naturally occurring rainfall rates, the results are presented in the logR domain. Similar to Fig. 5 for number concentration, for rainfall rate Fig. 9 shows interpolated contour lines to represent the number of simulation events in the contingency table, described in section 5 (step 9). The population values (logR) are represented on the y axis and the POSS estimates (logR̂) on the x axis for 960 spectra in the average. As an illustrative example, the axes are limited to the ranges in which p(logR̂|logR = 0) (bottom panel) and p(logR|logR̂ = 0) (right-hand panel) are nonnegligible. As before, the assumption is made when applying (4) that the pdf p(logR) is uniform. The inverse pdf shows in a qualitative sense the degree to which the sampling errors are normally distributed for a POSS rate estimate of 1 mm h−1.
Figure 10 gives the difference logR̂ minus the mean of the pdf p(logR|logR̂), referred to as the fractional bias when transformed by 10(logR̂-mean) − 1. The three curves represent the statistics for different numbers of spectra (k) in the averaging period. The estimation of the mean is insignificantly affected by k. The mean is overestimated by about 0.5% at 0.1 mm h−1 to about 1.5% at 100 mm h−1.
Also calculated is the standard deviation (SD) of the pdf p(logR|logR̂). This is transformed to a fractional value by 10(SD) − 1. As expected the standard deviation is approximately proportional to k−0.5. Also the higher the rate, the lower the standard deviation, because there are more drops present. The results can be compared to Fig. 13 of S90. At that time the processing speed was one spectrum per second. The standard deviation curve for 500 samples in the average decreased from about 8% to 4% over the range of rates from 0.1 to 100 mm h−1, respectively. In the present work the corresponding range for 500 samples (assuming the k−0.5 relationship) is about 11%–6%.
Using the lognormal pdf for logR described above, the inverse conditional pdf p(logR|logR̂) was recalculated from (4). The slope of the bias curves approximately doubled compared to those in Fig. 10. The effect on the standard deviation was insignificant.
As was done for the estimation of DSDs in Fig. 8, the fractional bias in the estimation of rainfall rate due to the inversion algorithm was calculated. The fractional bias increases from about 0.0% for estimates of 0.1 mm h−1 to 1.4% at 100 mm h−1. At 100 mm h−1, the inversion algorithm bias accounts for the entire sampling bias in Fig. 10.
c. Radar reflectivity factor
The radar reflectivity factor Z is defined by
and is calculated directly from the POSS-estimated DSD. Because of the large range of Z values, dBZ units are used by applying the transformation 10logZ. In this paper, unless stated otherwise, the parameter Z is in units of dBZ. Figure 11 gives the difference of Ẑ minus the mean of the pdf p(Z|Ẑ), as a function of Ẑ, where Z is the mean of the population, and Ẑ is the POSS estimate of Z. The simulation is shown for different numbers of Doppler spectra in the average. The underestimation of the mean decreases as the number of spectra in the averaging period increases. For the three averaging numbers analyzed, the maximum underestimation is about 0.25 dBZ at 10 dBZ and approaches agreement as Z increases.
Also calculated is the standard deviation in dBZ of p(Z|Ẑ). As for rate, the standard deviation of Z is approximately proportional to k−0.5. For the fastest processor speed, using a 1-min average, the SD ranges from about 0.9 dBZ at 10 dBZ to about 0.3 dBZ at 50 dBZ.
Figure 11a assumes that p(Z) is uniform in the simulation of section 5. Figure 11b uses the pdf p(Z) determined from the lognormal pdf p(logR) given above, again assuming the MP DSD. A comparison of Figs. 11a,b shows that their respective bias and standard deviation curves are similar until higher Z values are reached. Also the effect is greatest when the number of spectra in the average is least.
The inversion algorithm bias for Z is insignificant, fluctuating irregularly from −0.05 to +0.05 dBZ over the range from 10 to 50 dBZ.
d. Differential reflectivity
The differential reflectivity ZDR, in units of decibels, is defined by
where ZH and ZV are the radar reflectivity factors for horizontal and vertical polarization, respectively, and in this analysis are calculated using the T-matrix scattering cross section formulation for a raindrop axial ratio given in (2), and for C-band radar wavelengths. Figure 12 gives the difference of ẐDR minus the mean of the pdf p(ZDR|ẐDR), as a function of ẐDR, where ZDR is the mean of the population ZDR, and ẐDR is the POSS estimate of ZDR.
The underestimation of the population mean decreases as the number of spectra in the averaging period increases. For all averaging numbers the bias becomes more negative (underestimation) initially as ẐDR increases. At about 1 dB the negative bias is greatest and then becomes less negative as the model rate increases. This is because the major contribution to ZDR comes from the large diameter drops. At very low population model rainfall rates there are so few large drops that the value of ZDR is small, and undersampling is insignificant. As the model rate increases the number of large drops increases and the undersampling becomes important. At rates above 5 mm h−1 the number of large drops increases sufficiently so that undersampling is less significant and the underestimation of ZDR decreases.
Similarly, the standard deviation decreases as the number of spectra in the averaging period increases. For each averaging period, it increases with the model rate, reaching a maximum value at higher rates and then decreasing.
As described in section 6c for Z, the pdf p(ZDR) was determined from the lognormal pdf for rate. The use of this pdf in (4) has a significant effect on the biases and standard deviations as seen from a comparison of the results in Fig. 12b to those in Fig. 12a, which assumed a uniform pdf. The biases become more positive and standard deviation increases as ẐDR increases, if the lognormal pdf for rate is assumed.
The inversion algorithm bias for ZDR was relatively insignificant (−0.004 to +0.01) over the range presented in Fig. 12.
e. Specific differential phase
The specific differential phase KDP, in units of degrees per kilometer, is defined by Bringi et al. (1990) as
where λ is the wavelength, and Re[ fH(D) − fV(D)] is the difference of the real components of the forward-scattering amplitudes for horizontal and vertical polarized waves, calculated in this analysis using the T-matrix formulation for a raindrop of equivolumetric diameter D and axial ratio given in (2), and for C-band radar wavelengths.
Figure 13 gives the difference of K̂DP minus the mean of the pdf p(KDP|K̂DP), as a function of K̂DP, where KDP is the mean of the population KDP, and K̂DP is the POSS estimate of KDP. The underestimation of the population mean decreases as the number of spectra in the averaging period increases. For all averaging numbers the underestimation of the population mean increases as K̂DP increases.
Similarly, the standard deviation decreases as the number of spectra in the averaging period increases. For each averaging period, it increases as K̂DP increases.
Following the analysis of the other integral rainfall parameters, a comparison of Figs. 13a,b shows the effect of using a uniform pdf versus a lognormal pdf for rate to determine p(KDP). The standard deviation curves are similar but the biases become more positive, if the lognormal pdf for the rate is assumed. Also, the effect is greatest when the number of spectra in the average is least.
Below about 6° km−1 the inversion algorithm bias is insignificant, ranging from −0.005° km−1 to +0.005° km−1. Above this value the inversion bias increases as the contribution of the large diameter drops to KDP increases. Figure 8 shows that the inversion bias for the number concentrations of large drops causes the estimate to be about 40% (approximately antilog −0.4) of the population mean at model rainfall rates greater than 100 mm h−1. The overestimation of large drop number concentrations at the lowest model rainfall rates is not important because the population mean number concentrations are very small.
f. Example of sampling statistics applied to Z–R relationships
The POSS has been used to study the effect of instrumental uncertainties on Z–R relationships (e.g., Campos and Zawadzki 2000). The conventional relationship
transformed to dBZ units using
is used for the purpose of illustrating an application of the sampling errors presented here.
POSS 1-min-average Doppler spectra (380 spectra in each) were analyzed from May to August 1998, from Egbert. Minutely DSDs were estimated and rainfall rate and radar reflectivity factors calculated from (5) and (6), respectively. A bivariate regression analysis between Z, in dBZ units, and log R was used to determine the coefficients log A and b for each rain event during this period. An “event” was defined as a rain period with at least sixty 1-min-average spectra in which there were no contiguous 10-min periods without precipitation.
The software program (FITEXY) given in Numerical Recipes (Press et al. 1992) was used. Unlike standard least squares fitting routines where the independent variable x is assumed to be without error, in bivariate fitting both the x and y variables are assumed to have errors. Bivariate regression requires an estimate of the errors in both x and y, and these may be a function of x and y, respectively. The errors are assumed to have a Gaussian distribution. This is approximately true for precipitation rates greater than 0.1 mm h−1 and reflectivity factors greater than 10 dBZ; therefore, data outside this range was not included in the analysis. One advantage of bivariate regression is that the coefficients in (10) do not depend on the choice of independent variable.
Because the estimates of Z and R come from the same estimated DSD, their errors are correlated. It was necessary to modify the Numerical Recipes code to account for this correlation following Acton (1966). The standard deviations σR and σZ in units of log R and dBZ, respectively, are determined from the analysis given in Figs. 10 and 11, respectively. The correlation of these errors was also determined as a function of Ẑ in units of dBZ and is given in Fig. 14.
Figure 15 shows the log A versus b coefficients for each rain event. Also shown are values for log A and b given by MP48, and the convective and tropical Next Generation Weather Radar (NEXRAD) values (Fulton et al. 1998). The uncertainties in the coefficient estimates are represented by the error bars. The coefficients cover a wide range and the error bars indicate that they are significantly different from the MP and NEXRAD values, and different from each other with the exception of two events during this period.
The uncertainty in estimates of rainfall parameters by POSS has been characterized using a Monte Carlo simulation of falling raindrops at terminal velocity with vertical trajectories through the sensor’s measurement volume. An “inverse problem” analysis is used to calculate the statistics of the pdf of the underlying population for a given POSS estimate of a rainfall parameter. Results for several parameters are presented graphically. Previous POSS uncertainty analysis was derived from the forward model only, and this model only considered the sampling effects on the total scattered power by each diameter and not the effect on the Doppler spectral distribution of this power.
These results will allow researchers who use POSS estimates of rainfall parameters to evaluate their results knowing the uncertainties in the measurement due to sampling. This is necessary to determine if the differences in derived parameters are statistically significant.
The sampling error has been reduced by increasing the real-time processing speed so that more raindrops are sampled per averaging period. The previous Intel® 8086 processor calculates 380 Doppler velocity spectra per minute. The newest version of the sensor uses an Intel 386 processor that can calculate about 960 spectra per minute.
The motivation for this paper resulted from discussions with I. Zawadzki of McGill University, and N. Donaldson and P. Joe, both from Environment Canada.
Comparison of “Forward” versus “Inverse” Conditional PDFs
Not withstanding the specific sampling issues related to the antenna geometry of the POSS, there will always be sampling effects owing to the statistical nature of precipitation. If the precipitation process is stationary, and the samples are statistically independent, the number of drops in a diameter class, in a fixed measurement volume, can be assumed to be Poisson distributed.
If the mean number of drops in the measurement volume is N, then assuming there are no other instrumental sampling biases, the measured estimate N̂ will have a Poisson probability distribution function given by
Figure A1 presents interpolated contours of the conditional Poisson probability distribution function p(N̂|N) for a range of N from 0 to 30. Note that although the contours are shown as continuous, N̂ has a range of discrete integer values ≥0, and N has a range of continuous real values ≥ 0. An example of the forward conditional probability distribution function p(N̂|N) with N = 10 is given on the panel below the x axis. Choi (1994) has shown that the median of a Poisson distribution is ≥ (N − ln 2) and < (N + ⅓).
If it is assumed that p(N) is a uniform pdf in (4), then the inverse conditional pdf p(N|N̂) can be written as a gamma distribution:
where Γ(N̂ + 1) = N̂! is the standard gamma function. The mean of this pdf is N̂ + 1 and the median is approximately N̂ + ⅔ (Choi 1994). Therefore, N̂ underestimates the mean and median by 1 and about ⅔ counts, respectively. If the number of drops is large, the fractional difference in estimates of the mean or median resulting from use of the forward rather than inverse conditional pdf is insignificant.
Unlike disdrometers that are discrete counters, POSS estimates of the number of drops in its measurement volume are continuous. Poisson sampling variability is less important for POSS because of its larger measurement volume and it is ignored in the following example. However, because the Doppler signal’s amplitude and frequency depend on the raindrop’s location in the volume, there is a different sampling problem. For a single raindrop in the measurement volume, the pdf p(N̂|N = 1) is positively skewed. The lognormal distribution (Wilks 1995) has this property and was chosen to represent the pdf when N = 1:
where lnN is the mean of the natural logarithm of the population, and var(N) and var(lnN) are the variance and the variance of the natural logarithm of the population, respectively. In general, the relationships between N, lnN, var(N), and var(lnN) are defined by
In this example for a single drop, N = 1 is substituted in (A4), then
Substitution of (A6) and (A7) in (A3) defines p(N̂|1). The pdf p(N̂|N) for N > 1 was found by numerically autoconvolving the lognormal distribution p(N̂|1) N − 1 times. The resulting pdf is not lognormal. It approaches a Gaussian distribution as N increases in accordance with the central limit theorem.
For example, Fig. A2 shows interpolated contours of the forward conditional pdf p(N̂|N) when var(N) = 10 in (A7). The forward conditional pdf for N = 10 is given on the panel below the x axis. Again assuming a uniform pdf for p(N), the inverse conditional pdf for N̂ = 10 is shown on the panel beside the right-hand axis. The shape of these density functions is quite different. The mean of p(N̂|N = 10) is 10, and the median is about 7.5. The mean and median of the pdf p(N|N̂ = 10) are both about 15. Therefore, if the forward POSS pdf for a single drop were represented by the example given here, then an estimate of 10 drops would underestimate both the population mean and median by 5 drops.
Corresponding author address: B. E. Sheppard, Cloud Physics and Severe Weather Research Section, Environment Canada, 4905 Dufferin St., Toronto, ON M3H5T4, Canada. Email: email@example.com