Polarimetric Radar Quantitative Precipitation Estimation

Abstract

Radar quantitative precipitation estimation (QPE) is one of the primary tasks of weather radars. The QPE quality was substantially improved after polarimetric upgrade of the radars. This study provides an overview of existing polarimetric methodologies for rain and snow estimation and their operational implementation. The variability of drop size distributions (DSDs) is a primary factor affecting the quality of rainfall estimation and its impact on the performance of various radar rainfall relations at S, C, and X microwave frequency bands is one of the focuses of this review. The radar rainfall estimation algorithms based on the use of specific attenuation A and specific differential phase K_DP are the most efficient. Their brief description is presented and possible ways for their further optimization are discussed. Polarimetric techniques for the vertical profile of reflectivity (VPR) correction at longer distances from the radar are also summarized. Radar quantification of snow is particularly challenging and it is demonstrated that polarimetric methods for snow measurements show good promise. Finally, the article presents a summary of the latest operational radar QPE products available in the US by integration of the information from the WSR-88D radars via the Multi-Radar Multi-Sensor (MRMS) platform.

1. Introduction

Quantitative precipitation estimation (QPE) is one of the primary tasks of weather radars. It is hard to overestimate a societal impact of accurate precipitation measurements in the context of global warming with the rapidly increasing frequency of devastating flash flood events, particularly those associated with landfalling hurricanes and typhoons causing massive damage to property and infrastructure as well as human fatalities. While accurate precipitation measurements are critical for storm warnings and hydrology, they are also important for assimilation into the NWP models and optimization of their microphysical parameterization.

Traditional radar methods of rain and snow estimation were based on the use of a single radar variable, radar reflectivity factor Z. Significant progress in radar QPE was achieved during last 2–3 decades with introduction of polarimetric weather radars which are capable of complementing measurements of Z with additional radar variables such as differential reflectivity Z_DR, specific differential phase K_DP, and most recently, specific attenuation A. The advantages of polarimetric radar measurements for rainfall estimation were demonstrated in a number of research studies back in the nineties and early 2000s [1,2,3,4,5,6,7]. It was shown that the use of multiple radar variables helps to reduce the QPE uncertainty caused by the drop size distribution (DSD) variability whereas the utilization of a specific differential phase K_DP allows to address the issues of radar miscalibration, attenuation, and partial beam blockage.

In the 2000s, the polarimetric technology and methodology for rainfall estimation matured to the point to justify their operational implementation on the nationwide weather radar networks. The US spearheaded this effort after successful completion of the Joint Polarization Experiment (JPOLE) demonstration project [8] which culminated in the polarimetric upgrade of all WSR-88D radars in 2013. Another important step was integration of the information from the WSR-88D radars using the Multi-Radar Multi-Sensor (MRMS) platform that combines the polarimetric radar data with atmospheric environment data, satellite data, and lightning and rain gauge observations to generate a suite of severe weather and QPE products [9].

A comprehensive overview of the modern polarimetric radar methods for rain and snow estimation can be found in the recent monograph of Ryzhkov and Zrnic [10] (their Chapter 10). The purpose of this review paper is to briefly summarize the QPE methodologies described in [10] and to provide an update on the recent progress in polarimetric radar QPE and its operational implementation after the monograph was published. Although we attempted to provide a summary of all relevant research around the globe, the primary focus of this review is on the last decade developments in the US. In Section 2, we describe the evolution of polarimetric methods for rainfall estimation at the distances where the radar sampling volume is below the melting layer (ML) and is primarily filled with raindrops occasionally mixed with graupel and hail. The impact of the DSD variability on the performance of various rainfall relations is discussed along with the comparative analysis of their performance in rain of different microphysical types and intensity in the S, C, and X frequency bands. The emphasis is on the most recent rainfall estimation algorithm operationally implemented on the WSR-88D network based on the combined use of specific attenuation A, K_DP, and Z and its further optimization. At longer distances, where the radar samples the hydrometeors within the melting layer and above, the methods for the vertical reflectivity profile (VPR) correction are traditionally utilized. A description of the latest polarimetric VPR techniques is also presented in Section 2.

Radar measurements of snow are challenging due to tremendous variability of snow particle size distributions, their density, water content, shape, and orientation. Until recently, radar methods for estimation of snowfall rate or snow water equivalent (S or SWE) utilize multiple power-law relations between S and radar reflectivity Z. These relations exhibit roughly an order of magnitude difference in the estimates of snowfall rate for the same Z. Polarimetric methods for snow measurements have been introduced in the last few years and their description and performance results are provided and discussed in Section 3.

Section 4 contains a brief overview of the MRMS system for integration of the polarimetric radar data across the US and generation of various QPE products which are evaluated in real time with a multitude of surface rain gauges.

2. Rainfall Estimation

2.1. Polarimetric Rainfall Relations

The variability of drop size distributions is a primary source of uncertainty in the radar rainfall estimation. Because radar reflectivity Z is the 6th moment of DSD and rain rate R is approximately proportional to its 3.67th moment [10], the magnitude of Z is determined by the contribution of a few largest drops in the raindrop size spectrum whereas a bulk of smaller drops dominate rain rate. This dictated the use of multiple R(Z) relations for single-polarization radars. The choice of a particular relation depends on the type of rain (convective vs. stratiform, tropical vs. continental, etc.), season, and climate region. For example, the following five R(Z) relations are utilized on the US network of WSR-88Ds:

$$\begin{array}{l}Z=300R^{1.4} \mathrm{in} \mathrm{convective} \mathrm{rain},\\ Z=250R^{1.2} \mathrm{in} \mathrm{tropical} \mathrm{rain},\\ Z=200R^{1.6} \mathrm{in} \mathrm{summer} \mathrm{stratiform} \mathrm{rain},\\ Z=130R^{2.0} \mathrm{in} \mathrm{winter} \mathrm{stratiform} \mathrm{rain} (\mathrm{eastern} \mathrm{US}),\\ Z=75R^{2.0} \mathrm{in} \mathrm{winter} \mathrm{stratiform} \mathrm{rain} (\mathrm{western} \mathrm{US}).\end{array}$$

(1)

The choice of a particular R(Z) relation is at the discretion of a radar operator/meteorologist and there are no objective criteria for choosing an optimal one for a provided time and location.

Historically, the first polarimetric rainfall relations used a combination of Z and Z_DR [11]. Utilization of Z_DR, a polarimetric radar variable that depends on the median volume diameter of raindrops and is insensitive to their concentration, together with Z was expected to reduce the sensitivity of the rain rate estimate to the DSD variability. However, Z_DR is even more affected by a few largest raindrops than Z. Therefore, the benefit of using the R(Z, Z_DR) estimator instead of the R(Z) relation turned out to be relatively marginal. Utilization of the radar variables such as K_DP and A provided better improvement because both of these variables are proportional to the DSD moments that are much closer to the 3.67th moment than Z and are immune to the radar calibration errors, attenuation in rain, and partial beam blockage.

The sensitivity of the R(Z), R(Z, Z_DR), R(K_DP), and R(A) estimates to the DSD variability at S, C, and X bands is illustrated in Figure 1, Figure 2 and Figure 3. Scatterplots of rain rates versus their estimates from the four optimal relations were generated using a large disdrometer dataset obtained in Oklahoma [12]. It is obvious from all three figures that the R(Z, Z_DR) scatterplots are only slightly narrower than the R(Z) scatterplots and combining Z and Z_DR may not result in tangible QPE improvement. Moreover, additional measurement errors of Z_DR associated with its possible miscalibration and differential attenuation can make the performance of the R(Z, Z_DR) algorithm even worse than the R(Z) relation. Nevertheless, certain advantages of using R(Z, Z_DR) relations at S and C bands were reported in a number of studies [5,6,13,14,15,16,17,18,19,20]. Many of these relations can be found in Ryzhkov and Zrnic [10] (their table 10.3). It is not instrumental to utilize a combination of Z and Z_DR for rainfall estimation at X band because their biases caused by attenuation/differential attenuation are too large and difficult to account for.

The R(K_DP) scatterplots in Figure 1, Figure 2 and Figure 3 exhibit obvious benefit of the R(K_DP) relations for estimation of moderate-to-heavy rain at all three microwave frequency bands. A simple rule can be formulated as follows: the R(K_DP) relation is advantageous for measuring rain rate R if the value of R expressed in mm/h exceeds the value of radar wavelength λ in cm (e.g., [21]). This means that the R(K_DP) estimator is efficient for R > 10–11 mm/h at S band, for R > 5–6 mm/h at C band, and for R > 3 mm/h at X band. Utilization of R(K_DP) at lower rain rates is generally counterproductive because the R(K_DP) estimator becomes very sensitive to the DSD variability as can be seen from Figure 1, Figure 2 and Figure 3 and the estimates of K_DP itself become very noisy. Such a noisiness, however, can be reduced by averaging of K_DP over a large spatial/temporal domain and the R(K_DP) algorithm remains efficient even in light rain for areal rainfall estimation [22] or monthly rain accumulation [23].

A more quantitative assessment of relative performance of different rainfall relations can be obtained using Figure 4 where the dependencies of the fractional mean absolute error (fMAE) on rain rate for the R(Z), R(Z, Z_DR), R(K_DP), and R(A) relations optimized for the Oklahoma DSDs at S, C, and X bands are shown. The advantages of all three polarimetric relations compared with the R(Z) relation are obvious in a wide range of rain intensities at all three microwave bands. The R(A) relation (which will be discussed in more detail in the next Section 2.2) is an apparent winner at low-to-moderate rain rates. The R(A) advantage is particularly impressive at S band where the fMAE caused by the variability of DSD is less than 10% in a wide range of rain intensities. For higher rain rates exceeding 10 mm/h, fMAE tends to increase at C and X bands for all relations except R(K_DP) that is least sensitive to the effect of the resonance scattering on large raindrops which is primarily responsible for such an increase. Therefore, the K_DP-based relation is an obvious choice for heavy rain at C and X band.

Because K_DP is proportional to the radial derivative of differential phase Φ_DP, it is immune to radar miscalibration, attenuation, partial beam blockage, and presence of water or snow on the surface of a radome. The immunity of K_DP to attenuation makes it particularly beneficial for rainfall measurements at shorter wavelengths. The R(K_DP) relation is commonly used if rain is mixed with hail as was originally recommended by Balakrishnan and Zrnic [24]. Until recently, it was believed that hail does not contribute much to K_DP. This is not the case for melting graupel and water-coated hail of small size and high concentration. Moreover, small, and wet graupel/hail can produce anomalously high K_DP well-exceeding typical values of K_DP in rain. For example, Kumjian et al. [25] measured K_DP as high as 17 deg/km at S band in a storm producing large accumulations of small hail. If such a value of K_DP is converted to rain rate using a standard S-band relation R(K_DP) = 44.0|K_DP|^0.822, it would yield an unrealistic rain rate of 452 mm/h. In order to mitigate such an obvious overestimation, the R(K_DP) with a smaller multiplier is recommended if the cross-correlation coefficient ρ_hv drops below a certain level (0.97) in areas of high Z indicating the presence of hail [26], see Section 4. The K_DP-based algorithms for rainfall estimation were used in a large number of studies and operational applications (e.g., [1,2,3,5,6,7,13,14,15,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41]). A long list of various R(K_DP) relations can be found in [10] (their Table 10.3).

2.2. Attenuation-Based Rainfall Retrievals

A game-changer during the last decade was the introduction of the rainfall estimation methodology based on specific attenuation A [20,26,39,41,42,43,44,45,46,47,48]. Because this methodology is relatively new, herein we provide its description in more detail compared with other rainfall estimation techniques. A radial profile of A can be retrieved from the radial profile of attenuated Z_a and two-way path integrated attenuation PIA along the propagation path (r₁,r₂) in rain as [42]

$$A(r)=\frac{{[Z_a(r)]}^{b}C(b,PIA)}{I(r_1,r_2)+C(b,PIA)I(r,r_2)}$$

(2)

where

$$I(r_1,r_2)=0.46 b{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{[Z_a(s)]}^{b}ds}$$

(3)

$$I(r,r_2)=0.46 b{\displaystyle \underset{r}{\overset{r_2}{\int }}{[Z_a(s)]}^b}ds$$

(4)

$$C(b,PIA)=\exp (0.23bPIA)-1$$

(5)

In Equations (2)–(5), the factor b is a constant parameter (usually between 0.6 and 0.9 at microwave frequencies). The value of PIA cannot be estimated with a single-polarization radar but can be computed from the total span of a differential phase ΔΦ_DP measured by a dual-polarization radar:

$$PIA(r_1,r_2)=\alpha [{\Phi }_{DP}(r_2)-{\Phi }_{DP}(r_1)]=\alpha \Delta {\Phi }_{DP}$$

(6)

where α is equal to the average ratio of A and K_DP along the propagation path.

It follows from Figure 1, Figure 2, Figure 3 and Figure 4 that the R(A) estimate is much less affected by the DSD variability because the R(A) scatterplots are significantly narrower than the ones for the R(Z), R(Z, Z_DR), and R(K_DP) relations, particularly at S band. Therefore, the R(A) algorithm has a chance to work best of all these relations at S band where attenuation in rain is lowest. At first glance, this seems counterintuitive. However, even small path integrated attenuation PIA can be estimated reliably because of its strong correlation with ΔΦ_DP according to Equation (6). The R(A) relation is almost linear at S band:

$$R(A)=4120A^{1.03}$$

(7)

where λ = 11.0 cm and T = 20 °C ([42]). It deviates from linearity at C and X bands because of the more pronounced effects of resonance scattering on raindrops at shorter wavelengths. As a result, the effectiveness of utilizing specific attenuation at C and X bands for rainfall estimation is lower than at S band. It is evident from Figure 2, Figure 3 and Figure 4 that the R(A) relation has apparent advantages for lower rain rates whereas the R(K_DP) relation seems to be optimal for heavy rain. One of the advantages of the R(A) estimates compared with R(K_DP) is its better radial resolution which is similar to the one of Z and R(Z).

The radar retrieval formulas used for plotting scattergrams in Figure 1, Figure 2 and Figure 3 are optimized for a disdrometer dataset obtained in Oklahoma. The parameters of radar rainfall relations vary depending on rain type, therefore such relations should be adjusted using DSD data collected in a particular climate or geographical area. Microphysical characteristics of rain may vary in a wide range even in a particular geographical region such as Oklahoma. Different types of rain can be distinguished using the average dependencies of Z_DR on Z. Figure 5 shows Z–Z_DR dependencies corresponding to different percentiles of Z_DR for a given Z in rain simulated from 47114 DSDs measured in Oklahoma. The simulations are for S band at T = 20 °C. A tropical rain, a bulk of which is generated via a “warm rain” process below the melting layer, is commonly characterized by high concentration of small- and medium-size drops with a relative deficit of large drops which results in a lower value of Z_DR at a given Z as in the 0–20% range of Z_DR between two lowest curves in Figure 5. A continental rain, mostly originated from melting ice, is usually dominated by large raindrops with a lower concentration of small drops as opposed to the tropical rain and is characterized by higher values of Z_DR for a given Z. The two families of DSDs belonging to the 0–20% and 80–100% ranges can be designated as “very tropical” and “very continental” rain. In Figure 6, the log(R) vs. A scatterplots corresponding to the very tropical rain at S, C, and X bands are marked by red symbols, whereas the very continental scatterplots are presented in blue symbols. It is evident that the clusters associated with these two rain types almost coincide at S band and become more separated at C and X bands. The R(A) relations are especially beneficial for lower rain rates and for tropical rain of any intensity (red clusters). At C and X bands, the separation between red and blue clusters identifies that different R(A) relations must be utilized for tropical and continental rain. For a given rain rate, the value of A is larger for a continental rain which means that the multiplier in the R(A) relation should be lower for a continental rain. It can be shown that a similar rule is applied for the R(K_DP) and R(Z) relations.

Specific attenuation is a strong function of radar wavelength λ and temperature T; therefore, the parameters of the power-law R(A) relation are also wavelength- and temperature-dependent. At S band, this, however, does not preclude the use of a single R(A) relation specified by Equation (7). As is shown in [10], the value of A estimated from Equation (2) is approximately proportional to the magnitude of the parameter α which is proportional to A and depends on λ and T in a similar manner. This leads to almost complete cancellation of both dependencies in the final estimate of R which allows using a single R(A) relation provided by Equation (7). Note, that such cancellation is incomplete at C and X bands and a single R(A) relation is not as accurate as at S band.

A very important advantage of the K_DP- and A-based QPE algorithms is their immunity to the radar miscalibration, attenuation in rain, partial beam blockage, and impact of a wet radome. The immunity to partial beam blockage is illustrated in Figure 7 where the fields of 6-h rain total computed from R(Z) and R(A) are displayed. Large negative bias in the rain totals retrieved from R(Z) associated with partial beam blockage caused by trees near the radar is completely eliminated in the rain field estimated from R(A).

All discussed radar rainfall relations have advantages and disadvantages at different rain rates and radar wavelengths. Hence, it is beneficial to use composite algorithms by combining various relations. The operational polarimetric rainfall estimation algorithm currently implemented on the WSR-88D radar network in the US uses a combination of the R(Z), R(A), and R(K_DP) relations [26,46,47]. The R(K_DP) relations are used in areas of high Z where rain may be mixed with hail and the R(A) relation is utilized otherwise. The QPE algorithm resorts to various R(Z) relations if the radar sample is within the melting layer or above or if the total span of differential phase ΔΦ_DP in rain is too small or beam blockage is severe. The details of the operational WSR-88D QPE algorithm (labeled Q3DP) are in [26] and Section 4 of this review.

The Q3DP algorithm became operational in October 2020 and has demonstrated very good performance so far, particularly in heavy rain associated with flash floods. Figure 8 illustrates its performance for the most significant flash flood events of the 2021 warm season in the US. One of them occurred in Tennessee on August 21 and resulted in 27 fatalities and another one was associated with Hurricane Ida during the period from August 30 to September 2 with a death toll of 67 people. The Hurricane Ida had its landfall in Louisiana and transformed from a tropical depression to a post-tropical cyclone that rolled into the US northeast two days later. Color-coded maps in the upper images in Figure 8 show 24-h radar-estimated rain totals and the circles, their comparison with rain gauges. The size of the circle is proportional to the rain amount whereas different shades of red and blue indicate various levels of radar underestimation and overestimation, respectively. White or pale colors indicate either the absence of bias or a small bias. Indeed, the scatterplots of 24-h totals versus their radar estimates shown in the bottom panels indicate very good performance in all three cases with the average bias below 10% and the correlation coefficient well-above 0.9.

2.3. Further Optimization of the R(A) and R(K_DP) Relations

Although the polarimetric rainfall estimators R(A) and R(K_DP) have undeniable advantages compared with the old Z-based methods, they may require further optimization. It was mentioned before that the R(A) relation is practically insensitive to the DSD variability (at least at S band) but the estimate of A from Equations (2)–(6) depends on DSD because the factor α is a strong function of Z_DR (Figure 9). Such a dependence of α on Z_DR at S band and temperature 20 °C can be approximated by equation:

$$\alpha (Z_{DR})=0.008+\frac{0.009}{Z_{DR}-0.03}$$

(8)

where Z_DR > 0.2–0.3 dB. In Equation (8), Z_DR is expressed in dB. In the existing R(A) algorithm, the parameter α is optimized on the scan-to-scan basis and determined from the slope of the Z_DR(Z) dependence as described in [10,26,46]. The use of a single factor α for a whole radar coverage area generally leads to some underestimation of light rain and occasional overestimation of heavy rain. This is attributed to the fact that light rain is commonly characterized by larger α whereas significantly lower α might be more typical for heavier rain often associated with deep convection, and the use of a single “net” α (or <α>) for both rain types produces such biases (see Figure 9).

There are a couple of possible approaches to address this issue and to take into the account an azimuthal dependence of <α>. One of them is proposed in [49] using variational principles. According to [49], the state vector X consists of the “true” values of Z and Z_DR (not biased by attenuation/differential attenuation) along the radial in rain and the observation vector Y is composed of the radial profiles of the measured (attenuated) values of Z(Z^(m)) and Z_DR(Z_DR^(m)) and total differential phase Φ_DP^(m). The “forward operator” includes analytical relations for computing simulated K_DP^(s), A^(s), specific differential attenuation A_DP^(s), and simulated differential phase Φ_DP^(s) from the “true” values of Z and Z_DR. The components of the state vector are retrieved using the variational routine to match the simulated and measured vectors Y. After the “true” values of Z and Z_DR are determined_, the corresponding values of A and K_DP along each ray are computed using the forward operator, and the net <α> for a provided azimuth is estimated as the ratio ΣA⁽ⁱ⁾/ΣK_DP⁽ⁱ⁾ where summation is performed along the propagation path in rain. Such a routine is computationally expensive, requires well-calibrated values of the measured Z^(m) and Z_DR^(m) and its operational utilization must be demonstrated.

A possible alternative methodology implies the use of the radial profile of α(Z_DR) for estimation of <α> separately for each azimuth. The “net alpha” for a provided radial of data and propagation path in rain (r₁,r₂) is a ratio:

$$<\alpha >=\frac{PIA}{\Delta {\Phi }_{DP}(r_1,r_2)}=\frac{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}A(r)dr}}{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{K}_{DP}(r)dr}}=\frac{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}A(r)dr}}{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}\frac{A(r)dr}{\alpha [Z_{DR}(r)]}dr}}$$

(9)

Assuming A = aZ^b, Equation (9) can be rewritten as:

$$<\alpha >=\frac{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{Z}^b(r)dr}}{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}\frac{Z^b(r)}{\alpha [Z_{DR}(r)]}dr}}$$

(10)

where Z is the “true” reflectivity not biased by attenuation. The net α can be also estimated using a different formula by using the K_DP = cZ^d relation instead of the A(Z) relation.

$$PIA=2{\displaystyle \underset{r_1}{\overset{r_2}{\int }}A(r)dr=2{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{K}_{DP}(r) \alpha [Z_{DR}(r)]dr}}$$

(11)

Following [50], K_DP can be expanded as:

$$K_{DP}(r)=\frac{\Delta {\Phi }_{DP}(r_1,r_2)}{2{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{[Z(r)]}^{d}dr}}{[Z(r)]}^d$$

(12)

which after substituting into Equation (11) yields:

$$<\alpha >=\frac{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{[Z(r)]}^d\alpha [Z_{DR}(r)]dr}}{{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{[Z(r)]}^{d}dr}}$$

(13)

The values of Z and Z_DR in Equations (10) and (13) should be first corrected for attenuation/differential attenuation.

Azimuthal profiles of <α> estimated using Relations (10) and (13) for the case of a mesoscale convective system are presented in Figure 10 (upper right panel) using black and red dots respectively. It is obvious that Relations (10) and (13) provide very similar results. In the middle and bottom right panels, the azimuthal dependencies of Z and Z_DR averaged over the whole propagation path in rain are displayed. As expected, the azimuthal dependence of <α> exhibits strong anti-correlation with averaged Z_DR and the absence of correlation with averaged Z. The lowest value of <α> is observed in the azimuthal direction along the squall line in the NE sector.

This routine for <α> optimization, if proven efficient, can replace the existing rather cumbersome “Z_DR slope” methodology that does not allow capturing the azimuthal variability of <α>. Note that the value of <α> obtained using this technique can be used for polarimetric attenuation correction of Z that utilizes the relation ΔZ = <α> Φ_DP. Most of the existing attenuation correction techniques assume a fixed value of <α> whereas the actual value of <α> can vary quite significantly with time and azimuth as Figure 10 indicates. The value of <α> obtained from Equations (10) and (13) may also be utilized for optimization of the R(K_DP) relations because the parameters of these relations generally depend on DSDs and rain type. The multiplier a₁ in the $R=a_1{K}_{DP}^{b_1}$ power-law relation should be higher for tropical rain characterized by larger <α> and lower for continental rain with smaller <α>. It can be shown that the value of a₁ is equal to 35.8, 50.1, and 62.7 for b₁ = 0.82 for <α> = 0.015, 0.030, and 0.050 dB/deg respectively. Such tuning of the R(K_DP) relations based on the value of <α> is strongly recommended.

2.4. Polarimetric VPR

Once the radar beam intercepts the melting layer (ML) and snow/ice aloft, the relation between reflectivity (and other variables) and rain rate at the surface becomes complicated. Because Z in the ML is usually higher than in rain near the surface, the R(Z) relation valid in pure rain inevitably overestimates the surface rain rate at the distances where the radar samples ML. At longer distances where the radar resolution volume is in ice/snow above the ML, underestimation takes place. Various methods have been developed to address the melting layer problem (e.g., [51,52,53,54]). Common to these is the estimation of a vertical profile of reflectivity (VPR) and its use to retrieve rain rate near the ground. One of the ways to obtain the VPR is to scan the storm, reconstruct its structure near the radar, and extrapolate its effect to long range assuming horizontal homogeneity [53,55]. Such volume-scan VPR corrections are often part of operational algorithms for QPE. Andrieu and Creutin [56] and Andrieu et al. [57] added sophistication by applying discrete inverse theory based on hourly averaged VPRs at two elevation angles (1° and 3°). Local VPRs correction over small areas (20 km in range and 20° in azimuth) using an inverse filtering of the beam effects was proposed by Vignal et al. [58] and proven to be better in Switzerland [59] and US [60] than “global” corrections.

The concept of an “apparent” VPR or AVPR was introduced by Zhang and Qi [61] and Qi et al. [62]. The VPR correction is performed in the “bright band affected area” (BBA) within the stratiform part of the storm that is identified using radar reflectivity Z and cross-correlation coefficient ρ_hv. The average radial profile of Z at the lowest antenna tilt is determined within BBA and subtracted from the radial profiles of Z at each particular azimuth. Then the standard R(Z) relation is applied to estimate rain rate. Therefore, it is expected that the radial dependency of Z is entirely determined by the vertical profile of Z through the ML and the vertical structure of a storm is horizontally uniform within BBA. This methodology assumes that the BBA is relatively homogeneous azimuthally and is called “tilt-VPR”. Such an assumption can fail in the presence of frontal boundaries with large variations in the ML height which usually happens during cold season with lower ML. Hanft and Zhang [63] suggested the so-called dpVPR methodology to derive azimuthally dependent VPRs which results in the overall improvement of the VPR correction as illustrated in Figure 11. In the absence of the VPR correction, the 24-h radar derived rain totals are greatly overestimated (blue circles in Figure 11c,f). The positive QPE bias is dramatically reduced with tilt-VPR correction applied (Figure 11a,d) and the QPE results are further improved after the dpVPR correction scheme is utilized (Figure 11b,e). Note that, in addition to reducing the artificial enhancement of Z and R(Z) in BBA, the dpVPR algorithm also reduces the negative bias in R(Z) at longer distances from the radar where the radar beam is above the ML.

An alternative methodology for polarimetric VPR (PVPR) implies determination of the heights of the ML top and bottom over large areas using the recent melting layer detection algorithm (MLDA) described in [64]. According to this methodology, a multitude of the model radial profiles of Z and ρ_hv at lowest antenna tilts for different heights and “strengths” of the ML are simulated and stored in lookup tables. The “strength” of the ML characterizes both its thickness and the depth of the corresponding ρ_hv minimum that is usually well-correlated with the Z enhancement in the ML. The simulation routine explicitly takes into account radar beam broadening that has a strong impact on the ML signature at lower antenna elevations and large distances from the radar. Then an appropriate model profile of Z at every azimuth and elevation is determined by finding the best match between the model and observed profile of ρ_hv. Because all simulated model profiles of Z and ρ_hv are labeled by the intrinsic height and strength of the ML, such a matching routine allows for estimating the parameters of the ML such as the heights of its top and bottom. In addition, a model radial profile of Z associated with these parameters of the ML is identified and subtracted from the measured radial profile of Z at a particular azimuth. The difference between the PVPR and dpVPR routines is that PVPR is essentially coupled with the MLDA and uses the model radial profiles of Z whereas dpVPR is based on the piecewise linear approximation of the measured average radial profile of Z in a narrow azimuthal interval without quantification of the ML height and depth. Both VPR techniques have the same deficiency, an assumption of radial uniformity of Z and corresponding rain rate within the range interval affected by the “bright band contamination”. The tilt-VPR and dpVPR have been validated for a significant number of rain events [61,62,63] while massive validation of PVPR is in order. The performance of the PVPR algorithm is illustrated in Figure 12 and Figure 13 for the rain event with a low ML that is very well-identified in the QVPs of Z, Z_DR, ρ_hv, and K_DP (Figure 12). The maps of hourly rain total retrieved from R(Z) in the area affected by ML before and after PVPR correction of Z measured by the KLSX WSR-88D radar are shown in Figure 13. Artificial enhancement of the estimated rain accumulation due to the increase in radar reflectivity in the ML clearly visible in the left panel is completely eliminated after PVPR is utilized (right panel).

Concluding this section, we must mention recent studies where artificial intelligence (AI)/machine learning (ML) algorithms were used for polarimetric QPE. This is a rapidly growing area of research stimulated by accumulation of large datasets of polarimetric radar variables combined with surface rainfall measurements and facilitated by the availability of abundant modern computer resources. A detailed analysis of the AI/ML polarimetric QPE algorithms is beyond the scope of this review and we cited recent articles by Zhang et al. [65], Shin et al. [66], and Wolfensberger et al. [67] as examples for an interested reader.

3. Polarimetric Measurements of Snow

3.1. Historical Overview

Radar snow measurements present an enormous challenge due to immense variability of snow particle size distributions (PSDs), density, shapes, orientations, crystal habits, and terminal velocities. Historically, radar reflectivity factor at horizontal polarization Z was used to estimate snow water equivalent rates S in the form of numerous power-law relations (e.g., [68,69,70,71,72,73,74,75]). The majority of these relations assume that Z is proportional to S². Radar reflectivity is very close to the 4th PSD moment for low-density snow (e.g., aggregates) due to the inverse dependence of snow density on the snowflake size whereas S is proportional to significantly lower the order moment of PSD (2.1–2.2) [10]. This is one of the primary problems with existing S(Z) relations. The accuracy of S(Z) estimates is altitude dependent. At longer distances from the radar where the radar beam is much higher, the performance degrades due to the significant changes in Z with height, especially in aggregated snow. Despite these deficiencies, the conventional radar algorithms still utilize the standard S(Z) relations, which exhibit roughly an order of magnitude difference in the estimates of snowfall rate for the same reflectivity factor. This traditional approach is based on the usage of a multitude of climatological S(Z) relations optimized for different snow types and geographical areas. For example, the US National Weather Service (NWS) operationally uses several S(Z) relations for radar snow QPE (depending on the geographical region) as presented in

Table 1. Summary of Z(S) relations for dry snow listed in the literature and utilized by the WSR-88D network in the United States.

Source	Z(S) Relation for Dry Snow
Gunn and Marshall [68]	Z = 448 S2
Sekhon and Srivastava [69]	Z = 399 S2.21
Matrosov et al. [72]	Z = (100–130) S(1.3–1.55)
Huang et al. [76]	Z = (106–305) S(1.11–1.92)
Saltikoff et al. [77]	Z = 100 S2
Szyrmer and Zawadzki [73]	Z = 494 S1.44
Wolfe and Snider [78]	Z = 110 S2
Huang et al. [79]	Z = (130–209) S(1.44–1.81)
Von Lerber et al. [80]	Z = (53–782) S(1.19–1.61)
WSR-88D, Northeast	Z = 120 S2
WSR-88D, Great Lakes	Z = 180 S2
WSR-88D, North Plains/Upper Midwest	Z = 180 S2
WSR-88D, High Plains	Z = 130 S2
WSR-88D, Inter-mountain West	Z = 40 S2
WSR-88D, Sierra Nevada	Z = 222 S2

and it has not yet capitalized on the polarimetric capability of the WSR-88D radars.

Various alternatives to a Z-based approach for snow QPE were recently explored through the designated field campaigns. The GPM cold-season Precipitation Experiment (GCPEx, [81]), conducted by the National Aeronautics and Space Administration (NASA, Washington, DC, USA) and Environment Canada in Ontario, Canada, was one of such campaigns with a goal to investigate the utility of the multi-frequency (active and passive) microwave sensors for better classification and quantification of snow. Huang et al. [82] utilized a dual-wavelength ratio (DWR) of reflectivities at Ku and Ka bands available from the GCPEx data to develop the S(Z) relation parametrized by DWR. The difference between Zs at Ku and Ka bands increases with the characteristic size of snowflakes and the use of a multiplier in the S(Z) relations depending on DWR, helps to reduce the PSD-related uncertainty in the radar snow estimates.

Another campaign that significantly expands our knowledge of the snow microphysics was Biogenic Aerosols—Effects on Clouds and Climate (BAECC) 2014. Von Lerber et al. [80] utilized the snow data from BAECC to investigate the connection between the snow properties and radar observations using a video disdrometer and derived event-based power-law mass-size relations for ice/snow along with the corresponding S(Z) relations. They linked the changes in the S(Z) relation’s exponent to the variations in the mass-size relation’s exponent. The multiplier of the S(Z) relation is determined by the prefactors of the mass-size and velocity-size relations and by the intercept parameter N₀ of the particle size distribution. The changes in the S(Z) multiplier for a given N₀ are also associated with the changes in a liquid water path, a proxy for the degree of riming. In another study which partially used BAECC data in conjunction with the data from the University of Helsinki Hyytiälä station, Tiira and Moisseev [83] developed the unsupervised method for classification of vertical profiles of polarimetric radar variables using k-means clustering. They showed that the profiles of radar variables can be grouped into 16 snow classes which capture the most important snow growth and ice cloud processes. This can be used for the development of the snow class-based S(Z) relations that capture the variability of snow PSD, density, and fall velocity.

3.2. Radar Polarimetric Relations for Snow Estimation

Radar polarimetry opens new horizons for snow habit classification and quantification (e.g., [83]). Since the advance of polarimetry, only a handful of studies have explored polarimetric quantitative snow retrievals, mostly focused on ice water content (IWC) quantification. Vivekanandan et al. [84] and Lu et al. [85] utilized a specific differential phase K_DP for estimation of IWC. Aydin and Tang [86] investigated possible methods for IWC retrievals in the clouds composed of pristine ice crystals using the combination of differential reflectivity Z_DR and K_DP. Ryzhkov et al. [87,88] and Ryzhkov and Zrnic [10] used a combination of K_DP and Z_DR or K_DP alone to quantify IWC of pristine or lightly-to-moderately aggregated ice. Using the instrumented aircraft, Nguyen et al. [89] derived the empirical IWC(K_DP) and IWC(K_DP, Z_DR) relations from the airborne X-band polarimetric radar measurements and in situ microphysical probes. These empirical relations are very similar to the theoretical ones derived in [88]. Bukovčić et al. [90] introduced the first polarimetric algorithm (based on the K_DP and Z) for estimation of the extinction coefficient and visibility in snow which are traditionally linked to the snowfall rates (e.g., [91,92,93,94,95,96,97]). The algorithm, tested with S-band WSR-88D radar data, showed a good potential for quantification of extinction and visibility in snow.

Hassan et al. [98] proposed one of the first polarimetric relations for snowfall rate estimation using Z_DR and Z at C band. They showed that the addition of Z_DR to Z produces comparable results to standard S(Z) relations. Much more success in polarimetric snow QPE can be achieved if Z is combined with specific differential phase K_DP. A theoretical S(K_DP, Z) relation in a form of a bi-variate power-law was derived by Bukovčić et al. [99] assuming an exponential PSD, $N(D)=N_0\exp (-\Lambda D)$, and Rayleigh approximation:

$$S(K_{DP},Z)=\frac{c_s}{{(F_o{F}_s)}^{a_s}}{(K_{DP}\lambda )}^{a_s}{Z}^{b_s}$$

(14)

where λ is the radar wavelength in mm, and Z is in mm⁶m⁻³. It is also assumed that the snow density ρ_s and terminal velocity of snowflakes vs. are related to the particle equivolume diameter D as:

$${\rho }_s={\alpha }_s{D}^{-1}$$

(15)

and:

$$V_s={\tau }_s{D}^{{\delta }_s}$$

(16)

where in Equation (14), F_o and F_s are the factors characterizing orientation and shape of snowflakes and parameters a_s and b_s are the functions of δ_s whereas c_s is a function of α_s, δ_s, and the density of air ρ_a (or atmospheric pressure). The factor F_o is determined by the width of the canting angle distribution σ_s and the factor F_s is a function of the aspect ratio a/b of snowflakes specified in [10,99], see Figure 14. It can be shown that the theoretical Relation (14) is simplified as:

$$S(K_{DP},Z)=1.6K_{DP}^{0.62}{Z}^{0.38}$$

(17)

where λ = 11.0 cm, a/b = 0.65, σ_s = 0°, and α_s = 0.178 (as in [100]).

Relations (14) and (17) are derived for an exponential size distribution. A generalized empirical S(K_DP, Z) relation is derived from the extensive Oklahoma snow 2DVD measurements of PSDs and terminal velocities assuming wide ranges of particle aspect ratios (0.5–0.8) and the canting angle distribution width σ (0–40°) [99]. It has a form:

$$S(K_{DP},Z)=1.5K_{DP}^{0.62}{Z}^{0.33}$$

(18)

which is very close to the theoretical Relation (17). For practical purposes, an alternative S(K_DP, Z) relation can be used that explicitly takes into account the shape and orientation factors as well as air density (or atmospheric pressure p):

$$S(K_{DP},Z)=\frac{27.9\times {10}^{-3}}{{(F_o{F}_s)}^{0.62}}{\left(\frac{p_0}{p}\right)}^{0.5}{(K_{DP}\lambda )}^{0.62}{Z}^{0.33}$$

(19)

In Equation (19), p₀ = 1013 mb, K_DP is in deg/km, and λ is in mm.

The reason why the multivariate power-law S(K_DP, Z) relations outperform the ones based on a single Z is that the K_DP of snow is proportional to the first moment of its particle size distribution (PSD) if the bulk density of snow decreases with the size of snowflakes as specified by Equation (15) [10]. The product of K_DP (proportional to the first moment of PSD) and Z (proportional to its fourth moment) approximates snow rate S (which is close to the second PSD moment) better than the S(Z) power-law relation and is less sensitive to the PSD variability. This is well-illustrated in Figure 15 where the scatterplots of SWE measured by disdrometer during 16 Oklahoma snow events versus its estimate from the S(Z), S(K_DP, Z), and the S(K_DP, Z_DR) relations (discussed next) are displayed. The values of radar variables are simulated from the measured PSD assuming certain shape and orientation factors [99,101]. It is obvious that the impact of PSD variability is significantly reduced if polarimetric snow relations are utilized.

The biggest downside of the S(K_DP, Z) snow estimator is its sensitivity to the shape and orientation factors F_s and F_o. This emphasizes the need for realistic assumptions about snow particles’ shapes and orientations either from in situ microphysical measurements or from polarimetric radar retrievals. A possible alternative solution is to find radar relations that are not sensitive to the variability of shapes and orientations. Such relations based on the joint use of K_DP and Z_DR were first suggested by Ryzhkov et al. [87,88]). It can be shown that the ratio K_DP/(1 − Z_dr⁻¹) (where $Z_{dr}={10}^{0.1Z_{DR}(dB)}$) is practically immune to the changes in particle shapes and orientations. The Nguyen et al. [89] study confirms that the theoretical IWC(K_DP, Z_dr) relation from [87,88] is very close to their empirical relation derived from the side-looking airborne X-band radar measurements of K_DP and Z_DR and in situ observations onboard the Convair-580 aircraft during the High Altitude Ice Crystals—High Ice Water Content (HAIC—HIWC) campaign. This served as a motivation for derivation of S(K_DP, Z_dr) relation using the theoretical IWC(K_DP, Z_dr) relation as a starting point.

Because snow rate is proportional to the product of IWC and average terminal velocity vs. (or $D_m^{{\delta }_s}$) of snowflakes according to Equation (16):

$$S(K_{DP},Z_{dr})\sim IWC(K_{DP},Z_{dr})D_m^{{\delta }_s}$$

(20)

where D_m is the mean volume diameter of snowflakes and IWC and D_m are determined using formulas [10]:

$$IWC(K_{DP},Z_{dr})=4.0\times {10}^{-3}\frac{K_{DP}\lambda }{1-Z_{dr}^{-1}}$$

(21)

and:

$$D_m=-0.1+2.0{\left(\frac{Z(1-Z_{dr}^{-1})}{K_{DP}\lambda }\right)}^{1/2}$$

(22)

A formula for practical utilization has a form [99]:

$$S(K_{DP},Z_{dr})=10.8\times {10}^{-3}{\left(\frac{p_0}{p}\right)}^{1/2}\frac{K_{DP}\lambda }{1-Z_{dr}^{-1}}{D}_m^{0.15}$$

(23)

Although relation Equation (23) is immune to the variability of shapes and orientations, it is somewhat sensitive to the degree of riming quantified by the factor α_s in Equation (15). It is also quite vulnerable to possible biases of the Z_DR measurements, especially if Z_DR is low. Hence, relation Equation (23) is not recommended to be used if Z_DR < 0.2–0.3 dB.

Capozzi et al. [102] utilized an almost identical approach (using T-matrix scattering simulations instead of Rayleigh approximation) as Bukovčić et al. [101] to derive X-band polarimetric radar relations for snowfall rates using K_DP or a combination of K_DP and Z. They obtained comparable results to the wavelength-scaled S(K_DP, Z) relation from [101] using X-band radar measurements in southern Italy, confirming that combining K_DP with Z or Z_DR improves the accuracy of snowfall rate estimates compared with the Z-based relations, as demonstrated in the subsequent study of Bukovčić et al. [99].

Because K_DP in snow is usually quite low and noisy, especially at S band, its additional spatial averaging is required to obtain reliable and stable estimates of K_DP and S. Such averaging is implemented as part of the Quasi-Vertical Profiles (QVPs) technique introduced by Ryzhkov et al. [103] which prescribes azimuthal averaging of polarimetric radar variables measured at higher antenna elevations. Similar techniques such as range-defined QVPs (RD-QVP; [104]), enhanced vertical profiles (EVP; [105]), and columnar vertical profiles (CVPs; [106]) became recently available to address the noisiness of K_DP and Z_DR in snow.

3.3. Results of Observations

The performance of polarimetric algorithms for snow estimation is illustrated for a heavy snowfall event. This storm was a classical nor’easter named Stella, which occurred on 14–15 March 2017. The storm produced heavy snow throughout the inner mid-Atlantic and northeast US, causing hurricane-like wind gusts along the New England coast. Several daily snowfall records were broken (e.g., Binghamton, NY; Burlington, VT; Scranton, PA), with the total maximum of 42 in (~1 m) in West Winfield (NY). The storm ranked twelfth on the Regional Snowfall Index (RSI; [107]) and twenty-third on the Northeast Snowfall Impact Scale (NESIS; [108]), disrupting the daily routine of more than 60 million people.

Figure 16 shows an example of composite PPI of Z, K_DP, Z_DR, and cross-correlation coefficient ρ_hv measured by the KBGM WSR-88D radar at El = 0.5° during the storm. Specific differential phase K_DP exhibits inherent noisiness which, nevertheless, does not prevent to reveal a definite spatial structure with K_DP enhancement in the area that does not coincide with the region of increased Z. Aggressive azimuthal averaging implemented in the RD-QVP methodology helps to significantly reduce the noisiness of K_DP (see Figure 17).

The range-defined quasi-vertical profiles (RD-QVPs) of Z, K_DP, and Z_DR in a height-vs-time format obtained from the KBGM WSR-88D radar and vertical profiles of radar-retrieved snow water equivalent (SWE) accumulations are presented in Figure 17. The RD-QVP radar data were azimuthally averaged in a 15-km-radius circular area centered on the radar. Dashed lines indicate temperature isotherms obtained from the RAP model. Moderate to heavy snowfall rates persisted from 0700 to 1500 UTC, with a few intense snow bands later in the event. During these periods, the RD-QVP of Z exhibited elevated values greater than 30 dBZ (Figure 17a). Enhanced values of K_DP and especially Z_DR (Figure 17b,c) are noticeable in the dendritic growth layer (DGL) centered at about −15 °C. Such enhancements occasionally extended all the way to the ground (e.g., from 0700 to 1330 UTC). This marked a period of moderate concentration of smaller and denser dendrite-type crystals with a lower degree of aggregation at the surface. These can occasionally be mixed with the irregular-shape crystals, indicated by the enhanced values of K_DP aloft (~1030 to 1130 UTC). This K_DP increase, caused by high concentrations of irregular-shape crystals, comes as a precursor to the biggest low-level increase in Z, which shows a promising potential for snow nowcasting. In the later storm stages (1900 UTC onward) temperatures and cloud-top heights decrease, and the pristine dendrite-type crystals which do not aggregate, almost perfectly align with the −15 °C isotherm. These are also seen as Z_DR streaks reaching the ground at 2010 and 2210 UTC. Previous observational studies [109,110] saw similar signatures in their data.

Three algorithms were used to estimate total SWE accumulation profiles from the RD-QVP data. These include the S(Z) = 0.088Z^0.5 relation which is commonly used for the northeast US and polarimetric Relations (19) and (23) assuming that the particle aspect ratio is 0.6 and the canting angle spread σ is 20°. Actual SWE total accumulation at the surface of about 50 mm was measured by the Automatic Surface Observations Station (ASOS) in the close proximity to the KBGM radar. Both polarimetric estimates (red and green lines in Figure 17d) show more realistic representation of the vertical profile of accumulated SWE. The S(Z) relation produces apparent underestimation of snowfall (blue line in Figure 17d). Polarimetric relations are in accord with the notion that the bulk of precipitation (about 80 to 90%) is formed in the DGL ([111]). On the contrary, S(Z) produces about 40% of its total accumulation within the DGL (from ~2.0 to 4.5 km AGL) and moderately underestimates surface snow measurements (32 vs. 50 mm). Polarimetric relations at the lowest usable K_DP height level (about 450 m AGL) yield results which are very close to the surface ASOS measurements.

Instantaneous snowfall rates and accumulations show more detailed pictures about the temporal evolution of these estimates (Figure 18). The standard S(Z) relation moderately underestimates snowfall rates throughout the event, with peaks not greater than 3 mm h⁻¹ (gauge maximum is 5.5 mm/h). Polarimetric relations are in good agreement with the gauge, especially S(K_DP, Z), which more accurately reproduce the peaks in rates. Strong winds (>7 m/s) were reported after 2000 UTC, where the minor degradation is seen in polarimetric relations estimations. In the strong wind situations, particles have more random motions and appear quasi-spherical, negatively impacting K_DP (and Z_DR). Polarimetric-based accumulations are in good agreement with the gauge estimate whereas the S(Z) relation underestimates the SWE amounts by about 40%. Additional examples of polarimetric measurements of snow can be found in [99].

3.4. Discussion

Although the results of polarimetric estimation of snow in a limited number of cases are encouraging, they reveal obvious challenges as well. There is little doubt that the use of multiparameter polarimetric relations significantly reduces the impact of the PSD variability in snow. However, the S(K_DP, Z) estimates are sensitive to the variability of particle shapes and orientations which are not known apriori. This variability does not affect the S(K_DP, Z_DR) estimates but these are much more affected by the uncertainty in snow density or degree of riming. A common problem for both polarimetric relations is that they cannot be directly applied for heavily aggregated dry snow where K_DP and Z_DR are vanishingly small and noisy. If K_DP drops below 0.01–0.03 deg/km and Z_DR are below 0.2–0.3 dB, the polarimetric estimates of snow become very unstable. This may limit the efficiency of polarimetric radar estimates in “warm” snow with temperatures just slightly below 0 °C. One possibility to address this issue is to resort to the polarimetric estimates of snow rate at higher altitudes and lower temperatures where K_DP (or Z_DR) are sufficiently high and project these estimates down to the surface assuming either conservation of the snow flux or predicting its downward evolution using simple cloud models as discussed in [112]. An alternative option is to switch to the S(Z) relation optimized for a certain snow type in areas of unreliable K_DP and Z_DR [83]. Note that the polarimetric snow relations are expected to perform better at shorter radar wavelengths (e.g., at C and X bands) for which K_DP is higher and can be measured more reliably compared with S band.

Summarizing, we can conclude that, as opposed to the polarimetric methods for rainfall estimation which are quite mature and operationally usable, the algorithms for polarimetric snow measurements are still in their infancy and require substantial refinement and testing on a larger number of snow events.

4. Operational Implementation—MRMS QPE

Since the polarimetric upgrade of the US WSR-88D network, several QPE methodologies based on polarimetric radar variables have been developed for operations in the US National Weather Service (NWS). On the single radar applications, a hydrometeor classification algorithm (HCA, [113]) and an R(Z, Z_DR), R(K_DP) and R(Z) synthetic QPE [8] were developed and implemented in the operational WSR-88D Open Radar Product Generation (ORPG) system. The single radar synthetic QPE was later upgraded to an R(A) [42], R(Z, Z_DR), R(K_DP) and R(Z) synthetic QPE. On the multi-radar mosaic applications, a dual-pol (DP) radar quality control (dpQC) [114,115] to remove non-hydrometeor echoes and an R(A), R(K_DP) and R(Z) synthetic QPE [26,46,47] were implemented in the Multi-Radar Multi-Sensor (MRMS, [9]) system. This section provides an overview of the DP radar synthetic QPE in the MRMS system.

Figure 19 shows an overview flowchart of the MRMS DP synthetic QPE process. The input data include DP radar moments and environmental data such as the 3D temperature field, freezing level height, and surface wet bulb temperature. The dpQC includes a simple correlation coefficient (ρ_hv) filter to separate hydrometeor (high ρ_hv) and non-hydrometeor (low ρ_hv) areas and a set of heuristic rules to handle exceptions to the simple ρ_hv filter. Such exceptions include areas of hail, non-uniform beam filling (NBF) and melting layer (ML) that have low ρ_hv values. Another exception is random clutter and biological pixels with high ρ_hv values. The dpQC uses 3D reflectivity structure and environmental data to protect hail, NBF, and ML areas from being removed by the simple ρ_hv filter, and it uses spatial filters and vertical and horizontal consistency checks to remove random non-precipitation pixels that exhibit high ρ_hv values. The dpQC removes more than 99% of non-hydrometeor echoes with very high computational efficiency [114].

After the dpQC, the differential phase Φ_DP field is further processed for additional quality assurance. The Φ_DP processing include four steps. The first step includes a speckle filter and a ρ_hv screening. The speckle filter checks for each Φ_DP pixel and counts a number of non-missing Φ_DP pixels in a 4.5$°$ $\times$ 2.25 km box centered at the given pixel. If the number is less than 50% of the box total, then the given Φ_DP pixel at the center is removed. A ρ_hv > 0.8 filter is applied to remove noisy and unreliable Φ_DP data potentially associated with hail, severe NBF, or other residual contaminations after the dpQC. The second step unfolds Φ_DP field based on the gate-to-gate Φ_DP change and adjusts the unfolded Φ_DP field to assure a monotonic increase trend with range. The third step applies a 6.25 km radial running mean through all non-missing Φ_DP pixels to further reduce random errors and fluctuations. Finally, a linear interpolation along the radial is applied to fill in any gaps in the Φ_DP field.

The specific differential phase K_DP is calculated via a local linear fitting to Φ_DP along the radial direction in areas where hail might be present. Specific attenuation A along the propagation path (r₁, r₂) in rain is calculated following Equations (2)–(6). The distance r₁ corresponds to the nearest gate to the radar with precipitation and r₂, to the last precipitation gate or the gate just below the bottom of the apparent melting layer (ML), whichever is closer to the radar. The ML is delineated using a scheme similar to [116] but its bottom and top were further expanded based on ρ_HV and Z to encompass all pixels potentially impacted by the ML. Pixels in [r₁, r₂] with potential hail contamination are excluded from the PIA calculation.

In the scheme currently implemented on the WSR-88D network, the factor α in Equation (6) is estimated from the 0.5$°$ tilt data using the “Z_DR-slope” K:

$$\alpha =-0.75K+0.049$$

(24)

where K is a linear fit to the median Z_DR values of each 2 dBZ reflectivity bin within the 20–50 dBZ range (Figure 20). Equation (24) was obtained from disdrometer data and supported by an analysis of ~10 events with a wide spread of moderate to heavy precipitation events from different geographical regimes. Figure 20 illustrates the process of α calculation. For the mesoscale convective system (MCS) near the KHTX WSR-88D radar (Figure 20a–c) on 11 December 2021, there was an apparent increase in Z_DR with increasing Z, indicating the presence of large rain drops in the convective system. The Z_DR-slope was relatively large (K = 0.048 dB/dBZ) and α was relatively low (0.0126 dB/deg) and close to those for convective rain [42,46,47]. For the Hurricane Ida case near the KMOB WSR-88D radar on 30 August 2021 (Figure 20d–f), the Z_DR-slope was relatively flat (K = 0.017) and α (0.0361 dB/deg) was higher and closer to a tropical-type rain value.

Equation (24) is used to estimate α only if there are sufficient Z–Z_DR data pairs in each 2-dBZ bin between 20 and 50 dBZ. If not enough Z–Z_DR data pairs are found in the 20–50 dBZ interval, the number of Z–Z_DR pairs in each 2-dBZ bin between 10 and 30 dBZ are checked. If significant data samples were found in the 10–30 dBZ interval, then the precipitation is considered pure stratiform and a default stratiform α (0.035 dB/deg) is applied. Otherwise, a new linear fit may be applied to the median Z_DR values within the 10–40 dBZ range if sufficient Z–Z_DR pairs were found in that range. If insufficient data pairs were found in the 10–40 dBZ reflectivity range, the precipitation is considered sporadic. A default convective (0.015 dB/deg) or stratiform (0.035 dB/deg) α is applied depending on the reflectivity intensities.

Once α is obtained and A calculated for each radar gate, precipitation rates are estimated from A on the 0.5° tilt in areas where reflectivity is below 45 dBZ (an adaptable parameter) using Equation (7). The 45 dBZ constraint is used to prevent R(A) from potential hail contamination. Above 50 dBZ (an adaptable parameter), two R(K_DP) relationships are applied depending on the correlation coefficient, ρ_hv, field:

$$R(K_{DP})=29.0|K_{DP}{|}^{0.77} \mathrm{if} {\rho }_{hv} < 0.97$$

(25)

$$R(K_{DP})=44.0|K_{DP}{|}^{0.82} \mathrm{if} {\rho }_{hv} > 0.97$$

(26)

Equation (25) was derived using disdrometer measurements in the proximity of selected hailstorms in central Oklahoma, and Equation (26) provided optimal performance in heavy rain with large drops [6], also in central Oklahoma. Between 45 and 50 dBZ, a linear combination of R(A) and R(K_DP) is applied to create a smooth transition between the two rates. The R(A) and R(K_DP) are calculated from the 0.5° tilt currently with two additional constraints applied: (1) the beam blockage < 90% and (2) ΔΦ_DP $\ge$ 3°. If the blockage $\ge$ 90% or ΔΦ_DP < 3° at a given radial, the rate values for all the pixels in that radial are set to missing and R(Z), based rates are applied in the synthetic process at a later time (Figure 19).

Precipitation rate fields estimated via R(A) and R(K_DP) from individual radars are mosaicked through a physically based scheme shown in [117]. Since R(A) is only valid below the melting layer, the mosaic may have gaps in areas far away from the radars. Figure 21a shows such a gap between KHTX and KMRX radars. To fill in these gaps, a multiple R(Z), based QPE (“Q3RAD”, Figure 21b) described in [9] was applied. A flowchart of the Q3RAD QPE is shown in the right half of Figure 19. After the QC, a 2D hybrid scan reflectivity (HSR) field is derived from the lowest radar bins with less than 50% blockages, and the reflectivity bins with blockages greater than 10% are adjusted to compensate for the power loss. The apparent melting layer is delineated for all the tilts based on a scheme similar to [118] and a vertical profile of reflectivity (VPR) correction was applied to the HSR field within and above the melting layer [61]. The HSR fields from individual radars are then mosaicked based on the methodology developed in [117]. A precipitation type field of seven categories (warm stratiform rain, cool stratiform rain, convective rain, tropical-stratiform rain mix, tropical-convective rain mix, hail, and snow) is generated using an automated precipitation classification described in [9]. Finally, multiple R(Z) relationships are applied to the HSR field according to the precipitation types to obtain the Q3RAD rate. The main Q3RAD R(Z) relationships applied in the warm season are as follows:

Stratiform rain: R_stra = max (0.0365 Z^0.625, 0.1155 Z^0.5);

(27)

Convective rain: R_conv = 0.017 Z^0.714;

(28)

Tropical rain: R_trop = 0.010 β Z^0.833

(29)

where β is a rate multiplier ranging from 1.0 to 1.5 depending on the month of the year and the proximity to frequent hurricane zones [9,118]. For the coastal areas of the Gulf of Mexico and Atlantics and southeastern US, β can be as high as 1.45 for September and 1.5 for October. The convective R(Z) was capped at a value that ranges from 53.8 mm/h to 150 mm/h depending on the climatological region and real-time hail size estimation. The cap is lower in areas of continental rain regime with severe hail and higher in areas of maritime warm rain and hurricanes.

Figure 21c,d show the example of the HSR mosaic and precipitation type fields corresponding to the Q3RAD rate in Figure 21b. The final dual-pol radar synthetic QPE (“Q3DP”) is shown in Figure 21e,f which shows the contributing sources of the rates to the Q3DP from different relationships. The dark blue areas in Figure 21f indicate where Q3RAD fills in when the radar beam intersects or overshoots the melting layer and light blue when ΔΦ_DP < 3°. The yellow areas indicate a transition zone (default width: 50 km) near the outer boundary of R(A) where a weighted mean of the R(A) and R(Z) rates is applied to prevent discontinuities between the two. The red areas indicate where R(K_DP) is applied and orange the linear transition between R(A) and R(K_DP). Finally, an evaporation correction is applied to the mosaicked precipitation rate field based on the methodology provided by Martinaitis et al. [119]. The correction factor is a function of the rate, the height of the radar estimate, and atmospheric properties between the radar data height and the ground. The correction helps reducing false light precipitation in the radar QPE due to virga. Figure 21g shows the same Q3DP rate field as in Figure 21e but with the evaporation correction.

Q3DP was evaluated across CONUS during the period of 15 September–31 October 2018. The Community Collaborative Rain, Hail and Snow Network, www.cocorahs.org/, (accessed on 1 February 2022) (CoCoRaHS) rain gauges were used to compare the performances of the Q3DP and Q3RAD 24-h QPEs. The statistic scores for both QPE techniques for different rainfall categories are summarized in

Table 2. Mean bias ratio (MBR), correlation coefficient (CC), mean absolute error (MAE) and fractional MAE of the Q3RAD and Q3DP 24 h QPEs for predefined 24 h rainfall categories. The bold numbers indicate the better scores among the two products.

Category		VL	L	M	H	VH
QPE-G pairs		53,478	13,796	9834	3404	710
G-mean (in)		0.17	0.71	1.37	2.66	5.45
Q3RAD	MBR	1.19	0.91	0.82	0.80	0.81
	CC	0.69	0.39	0.43	0.52	0.57
	MAE (in)	0.09	0.21	0.38	0.78	1.73
	fMAE (%)	52.94	29.58	27.74	29.32	31.74
Q3DP	MBR	1.04	0.91	0.88	0.92	0.95
	CC	0.68	0.39	0.49	0.59	0.71
	MAE (in)	0.08	0.22	0.36	0.59	1.14
	fMAE (%)	47.06	30.99	26.28	22.18	20.92

. The categories are subjectively defined as very light (VL: G < 0.5 in.), light (L: 0.5 $\le$ G $<$ 1 in.), moderate (M: 1 $\le$ G < 2 in.), heavy (H: 2 $\le$ G < 4 in.), and very heavy (VH: G $\ge$ 4 in.).

Q3DP showed a reduction in systematic bias and random errors over Q3RAD, especially for the heavy amounts higher than 2 in., Q3RAD had a 19% overestimation bias in the VL rain and 18%, 20%, and 19% underestimation biases in the M, H, and VH rain, respectively. Q3DP reduced the biases by 15%, 6%, 12%, and 16% for VL, M, H, and VH categories, respectively, and remained the same for the L category. A similar trend was found in the MAE score, where Q3DP reduced the errors for VL, M, H, and VH categories but slightly increased the error in the L category. For the correlation coefficient, Q3DP was about the same as Q3RAD in the VL and L and was higher in the M–VH categories. Q3DP showed a significantly reduced uncertainty with respect to Q3RAD for rainfall amounts. The reduced uncertainty in Q3DP was also reflected in its lower MAEs/fMAEs for the H and VH categories. Q3DP reduced the fMAE by 6%, 7%, and 11% for the VL, H, and VH categories, respectively, while the fMAEs of Q3DP and Q3RAD for the L and M categories remained similar.

5. Conclusions

A primary outcome of the radar rainfall estimation research during last two decades is the realization that the polarimetric methodologies based on the use of specific attenuation A and specific differential phase K_DP have undeniable benefits compared with the techniques that utilize either a single radar reflectivity factor Z or a combination of Z with differential reflectivity Z_DR. This is attributed to lower sensitivity of the R(A) and R(K_DP) relations to the DSD variability and immunity of the estimates of A and K_DP to the radar calibration errors, attenuation in rain, partial beam blockage, and impact of wet radome. These benefits are well-proven for all three S-, C-, and X-microwave frequency bands at which polarimetric weather surveillance radars operate. It is demonstrated that the composite rainfall algorithms utilizing R(A) relations for lower rain rates and R(K_DP) relations for heavier rain show the best performance. The A-based methods are particularly advantageous at S band because the S-band R(A) relation is almost linear and can be utilized in a wider range of rain intensities compared with C and X bands.

There is still room for additional optimization of the R(A) and R(K_DP) algorithms to further mitigate their residual sensitivity to the DSD variability. At S band, this must be performed by optimizing the parameter α = A/K_DP in the estimation of A and the multiplier in the power-law R(K_DP) relation.

Notable progress was achieved in mitigation of the impact of the “bright band” contamination on the surface rainfall estimates at longer distances from the radar using polarimetric VPR correction (PVPR) although these techniques require more testing and refinement. Taking into the account the vertical profiles of reflectivity below the melting layer and above, is a subject of ongoing research.

The algorithms for polarimetric measurements of snow were recently introduced and show good promise as was demonstrated at S band. They utilize combinations of Z and K_DP or K_DP and Z_DR and may produce even better results at C and X bands due to higher value of K_DP at shorter radar wavelengths. A common problem for this methodology is a reduced value of K_DP in dry aggregated snow near the surface and the inherent uncertainty related to the K_DP dependency on the shapes and orientations of snowflakes.

The advanced polarimetric radar QPE methods utilizing the R(A), R(K_DP), and R(Z) relations were implemented on the operational network of the WSR-88D radars in the US using the Multi-Radar Multi-Sensor (MRMS) platform that integrates rainfall estimates from all operational radars and provides comparison of radar rainfall products with dense networks of rain gauges.

With all the progress achieved during last two decades, a number of challenges and opportunities remains. These include the need to improve the rainfall estimation performance at longer distances from the radar, further development of the algorithms for polarimetric quantification of snow, and more aggressive utilization of the methodologies involving artificial intelligence/machine learning.