Estimating raindrop size distributions using microwave link measurements: potential and limitations
- 1Hydrology and Quantitative Water Management Group, Wageningen University, Wageningen, the Netherlands
- 2R&D Observations and Data Technology, Royal Netherlands Meteorological Institute (KNMI), De Bilt, the Netherlands
Correspondence: Remko Uijlenhoet (firstname.lastname@example.org)
We present a novel method of using two or three collocated microwave link instruments to estimate the three parameters of a gamma raindrop size distribution (DSD) model. This allows path-average DSD measurements over a path length of several kilometers as opposed to the point measurements of conventional disdrometers. Our method is validated in a round-trip manner using simulated DSD fields as well as five laser disdrometers installed along a path. Different potential link combinations of frequency and polarization are investigated. We also present preliminary results from the application of this method to an experimental setup of collocated microwave links measuring at 26 and 38 GHz along a 2.2 km path. Simulations show that a DSD retrieval on the basis of microwave links can be accurate under idealized conditions. We found that a triple-link retrieval provides little added benefit over a dual-link retrieval in terms of accuracy or precision. In practice, the accuracy and success rate of any retrieval is highly dependent on the stability of the base power level as well as the precision of the instruments and in particular the quantization applied to the recorded power level.
The use of microwave links to measure rainfall intensity has received significant attention in the last decade. The main driver for this has been the insight that the backhaul links of mobile communication networks are suitable for such measurements and are available in greater numbers than dedicated rainfall measurement stations in many countries (Messer et al., 2006; Leijnse et al., 2007b; Gosset et al., 2016). The usage of microwave uplinks to geostationary satellites has been investigated as well (Giannetti et al., 2017). In addition, efforts have been taken to expand the range of atmospheric phenomena that can be measured with microwave links, including fog (David et al., 2013), solid precipitation (Cherkassky et al., 2014) and evapotranspiration (Leijnse et al., 2007a). Another enticing possibility is the use of multiple link instruments along the same path to measure not just the bulk rainfall intensity but also the path-average raindrop size distribution (DSD).
DSD estimates can be used to derive all other bulk rainfall variables and are therefore valuable for a variety of purposes (see, e.g., Uijlenhoet and Stricker, 1999, for an overview of relevant statistical moments). Examples include precipitation microphysics (e.g., Uijlenhoet et al., 2003), soil erosion by rain (e.g., Angulo-Martínez and Barros, 2015; Salles et al., 1999; Salles and Poesen, 1999, 2000) and radar validation (e.g., Hazenberg et al., 2014). The use of microwave links for such purposes promises measurements that are more spatially representative of radar or satellite pixels than what is offered by the usually sparse networks of impact, laser and video disdrometers that are currently the most common way DSDs are measured.
In order to estimate the DSD from a limited number of statistical moments, a parameterization has to be used. The gamma distribution with three parameters provides a good estimate for a wide variety of rainfall types (Ulbrich, 1983). Rincon and Lang (2002) were the first to attempt a gamma DSD parameter retrieval using two microwave links, with promising results. However, their methods require an a priori parameter estimate in addition to the parameters derived from the two measurements. To the best of our knowledge, no further research has been published regarding DSD estimations using microwave links. However, the adjacent field of polarimetric radar measurement has seen much development (see, e.g., Fabry, 2015), and many different methods to estimate DSDs from polarimetric radar have now been developed and tested.
A handful of different techniques can be identified that have been developed to retrieve DSD parameters from a limited number of polarimetric radar moments: the constrained-gamma method developed by Zhang et al. (2001), the β method proposed by Gorgucci et al. (2002) and double-moment DSD normalization (Raupach and Berne, 2017). Several machine learning approaches have also been used: neural networks (Vulpiani et al., 2006), Bayesian regression (Cao et al., 2010) and tree-based genetic programming (Islam et al., 2012).
In this paper we explore the potential of a numerical retrieval from microwave attenuation and/or differential propagation phase shift. We consider two different techniques: the first method uses three measured microwave link variables to derive the three parameters of the gamma distribution. Here, the gamma parameters are weakly constrained (i.e., only a limited range of parameter values is allowed). The other method uses two measured microwave link variables to derive two parameters of the gamma distribution. Here, the third parameter is completely constrained by the other two, similar to the method used by Zhang et al. (2001) for radar moments. We apply these techniques first using a simulated dataset based on radar and disdrometer measurements from the Ardèche region of France. Next, we apply the method to microwave link variables derived from a set of path-averaged measurements from five laser disdrometers in Wageningen, the Netherlands. Finally, we test the methods on measurements from microwave link instruments along the same path as the disdrometers in Wageningen (van Leth et al., 2018a) and compare the resulting DSDs with the DSDs measured from the disdrometers.
The paper is organized as follows: in Sect. 2 we present in more detail the datasets used in this paper. In Sect. 3 we present the theory and the methods used to retrieve DSDs from microwave link variables. We also describe the validation methods employed. In Sect. 4 we discuss the results retrieved from the Ardèche dataset. In Sect. 5 we discuss the results of several tests using simulated attenuations from the Wageningen disdrometer dataset. We also consider the efficacy of different frequency combinations and the robustness of the retrieval to measurement uncertainty. In Sect. 6 we apply the developed methods to actual link measurements that were obtained in Wageningen. In Sect. 7 we present our thoughts on the feasibility of the techniques in practice and the choices made in this paper. Finally, in Sect. 8 we come to conclusions and give recommendations for further study.
2.1 Ardèche DSD reanalyses
Our first test case consists of a two-dimensional interpolated DSD field based on polarimetric radar data measured in the Ardèche, France, as part of the HyMeX campaign (Raupach and Berne, 2016, 2017). This field was generated using an advection-based temporal interpolation technique (De Vos et al., 2018), where disdrometer data were used to train the technique. This technique does not produce μ and Λ parameters, but instead produces a full DSD histogram per pixel per time step. The field has spatial dimensions of 20×20 km and a resolution of 100 m. It covers two distinct events on 27 November 2012 and 27 October 2013. Their durations are on the order of several hours, and the time step is 30 s. Both events can be classified as orographic or convective events. The second event is more spatially heterogeneous than the first, with a decorrelation distance of 2.8 km vs. 11 km at 30 s accumulation intervals. The DSD is divided into 20 diameter bins with unequal bin width based on the detection bins of the OTT Parsivel2 laser disdrometer measurements, upon which the reanalysis is partially based. The smallest drop diameter class is 0.3 mm, while the largest diameter is 6.5 mm. A more detailed description of the dataset and the technique used to generate it can be found in De Vos et al. (2018).
We take a transect through this field over the entire length of the field (Fig. 1a) and average the DSD over this transect to approximate the footprint of a microwave link. We calculate the attenuations and differential phase shifts for a number of frequencies: 15, 26, 32 and 38 GHz. This includes the frequencies employed in the Wageningen experiment (see below) and covers the range of most commonly employed carrier frequencies in mobile phone networks.
In order to calculate the microwave link variables from the DSD data, we first need to interpolate the data in the diameter dimension from the irregular bins to a regular diameter grid with a resolution of 0.1 mm. The main reason for this is to accurately reproduce the shape of the scattering cross sections as a function of diameter, as simply assuming the scattering cross section to be constant for the entire diameter interval would introduce too much error. We have used a simple linear interpolation method for this purpose.
2.2 Wageningen link experiment
Our second test case is a microwave link experimental setup in Wageningen, the Netherlands (van Leth et al., 2018a). The setup consists of three collocated microwave links arranged between two buildings 2.2 km apart and covering mostly built-up terrain. The setup contains one dual-polarization 38 GHz link, which also measures the phase difference between the two polarizations, one additional single-polarization 38 GHz link (not used for this paper) and a single-polarization 26 GHz link (all with sampling frequencies of 20 Hz; down-sampled – i.e., averaged – to 30 s intervals). In addition, five OTT Parsivel laser disdrometers (providing DSDs at 30 s intervals) are positioned at four locations beneath the link path, including the sites of the transmitting and receiving antennas, as shown in Fig. 1b. The data used in the following analyses are all taken from the period between 1 April 2015 and 1 January 2016. We also specifically focus on one event on 27 July 2015 for illustrative purposes.
In order to use the disdrometer data to represent the DSD of the link path, we take a weighted spatial average over the disdrometer data. As with the Ardèche data, we interpolate the DSD data in the diameter dimension from the irregular bins to regular intervals before calculating the microwave link variables. In order to improve the robustness of the results, we only consider measurement intervals where each one of the five disdrometers counted at least 50 drops (see Uijlenhoet et al., 2006, for the effect of drop sample size on the robustness of the DSD estimate). We also apply the correction method of Raupach and Berne (2015) to the disdrometer DSDs.
3.1 Basic procedure
To determine the underlying DSD from a limited number of statistical moments, we need an approximation with a limited number of parameters. One of the most widely used approximations for rain DSD is the gamma distribution suggested by Ulbrich (1983),
where D is the drop diameter in millimeters, and N0, μ and Λ are the parameters determining the shape and drop concentration. The parameter μ is a dimensionless shape parameter, and Λ is a slope parameter with a unit of per millimeter (mm−1). Note that the dimension of N0 is dependent on the value of μ. Therefore it is convenient to also use a derived parameter,
where Γ is the gamma function. NT has the unit per cubic meter (m−3) and is equal to the total drop concentration (assuming integration limits of 0 to ∞), resulting in
The specific attenuation of a link signal in decibels per kilometer (dB km−1) at a given frequency can be described in terms of the DSD,
where λ is the wavelength of the incoming and outgoing waves in millimeters (mm), ℑ is an operator indicating the imaginary part of its argument and is a unit conversion factor. SHH and SVV (dimensionless) are the diagonal components of the forward-scattering amplitude matrix S defined by the relationship
where E0 is the electric field strength of the incoming electromagnetic (EM) wave and E is the electric field strength of the outgoing EM wave. The specific phase shift (in rad km−1) between the horizontal and vertical components of an outgoing EM wave of a given frequency due to forward scattering can also be described as an integral over the DSD,
where ℜ is an operator indicating the real parts instead of the imaginary parts of the diagonal components of the scattering amplitude matrix. is another unit conversion factor.
The corresponding dimensionless scattering efficiencies are defined as
It is the subtle differences in the shapes of the scattering efficiencies as functions of diameter for different frequencies and polarizations that contain the information necessary to retrieve the DSD. We calculate the scattering amplitude matrix using the T-matrix method (Waterman, 1965) following the approach of Mishchenko et al. (1996), Mishchenko and Travis (1998) and Mishchenko (2000). We assume the drops to be oblate spheroids with axis ratios dependent on the drop diameters following the relationship of Thurai et al. (2007) and averaging over canting angles following a Gaussian distribution with a mean of 0∘ (vertical) and a standard deviation of 2∘. Resulting scattering efficiencies at some relevant frequencies are shown in Fig. 2 as a function of volume-equivalent diameter.
The variables used as input in the retrieval could be attenuation of horizontally or vertically polarized radiation or phase differences between horizontally and vertically polarized radiation at one or several frequencies. In order to be able to use attenuations and phase differences interchangeably in the retrieval algorithm, we rearrange Eqs. (4), (5) and (7) to have the same form:
where , , and the scattering cross sections (σX) are the scattering efficiencies (Eqs. 8–10) multiplied by the drop cross section (): , and . We will from here on refer to an arbitrary input microwave link variable as , where , such that
3.2 Three-parameter method
Note that these ratios do not depend on the N0 or NT parameter. If we now replace the integrals by discrete summations, we can use an iterative nonlinear root-finding technique to find values for μ and Λ from two ratios of microwave moments. Knowing μ and Λ, we can directly solve the discretized equation of one of the microwave link variables for NT.
Figure 3 shows several possible ratios of microwave link variables as a function of μ and Λ. These observables are calculated using the forward model of Eq. (15). Because it is possible to detect the attenuation of the horizontally polarized signal and the vertically polarized signal and the phase difference at a given frequency using a single set of antennae, we prefer this combination of variables over combinations of attenuations at multiple frequencies. We can see in Fig. 3 that the ratios and provide complementary information and are therefore in principle suitable.
We use the Powell hybrid method (Powell, 1970) to solve the system of equations of Eq. (15). We also tested several other gradient-based root-finding methods such as Levenberg–Marquardt and found that this makes little difference in the stability of the retrieval. The convergence of the retrieval is highly dependent on the initial guess. This is problematic, as we want to automatically retrieve a large number of DSDs without manual input. In many retrieval attempts the gradient-based root-finding algorithm diverges towards infinity. In order to prevent this, we restrict the root-finding algorithm to a limited range of parameter values. If the estimates reach the edge of the range, we reset the parameters to a new initial guess which is taken from within the allowed range but offset from the first guess by or Δμ=0.2. We systematically traverse the μ–Λ space in this way starting from (−1, 0) until convergence is reached.
Even when the algorithm converges, the solution of the system of equations is not necessarily unique. We need an extra set of constraints to make sure we retrieve the parameters that are the most plausible. In order to do so, we use the parameters retrieved by the analytical method of moments of Tokay and Short (1996) (TS96) as a rough indication of plausible combinations of μ and Λ values. We have calculated the parameters using this method for all individual time steps in the 9-month record of the spatially weighted average of the disdrometers placed in Wageningen. Figure 4b shows all these individual retrievals in the (Λ,μ) space. We then apply a kernel density estimator to this dataset (shown in Fig. 4a) and calculate the total fraction of data points contained within the contour lines. We then choose the contour corresponding to the 0.95 quantile as a mask. This mask is also shown in Fig. 4b.
The mask is then applied to all attenuation/phase-based numerical retrievals to define the range of allowed parameter values. If the estimate is outside this contour, we reset the root-finding algorithm with a new initial guess that is within the contour but slightly perturbed from the previous initial guess as described above. This is continued until we find convergence within the contour or we reach the maximum number of iterations without a solution.
3.3 Two-parameter method
When only two microwave link variables are available, the system of equations of Eq. (15) is reduced to just one equation:
In order to still solve for the two parameters, an additional equation is required for the relationship between μ and Λ. We obtain this relationship empirically. We obtain gamma DSD parameter values via the analytical method of moments from the Wageningen disdrometer dataset. We then fit a second-order polynomial function (as shown in Fig. 4b) to the values of μ and Λ with a linear least-squares method. All 9 months of data were used to fit this function. The resulting relationship is
which is similar in magnitude but still somewhat different to the relationship found by Zhang et al. (2003). We can substitute this equation for Λ in Eq. (16) and then solve for μ using Brent's root-finding method (Brent, 1973). We prefer this method because it is not based on gradients and therefore guaranteed to find a solution if it exists. It is, however, not applicable to multivariate problems.
Figure 5 shows the ratio of Eq. (16) as a function of μ for two combinations of microwave attenuations that can be measured with the Wageningen link setup: two polarizations at 38 GHz and two frequencies (26 and 38 GHz) at horizontal polarization. From this we can see that the dual-polarization configuration is more suitable for retrieving the DSD; the dual-frequency configuration has non-unique solutions for high underlying values of μ, whereas the dual-polarization model is monotonously decreasing over the entire range of valid μ values. On the other hand, the dual-frequency ratio of attenuations is much more sensitive to changes in μ for μ between −2 and 8, potentially yielding more accurate estimates of this parameter. For the dual-polarization ratio (at 38 GHz) it can also be seen that ratios lower than 1 and higher than 1.25 are not valid. Such ratios would yield no solution.
3.4 Validation methods
We test the capability of the methods to accurately retrieve DSDs and their associated statistical moments with two different datasets of drop size distributions: one simulated and one measured. We use Eqs. (11), (12) and (13) to calculate the microwave link variables from the known DSDs. We then use those variables as input to the retrieval algorithms described in the previous sections to retrieve the parameters of the gamma distribution approximating the DSD. From those parameters the complete DSD and its statistical moments are reconstructed. We then compare these to the original DSD and its moments to assess the systematic bias and random error in the retrieval. To be able to distinguish between cases where the gamma distribution is simply not a good fit for the measured DSD and cases where the retrieval itself is the cause for inaccuracies, we also calculate the gamma parameters using the TS96 method. If the gamma distribution is not a good fit for the measured DSD, then the results from the TS96 method would deviate appreciably from the measured DSD. If this is not the case but the retrieved DSD does deviate considerably from the measured DSD, then the retrieval itself is the main cause of the inaccuracy in the DSD. The TS96 method can be used to directly evaluate the parameters of the distribution themselves as well if we assume that the parameters yielded by the TS96 approach are the “true” parameters.
In order to assess the performance of the retrieval methods, we will use a number of statistical measures throughout the Results section. As a measure of the accuracy of the retrieval we use the median of the residuals (MOR). As a measure of the precision of the retrieval we use the median absolute deviation (MAD) of the residuals with respect to the median of the residuals. We chose MAD and MOR over the use of standard deviation and means, because there are a relatively small number of extreme deviations which would otherwise have too much influence and thus would not give much information about the typical precision. The statistical metrics employed here are less influenced by non-normality and outliers. Because the number and severity of the extreme deviations are also an important part of the performance assessment of the retrieval, we also compute the 95th-percentile absolute deviation with respect to the median of the residuals (95AD). Together with the MAD this gives a more complete picture of the distribution of the errors, while still being insensitive to the true outliers. All metrics are normalized with respect to the median of the original quantities; hence they are dimensionless. Furthermore, we also compute the fraction of non-convergent retrievals compared to the total number of retrievals (taking into account the filtering described in Sect. 2.2). This “failure ratio” (also dimensionless) is necessary for a complete picture of the robustness of the method since the other metrics naturally exclude these intervals.
4.1 Single retrieval
Figure 6a shows a typical three-parameter retrieval result from the Ardèche dataset using horizontal attenuation, vertical attenuation and phase difference at a single frequency (38 or 26 GHz in this case). The normalized difference between the retrieved DSD and the original simulation procedure,
where Nr is the retrieved DSD, Nm is the original DSD and NT is the integral over all diameters of the original DSD – is illustrated in Fig. 6c. ΔN* is within 10−4 for particles larger than 1 mm. However, at the smallest sizes the results tend to diverge up to . In practice we are often more interested in quantities that scale with the statistical moments of the DSD rather than in the DSD itself. Important quantities are, for example, liquid water content, W (third-order moment); rain intensity, R (close to fourth order); and radar reflectivity, Z (sixth order). For higher-order moments, the contribution of the smallest drop sizes decreases. To illustrate this, we also show the specific rain intensity as a function of drop diameter,
where Ro is the total rain intensity based on the original DSD, rr(D) is the retrieved diameter-specific rain intensity and ro(D) is the original diameter-specific rain intensity. From Fig. 6d it can be seen that for all drop sizes. The difference in the total drop concentration is in the first case, while the difference in the total rain intensity is . The assessment of the accuracy of such a retrieval must therefore take into account its most likely application. We will focus our attention in this section on the rainfall intensity and to a lesser extent on the individual gamma DSD parameters. In further sections, when more detailed comparison is required, we will also look at a range of integer moments.
4.2 Complete events
Figure 7 shows the results of the three-parameter retrieval for the first Ardèche event, which took place from 26 November 2012 at 04:54 UTC to 27 November at 04:36 UTC. The total duration is almost 24 h, and the precipitation intensity averages at 1.77 mm h−1 with a maximum of 10.29 mm h−1. The gamma DSD parameters of the retrieval are mostly very close to those of the TS96 procedure, with a few rather large exceptions, particularly noticeable in the NT parameter. These outliers do not seem to correspond with any particularly high or low precipitation intensity, but they do correspond with high drop concentrations. The temporal evolution of μ is very close to the temporal evolution of Λ, with a correlation coefficient of 0.86 (not shown). Similar results can be observed for the second event (see Fig. 8), which took place on 27 October 2013 from 03:22 UTC to 08:48 UTC, with a total duration of 3.5 h, an average intensity of 1.58 mm h−1 and a maximum intensity of 27.05 mm h−1. The correlation coefficient between μ and Λ is ρ=0.87. Outliers for this event also occur at lower drop concentrations.
Looking at the rainfall intensity for both events (Figs. 7d and 8d), we see that the retrieval corresponds closely not only to the TS96 approach but also to the rainfall intensity derived from the original DSD. There are only a handful of outliers here, the exact timing of which is dependent on the carrier frequencies for which the retrieval is attempted. The MOR, MAD and 95AD of the rain intensity are given in Table 1 for retrievals based on several different carrier frequencies. Another way of assessing the performance of the retrievals is to consider the temporal mean of the DSD. The retrieval (not shown) gives an overestimation for small diameters (with the exception of very low carrier frequencies), while becoming more accurate for larger diameters. For diameters larger than 1 mm the residuals are lower than 0.1 .
We apply both the two-parameter and three-parameter retrieval algorithms to the Wageningen disdrometer dataset. The results are shown in Fig. 9. We selected a single rain event on 27 July 2015 starting at 09:40 UTC and ending at 12:00 UTC. The first notable difference between the Ardèche and the Wageningen datasets is that in the latter the values for μ and Λ are generally much higher. For both the TS96 approach and the numerical retrieval, values higher than 20 occur regularly for both the μ and the Λ parameter, which is not consistent with what is typically found in the literature (e.g., Raupach and Berne, 2016, for the Ardèche; Zhang et al., 2003, for Florida; and Atlas and Ulbrich, 2006, for Kapingamarangi Atoll (western tropical Pacific)). We can also see that the two-parameter retrieval yields an overestimation of the total drop concentration, while the μ parameter gets underestimated. Because these biases compensate for one another, this results in a rainfall intensity that is not significantly biased from the original. The three-parameter retrieval method yields a significant overestimation of μ and Λ, including a large number of outliers. The three-parameter retrieval produces a less accurate result in terms of rainfall intensity as well. Furthermore, for both retrieval methods, when averaged over an entire event, the retrieval (not shown) overestimates the number of drops with a diameter smaller than 1 mm and underestimates the number of drops with a diameter between 1 and 7.5 mm when compared to the measured DSD (except for very low carrier frequencies). The correlation between the TS96-derived μ and Λ parameters is very high at ρ=0.96 for the two-parameter retrieval. The MOR, MAD and 95AD are given in Table 2. In general, the precision of the retrieval is an order of magnitude lower (higher MAD and 95AD), while the accuracy is actually higher (lower MOR) when compared with the Ardèche dataset.
5.1 Differences between two- and three-parameter retrievals
We summarize the differences in accuracy and precision between the two-parameter and the three-parameter method in Table 3. The analyses are all performed with respect to a 38 GHz dual-polarization retrieval (the third microwave link variable is the differential propagation phase). In addition, the number of failed retrievals (no solution at all) was 1.7 % when using only two moments, whereas there where no failed retrievals using three parameters within the filtered dataset. This indicates that there is at least some advantage to using three moments.
However, the differences in accuracy and precision of the retrieval between the three-parameter and two-parameter retrieval as measured by a range of integer moments is small, and in many cases the two-parameter retrieval proved to be more reliable. Especially the number of sub-millimeter raindrops is severely overestimated by using the three-parameter method, as shown in Fig. 10a. Figure 10b also demonstrates that there is little difference in precision between the two-parameter and three-parameter method for any diameter class and even a slight advantage for the two-parameter method. Therefore, the addition of a third microwave link variable does not improve the retrieval (in many cases it actually harms the retrieval) and is unnecessary. Aside from needing one fewer measured moment, the two-parameter retrieval is also orders of magnitude faster. This is because the numerical part is univariate and therefore the dimensionality of the problem is reduced and also more efficient root-finding methods other than gradient-based methods can be used (such as Brent's method). We do not need the workaround for local minima either, which is computationally very inefficient.
Because these results show that a three-parameter retrieval provides little added value above a two-parameter retrieval and because the two-parameter retrievals are far less computationally intensive than three-parameter retrievals, we will restrict ourselves to two-parameter retrievals in the remainder of this paper.
5.2 Dependence on link frequency
In order to determine the effect of the carrier frequencies of the links on the accuracy and precision of the retrieval, we perform two-parameter DSD retrievals at many different frequencies and calculate MOR, MAD and 95AD for the third-order moment of the retrieval compared with the third-order moment directly calculated from the measured DSD. We consider both dual-polarization retrievals with frequencies ranging from 10 to 45 GHz (with steps of 1 GHz) and dual-frequency retrievals using every combination between 10 and 45 GHz. This range contains the bulk of microwave link frequencies found in typical communication networks. The results are shown for dual-polarization in Fig. 11 and for dual-frequency in Fig. 12. We can see in Fig. 11 that the accuracy and precision is high for all frequencies. However, the accuracy is pessimal for frequencies between 22 and 34 GHz. Similarly, MAD is largest around 25 GHz, and 95AD is largest around 17 GHz. Therefore, for an optimal retrieval those intermediate frequencies should be avoided. Figure 12 shows that the accuracy and precision of dual-frequency retrievals is highest when both frequencies are high. The difference between the two frequencies does not seem to matter much; when the two frequencies are far apart, the precision and frequency are actually slightly lower. Predictably, there are no solutions found when the frequencies are exactly the same. It should be noted that for this simulation the effect of noise is not taken into account. It is expected that this would influence the retrieval the most when the frequencies are close together.
5.3 Sensitivity to attenuation bias
Because our retrieval algorithm uses ratios of attenuations as input, it is important that a reliable baseline power level is established from which to calculate the attenuations. To assess the sensitivity of the retrieval technique to inaccuracies in the baseline (dry) power level, we perform the two-moment retrievals based on the simulated attenuations from the disdrometer measurements in Wageningen but with an (equal) offset added to all input attenuations. Figure 13 shows the resulting DSDs averaged over 9 months for attenuation offsets between 0 and 5 dB and a path length of 2.2 km. For this analysis we chose two combinations of links that we will also later use for the actual link measurements: a dual-polarization retrieval at 38 GHz and a dual-frequency retrieval at 26 and 38 GHz. We can see in Fig. 13a that the addition of an offset to the attenuation leads to an overestimation below 2 mm and an underestimation above 2 mm for the dual-polarization retrieval. The introduction of an offset has no effect at 2 mm, and the effects are largest below 1 mm. The effects for a dual-frequency retrieval are quite different, as can be seen in Fig. 13b. There is an overestimation at all diameters. The overestimation is smallest around 1 mm and increasing towards higher and lower diameters. Overall, the mean bias is larger than for dual-polarization retrievals, and the shape of the DSD is especially sensitive to small offsets in the base power level.
5.4 Sensitivity to power quantization error
Aside from systematic measurement bias as discussed in Sect. 5.3, there can also be measurement limitations that affect the precision of the attenuation measurements. Because the retrieval method relies on ratios of attenuations, we expect that the retrieval is highly sensitive to such limitations as well. In practice, when processing link attenuation data from operational telecommunication networks, the power quantization error of the analog-to-digital conversion completely overshadows any instrumental error in the analog detector. Therefore, we will focus our analysis here exclusively on such quantization errors. Consequently, we do not have to make any assumptions about the instrumental precision of any particular transceiver model.
We have applied several different magnitudes of rounding to the attenuations calculated from the disdrometer data and performed the retrieval on the rounded attenuations for the complete 9-month dataset. The results are shown in Table 4. Common attenuation levels in operational networks are 0.1, 0.5 and in some cases 1 dB. It is clear, then, that the effect of quantization on the performance of the retrieval method is significant. Especially the number of non-converging retrievals (from 1.7 % to 66 %) limits the prospective of successful application to current networks. The MAD calculated from the remaining successful retrievals is an order of magnitude higher than without quantization applied (see Table 2). The systematic bias is more than 2 orders of magnitude higher, but still limited (MOR=4.3 % at 0.1 dB). To achieve non-convergence ratios of less than 10 %, a quantization of 0.001 dB or less is required, which makes this unachievable with current-generation operational networks. It should also be noted that taking into account the quantization error in the analysis favors the dual-frequency method over the dual-polarization method. This can be attributed to the steeper slope of the attenuation-ratio–μ relationship within the band of common DSD shapes as shown in Fig. 5.
Using the double-moment retrieval method, we estimated the DSDs from actual link measurements of the Wageningen setup. The baseline power level of the links showed considerable fluctuations over the course of the measurement period. Therefore, it was not feasible to perform retrievals for the entire 9-month dataset. We selected the event of 27 July 2015 (see Fig. 14), because the power levels of the links in the period surrounding this event showed relatively little fluctuations. We determined a suitable constant baseline power level calibrated for this event. The measured attenuation after subtracting a baseline for this event is shown in Fig. 14 and compared with the path-average attenuation derived from the disdrometer measurements; the retrieval results are given in Fig. 15. It should be noted that for this experiment we used the analog detector voltage directly fed into a data logger with a 13 bit analog-to-digital converter (ADC) and no further quantization applied, so – unlike for operational networks – we expect quantization effects to be small compared to instrumental errors.
The resulting DSD is very similar in shape to that obtained in the simulations, with overestimations especially at smaller diameters, but with the general shape of the DSD preserved (not shown). Closer inspection reveals that the bias and scatter compared to the original DSD are actually up to 2 orders of magnitude higher than in the simulations, as can be seen in Table 5 when compared to Table 3. The scatter as indicated by the MAD is about 1 order of magnitude higher than the simulations across all moments. However, when taking the 95AD as a measure of scatter, the order of magnitude is the same. There are more intervals with no solution at all when compared to the simulations. These correspond with ratios of observables that are outside the range of the forward model (see Fig. 5). The ratios for the dual-polarization retrieval are illustrated in Fig. 14b. These time intervals with extremely low or high ratios between attenuations mostly (but not always) occur when the rain intensity is low and thus other sources of signal variability are more dominant. Overall, the dual-polarization and dual-frequency retrievals have a similar performance in this case. However, using two frequencies instead of two polarizations leads to a higher accuracy for low-order moments (up to fourth order), but lower accuracy for higher-order moments. Similar to the simulated retrievals based on the disdrometer dataset, there is an overestimation of the NT parameter and underestimation of the μ parameter. There is also an overestimation of the Λ parameter. Especially the dual-polarization retrieval yields some outliers which overestimate μ but underestimate Λ. No plausible solutions for the retrieval were obtained when using the phase difference instead of one of the attenuations or with a triple variable retrieval. Very few intervals showed convergence at all with this configuration.
7.1 Feasibility in practice
Constraints on the feasibility of the proposed methods in practice fall into three broad categories: availability of multiple link signals on the same path, quality of the available signals and real-time processing speed.
The use of a three-moment retrieval means that three moments on the same path need to be available. This is rare in commercial networks; therefore this method is most readily applicable to dedicated research networks. There are several different combinations of moments to choose from. However, in our approach we focused on the combination of a horizontal attenuation, vertical attenuation and phase difference at the same carrier frequency. This allows the use of a single set of antennae for all three moments, allowing for the use of a more compact and less expensive device (such as the device that was used in our test setup).
The second concern is with regard to the quality and reliability of the signal. In order to apply the method in practice, it is essential that a baseline (no rain) signal is accurately determined. Because the method relies on the ratios of attenuations, small deviations in the baseline determination can result in large deviations in the retrieval (and even non-convergence). No such problem exists in principle with regard to the phase difference; it is independent of any power baseline. However, phase difference on its own is not sufficient for the retrieval. To have a chance at a successful retrieval, the baseline needs to be as invariant as possible. Where it is not, the variability should be accurately modeled and predicted from auxiliary measurements. In our own preliminary attempts we found our instruments lacking in stability. In particular the clinging of drops to the antenna cover (as described in van Leth et al., 2018a) seems an intractable problem. However, as was also described in that paper, we found that a former commercial microwave link had a much stabler baseline and furthermore the effect of wet antennas was much more manageable for that particular device. This provides a hopeful perspective for the application of this method to commercial networks, in particular if data could be logged with high precision (Chwala et al., 2012). However, when using data from such commercial networks, the quantization of the signal plays a much bigger role. We have seen that the typical quantization levels applied to the data collected by the network operators are insufficient to allow reliable DSD retrievals. The only way to apply such retrievals to currently operational unmodified link networks consistently is to install dedicated data loggers at selected link locations to read out the analog signal directly. This may not be feasible. Therefore, a higher-resolution default storage of link-received power level by future models of transceiver units would be needed before wide-scale application could take place.
The third constraint is only relevant when real-time processing is required. The three-moment method is relatively wasteful with computing cycles because of the repeated re-initialization of the root-finding process. We might expect this to become a bottleneck. However, in its current implementation the processing of 9 months of data from one link requires roughly 3 h of wall time on a high-end desktop workstation. This is while utilizing 10 CPU cores simultaneously. Such a setup is therefore expected to be able to process a network of more than 2000 links in real time. Nevertheless, reducing this computational load would make real-time retrievals more feasible on low-end (embedded) hardware as well. We already found that the specific programmatic implementation of the root-finding algorithm can make an order-of-magnitude difference in computation time. The use of pre-computed lookup tables may help to bring down computation time in a real-time setting. The two-moment method requires far less computation (e.g., processing the same 9 months of data with the same workstation requires only 5 min), which makes it far more suitable for real-time processing.
There are a number of caveats to our methods which could influence the interpretation of the results: firstly, we use a threshold of 50 drops per disdrometer to filter out low quality measurements before calculating the mask or μ–Λ fit. How high this number should be is debatable (e.g., Uijlenhoet et al., 2006). A threshold that is too low might allow for many erroneous measurements that are not representative of rain. On the other hand, a threshold that is too high might result in too many reasonable measurements being rejected. This might mean the results are less statistically representative and possibly biased towards high rain intensities. The number used here is therefore a compromise, and the resulting relationships do change somewhat depending on the threshold chosen. A similar consideration applies to the choice to filter on a per-instrument basis instead of on the basis of the total drop count. Filtering on a per-instrument basis makes it more likely that all instruments were measuring correctly. However, it does bias the measurements towards more homogeneous rainfall. Considering that the employed disdrometers were located relatively close together and therefore that their measurements are strongly correlated, we accepted this potential bias.
Another consideration is the use of the mask itself. The mask is determined on the basis of the measured disdrometer data. We then use this mask in, among others, the retrieval of the DSD from the disdrometer-derived simulated variables. This could potentially lead to a retrieval procedure that is biased towards these particular circumstances and therefore yields more accurate retrievals in this simulation than would be representative for a general application. Nevertheless, these data are never used as input for the root-finding procedure itself. They are only used a posteriori to assess whether the results fall into a plausible range of values. Another potential issue is the use of the predetermined μ–Λ relationship in the case of the two-parameter retrieval. In this case, the relationship determined on the basis of the disdrometer measurements is used directly as an assumption in the retrieval of these parameters from the disdrometer measurements. However, this is justified by the fact that both the μ–Λ relationship and the mask are determined from the total of all 9 months of disdrometer measurements, not from the specific event in question. It should also be noted that for the retrievals from the Ardèche dataset we have also used both the mask and the μ–Λ relationship determined from the Wageningen dataset.
The third consideration is the underlying assumption that the gamma distribution is an adequate representation of the actual DSD and that the untruncated gamma distribution is applicable throughout the diameter domain. It is on the basis of this assumption that we treat the values of μ and Λ derived from the TS96 method of moments as the correct parameter values and the DSD resulting from that as representative of the actual DSD. The gamma distribution has a nonzero value up to positive infinity. Meanwhile, there are both physical and instrumental cutoffs to the maximum drop size that can occur. This would suggest that a truncation should be included in the expression of the gamma distribution. However, the truncation of the gamma distribution at large diameters is not relevant in this case because we can see in Figs. 6a, 10a and 13 that the gamma distribution corresponding to the DSDs under consideration tends to zero (or ) before the instrumental cutoff. We can also observe from those figures that the systematic deviation between the measured (interpolated) DSD and the DSD obtained from the method of moments is small, which suggests that the gamma distribution is a valid approximation for the 30 s aggregation interval that has been used here. Instrumental cutoff at the small end of the diameter scale is relevant in this case, but the effects of this on higher-order moments is minute. Regardless, since the attenuation/phase-based retrieval in this case is limited to three parameters, it is not possible to include a cutoff there. A retrieval using even more signals might make this possible, but this would further limit the practical applicability of this method.
Using simulated link data, we have shown that a DSD retrieval on the basis of multiple microwave link variables can be successful and accurate, but only when precise high-resolution records of rain-induced attenuation are available. This was confirmed when applied on actual link data, where baseline variations prohibited accurate DSD retrievals. The use of both dual-polarization and dual-frequency retrievals is feasible. However, the use of dual-polarization is less sensitive to systematic inaccuracies in the base power level while being more sensitive to quantization errors than the use of dual-frequency links. Simulated retrievals using a variety of frequencies show that, at least between 10 and 45 GHz, the accuracy and precision of the retrievals is very high for all frequency combinations. Therefore, the frequency chosen for a dual-polarization retrieval is not very important. Nevertheless, when a choice is available, our simulations indicate that frequencies at the lowest end or the highest end of the range are preferred. Furthermore, at the lowest end of the frequency range there is less attenuation in general and therefore a smaller difference in attenuation. This could more easily be obscured by noise or quantization effects. Therefore, the higher end of the tested frequencies is optimal. For dual-frequency retrievals, bias and random error are an order of magnitude higher than for the dual-polarization retrievals when no input error is assumed. If a dual-frequency retrieval is attempted, both frequencies should be as high as possible to minimize the bias and random error.
In our field experiment we tested a dual-frequency retrieval using 26 and 38 GHz as well as a dual-polarization retrieval at 38 GHz. Both retrievals produced some reasonable results for a selected summer event where other attenuating atmospheric phenomena were not present, but there were many intervals within this event where no solution was obtained at all. The feasibility of the retrieval depends very much on a stable base power level, which was not guaranteed in our experiment. The Nokia link (former commercial link) is promising in this respect, because it was much stabler and less sensitive to, for example, temperature fluctuations and antenna wetting. A follow-up experiment using two commercial links is planned for the near future.
Using phase differences in addition to attenuations is feasible in the simulations. However, in practice these measurements are not accurate enough to yield meaningful solutions. In most instances no convergence was obtained. We have also shown that using three microwave link variables yields no improvements over a retrieval using only two variables, which is also computationally faster and more readily applicable in operational settings. At least in comparable climatologies to those treated here, a predetermined μ–Λ relation suffices to determine the gamma DSD parameters from two attenuations.
A follow-up experiment using different microwave links of similar frequencies (preferably commercially available ones) is needed to determine if the base power level of commercial links is sufficiently stable for reliable continuous observations. A tally should also be done on the number of dual-polarized links in cellular communications networks to determine if it is feasible to retrieve spatial DSD information from such networks or whether this technique is only applicable to some individual link paths, either from commercial or research networks.
Another concern is the quantization of data from commercial link networks. As the difference between the attenuation of the horizontally polarized signal and the vertically polarized signal is often a fraction of a decibel and data available from such networks are often rounded to 0.1, 0.5 or 1 dB, this limits the applicability of this method severely. Therefore, the possibility of collecting higher resolution data from such instruments should be looked into.
The underlying radar and disdrometer data employed to simulate the DSD dataset described in Sect. 2.1 are available through the HyMeX database (https://doi.org/10.14768/MISTRALS-HYMEX.721, Raupach and Berne, 2017). The link and disdrometer data from the Wageningen link experiment (described in Sect. 2.2) are publicly available through the 4TU data repository (https://doi.org/10.4121/uuid:1dd45123-c732-4390-9fe4-6e09b578d4ff, van Leth et al., 2018b).
This paper was conceived by TCvL, with some suggestions from RU. TCvL devised the algorithm, wrote the necessary software to perform the retrievals and analyzed the results, with feedback by HL and RU. TCvL wrote the first draft and, after critical revision by all co-authors, wrote the final version.
The authors declare no competing interests.
The Ardèche dataset was kindly provided by Timothy H. Raupach.
This research has been supported by the Netherlands Organization for Scientific Research (NWO-TTW; project number 11944).
This paper was edited by Saverio Mori and reviewed by four anonymous referees.
Angulo-Martínez, M. and Barros, A.: Measurement uncertainty in rainfall kinetic energy and intensity relationships for soil erosion studies: An evaluation using PARSIVEL disdrometers in the southern Appalachian mountains, Geomorphology, 228, 28–40, https://doi.org/10.1016/j.geomorph.2014.07.036, 2015. a
Atlas, D. and Ulbrich, C. W.: Drop size spectra and integral remote sensing parameters in the transition from convective to stratiform rain, Geophys. Res. Lett., 33, https://doi.org/10.1029/2006GL026824, 2006. a
Brent, R.: Algorithms for minimization without derivatives, Prentice-Hall, Englewood Cliffs, New Jersey, USA, 1973. a
Cao, Q., Zhang, G., Brandes, E. A., and Schuur, T. J.: Polarimetric radar rain estimation through retrieval of drop size distribution using a Bayesian approach, J. Appl. Meteorol. Climatol., 49, 973–990, https://doi.org/10.1175/2009JAMC2227.1, 2010. a
Cherkassky, D., Ostrometzky, J., and Messer, H.: Precipitation classification using measurements from commercial microwave links, IEEE T. Geosci. Remote, 52, 2350–2356, https://doi.org/10.1109/TGRS.2013.2259832, 2014. a
Chwala, C., Gmeiner, A., Qiu, W., Hipp, S., Nienaber, D., Siart, U., Eibert, T., Pohl, M., Seltmann, J., Fritz, J., and Kunstmann, H.: Precipitation observation using microwave backhaul links in the alpine and pre-alpine region of Southern Germany, Hydrol. Earth Syst. Sci., 16, 2647–2661, https://doi.org/10.5194/hess-16-2647-2012, 2012. a
David, N., Alpert, P., and Messer, H.: The potential of commercial microwave networks to monitor dense fog-feasibility study, J. Geophys. Res.-Atmos., 118, 750–761, https://doi.org/10.1002/2013JD020346, 2013. a
De Vos, L. W., Raupach, T. H., Leijnse, H., Overeem, A., Berne, A., and Uijlenhoet, R.: High-resolution simulation study exploring the potential of radars, crowdsourced personal weather stations, and commercial microwave links to monitor small-scale urban rainfall, Water Resour. Res., 54, 10293–10312, https://doi.org/10.1029/2018WR023393, 2018. a, b
Giannetti, F., Reggiannini, R., Moretti, M., Adirosi, E., Baldini, L., Facheris, L., Antonini, A., Melani, S., Bacci, G., Petrolino, A., and Vaccaro, A.: Real-time rain rate evaluation via satellite downlink signal attenuation measurement, Sensors, 17, 1864, https://doi.org/10.3390/s17081864, 2017. a
Gorgucci, E., Chandrasekar, V., Bringi, V. N., and Scarchilli, G.: Estimation of raindrop size distribution parameters from polarimetric radar measurements, J. Atmos. Sci., 59, 2373–2384, https://doi.org/10.1175/1520-0469(2002)059<2373:EORSDP>2.0.CO;2, 2002. a
Gosset, M., Kunstmann, H., Zougmore, F., Cazenave, F., Leijnse, H., Uijlenhoet, R., Chwala, C., Keis, F., Doumounia, A., Boubacar, B., Kacou, M., Alpert, P., Messer, H., Rieckermann, J., and Hoedjes, J.: Improving rainfall measurement in gauge poor regions thanks to mobile telecommunication networks, B. Am. Meteorol. Soc., 97, ES49–ES51, https://doi.org/10.1175/BAMS-D-15-00164.1, 2016. a
Hazenberg, P., Leijnse, H., and Uijlenhoet, R.: The impact of reflectivity correction and accounting for raindrop size distribution variability to improve precipitation estimation by weather radar for an extreme low-land mesoscale convective system, J. Hydrol., 519, 3410–3425, https://doi.org/10.1016/j.hydrol.2014.09.057, 2014. a
Islam, T., Rico-Ramirez, M. A., and Han, D.: Tree-based genetic programming approach to infer microphysical paramters of the DSDs from the polarization diversity measurements, Comput. Geosci., 48, 20–30, https://doi.org/10.1016/j.cageo.2012.05.028, 2012. a
Leijnse, H., Uijlenhoet, R., and Stricker, J.: Rainfall measurement using radio links from cellular communication networks, Water Resour. Res., 43, W03201, https://doi.org/10.1029/2006WR005631, 2007b. a
Mishchenko, M. I. and Travis, L. D.: Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers, J. Quant. Spectrosc. Ra., 60, 309–324, https://doi.org/10.1016/S0022-4073(98)00008-9, 1998. a
Mishchenko, M. I., Travis, L. D., and Mackowski, D. W.: T-matrix computations of light scattering by nonspherical particles: a review, J. Quant. Spectrosc. Ra., 55, 535–575, https://doi.org/10.1016/0022-4073(96)00002-7, 1996. a
Powell, M. J.: A hybrid method for nonlinear equations, in: Numerical methods for nonlinear algebraic equations, edited by: Rabinowitz, P., pp. 87–114, Gordon and Breach Science Publishers, New York, USA, 1970. a
Raupach, T. H. and Berne, A.: Correction of raindrop size distributions measured by Parsivel disdrometers, using a two-dimensional video disdrometer as a reference, Atmos. Meas. Tech., 8, 343–365, https://doi.org/10.5194/amt-8-343-2015, 2015. a
Raupach, T. H. and Berne, A.: Retrieval of the raindrop size distribution from polarimetric radar data using double-moment normalisation, Atmos. Meas. Tech., 10, 2573–2594, https://doi.org/10.5194/amt-10-2573-2017, 2017 (data available at: https://doi.org/10.14768/MISTRALS-HYMEX.721). a, b, c
Rincon, R. F. and Lang, R. H.: Microwave link dual-wavelength measurements of path-average attenuation for the estimation of drop size distributions and rainfall, IEEE T. Geosci. Remote, 40, 760–770, https://doi.org/10.1109/TGRS.2002.1006324, 2002. a
Salles, C. and Poesen, J.: Performance of an optical spesctro pluviometer in measuring basic rain erosivity characteristics, J. Hydrol., 218, 142–156, https://doi.org/10.1016/S0022-1694(99)00031-1, 1999. a
Salles, C. and Poesen, J.: Rain properties controlling soil splash detachment, Hydrol. Process., 14, 271–282, https://doi.org/10.1002/(SICI)1099-1085(20000215)14:2<271::AID-HYP925>3.0.CO;2-J, 2000. a
Salles, C., Poesen, J., and Borselli, L.: Measurement of simulated drop size distribution with an optical spectro pluviometer: Sample size considerations, Earth Surf. Process. Landf., 24, 545–556, https://doi.org/10.1002/(SICI)1096-9837(199906)24:6<545::AID-ESP3>3.0.CO;2-D, 1999. a
Thurai, M., Huang, G., Bringi, V., Randeu, W., and Schönhuber, M.: Drop shapes, model comparisons, and calculations of polarimetric radar parameters in rain, J. Atmos. Ocean. Tech., 24, 1019–1032, 2007. a
Tokay, A. and Short, D. A.: Evidence from tropical raindrop spectra of the origin of rain from stratiform versus convective clouds, J. Appl. Meteorol., 35, 355–371, https://doi.org/10.1175/1520-0450(1996)035<0355:EFTRSO>2.0.CO;2, 1996. a
Uijlenhoet, R. and Stricker, J.: A consistent rainfall parameterization based on the exponential raindrop size distribution, J. Hydrol., 218, 101–127, https://doi.org/10.1016/S0022-1694(99)00032-3, 1999. a
Uijlenhoet, R., Smith, J. A., and Steiner, M.: The microphysical structure of extreme precipitation as inferred from ground-based raindrop spectra, J. Atmos. Sci., 60, 1220–1238, https://doi.org/10.1175/1520-0469(2003)60<1220:TMSOEP>2.0.CO;2, 2003. a
Uijlenhoet, R., Porrà, J. M., Sempere Torres, D., and Creutin, J.-D.: Analytical solutions to sampling effects in drop size distribution measurements during stationary rainfall: Estimation of bulk rainfall variables, J. Hydrol., 328, 65–82, https://doi.org/10.1016/j.jhydrol.2005.11.043, 2006. a, b
Ulbrich, C. W.: Natural variations in the analytical form of the raindrop size distribution, J. Climate Appl. Meteorol., 22, 1764–1775, https://doi.org/10.1175/1520-0450(1983)022<1764:NVITAF>2.0.CO;2, 1983. a, b
van Leth, T. C., Overeem, A., Leijnse, H., and Uijlenhoet, R.: A measurement campaign to assess sources of error in microwave link rainfall estimation, Atmos. Meas. Tech., 11, 4645–4669, https://doi.org/10.5194/amt-11-4645-2018, 2018a. a, b, c, d
van Leth, T. C., Overeem, A., Uijlenhoet, R., and Leijnse, H.: Wageningen Urban Rainfall experiment, 4TU.Centre for Research Data, https://doi.org/10.4121/uuid:1dd45123-c732-4390-9fe4-6e09b578d4ff, 2018b. a
Vulpiani, G., Marzano, F. S., Chandrasekar, V., Berne, A., and Uijlenhoet, R.: Polarimetric weather radar retrieval of raindrop size distribution by means of a regularized artificial neural network, IEEE T. Geosci. Remote, 44, 3262–3275, https://doi.org/10.1109/TGRS.2006.878438, 2006. a
Zhang, G., Vivekanandan, J., and Brandes, E.: A method for estimating rain rate and drop size distribution from polarimetric radar measurements, IEEE T. Geosci. Remote, 39, 830–841, https://doi.org/10.1109/36.917906, 2001. a, b
Zhang, G., Vivekanandan, J., Brandes, E. A., Meneghini, R., and Kozu, T.: The shape-slope relation in observed gamma raindrop size distributions: statistical error of useful information?, J. Atmos. Ocean. Tech., 20, 1106–1119, https://doi.org/10.1175/1520-0426(2003)020<1106:TSRIOG>2.0.CO;2, 2003. a, b