The Re-Foundation of Radar Meteorology

The retrieval of the rainfall rate from dual polarimetric radar measurements assumes that there is a functional dependence (map) $\mathcal{M}$ of the rainfall rate $R$ on the radar variables [1, 2]. This approach extends traditional single-polarization radar methods for estimating rainfall rate from reflectivity, built on early work on raindrop size distributions by Marshall and Palmer [3]. The radar variables most commonly used are the horizontal polarization reflectivity $Z_h$, the ratio between the horizontal and vertical polarization reflectivities (differential reflectivity) $Z_{DR}$, and the differential phase $K_{DP}$. Almost all conventional rainfall estimators use scaling relations involving one or two radar variables [4]. Other approaches use all three observables ($Z_h$, $Z_{DR}$, $K_{DP}$) to retrieve DSD parameters and subsequently estimate rainfall rate [5]. In this manuscript, we focus our analysis on the commonly used couple ($Z_h$, $Z_{DR}$) [6, 7] although the methodology can be extended to any couple of variables. In this case, the retrieval algorithm is the so called ZR–relationship

Mathematically, the above equation implies that there is a unique value of the rainfall rate given the values of the horizontal and differential reflectivity. In practice, it means that the couple ($Z_h$, $Z_{DR}$) "reasonably" fixes the value of the rainfall rate. In topological terms, it means that, for all practical purposes, the rainfall phenomenon has two degrees of freedom.

The values of the rainfall rate and the radar variables are a function of the drop concentration per unit volume and the probability density function of drop diameters: the product of these two quantities is commonly referred to as the Drop Size Distribution (DSD): rain drops are not perfectly spherical, hence "size" is the diameter of the sphere with the same volume as the rain drop. The validity of Eq. (1) implies that two independent parameters are enough to quantify the statistical variability of the DSD. The most used DSD functional forms are the Gamma distribution (three parameters) [8] and the generalized Gamma model (four parameters) [9]. Studies such as [10, 11, 12], and the references therein, discuss methods for estimating DSD parameters and their interrelationships, while [13] examines their variability across climatic regimes and between convective and stratiform precipitation.

One minute disdrometer counts are commonly used to derive enough samples to estimate the triplet of parameters $(a, b, c)$. While the calculation of the rainfall rate is straightforward, the calculation of radar variables depends on the drops aspect ratio, drops orientation, drops velocity, air temperature, and angle of the incident radar beam. The introduction of the T-matrix algorithm [14, 15] has enabled a much more sophisticated and accurate estimation of the radar variables given a DSD, and thus the estimation of rainfall rate from dual polarimetric radar measurements.

In this manuscript, we build a statistical retrieval algorithm which outperforms the scaling functional relation of Eq. (1) even when the scaling condition is relaxed in favor of an empirically estimated one. The rationale is that no functional relationship exists, the relationship is statistical: the couple $(Z_h, Z_{DR})$ does not "reasonably" fix the value of the rainfall rate. To a point $z$ in the $Z_hZ_{DR}$-space is associated not an unique value of the rainfall rate but a distribution of possible values $f_z(R)$. Thus, building a retrieval algorithm based on a "mapping" of $R$ in the $Z_hZ_{DR}$-space requires the estimation $f_z(R)$ for all possible points $z$, with the possible complication that $f_z(R)$ might also depend on the synoptic condition. Over the years, more than 100 different ZR-relationships have been reported in Literature [16], and references therein, while Table 1 reports the variability observed for the 149 datasets used in our analysis. This "variability" of the ZR-relationship, in Literature, is usually solely attributed to synoptic effects. We think, that this interpretation is erroneous, and that while some synoptic effect might be present the "variability" of the ZR-relationship is an indication of the statistical nature of the relation between $R$ and the couple $(Z_h, Z_{DR})$.

Parameter	Mean	Std Dev	Min	Q1	Median	Q3	Max
a	0.0235	0.0064	0.0091	0.02078	0.0259	0.0284	0.0328
b	0.807	0.053	0.672	0.782	0.815	0.845	0.939
c	-6.4	2.4	-11.6	-8.4	-6.4	-4.4	-2.1

A very strong evidence supporting our thesis is that the degrees of freedom of the drop size distribution are at minimum three [17]. The Authors of [17] show that given a disdrometer drop count and not depending on synoptic conditions, the number of drops $N$, the mean $\mu$, and the skewness $\gamma$ of the probability density function of drop diameters $p(D)$ are three independent parameters that resolve $\ge95\%$ (via Principal Component Analysis) of the observed variability of a disdrometer count [18, 19, 20]. Because the drop size distribution has a minimum of three degrees of freedom, there cannot be any functional dependence of the rainfall rate from a couple of radar variables. Any determination of the rainfall rate via a couple of radar variables has an accuracy limit which can’t be surpassed due to the extra unresolved degree of freedom. Conversely, no retrieval algorithm of the rainfall rate based only on two radar variables can have an average error below a given threshold. This situation constitutes a sort of “rainfall uncertainty principle”, similar to the one introduced by Heisenberg for conjugate variables in quantum mechanics.

The statistical retrieval algorithm hereby present uses the DSD parametrization of [17] to build a better retrieval model than the $ZR$-relationship. To demonstrate it, we adopt a classical data science setup. Useful models need to reasonably fit the dataset used to build them (the training dataset) without a sensible reduction of accuracy on a dataset which was not used for training (the validation dataset). Thus, the winning model among a group of models trained on the same dataset is the one with the best validation performance. Hereby, we use 149 disdrometer datasets from 76 locations on Earth’s surface (Figure 1). In terms of the disdrometer type, our database is composed by 27 2dvd, 36 Parsivel, 30 Joss-Waldvogel RD69, and 56 Joss-Waldvogel RD80 datasets. When a dataset is used to build a model, the remaining 148 are considered as validation. In this way, for each model (rainfall rate retrieval algorithm) we have 149 training accuracies and 22,052 $(149 \times148)$ validation accuracies.

The three models adopted in this manuscript are:

1. Scaling Law Model (SLM): this is the scaling relationship of Eq. (1). Each data set is used to derive the parameters $(a, b, c)$ via least square fit in a log-log plot.

2. Direct Fit Model (DFM): each dataset is used to derive the empirical fit of the map $\mathcal{M}$ such that $R = \mathcal{M}(Z_h, Z_{DR})$. This is a generalization of Eq. (1) as we relax the scaling hypothesis. We do this to further test the fact that there is no functional dependency of the rainfall rate on the couple of radar variables considered: this is so because the degrees of freedom of the DSD are at minimum 3.

3. Statistical Retrieval Model (SRM): this model uses the fact that, given a disdrometer count, for any rainfall bulk variable $X$

where $I$ is an instrument factor, $N$ the drop count, and $\mathcal{F}_X$ is a shape factor depending only on the particular shape of $p(D)$. The first equality is nothing else but a useful factorization, e.g. given a disdrometer count the shape factor for the rainfall rate is just the third moment of the pdf $p(D)$. It is important to notice that $\mathcal{F}_X(p(D))$ is a functional dependence in the mathematical sense: there is a unique value associated to a particular shape $p(D)$. However, it is the second equality (demonstrated in [17]) that is fundamental for building the SRM as the shape factor of any bulk variable is just a function defined in the $\mu\gamma$-plane. The couple $(\mu, \gamma)$ being the mean and skewness of $p(D)$ accounts for $\gtrsim95\%$ (via Principal Component Analysis) of its variability [17]. Given a disdrometer dataset, the SRM uses Eq. (2) to build the shape factors of the horizontal reflectivity $F_{Z_h}(\mu,\gamma)$, differential reflectivity $F_{Z_{DR}}(\mu,\gamma)$, and rainfall rate $F_R(\mu,\gamma)$. Note that $Z_{DR}$ is independent from the drop count $N$ because the differential reflectivity is the ratio of two bulk variables. As matter of fact $Z_{DR} \equiv F_{Z_{DR}} = F_{Z_h}/F_{Z_v}$. The retrieval algorithm is as follows

• The observed value $z_{DR}$ of the differential reflectivity, and the shape factor $F_{Z_{DR}}(\mu,\gamma)$ are used to determine the set $S_{\mu\gamma}$ of couples $(\mu,\gamma)$ compatible with the observed value.

• The values of each couple $(\mu,\gamma)$ in $S_{\mu\gamma}$ are plugged into Eq. (2) for $X = Z_h$, and retrieve the set $S_N$ of drop counts compatible with the observed value of horizontal reflectivity $z_h$.

• The results of the two previous steps are used to derive the set $S_{N\mu\gamma}$ of all possible triplets $(N,\mu,\gamma)$ compatible with the observed values $(z_h,z_{DR})$.

• Finally, we use Eq. (2) with $X = R$ to obtain all possible values of the rainfall rate associated with a given dual polarimetric radar return echo: the estimate is the median.

To measure the train and validation accuracies of all the retrieval models, we use the bias $\mathcal{B}$ and dispersion $\mathcal{D}$. The bias of a model when predicting the values of the target variable is the median of all relative errors, while the dispersion is the interquartile range of all relative errors. The ideal model would have zero bias and zero dispersion: $(\mathcal{B}, \mathcal{D}) = (0,0)$. The train accuracy of a retrieval model built on the dataset $i$ is measured by the bias and dispersion when the model fits the rainfall rates of the dataset $i$ itself. We speak of validation accuracy of a retrieval model, when we compute the couple $(\mathcal{B},\mathcal{D})$ for the predictions of the rainfall rates on a dataset $j$ different from the training dataset $i$.

Overall, to compare the accuracies of the three rainfall rate retrieval models considered in this manuscript, we compute for each model all the possible couples $(B_{ij},D_{ij})$ where $i,j = 1,2,\ldots,K$ ($K$ being the number of datasets considered). Training accuracies are all the couples $(B_{ij},D_{ij})$ where $i = j$, validation accuracies are all the couples $(B_{ij},D_{ij})$ where $i \ne j$. We investigate the probability density function of these accuracies to determine the best model.

Results

Figure 2 shows how the couples $(B_{ij},D_{ij})$, as expressed in percent values, are distributed in the bias-dispersion plane within a regular grid of 1% size. Every grid square is colored according to the presence in the grid of bias-dispersion couples from each of the three models adopted: white if no couples are present in the grid, light blue if only couples relative to the SRM are present, light green if only couples relative to the DFM are present, light red if only couples relative to the SLM are present, all other colors if there is a mix of couples relative to different models (in which case the color reflects the proportion of each model). The left panel shows training ($i = j$) results, and the right panel validation results ($i \ne j$). The results reported are limited to couples of datasets relative to the same type of disdrometer: one row for each type of instrument. We clearly see, both in training and in validation and not depending on the type of instrument, how the SRM model has the smaller dispersions values, and that the bias-dispersion SRM couples are clearly separated from the other two models. The SLM model, which is the current state of the art retrieval methodology, is clearly the worst of the three models: the ZR-relationship fitting a dataset is generally a poor fit when applied to a different dataset. The DFM model is better than the SLM model as the scaling condition of Eq. (1) is relaxed.

Vector (PDF)

In Table 2, we examine the cross-instrument effect on the validation dispersion values for all three models. We report the median, the 5-th and 95-th percentiles of validation dispersions as organized by fitting instrument (the instrument of the training dataset), fitted instrument (the instrument of the dataset used for validation), and the model. We see how: 1) The SRM is still the best of the three models 2) The SLM is still the worst of the three models 3) The performance of each model degrades when validated on a dataset relative to a different instrument. 4) The largest cross-instrument degradation occurs between the 2DVD and the remaining instruments, while the degradation within the PARSIVEL, RD69, and RD80 group are more moderate. These cross-instrument errors may reflect differences in diameter binning, sensitivity to small drops, and instrument response. Intercomparisons have documented differences among Joss–Waldvogel, PARSIVEL, and 2DVD measurements [21, 22], and small-drop measurement limitations have also been reported for the 2DVD [9].

Fitting	Fitted	Model	50% $D_{ij}$	5% $D_{ij}$	95% $D_{ij}$	Samples
2DVD	2DVD	SRM	0.124	0.088	0.158	695
2DVD	2DVD	DFM	0.212	0.165	0.331	702
2DVD	2DVD	SLM	0.259	0.160	0.325	702
2DVD	PARSIVEL	SRM	0.247	0.176	0.360	950
2DVD	PARSIVEL	DFM	0.431	0.335	0.658	958
2DVD	PARSIVEL	SLM	0.477	0.339	0.613	972
2DVD	RD69	SRM	0.353	0.264	0.494	781
2DVD	RD69	DFM	0.413	0.325	0.642	809
2DVD	RD69	SLM	0.619	0.311	0.815	810
2DVD	RD80	SRM	0.354	0.252	0.528	1475
2DVD	RD80	DFM	0.405	0.310	0.685	1512
2DVD	RD80	SLM	0.542	0.306	0.775	1512
PARSIVEL	2DVD	SRM	0.147	0.098	0.216	972
PARSIVEL	2DVD	DFM	0.217	0.161	0.306	972
PARSIVEL	2DVD	SLM	0.297	0.212	0.404	972
PARSIVEL	PARSIVEL	SRM	0.127	0.102	0.174	1260
PARSIVEL	PARSIVEL	DFM	0.254	0.204	0.340	1260
PARSIVEL	PARSIVEL	SLM	0.352	0.273	0.457	1260
PARSIVEL	RD69	SRM	0.161	0.129	0.206	1078
PARSIVEL	RD69	DFM	0.227	0.194	0.258	1080
PARSIVEL	RD69	SLM	0.380	0.279	0.615	1080
PARSIVEL	RD80	SRM	0.154	0.118	0.204	2015
PARSIVEL	RD80	DFM	0.229	0.187	0.310	2014
PARSIVEL	RD80	SLM	0.372	0.259	0.614	2016
RD69	2DVD	SRM	0.200	0.140	0.280	810
RD69	2DVD	DFM	0.270	0.196	0.346	810
RD69	2DVD	SLM	0.406	0.267	0.458	810
RD69	PARSIVEL	SRM	0.168	0.117	0.250	1078
RD69	PARSIVEL	DFM	0.281	0.224	0.355	1080
RD69	PARSIVEL	SLM	0.404	0.320	0.485	1080
RD69	RD69	SRM	0.148	0.122	0.192	870
RD69	RD69	DFM	0.229	0.198	0.266	870
RD69	RD69	SLM	0.356	0.303	0.493	870
RD69	RD80	SRM	0.151	0.117	0.295	1680
RD69	RD80	DFM	0.229	0.198	0.287	1680
RD69	RD80	SLM	0.391	0.311	0.485	1680
RD80	2DVD	SRM	0.192	0.123	0.295	1511
RD80	2DVD	DFM	0.259	0.178	0.356	1512
RD80	2DVD	SLM	0.352	0.233	0.473	1512
RD80	PARSIVEL	SRM	0.157	0.112	0.251	1975
RD80	PARSIVEL	DFM	0.300	0.222	0.403	2008
RD80	PARSIVEL	SLM	0.391	0.303	0.511	2016
RD80	RD69	SRM	0.144	0.113	0.188	1679
RD80	RD69	DFM	0.230	0.191	0.305	1679
RD80	RD69	SLM	0.348	0.288	0.450	1680
RD80	RD80	SRM	0.143	0.109	0.188	3080
RD80	RD80	DFM	0.233	0.191	0.322	3080
RD80	RD80	SLM	0.355	0.270	0.460	3080

Conclusions

The Statistical Retrieval Model presented in this manuscript outperforms the standard Scaling Law retrieval Model because it does NOT attempt to fit any functional form to the relationship between the rainfall rate and the horizontal and differential reflectivity. The non functional (statistical) nature of this relation is responsible for the strong dependence of the parameters of Eq. (1) on the dataset adopted: there is not a fixed manifold where the triplets $(R, Z_h, Z_{DR})$ reside. Conversely, the shape factors of rainfall bulk variables have a real functional dependence on the shape of $p(D)$ ($\mathcal{F}_X(p(D))$) which can be estimated using any disdrometer dataset, no matter its synoptic mixture: there is a fixed manifold where the couples $(p(D), \mathcal{F}_X)$ reside. In practice, the estimation of the shape factors is possible because the shape of $p(D)$ is uniquely determined [17] by the mean and skewness of the distribution: $\mathcal{F}_X(p(D)) \to\mathcal{F}_X(\mu,\gamma)$.

In principle one could use the triplet of radar variables $(Z_h, Z_{DR}, K_{DP})$ to fully resolve all three degrees of freedom of the drop size distribution, but this is not the scope of this manuscript for two main reasons. 1) While given a disdrometer record, one can estimate any radar variable, the reality of radar measurements is that the three variables $(Z_h, Z_{DR}, K_{DP})$ cannot always be measured with an adequate accuracy [4]. Thus, there is a need to estimate the rainfall rate only with a couple of radar variables. 2) We want to prove that without imposing an analytical functional dependence to the Drop Size Distribution (Gamma distribution or generalized Gamma distribution), or adopting a scaling law among the rainfall rate and a couple of radar variables, it is possible to build an efficient and universal (reasonably good validation accuracy) rainfall retrieval algorithm.

We think that the Gamma distribution and the scaling law of Eq. (1) are vestiges of the historical development of radar meteorology offering an analytical treatment of the rainfall phenomenon: hence the wide utilization. However, as shown in [23, 16, 17], and in the present manuscript they do not properly capture its statistical properties. The SRM algorithm offers a proper statistical foundation to the remote sensing retrieval of rainfall rates, making the question "which ZR relationship should I use?" a question of the past.

Data and Methods

The datasets used in this manuscript are drawn from the collection described in [17]. They undergo the exact same processing methodology detailed in that article and its online supplement.

Building a SLM

Given a set of disdrometer observations the associated SLM model is built by calculating, via a least square fit, the parameters of Eq. (1) in the $(\log R, \log Z_h, \log Z_{DR})$ space. The predicted rainfall rate is obtained by plugging the observed values of the couple $(z_h, z_{DR})$ into Eq. (1).

Building a DFM

Given a set of disdrometer observations the associated DFM model is built by using the Local Adaptive Fit (LAF) methodology introduced in [17] for the rainfall rate $R$ as a function of the couple $(Z_h, Z_{DR})$. The LAF fit adopts a grid on the space of the independent variables. Therefore, the predicted rainfall rate is the value of the mapping $\mathcal{M}(Z_h, Z_{DR})$ for the closest grid point to the observed values of the couple $(z_h, z_{DR})$.

Building a SRM

Given a set of disdrometer observations we build the associated SRM model as follows. For each disdrometer count we calculate the octuplet $(N,\mu,\gamma,R,Z_h,Z_{DR},\mathcal{F}_{Zh},\mathcal{F}_R)$: drop count, mean and skewness of the pdf $p(D)$, rainfall rate, horizontal reflectivity, differential reflectivity, and the rainfall rate and horizontal reflectivity shape factors (obtained by inverting Eq. (2)). The values of $Z_h$ and $Z_{DR}$ are obtained using the drop terminal velocity–diameter relation of [24], and the pytmatrix Python package [25] with the following setting= radar wavelength=S-band, complex refractive indices of water at temperature=10°C, drop orientation distribution=gaussian with a standard deviation of 10° degrees, drop aspect ratio functional form as in [26], and all remaining settings = default settings.

The SRM model is built by using the octuplets to derive empirical curves for $Z_{DR}$, $\mathcal{F}_{Zh}$, and $\mathcal{F}_R$ as a function of the couple $(\mu,\gamma)$. This task is done using the Local Adaptive Fit (LAF) methodology introduced in [17]. These three empirical mappings are all that is needed to retrieve the rainfall rate given an observed horizontal and differential reflectivity couple:

• Given a differential reflectivity value $z_{DR}$, we use $Z_{DR}(\mu,\gamma)$ to obtain the set $S_{\mu\gamma}$ of couples $(\mu,\gamma)$ such that $(Z_{DR}(\mu,\gamma)-z_{DR})=0$. In practice this is done by retrieving the couples $(\mu_j,\gamma_j)$ such that the empirical fit of $Z_{DR}(\mu,\gamma)$ is within a small fraction $\epsilon$ of $z_{DR}$. Therefore, $S_{\mu\gamma}=\{(\mu_j,\gamma_j),\ j=1,\ldots,M\}$ such that $\lvert Z_{DR}(\mu_j,\gamma_j)-z_{DR}\rvert\leq\epsilon z_{DR}$. We set $\epsilon=0.01.

• Given a horizontal reflectivity value $z_h$, we use Eq. (2) to build the set $S_N = \{N_j,\ j = 1,\ldots,M$ such that $N_j = z_h/(I \times F_{Z_h}(\mu_j,\gamma_j))\ \forall\ (\mu_j,\gamma_j) \in S_{\mu\gamma}\}$.

• We combine the sets $S_{\mu\gamma}$ and $S_N$ to obtain the set $S$ of all the triplets $(N_j,\mu_j,\gamma_j)$ $(j = 1,\ldots,M)$ which are compatible with the observed values of horizontal reflectivity $z_h$ and differential reflectivity $z_{DR}$.

• We use Eq. (2) to obtain $M$ possible values of rainfall rate, and take the median value $r$ to be the predicted rainfall rate given $(z_h,z_{DR})$.

The SRM model finds the values of drop count and the shapes of the drop size distribution (and thus the possible rainfall rate values) that are statistically compatible with the observed reflectivities: the final predicted rate being the median value of the compatible rainfall rates. Note that $M$ is not a fixed constant as it depends on the particular value of the couple $(z_h,z_{DR})$.

Assessing model accuracy

A model $\mathcal{M}$ built on the dataset $i$ and used to fit the dataset $j$ has a bias $B_{ij}$ and dispersion $D_{ij}$ which are defined, respectively, as the median and the interquartile range of all relative errors. The training accuracy of a model is defined by the couple $(B_{ii},D_{ii})$. The case $i \ne j$ constitutes a validation accuracy. For each rainfall rate retrieval algorithm we can build 149 models (one for each dataset) leading to 149 possible values of training accuracy and $149 \times148 = 22,052$ possible values of validation accuracies.

Acknowledgment

We wish to sincerely thank colleagues and institutions who have generously provided disdrometer data: 1) Dr. C. R. Williams, the National Oceanic and Atmospheric Administration (NOAA) and Physical Sciences Laboratory (PSL) for the public availability of the datasets recorded at Alta (16), Bodega Bay (DSN 17), Erie (DSN 52), Cazadero (DSN 27 to DSN 38), Colfax (DSN 18 to DSN 26), Darwin (DSN 145), Davis (DSN 39), Estacion Obispo (DSN 50 and 51), Fort Ross (DSN 40), Grass Valley (DSN 41), Point Sur (DSN 48 to DSN 49), Rocky Mtn. Arsenal (DSN

54), Spanaway (DSN 10), Westport (DSN 11 to DSN 15); 2) Dr. Ali Tokay and National Aeronautics and Space Administration (NASA) for the public availability of the following datasets: Wallops Precipitation Research Facility (DSN 91 to DSN 105), IFloodS Alabama disdrometer data (DSN 78 and 79), Crystal Face Florida disdrometer data (DSN 148), Kwajalein Disdrometer Data (DSN 1 to DSN 9), LPVex Disdrometer Data (DSN 141 to DSN 143), and as part of the Global Precipitation Measurement (GPM) mission Ground Validation campaigns: HYMEX (DSN 133 to DSN 140), IPHEX (DSN 74 to DSN 77, and DSN 80 to DSN 90), MC3E (DSN 55 to DSN 73); 3) the National Institute of Information and Communications Technology, Japan for the Kashima data set (DSN 147); 4) the Institute of Observational Research for Global Change together with the Japan Agency for Marine-Earth Science and Technology for the Bukit Koto Tabang data set (DSN 144); 5) Dr. Dan Brawn for the Hassel data set (DSN 146); 6) Dr. Martin Hagen of the Institut für Physik der Atmosphäre Deutsches Zentrum für Luft- und Raumfahrt, Wessling, Germany for the Macunaga data set (DSN 149); 7) Dr. Charles Wrench and the British Atmospheric Data Centre (BADC) for Chilbolton (DSN 106 to DSN 121) and Sparsholt (DSN 122 to DSN 131) data archives.