What Polarimetric Weather Radars O ﬀ er to Cloud Modelers: Forward Radar Operators and Microphysical / Thermodynamic Retrievals

: The utilization of polarimetric weather radars for optimizing cloud models is a next frontier of research. It is widely understood that inadequacies in microphysical parameterization schemes in numerical weather prediction (NWP) models is a primary cause of forecast uncertainties. Due to its ability to distinguish between hydrometeors with di ﬀ erent microphysical habits and to identify “polarimetric ﬁngerprints” of various microphysical processes, polarimetric radar emerges as a primary source of needed information. There are two approaches to leverage this information for NWP models: (1) radar microphysical and thermodynamic retrievals and (2) forward radar operators for converting the model outputs into the ﬁelds of polarimetric radar variables. In this paper, we will provide an overview of both. Polarimetric measurements can be combined with cloud models of varying complexity, including ones with bulk and spectral bin microphysics, as well as simpliﬁed Lagrangian models focused on a particular microphysical process. Combining polarimetric measurements with cloud modeling can reveal the impact of important microphysical agents such as aerosols or supercooled cloud water invisible to the radar on cloud and precipitation formation. Some pertinent results obtained from models with spectral bin microphysics, including the Hebrew University cloud model (HUCM) and 1D models of melting hail and snow coupled with the NSSL forward radar operator, are illustrated in the paper.


Introduction
Doppler polarimetric radars have become the standard for operational weather radars around the world. Compared to single-polarization radars, dual-polarization radars substantially improve data quality, precipitation estimation, hydrometeor classification, and severe weather warnings. Comprehensive overviews of radar polarimetry and its utilization for weather observations are described in the monographs of Bringi and Chandrasekar [1], Zhang [2] and Ryzhkov and Zrnic [3].
The utilization of polarimetric weather radars to optimize numerical weather prediction (NWP) models is a new frontier of research. It is widely understood that inadequate microphysical parameterization schemes in NWP models are a primary source of forecast uncertainties (e.g., Morrison et al. [4] and Fan et al. [5]). Observational data provided by ground-based weather surveillance and cloud radars, lidars, radiometers, and airborne or spaceborne remote sensors are crucial for constraining model microphysical schemes. Due to its ability to distinguish between hydrometeors with different microphysical habits and to identify the "polarimetric fingerprints" of various microphysical processes, polarimetric radar has emerged as an important source of needed information.
The purpose of this paper is to provide a review of some of the important information that polarimetric weather radars can provide that may be of use to the cloud modeling community. We hope that through increased awareness of polarimetric radar forward operators, microphysical and thermodynamic retrievals and other polarimetric products and data can be integrated into cloud models through improved data assimilation and leveraged to assess and validate numerical models and/or microphysics schemes.
In this review, we will focus on the applications of longer-wavelength polarimetric radars operating at centimeter radar wavelengths at S, C, and X microwave frequency bands which commonly comprise operational weather radar networks around the world. These radars are used for surveillance of large areas with extended radar coverage as opposed to cloud radars operating at millimeter wavelengths at Ka-or W-band frequencies. Weather surveillance radars do not have sufficient sensitivity to detect echoes from atmospheric aerosols or cloud particles (liquid or solid) with sizes less than 50-100 microns but have great advantage for remote sensing of deep convective storms containing large hydrometeors producing prohibitively strong attenuation at cloud radar frequencies.
Although centimeter-wavelength radars do not see small cloud droplets, they are often capable of detecting non-precipitating clouds containing ice or clear-air turbulence associated with Bragg scattering. Melnikov et al. [6] claim that longer-wavelength operational radars such as WSR-88D in the US have the potential to be used as "cloud radars".
In addition to the three radar variables commonly measured by single-polarization Doppler radars (i.e., radar reflectivity Z, mean Doppler velocity, and Doppler spectrum width), polarimetric radars can provide estimates of differential reflectivity Z DR , specific differential phase K DP , linear depolarization ratio (LDR), cross-correlation coefficient ρ hv , specific attenuation A, and specific differential attenuation A DP . For definitions and the physical meanings of these polarimetric variables, a reader can consult basic monographs [1][2][3]. Using such multi-parameter measurements, one can discriminate between hydrometeors with different shapes, orientations, and phase compositions; reduce uncertainties in the retrievals of size distributions of atmospheric particles that are generally characterized by at least three parameters; and identify polarimetric signatures attributed to different microphysical processes. These include the evaporation, coalescence, and breakup of raindrops; depositional growth, sublimation, and aggregation of ice/snow; and the processes involving mixed-phase hydrometeors such as melting, freezing, refreezing, and riming. Although the processes of cloud and ice nucleation and the initial condensation growth of small particles are invisible to the radar, their impacts can be indirectly identified through their influence on the polarimetric characteristics of clouds at later stages of cloud development when larger particles are formed. It is noteworthy that multiple polarimetric variables are measured in the same radar resolution volume as opposed to multi-frequency measurements from multiple radars, for which the matching of the sampling volumes is always a challenge. This is a reason why multi-frequency retrievals of microphysical and thermodynamic characteristics of clouds with ground-based radars are primarily confined in the research domain whereas polarimetric radar retrievals can be easily implemented on existing operational radar networks. Of course, the use of polarimetric radars in conjunction with cloud models is still a subject of exploratory research, but we believe that operational polarimetric radars are a primary resource for optimizing NWP models on a national scale in the foreseeable future.
There are two approaches for incorporating information from dual-polarization radars into NWP models: (1) microphysical and thermodynamic retrievals from radar (which is a way of converting radar information into quantities that the model prognoses) and (2) forward operators that convert the model output into polarimetric radar fields. In this paper, we will provide an overview of both. The idea of using forward operators for assimilating polarimetric radar information into NWP models is also explored in the recent article of Zhang et al. [7]. The Zhang et al. paper focuses on a unified statistical Atmosphere 2020, 11, 362 3 of 34 approach for model performance optimization based on synthesizing observation-based retrievals and model analyses. Such an approach requires knowledge of the background error covariance and observation error covariance stipulated by the theory of optimal estimation. In our article, we focus exclusively on the direct, observation-based retrievals (microphysical and thermodynamic) using polarimetric radar data.
The paper is organized as follows. We begin with a general overview of the electromagnetic model of scattering and forward operators in Section 2. The microphysical and thermodynamic radar retrievals are described in Sections 3 and 4, respectively, followed by conclusions in Section 5. Formulas for computing different polarimetric radar variables from the output of cloud models are in the Appendix A.

Scattering by Individual Hydrometeor
The polarimetric properties of most meteorological scatterers are generally well-quantified if they are modeled as spheroids with oblate or prolate shapes. The characteristics of electromagnetic scattering by spheroidal particles depend on their size, shape, orientation, and physical composition and are described in terms of the scattering matrix S, which has the following form for a single spheroidal hydrometeor with an axis of symmetry a and transverse axis b (Holt [8], Ryzhkov [9], and Ryzhkov and Zrnic [3]): In Equation (1), s a (π) is the backscattering amplitude if the electric field vector of the incident electromagnetic wave is parallel to the symmetry axis of the hydrometeor, and s b (π) stands for the backscattering amplitude if the electric field vector is perpendicular to the symmetry axis. Angles α and Ψ characterize the orientation of the spheroid and are illustrated in Figure 1. Strictly speaking, expression (1) for S is valid for particles such as raindrops at radar frequencies up to 35 GHz provided that Ψ > 80 • (Holt and Shepherd [10]), but it can be universally applied at lower frequencies for the majority of practical applications. For more rigorous modeling of scattering by spheroids, the reader is referred to Vivekanandan et al. [11].
synthesizing observation-based retrievals and model analyses. Such an approach requires knowledge of the background error covariance and observation error covariance stipulated by the theory of optimal estimation. In our article, we focus exclusively on the direct, observation-based retrievals (microphysical and thermodynamic) using polarimetric radar data. The paper is organized as follows. We begin with a general overview of the electromagnetic model of scattering and forward operators in Section 2. The microphysical and thermodynamic radar retrievals are described in Sections 3 and 4, respectively, followed by conclusions in Section 5. Formulas for computing different polarimetric radar variables from the output of cloud models are in the Appendix.

Scattering by Individual Hydrometeor
The polarimetric properties of most meteorological scatterers are generally well-quantified if they are modeled as spheroids with oblate or prolate shapes. The characteristics of electromagnetic scattering by spheroidal particles depend on their size, shape, orientation, and physical composition and are described in terms of the scattering matrix S, which has the following form for a single spheroidal hydrometeor with an axis of symmetry a and transverse axis b (Holt [8], Ryzhkov [9], and Ryzhkov and Zrnic [3]): In Equation (1), sa (π) is the backscattering amplitude if the electric field vector of the incident electromagnetic wave is parallel to the symmetry axis of the hydrometeor, and sb (π) stands for the backscattering amplitude if the electric field vector is perpendicular to the symmetry axis. Angles α and Ψ characterize the orientation of the spheroid and are illustrated in Figure 1. Strictly speaking, expression (1) for S is valid for particles such as raindrops at radar frequencies up to 35 GHz provided that Ψ > 80° (Holt and Shepherd [10]), but it can be universally applied at lower frequencies for the majority of practical applications. For more rigorous modeling of scattering by spheroids, the reader is referred to Vivekanandan et al. [11].
If the electromagnetic waves transmitted and received by the radar propagate in a medium filled with nonspherical hydrometeors, the complex voltages of the radar return at orthogonal polarizations depend not only on the elements of the scattering matrix S of the hydrometeor but also on the polarization properties of the propagation medium characterized by the transmission matrix T. The resulting voltage vector V r of the received signal is given by (Ryzhkov and Zrnic [3]). r t V TSTV = (2) Figure 1. Scattering geometry. Direction N denotes the orientation of the symmetry axis of the particle, and k represents the direction of wave propagation (which is perpendicular to the polarization plane and lies in the x, z plane). The polarization plane is depicted by the grey ellipse. The x axis is true vertical, and y and z are the orthogonal horizontal directions. The canting angle α is the angle between the projections of vector N and true vertical x on the polarization plane, ψ is the angle between N and k, and the antenna elevation is β. If the electromagnetic waves transmitted and received by the radar propagate in a medium filled with nonspherical hydrometeors, the complex voltages of the radar return at orthogonal polarizations depend not only on the elements of the scattering matrix S of the hydrometeor but also on the polarization properties of the propagation medium characterized by the transmission matrix T. The resulting voltage vector V r of the received signal is given by (Ryzhkov and Zrnic [3]).
where V t is the voltage vector characterizing a transmitted wave. Whereas the scattering matrix S is characterized by the "backscattering" amplitudes s a (π) and s b (π) for the scatterer in the radar resolution volume, the transmission matrix T is determined by the "forward scattering" amplitudes s a (0) and s b (0) of the scatterers located along the propagation path (e.g., Bringi and Chandrasekar [1], Zhang [2], and Ryzhkov and Zrnic [3]). If the hydrometeors are relatively small compared to the radar wavelength λ and are modeled as oblate or prolate spheroids, then simple analytical formulas for the forward scattering and backscattering amplitudes s can be obtained in the Rayleigh approximation (Van de Hulst [12]): where D = (ab 2 ) 1/3 is the equivolume diameter of the particle, ε is its dielectric constant, and L a and L b are the shape parameters defined as for oblate spheroids (a < b) and for prolate spheroids (a > b). The dielectric constant ε is a function of the radar wavelength; phase composition (e.g., water, ice, air, or some combination thereof); temperature; and density of the scattering material. For mixed-phase particles, ε depends on the relative fractions of water, solid ice, and air, with various mixing formulas used for computing ε of such hydrometeors, the most popular of which were suggested by Maxwell Garnett [13]. Uncertainties arise when using the Maxwell Garnett mixing formulas, as they do not yield a unique value of ε for a given mass water fraction but, rather, values dependent on the distribution of water within the mixed-phase particle. A model of a two-layer spheroid is a reasonable choice to treat water-coated graupel and hail or ice-coated supercooled freezing raindrops. A closed-form analytical solution also exists for two-layer spheroids within the Rayleigh scattering regime (Bohren and Huffman [14]).
The Rayleigh approximation (3) is valid if the ratio (or "resonance parameter") RP = D ε 1/2 /λ is less than 0.3-0.4 (Ryzhkov et al. [15]). It provides a very simple and economic way to compute scattering amplitudes for a wide range of practical problems. Moreover, the closed-form solution is instrumental for straightforward physical interpretation of the results. If the scatterers are sufficiently large with respect to the radar wavelength, then numerical methods such as the "T-matrix" are used to compute the scattering amplitudes (e.g., Waterman [16], Barber and Yeh [17], Mishchenko [18], and Bringi and Chandrasekar [1]). For large, uniformly filled spheroids, the Mishchenko [18] T-matrix code is a standard solution. A two-layer version of the T-matrix code developed by Bringi and Seliga [19] can be used for spheroidal ice-coated freezing drops and water-coated graupel particles and hailstones (e.g., Aydin and Zhao [20], Depue et al. [21], and Ryzhkov et al. [15,22]).
Certain types of ice hydrometeors, such as graupel, hail, and dendritic ice crystals, may have shapes that can differ significantly from spheroids. A relatively simple Rayleigh-Gans approximation can be utilized for computations of the backscattering cross-sections for snowflakes having arbitrary shapes and very low density (Lu et al. [23]). While the Rayleigh-Gans methodology is capable of capturing the effects of resonance scattering in snow at millimeter radar wavelengths, it cannot be used for quantification of its polarimetric properties which are significantly "muffled" for low-density snow.
More advanced techniques exist for calculating the scattering from nonspheroidal particles. Aydin and Seliga [24] used Waterman's [25] extended boundary condition technique to compute cross-sections at orthogonal polarizations for conical graupel. Lu et al. [26] created a database of single-scattering characteristics of select ice particles with complex shapes (such as plates, columns, branched planar crystals, and conical graupel) using numerical methods, including the Generalized Multi-particle Mie (GMM) method (Xu [27]) and the discrete dipole approximation (DDA) method (Yurkin and Hoekstra [28]). GMM models complex particles as clusters of spheres, whereas DDA models them as cluster of dipoles. Mirkovic [29] computed the elements of the scattering matrix for hail of irregular shapes using the Method of Moments (MoM) (Harrington [30]) as part of the WIre-PLate-Dielectric (WIPL-D) solver package (WIPL-D [31]), which can be applied for scatterers with arbitrary shapes. Although the spheroidal model of hydrometeors is widely accepted in the majority of scattering computations and interpretation of observations, some recent simulation studies (e.g., Botta et al. [32], Tyynela et al. [33], Leinonen et al. [34], Tyynela and Chandrasekar [35], Ori et al. [36], Leinonen and Moisseev [37], and Schrom and Kumjian [38]) indicate that this model has to be used with caution for branched planar ice crystals and aggregated snowflakes, particularly for radar frequencies higher than 10 GHz. Botta et al. [32] utilized the GMM technique for computing backscattering cross-sections σ b of dry and melting snow aggregates and found out that, for sizes larger than 3 mm, the GMM yield values of σ b several dB higher than the traditional T-matrix method that assumes uniformly filled spheroidal particles. The difference is well over 7 dB at the Ka band (35.6 GHz). Tyynela et al. [33] and Ori et al. [36] reported even larger discrepancies at Ka and W bands using the DDA computations. Schrom and Kumjian [38] also demonstrated significant differences in the scattering amplitudes computed by the DDA method for branched planar crystals and their lower-density spheroidal equivalents with the same mass, aspect ratio, and maximal dimensions.
All these findings indicate three primary limitations of the T-matrix method and uniform spheroidal model for computations of the scattering characteristics of ice and snow particles. Firstly, it may not be applicable at radar frequencies higher than 10 GHz. Secondly, T-matrix computations yield unacceptable errors for very anisotropic particles with aspect ratios of less than 0.2 (for oblates) and greater than 5 (for prolates). Finally, the intrinsic uncertainty in the determination of the effective dielectric constant of a uniformly filled spheroid composed of ice, water, and air can have a significant impact on the scattering parameters of hydrometeors regardless of the radar wavelength (Fabry and Szyrmer [39]).
The GMM, DDA, and MoM techniques do not have such limitations. They can be applied at any radar wavelength for scatterers of arbitrary shape, and there is no need to use mixing formulas for the effective dielectric constant of a mixture. It is important that these models take into account electromagnetic interactions between elementary blocks (dipoles or spheres) comprising the modeled hydrometeors (e.g., snow particles). These methods substantially reduce the constructive interference that leads to strong resonance effects typical of spheroidal models with rigid boundaries. However, such benefits come at high computational cost. The pros and cons for the choice of a particular scattering technique (GMM, DDA, or MoM) are not well-established, although Lu et al. [26] consider DDA more accurate than GMM for computing the scattering properties of pristine ice crystals, and Chobanyan et al. [40] claim that the MoM algorithm works much faster than DDA for relatively simple hydrometeor shapes. While more accurate, these methods are computationally expensive and cannot be used for online computations. Therefore, in order to widely use these approaches, precalculated lookup tables or databases of the scattering characteristics have to be created and made available to the research community. Examples of such databases and relevant references can be found in Lu et al. [26].

Polarimetric Radar Forward Operators
Polarimetric radar forward operators (PRFOs) are used to convert the output from NWP models into commonly observed polarimetric variables: radar reflectivity factor at horizontal polarization Z H , differential reflectivity Z DR , specific differential phase K DP , linear depolarization ratio LDR, cross-correlation coefficient ρ hv , specific attenuation A h at horizontal polarization, and specific differential attenuation A DP . All of these radar variables can be expressed via second moments of the scattering matrix S of an ensemble of hydrometeors in the radar resolution volume (Bringi and Chandrasekar [1], Zhang [2], and Ryzhkov and Zrnic [3]).
The scattering matrix of an ensemble of hydrometeors depends on the distributions of their sizes, shapes, and orientations. Cloud models with spectral bin microphysics provide size distributions (SDs) of hydrometeors explicitly over discrete size bins without assuming a certain functional form for each SD. Cloud models with bulk microphysics solve prognostic equations for one or more quantities proportional to one or more bulk moments of the SD, such as mass mixing ratio (proportional to the 3rd moment and commonly used in single-moment schemes) and total number concentration (proportional to the 0th moment and commonly used with the mass mixing ratio in double-moment schemes). The triple-moment scheme of Milbrandt and Yau [41] additionally prognoses reflectivity, which is proportional to the 6th moment of the SD. Commonly, either exponential or gamma size distributions are assumed with the parameters derived from the prognostic moments of the SDs. The degrees of freedom for each SD and, thus, the range of SDs that can be modeled depend upon the number of moments in the scheme. A triple-moment scheme, for example, effectively allows the three parameters of the gamma distribution to vary independently, whereas a single-moment scheme typically only allows the mass mixing ratio to vary with fixed number concentration.
The shapes or aspect ratios of hydrometeors are usually not predicted in cloud models but are assumed, often very simply, within the microphysical parameterization scheme. An exception to this is the bulk adaptive habit model (AHM; Harrington et al. [42] and Sulia and Kumjian [43,44]), which explicitly predicts the evolution of ice particle aspect ratio due to depositional growth and sublimation. For example, it is known that nonspherical ice becomes more nonspherical during depositional growth and less nonspherical during sublimation. The aspect ratio of ice particles also changes during riming and aggregation (Li et al. [45]). The aspect ratios of different habits of ice crystals as a function of size are summarized by Matrosov et al. [46]. Snow aggregates and irregularly shaped ice particles are more spherical than pristine crystals, with aspect ratios usually varying between 0.5 and 0.7 (Straka et al. [47] and Korolev and Isaac [48]). Reasonable estimates of the aspect ratios of ice/snow particles can be made using polarimetric radar measurements, as suggested by Matrosov et al. [49] and Matrosov [50]. The shapes and aspect ratios of melting crystals and snow are not as well-documented as those for dry snow. It is reasonable to assume that the shape of the melting ice particle gradually changes with increasing mass water fraction f w so that it eventually acquires the shape of the raindrop with the same mass. The forward operators put forth in Jung et al. [51] and Ryzhkov et al. [15] suggest varying the aspect ratio of melting snowflakes and hailstones as a linear function of f w . For raindrops, there is a wide consensus that the relation suggested by Brandes et al. [52] is a good approximation for the dependence of the axis ratios of raindrops on their equivolume diameters.
The distributions of hydrometeor orientations also are not specified by microphysics schemes in NWP models and have to be defined within the PRFO. It is reasonable to assume Gaussian canting angle distributions with a mean canting angle of 0 • for raindrops and ice particles of oblate shapes (Ryzhkov et al. [15]) and some non-zero width σ α . For rain, assuming σ α = 10 • is a reasonable choice. In ice and snow, σ α may vary from 10 • for pristine crystals and lightly aggregated snowflakes to Atmosphere 2020, 11, 362 7 of 34 40 • for heavily aggregated dry snowflakes (see discussion in Ryzhkov and Zrnic [3]), Section 4.3.6). For prolate hydrometeors, it is more appropriate to assume random orientation distributions with some flutter angle (Melnikov and Straka [53]). As with the aspect ratio, the width of the canting angle distribution σ α for mixed-phase particles is commonly assumed to be a linear function of f w so that the width of the distribution decreases as f w increases, as it does during the melting process (Ryzhkov et al. [15]). Ryzhkov et al. [54] mention the possibility to retrieve the orientation parameter σ α using a combination of LDR, Z DR , and hv (their Equation (12)).
Our theoretical simulations show that the change of the aspect ratio of ice particles from 0.5 to 0.7 leads to the decrease of K DP (in deg km −1 ) and Z DR (in dB) by a factor of two. A similar reduction of K DP and Z DR occurs if the parameter σ α increases from 10 • to 30 • .
Formulas for the dielectric constant ε of fresh water and solid ice as functions of radar wavelength and temperature can be found in Ray [55]. The mixing formulas for ε for dry and wet snow, graupel, and hail composed of air, ice, and water are summarized in Ryzhkov et al. [15] and Ryzhkov and Zrnic [3]. As mentioned before, there is an inherent uncertainty associated with the use of Maxwell Garnett mixing formulas, because the way water is distributed within the melting particles affects the value of their dielectric constant. Fabry and Szyrmer [39] describe at least six different models of such distribution.
The angular moments characterizing the orientations of atmospheric particles are determined by the antenna elevation angle β (see Figure 1) and the width of the canting angle distribution σ α . The corresponding formulas for different particle orientation types can be found in Ryzhkov [9], Ryzhkov et al. [15], and Ryzhkov and Zrnic [3] (their appendix B).
If the distribution of particle orientations is independent of the distributions of size and shape, then the polarimetric radar variables (e.g., Z H,V , Z DR , K DP , LDR, A h , A DP , and hv ) can be estimated by integrating (or averaging) the products of the scattering amplitudes s a,b (0,π) and angular moments A 1 -A 5 over the SD, as specified in the Appendix A (Equations (A1)-(A8)) where the angular moments are defined via (A10). In the Rayleigh approximation, the scattering amplitudes are given by (3)-(5), whereas lookup tables or scattering databases should be utilized for resonance-size particles. The effects of attenuation and differential attenuation on Z H , Z DR , and LDR are quantified by Equations (A11)-(A13) in the Appendix A. The observed radar variables can also be affected by beam broadening and depolarization on propagation. The latter is caused by the cross-coupling of horizontally and vertically polarized waves when they are simultaneously transmitted (and received). The associated biases can be significant and difficult to quantify (Ryzhkov [56], Ryzhkov and Zrnic [57], and Hubbert et al. [58]). None of the existing PRFOs fully account for such artifacts in the radar observations, although some of them attempt to capture the impact of beam broadening on Z H .
Over the past decade, several PRFOs have been developed. One of the first ones was designed by Pfeifer et al. [59] for the European COnsortium for Small-scale MOdeling (COSMO) model (Baldauf et al. [60]). It was able to simulate Z H , Z DR , and LDR for rain, dry and wet snow, and graupel using T-matrix computations. Several assumptions were made regarding "free parameters" characterizing the shape and orientation of hydrometeors which are not predicted by COSMO. The shape of raindrops was assumed to follow the dependence on diameters specified by Andsager et al. [61] and the aspect ratio of ice particles varied from 0.3 to 1.0. A problematic feature of the Pfeifer et al. PRFO is its unrealistic assumptions about the mean canting angle of hydrometeors, which were allowed to vary from 0 to 40 • to produce sufficiently large values of LDR. Such an assumption is apparently inconsistent with polarimetric observations (e.g., Ryzhkov et al. [62]), and all subsequent PRFOs imply a zero mean canting angle, with the only exception in the case of ice crystals oriented by strong electric fields (Ryzhkov and Zrnic [57]).
In the United States, Jung et al. [51] proposed another PRFO coupled with the advanced regional prediction system (ARPS) (Xue et al. [63]) that additionally computes K DP and hv . This PRFO was further refined by Dawson et al. [64] and Snyder et al. [65,66]. This PRFO makes more realistic assumptions about particle orientations. The mean canting angle of hydrometeors is set to 0 • , raindrops are assumed equioriented, and the width of the canting angle distribution σ α for ice decreases from 20 • for dry snow and 60 • for dry hail to 0 • depending on the fractional water content f w for melting snow and hail. The Brandes et al. [52] relation for the aspect ratio of raindrops is used, and a fixed aspect ratio of 0.75 is assumed for snow and hail. Since the fractional water content f w is not determined/predicted in models such as ARPS that use popular bulk microphysics schemes, Jung et al. [51] and Dawson et al. [64] suggested creating artificial classes of "melting snow" and "melting hail" by "borrowing" some amount of water and ice from the pure liquid and ice categories in the model and mixing them so that the value of f w is determined from the ratio of liquid and ice water contents. In the Jung et al. [51] PRFO, a fixed density of 0.1 g cm −3 is assumed for dry snow, and the density of wet snow increases from 0.1 to 1.0 g cm −3 proportional to f w as snowflakes melt.
Ryzhkov et al. [15] developed an advanced PRFO for the Hebrew University cloud model (HUCM) with spectral microphysics (Khain and Pinsky [67]). This operator can also be used with cloud models employing bulk microphysics. It parameterizes the aspect ratio and the canting angle distribution width of mixed-phase hydrometeors as functions of f w, which is explicitly determined by the spectral microphysics scheme. Additionally, the operator utilizes a two-layer version of the T-matrix code and modified Rayleigh formulas for ice-coated freezing drops and water-coated snow, graupel, and hail, though the more common T-matrix code that assumes a homogeneous mixture can be used for any of these species as well. The radar variables are expressed via multiple angular moments characterizing particle orientations, and the effects of attenuation/differential attenuation are taken into account (see Appendix A).
Efforts to develop additional PRFOs have renewed over the past several years (Augros et al. [68], Wolfensberger and Berne [69], Matsui et al. [70], Mendrok et al. [71], and Oue et al. [72]). Augros et al. [68] created a PRFO for the French nonhydrostatic mesoscale research NWP model Meso-NH (Lafore et al. [73]). This operator simulates the whole set of polarimetric radar variables for rain, dry snow, and dry and wet graupel at the S, C, and X bands. Beam broadening and bending are accounted for along with attenuation in rain. Wolfensberger and Berne [69] developed an advanced PRFO for COSMO, and Mendrok et al. [71] reported on the progress to add polarimetric capabilities to the Efficient Modular Volume scanning RADar Operator (EMVORADO; Zheng et al. [74]) in COSMO and the ICOsahedral Nonhydrostatic model (ICON). These PRFOs include the effects of beam broadening, with the latter two computing Doppler radar variables (including the Doppler spectrum in Wolfensberger and Berne [69]). Similar to Jung et al. [51] and Dawson et al. [64], the Wolfensberger and Berne forward operator generates an artificial class of melting particles. The effects of beam broadening, attenuation, and atmospheric refraction are explicitly treated. A distinctive feature of this operator is that the aspect ratios and orientation parameters of snow and graupel particles are formulated as functions of particle size based on the results of massive observations with a multi-angle snowflake camera (MASC).
Matsui et al. [70] introduced a synthetic polarimetric radar simulator and retrieval package, POLArimetric Radar Retrieval and Instrument Simulator (POLARRIS), which generates polarimetric radar variables as well as vertical Doppler velocity, rain rate, and a synthetic hydrometeor classification product retrieved from cloud-resolving model output. The Cloud Resolving Model Radar Simulator (CR-SIM) (Oue et al. [72]) produces polarimetric and Doppler radar moment variables and lidar observables for a wide range of radar and optical frequencies. In its polarimetric radar component, CR-SIM utilizes the set of formulas and assumptions from the PRFO of Ryzhkov et al. [15].

Spectral Bin Models
The coupling of PRFOs with cloud models of varying complexity that use spectral bin microphysics, including the HUCM and a number of Lagrangian 1D or 2D cloud models, has aided the identification of "polarimetric fingerprints" of various microphysical processes and has been useful for finding corresponding deficiencies in microphysical parameterization schemes. One of the most ubiquitous polarimetric signatures observed in deep moist convection is the Z DR column. Z DR columns are vertical Atmosphere 2020, 11, 362 9 of 34 protrusions of nontrivially positive Z DR (i.e., exceeding 1 dB) above the environmental 0 • C level and indicate the presence of wet ice particles and large supercooled raindrops lofted in updrafts. In Figure 2, a Z DR column is displayed for a tornadic storm in Oklahoma that also produced large hail, which is marked by very high Z and near-zero Z DR down to the surface (at a range of 110 km from the radar). The Z DR column stretching up to 8 km in height is observed in near proximity of the hail core and is associated with a strong updraft and bounded weak echo region (BWER).

Spectral Bin Models
The coupling of PRFOs with cloud models of varying complexity that use spectral bin microphysics, including the HUCM and a number of Lagrangian 1D or 2D cloud models, has aided the identification of "polarimetric fingerprints" of various microphysical processes and has been useful for finding corresponding deficiencies in microphysical parameterization schemes. One of the most ubiquitous polarimetric signatures observed in deep moist convection is the ZDR column. ZDR columns are vertical protrusions of nontrivially positive ZDR (i.e., exceeding 1 dB) above the environmental 0 °C level and indicate the presence of wet ice particles and large supercooled raindrops lofted in updrafts. In Figure 2, a ZDR column is displayed for a tornadic storm in Oklahoma that also produced large hail, which is marked by very high Z and near-zero ZDR down to the surface (at a range of 110 km from the radar). The ZDR column stretching up to 8 km in height is observed in near proximity of the hail core and is associated with a strong updraft and bounded weak echo region (BWER). Initially, the HUCM was not able to simulate realistic ZDR columns when applying the PRFO of Ryzhkov et al. [15], because the mechanism of slow freezing of lofted raindrops was missing. This is illustrated in Figure 3 where simulated and observed fields of Z and ZDR are shown in a vertical cross-section of an Oklahoma hailstorm. Indeed, the simulated ZDR column is 2 km shorter than the observed one. Initially, the HUCM was not able to simulate realistic Z DR columns when applying the PRFO of Ryzhkov et al. [15], because the mechanism of slow freezing of lofted raindrops was missing. This is illustrated in Figure 3 where simulated and observed fields of Z and Z DR are shown in a vertical cross-section of an Oklahoma hailstorm. Indeed, the simulated Z DR column is 2 km shorter than the observed one.
After a slow-freezing process was implemented in the HUCM model, much more realistic and taller Z DR columns were reproduced (Kumjian et al. [75], Snyder et al. [76], and Ilotoviz et al. [77]). The evolution of vertical cross-sections of Z H , Z DR , and vertical air velocity w during the mature stage of a simulated hail-bearing storm after this change was made are illustrated in Figure 4. The evolution and cyclic growth and demise of the Z DR columns, and their direct association with vertical air velocities, are consistent with radar observations. Cumulative histograms of the height of Z DR columns above the 0 • C isotherm from idealized HUCM simulations and observations from tornadic storms in Oklahoma on 19 May 2013 are shown in Figure 5. While differences between the observed and simulated cases preclude a more direct comparison, the overall climatology and distribution of simulated Z DR columns appear to be in good agreement with those observed, lending confidence to the HUCM's underlying microphysics responsible for Z DR column development and evolution. After a slow-freezing process was implemented in the HUCM model, much more realistic and taller ZDR columns were reproduced (Kumjian et al. [75], Snyder et al. [76], and Ilotoviz et al. [77]). The evolution of vertical cross-sections of ZH, ZDR, and vertical air velocity w during the mature stage of a simulated hail-bearing storm after this change was made are illustrated in Figure 4. The evolution and cyclic growth and demise of the ZDR columns, and their direct association with vertical air velocities, are consistent with radar observations. Cumulative histograms of the height of ZDR columns above the 0ºC isotherm from idealized HUCM simulations and observations from tornadic storms in Oklahoma on 19 May 2013 are shown in Figure 5. While differences between the observed and simulated cases preclude a more direct comparison, the overall climatology and distribution of simulated ZDR columns appear to be in good agreement with those observed, lending confidence to the HUCM's underlying microphysics responsible for ZDR column development and evolution. The anatomy of Z DR columns simulated with the HUCM is described in the study of Kumjian et al. [75], whereas Ilotoviz et al. [77] demonstrated a strong correlation between the height of the Z DR column above the freezing level and the updraft speed ( Figure 6). Moreover, Ilotoviz et al. [77] also showed that the height and intensity of Z DR columns depend on the cloud condensation nuclei (CCN) concentration near the surface.
Another example of how PRFOs helped improve the microphysical parameterization in the HUCM was the implementation of spontaneous raindrop breakup to supplement collision-induced breakup, which resulted in a reduction of the number of extremely large raindrops and an associated decrease in simulated Z DR values in much better agreement with the observations (Ilotoviz et al. [77]). In Figure 7, the results of simulations of Z and Z DR in a hailstorm are compared if spontaneous breakup is turned on (left panels) and off (right panels). Abnormally high values of Z and Z DR are obtained in the absence of spontaneous breakup. Ilotoviz et al. [77] also checked the consistency of the model values of Z and Z DR at different altitudes in a simulated hailstorm with a large statistical dataset containing polarimetric radar data collected in a multitude of hailstorms across the US with more than 3000 hail reports (Ortega et al. [78]). Distributions of Z and Z DR in different height intervals for different hail size categories at the surface simulated by HUCM (for a single storm) and reported from observations are displayed in Figure 8. It is obvious that the HUCM is able to reproduce quite realistic values of Z and Z DR at various height levels in a hailstorm. Atmosphere 2020, 11, x FOR PEER REVIEW 11 of 35 The polarimetric emulator used to produce these images is described in Ryzhkov et al. [15]. White letters "A", "B", and "C" highlight three ZDR columns. From Ryzhkov and Zrnic [3]. Used with permission from Springer. The polarimetric emulator used to produce these images is described in Ryzhkov et al. [15]. White letters "A", "B", and "C" highlight three Z DR columns. From Ryzhkov and Zrnic [3]. Used with permission from Springer.  The anatomy of ZDR columns simulated with the HUCM is described in the study of Kumjian et al. [75], whereas Ilotoviz et al. [77] demonstrated a strong correlation between the height of the ZDR column above the freezing level and the updraft speed ( Figure 6). Moreover, Ilotoviz et al. [77] also showed that the height and intensity of ZDR columns depend on the cloud condensation nuclei (CCN) concentration near the surface. Another example of how PRFOs helped improve the microphysical parameterization in the HUCM was the implementation of spontaneous raindrop breakup to supplement collision-induced breakup, which resulted in a reduction of the number of extremely large raindrops and an associated decrease in simulated ZDR values in much better agreement with the observations (Ilotoviz et al. [77]). In Figure 7, the results of simulations of Z and ZDR in a hailstorm are compared if spontaneous breakup is turned on (left panels) and off (right panels). Abnormally high values of Z and ZDR are obtained in the absence of spontaneous breakup. Ilotoviz et al. [77] also checked the consistency of the model values of Z and ZDR at different altitudes in a simulated hailstorm with a large statistical dataset containing polarimetric radar data collected in a multitude of hailstorms The anatomy of ZDR columns simulated with the HUCM is described in the study of Kumjian et al. [75], whereas Ilotoviz et al. [77] demonstrated a strong correlation between the height of the ZDR column above the freezing level and the updraft speed ( Figure 6). Moreover, Ilotoviz et al. [77] also showed that the height and intensity of ZDR columns depend on the cloud condensation nuclei (CCN) concentration near the surface. Another example of how PRFOs helped improve the microphysical parameterization in the HUCM was the implementation of spontaneous raindrop breakup to supplement collision-induced breakup, which resulted in a reduction of the number of extremely large raindrops and an associated decrease in simulated ZDR values in much better agreement with the observations (Ilotoviz et al. [77]). In Figure 7, the results of simulations of Z and ZDR in a hailstorm are compared if spontaneous breakup is turned on (left panels) and off (right panels). Abnormally high values of Z and ZDR are obtained in the absence of spontaneous breakup. Ilotoviz et al. [77] also checked the consistency of the model values of Z and ZDR at different altitudes in a simulated hailstorm with a large statistical dataset containing polarimetric radar data collected in a multitude of hailstorms across the US with more than 3000 hail reports (Ortega et al. [78]). Distributions of Z and ZDR in different height intervals for different hail size categories at the surface simulated by HUCM (for a single storm) and reported from observations are displayed in Figure 8. It is obvious that the HUCM is able to reproduce quite realistic values of Z and ZDR at various height levels in a hailstorm.
In a similar study, Kumjian and Prat [79] used the column spectral bin model of warm rain from Prat and Barros [80] and the PRFO from Ryzhkov et al. [15] to reveal problems in the model treatment of raindrop coalescence and breakup through comparison with real polarimetric radar observations. Similar approaches employing the PRFO of Jung et al. [51] for the evaluation of existing microphysical parameterization schemes were used by Johnson et al. [81,82] for idealized supercells and an observed typhoon, respectively.   [77,78]. From Ilotoviz et al. [77]. Used with permission from American Meteorological Society.
The use of simple Lagrangian spectral bin models combined with PRFOs has also proven to be very efficient for studying the mechanisms behind the "polarimetric fingerprints" of individual microphysical processes. These include models for size sorting (Kumjian and Ryzhkov [83,84]), evaporation (Kumjian and Ryzhkov [85]), freezing of raindrops in convective updrafts (Kumjian et al. [86]), melting of hail (Ryzhkov et al. [22]), and melting of snow (Carlin and Ryzhkov [87]). As an example, Kumjian and Ryzhkov [84] adequately reproduced the ubiquitous polarimetric signatures of size sorting, manifested as enhanced ZDR combined with low ZH and KDP in rain, using an idealized 2D bin model of raindrop fallout. The results from a model of a 2D rain shaft in the presence of vertical wind shear are illustrated in Figure 9, where a displacement of the ZDR maximum near the surface with respect to the maxima of ZH and KDP is reproduced and compares well with the results of observations shown in Figure 10. The observed fields of Z and ZDR beneath the melting layer are qualitatively similar to the model rain shaft in Figure 9. Namely, the highest ZDR is located at the leading edge of the shaft along the gradient in Z (located at about 29 km in range and enclosed in the black oval), whereas the higher Z is offset and coincident with lower ZDR (30-31 km in range). The definition of height intervals/classes is in [77,78]. From Ilotoviz et al. [77]. Used with permission from American Meteorological Society.
In a similar study, Kumjian and Prat [79] used the column spectral bin model of warm rain from Prat and Barros [80] and the PRFO from Ryzhkov et al. [15] to reveal problems in the model treatment of raindrop coalescence and breakup through comparison with real polarimetric radar observations. Similar approaches employing the PRFO of Jung et al. [51] for the evaluation of existing microphysical parameterization schemes were used by Johnson et al. [81,82] for idealized supercells and an observed typhoon, respectively.
The use of simple Lagrangian spectral bin models combined with PRFOs has also proven to be very efficient for studying the mechanisms behind the "polarimetric fingerprints" of individual microphysical processes. These include models for size sorting (Kumjian and Ryzhkov [83,84]), evaporation (Kumjian and Ryzhkov [85]), freezing of raindrops in convective updrafts (Kumjian et al. [86]), melting of hail (Ryzhkov et al. [22]), and melting of snow (Carlin and Ryzhkov [87]). As an example, Kumjian and Ryzhkov [84] adequately reproduced the ubiquitous polarimetric signatures of size sorting, manifested as enhanced Z DR combined with low Z H and K DP in rain, using an idealized 2D bin model of raindrop fallout. The results from a model of a 2D rain shaft in the presence of vertical wind shear are illustrated in Figure 9, where a displacement of the Z DR maximum near the surface with respect to the maxima of Z H and K DP is reproduced and compares well with the results of observations shown in Figure 10.
The observed fields of Z and Z DR beneath the melting layer are qualitatively similar to the model rain shaft in Figure 9. Namely, the highest Z DR is located at the leading edge of the shaft along the gradient in Z (located at about 29 km in range and enclosed in the black oval), whereas the higher Z is offset and coincident with lower Z DR (30-31 km in range).   [84]. Used with permission from American Meteorological Society.

Bulk Models
Most operational NWP models utilize bulk microphysical parameterizations, which commonly prescribe a priori exponential or gamma-size distributions of particles with different habits; simulate mass mixing ratios, total number concentrations, and/or reflectivity for pure liquid or pure ice hydrometeors; and do not explicitly treat mixed-phase hydrometeors such as melting hail or snow that produce the most notable and pronounced polarimetric radar signatures. As mentioned in Section 2.2, Jung et al. [51], Dawson et al. [64], and Wolfensberer and Berne [69] suggested creating artificial classes of mixed-phase particles in order to address this deficiency.
With such modifications, several authors have had success in reproducing realistic polarimetric signatures in areas of mixed-phase precipitation using outputs of cloud models with double-or triple-moment bulk microphysics. Jung et al. [88] simulated a supercell storm using the ARPS model with single-and double-moment microphysics. The PRFO described by Jung et al. [51]

Bulk Models
Most operational NWP models utilize bulk microphysical parameterizations, which commonly prescribe a priori exponential or gamma-size distributions of particles with different habits; simulate mass mixing ratios, total number concentrations, and/or reflectivity for pure liquid or pure ice hydrometeors; and do not explicitly treat mixed-phase hydrometeors such as melting hail or snow that produce the most notable and pronounced polarimetric radar signatures. As mentioned in Section 2.2, Jung et al. [51], Dawson et al. [64], and Wolfensberer and Berne [69] suggested creating artificial classes of mixed-phase particles in order to address this deficiency.
With such modifications, several authors have had success in reproducing realistic polarimetric signatures in areas of mixed-phase precipitation using outputs of cloud models with double-or triple-moment bulk microphysics. Jung et al. [88] simulated a supercell storm using the ARPS model with single-and double-moment microphysics. The PRFO described by Jung et al. [51]  Adapted from Kumjian and Ryzhkov [84]. Used with permission from American Meteorological Society.

Bulk Models
Most operational NWP models utilize bulk microphysical parameterizations, which commonly prescribe a priori exponential or gamma-size distributions of particles with different habits; simulate mass mixing ratios, total number concentrations, and/or reflectivity for pure liquid or pure ice hydrometeors; and do not explicitly treat mixed-phase hydrometeors such as melting hail or snow that produce the most notable and pronounced polarimetric radar signatures. As mentioned in Section 2.2, Jung et al. [51], Dawson et al. [64], and Wolfensberer and Berne [69] suggested creating artificial classes of mixed-phase particles in order to address this deficiency.
With such modifications, several authors have had success in reproducing realistic polarimetric signatures in areas of mixed-phase precipitation using outputs of cloud models with double-or triple-moment bulk microphysics. Jung et al. [88] simulated a supercell storm using the ARPS model with single-and double-moment microphysics. The PRFO described by Jung et al. [51] was applied to examine its ability to simulate polarimetric signatures reported in observational studies (e.g., Kumjian and Ryzhkov [89]). The authors compared the simulations and found that certain signatures, such as the Z DR arc and midlevel Z DR and ρ hv rings, cannot be reproduced when using single-moment microphysics. Snyder et al. [65,66] used the same model with the triple-moment microphysics scheme of Milbrandt and Yau [41] and a similar polarimetric operator to reproduce realistic-looking Z DR and K DP columns, as well as the Z DR and ρ hv midlevel rings enclosing convective updrafts in simulated supercells. The results of simulations of Z H , Z DR , ρ hv , and w at a height of 4.6 km illustrated in Figure 11 are generally consistent with radar observations of a typical supercell storm (panel c in Figure 11 [89]). The Z DR arc signature, an important attribute of supercell storms resulting from size sorting, was also successfully reproduced in the study of Dawson et al. [64] (see Figure 12).
certain signatures, such as the ZDR arc and midlevel ZDR and ρhv rings, cannot be reproduced when using single-moment microphysics. Snyder et al. [65,66] used the same model with the triplemoment microphysics scheme of Milbrandt and Yau [41] and a similar polarimetric operator to reproduce realistic-looking ZDR and KDP columns, as well as the ZDR and ρhv midlevel rings enclosing convective updrafts in simulated supercells. The results of simulations of ZH, ZDR, ρhv, and w at a height of 4.6 km illustrated in Figure 11 are generally consistent with radar observations of a typical supercell storm (panel c in Figure 11 [89]). The ZDR arc signature, an important attribute of supercell storms resulting from size sorting, was also successfully reproduced in the study of Dawson et al. [64] (see Figure 12).

Microphysical Retrievals
Multi-parameter measurements available from dual-polarization radars provide ample opportunities for the retrieval of key microphysical variables: rain or snow rates, liquid or ice water contents, and particle size distribution parameters such as the mean volume diameter or number concentration. The latest methodologies for polarimetric microphysical retrievals are described in the monograph of Ryzhkov and Zrnic [3] (Chapter 11). A brief summary of those is presented herein. Figure 11. The ρhv and ZDR rings simulated in a supercell storm at ~4.6 km above ground level (AGL), as shown by (a) ρhv (colors) and w (contoured every 10 m s −1 ), (b) Z (dBZ; colors) and ZDR (dB; contoured every 1 dB), and a typical midlevel ρhv ring observed in the supercell (panel c). The polarimetric variables are calculated for the X band. Adapted from Snyder et al. [65,66] and Kumjian and Ryzhkov [89]. Used with permission from American Meteorological Society. Figure 11. The ρ hv and Z DR rings simulated in a supercell storm at~4.6 km above ground level (AGL), as shown by (a) ρ hv (colors) and w (contoured every 10 m s −1 ), (b) Z (dBZ; colors) and Z DR (dB; contoured every 1 dB), and (c) a typical midlevel ρ hv ring observed in the supercell. The polarimetric variables are calculated for the X band. Adapted from Snyder et al. [65,66] and Kumjian and Ryzhkov [89]. Used with permission from American Meteorological Society.

Radar Microphysical Retrievals in Rain
Rain microphysical retrievals include the determination of the rain rate ®, liquid water content (LWC), mean volume diameter (Dm), and total concentration of raindrops (Nt). For the rain rate estimation, a novel method based on the combination of Ah and KDP was suggested by Ryzhkov et al. [90] and proved to be the best rainfall estimation algorithm at the S band based on the results of a large-scale evaluation of different QPE techniques using the operational network of WSR-88D radars in the United States (Wang et al. [91], Cocks et al. [92], and Zhang et al. [93]). The new R(Ah, KDP) algorithm outperforms the one of Giangrande and Ryzhkov [94] currently used on polarimetric WSR-88Ds and the Q3RAD technique previously implemented on the multi-radar multi-sensor (MRMS) platform (Zhang et al. [95]).
Prior to the introduction of polarimetric radars, LWC was estimated using only ZH. Such estimates are notoriously inaccurate and sensitive to the variability of the raindrop size distribution (DSD; see Figure 13a). Much better accuracy is achieved if either a combination of ZH and ZDR or Ah is used (Figure 13b,c). The scatterplots in Figure 13 were obtained using the following radar relations optimized for DSDs measured in Oklahoma at the S band:  is in mm 6 m −3 , and ZDR is in dB. The advantage of using the LWC(Ah) relation is that it is least-sensitive to DSD variability and is immune to the biases of ZH and ZDR measurements. It can be shown that, similar to the R(Ah) algorithm, the dependence of Ah on air temperature in the LWC(Ah) relation can be ignored to a first

Microphysical Retrievals
Multi-parameter measurements available from dual-polarization radars provide ample opportunities for the retrieval of key microphysical variables: rain or snow rates, liquid or ice water contents, and particle size distribution parameters such as the mean volume diameter or number concentration. The latest methodologies for polarimetric microphysical retrievals are described in the monograph of Ryzhkov and Zrnic [3] (Chapter 11). A brief summary of those is presented herein.

Radar Microphysical Retrievals in Rain
Rain microphysical retrievals include the determination of the rain rate R, liquid water content (LWC), mean volume diameter (D m ), and total concentration of raindrops (N t ). For the rain rate estimation, a novel method based on the combination of A h and K DP was suggested by Ryzhkov et al. [90] and proved to be the best rainfall estimation algorithm at the S band based on the results of a large-scale evaluation of different QPE techniques using the operational network of WSR-88D radars in the United States (Wang et al. [91], Cocks et al. [92], and Zhang et al. [93]). The new R(A h , K DP ) algorithm outperforms the one of Giangrande and Ryzhkov [94] currently used on polarimetric WSR-88Ds and the Q3RAD technique previously implemented on the multi-radar multi-sensor (MRMS) platform (Zhang et al. [95]).
Prior to the introduction of polarimetric radars, LWC was estimated using only Z H . Such estimates are notoriously inaccurate and sensitive to the variability of the raindrop size distribution (DSD; see Figure 13a). Much better accuracy is achieved if either a combination of Z H and Z DR or A h is used (Figure 13b,c). The scatterplots in Figure 13 were obtained using the following radar relations optimized for DSDs measured in Oklahoma at the S band: where LWC is expressed in g m −3 , Z h = 10 0.1Z H (dBZ) is in mm 6 m −3 , and Z DR is in dB.  Owing to the monotonic relation between raindrop size and oblateness, the mean volume diameter Dm of raindrops, defined as the ratio of the 4th and 3rd moments of the DSD, or median volume diameter D0 (which is very close to Dm) are traditionally estimated from ZDR. A summary of various proposed D0(ZDR) relations is provided in Chapter 11 of Ryzhkov and Zrnic [3] for different radar wavelengths. Herein, we present two popular D0(ZDR) relations often used at the S band: from Cao et al. [97]. In Equations (9)-(10), D0 is in mm, and ZDR is in dB. The relations (9)-(10) were obtained using disdrometer measurements in Florida and Oklahoma, respectively. These two estimates are relatively close, and any of them can be utilized interchangeably. However, all rain retrieval relations may need some minor modifications, depending on the climate region. It would be a good practice to check the validity of the suggested relations in a concrete geographical region where sizeable DSD datasets exist and modify them if needed. The Nt (in units of number per L) can be estimated from the combination of Z and ZDR as Owing to the monotonic relation between raindrop size and oblateness, the mean volume diameter D m of raindrops, defined as the ratio of the 4th and 3rd moments of the DSD, or median volume diameter D 0 (which is very close to D m ) are traditionally estimated from Z DR . A summary of various proposed D 0 (Z DR ) relations is provided in Chapter 11 of Ryzhkov and Zrnic [3] for different radar wavelengths. Herein, we present two popular D 0 (Z DR ) relations often used at the S band: suggested by Brandes et al. [96] and from Cao et al. [97]. In Equations (9) and (10), D 0 is in mm, and Z DR is in dB. The relations (9)-(10) were obtained using disdrometer measurements in Florida and Oklahoma, respectively. These two estimates are relatively close, and any of them can be utilized interchangeably. However, all rain retrieval relations may need some minor modifications, depending on the climate region. It would be a good practice to check the validity of the suggested relations in a concrete geographical region where sizeable DSD datasets exist and modify them if needed. The N t (in units of number per L) can be estimated from the combination of Z and Z DR as where Z is expressed in dBZ, and Z DR is in dB. The standard deviation of the estimates D 0 (Z DR ) in (9) and (10) increases with D 0 , but the fractional standard deviation (FSD) is constant at about 10%. The standard deviation of the log(N t ) estimate is about 0.3 for the majority of DSDs and tends to be larger for very high (log(N t ) > 0) and very low (log(N t ) < −1) raindrop concentrations.
The described DSD retrieval is deterministic and does not consider the uncertainty caused by measurement/model errors. A statistical approach has been suggested in several recent studies based on the Bayesian theory (e.g., Cao et al. [98]) or an inverse model (Wen et al. [99]).

Radar Microphysical Retrievals in Ice and Snow
Recently, substantial progress has been made in microphysical retrievals of ice and snow from polarimetric radar. The following polarimetric relations for the estimation of the ice water content (IWC), D m , and N t for ice particles have been suggested by Ryzhkov et al. [54] and Ryzhkov and Zrnic [3]: (14) where Z DP = Z h − Z v is the reflectivity difference in mm 6 m −3 , Z dr is the differential reflectivity in linear scale, and λ is the radar wavelength in mm. In (12)- (14), IWC is expressed in g m −3 , D m is in mm, and N t is in 1/L. Equations (12)- (14) were derived assuming that the bulk density of ice or snow ρ s is inversely proportional to the particle's equivolume diameter (e.g., Brandes et al. [100]): where ρ s is expressed in g cm −3 , and D is in mm. The factor α is proportional to the degree of riming of the snowflake. Twelve different types of ice habits have been considered in the simulations, including snow aggregates, dendrites, hexagonal plates, solid thick plates, elementary needles, solid bullets, hollow bullets, solid columns, and hollow columns of various shapes. The aspect ratios of the aggregates were allowed to vary between 0.4 and 0.9, and those of pristine ice were determined as in Matrosov et al. [46]. Computations of the microphysical parameters and polarimetric radar variables were performed for a multitude of gamma size distributions: (16) and the width of the canting angle distribution σ α varying from 10 • to 40 • . The ratio Z DP /K DP in Equations (12)- (14) is very robust with respect to the variability of the particles' aspect ratios and orientations, because these two factors affect Z DP and K DP similarly, and, therefore, the ratio remains almost intact (Ryzhkov et al. [101]). Some results of simulations are presented in Figure 14 as D m versus [Z DP /(K DP λ)] 1/2 dependencies for all ice habits. Remarkably, the curves for all 12 habits corresponding to exponential size distributions are practically indistinguishable. Although the estimates from Equations (12)- (14) are almost insensitive to the variability of the shapes and orientations of ice particles, they are somewhat sensitive to the variability of the degree of riming of snow and changes of the shape parameter µ of the size distributions (Ryzhkov et al. [54,101] and Ryzhkov and Zrnic [3]). Theoretical simulations show that the fractional standard deviation (FSD) of the IWC estimates is within 20% if −1 < µ < 1, and IWC tends to be overestimated for µ < −1. The estimate (12) for the IWC is almost insensitive to the variability of α (or degree of riming). Equations (13) and (14) for D m and N t were optimized for α = 0.2. For a given α, the accuracy of the D m estimate is within 20%, but some adjustment of the Eq (13) is needed depending on the degree of riming. The accuracy of the log(N t ) estimate (14) is dependent on the accuracy of the D m retrieval and may vary from 0.7 to 1.0 (in log units) if D m is estimated with an accuracy of 20%. One of the primary advantages of using K DP instead of Z H for ice retrievals is that K DP is proportional to the 1st moment of the SD of ice particles, whereas Z H is proportional to its 4th moment (assuming an inverse dependence of the density of snowflakes on their size) and is therefore heavily weighted by the largest ice particles.
Atmosphere 2020, 11, x FOR PEER REVIEW 20 of 35 on the accuracy of the Dm retrieval and may vary from 0.7 to 1.0 (in log units) if Dm is estimated with an accuracy of 20%. One of the primary advantages of using KDP instead of ZH for ice retrievals is that KDP is proportional to the 1 st moment of the SD of ice particles, whereas ZH is proportional to its 4 th moment (assuming an inverse dependence of the density of snowflakes on their size) and is therefore heavily weighted by the largest ice particles. Reliable retrievals of size distribution parameters can be made in areas where KDP and ZDR are not close to zero. Therefore, the most accurate retrievals in the cold parts of stratiform clouds are usually produced within the dendritic growth layer (DGL) centered at the −15 °C isotherm where KDP and ZDR reach their maxima. This is important, because the bulk of snow precipitation is formed within the DGL (e.g., Hobbs and Rangno [102]). To ensure robust estimates of the microphysical variables from polarimetric measurements, KDP and ZDR have to be measured with sufficient accuracy. This dictates the need to do quite aggressive spatial averaging of KDP to reduce statistical errors of its estimate and to remove the possible miscalibration bias of ZDR [3].
The validity of Equation (12) was evaluated in the observational study of Nguyen et al. [103] using X-band polarimetric radar measurements and in situ observations onboard research aircraft flying through tropical clouds. A summary of seven research flights shows the average bias of the polarimetrically retrieved IWC to be 0.045 g m −3 and the rms error to be 0.52 g m −3 for IWC varying from 0.3 to 3 g m −3 .
Note that Equations (12)- (14) have been derived in the Rayleigh approximation and are not valid for graupel and hail, because such particles tend to be larger in size, and their density is not inversely proportional to their equivolume diameter. These formulas may also produce erroneous results for mixtures of different habits of small ice with comparable contributions to ZH and KDP. This may happen in the DGL usually centered at −15 °C where dendrites or hexagonal plates locally grow or in the temperature interval centered at −6° with favorable conditions for needle growth. In the DGL, these equations can be safely used if KDP is high and ZDR is low which means that a mixture is dominated by quasi-spherical particles (with aspect ratios 0.5-0.7) falling from higher altitudes in the cloud with a cold cloud top rather than by strongly anisotropic dendrites and plates with much lower aspect ratios which rapidly grow in the DGL (Griffin et al. [104]).
Retrievals for LWC and IWC in rain mixed with graupel/hail present additional complications, as ZH and ZDR can be completely dominated by even small IWC. Some ideas for LWC and IWC retrievals in mixtures of rain and graupel/hail using ZDR and KDP were proposed by Carlin et al. [105] based on simulations of deep convective storms using a one-dimensional model of melting Figure 14. Dependence of the mean volume diameter (D m ) of ice particles on [Z DP /(K DP λ)] 1/2 for 12 various ice habits for exponential size distribution (µ = 0) and the factor α in (15) equal to 0.2. From Ryzhkov and Zrnic [3]. Used with permission from Springer.
Reliable retrievals of size distribution parameters can be made in areas where K DP and Z DR are not close to zero. Therefore, the most accurate retrievals in the cold parts of stratiform clouds are usually produced within the dendritic growth layer (DGL) centered at the −15 • C isotherm where K DP and Z DR reach their maxima. This is important, because the bulk of snow precipitation is formed within the DGL (e.g., Hobbs and Rangno [102]). To ensure robust estimates of the microphysical variables from polarimetric measurements, K DP and Z DR have to be measured with sufficient accuracy. This dictates the need to do quite aggressive spatial averaging of K DP to reduce statistical errors of its estimate and to remove the possible miscalibration bias of Z DR [3].
The validity of Equation (12) was evaluated in the observational study of Nguyen et al. [103] using X-band polarimetric radar measurements and in situ observations onboard research aircraft flying through tropical clouds. A summary of seven research flights shows the average bias of the polarimetrically retrieved IWC to be 0.045 g m −3 and the rms error to be 0.52 g m −3 for IWC varying from 0.3 to 3 g m −3 .
Note that Equations (12)- (14) have been derived in the Rayleigh approximation and are not valid for graupel and hail, because such particles tend to be larger in size, and their density is not inversely proportional to their equivolume diameter. These formulas may also produce erroneous results for mixtures of different habits of small ice with comparable contributions to Z H and K DP . This may happen in the DGL usually centered at −15 • C where dendrites or hexagonal plates locally grow or in the temperature interval centered at −6 • with favorable conditions for needle growth. In the DGL, these equations can be safely used if K DP is high and Z DR is low which means that a mixture is dominated by quasi-spherical particles (with aspect ratios 0.5-0.7) falling from higher altitudes in the cloud with a cold cloud top rather than by strongly anisotropic dendrites and plates with much lower aspect ratios which rapidly grow in the DGL (Griffin et al. [104]).
Retrievals for LWC and IWC in rain mixed with graupel/hail present additional complications, as Z H and Z DR can be completely dominated by even small IWC. Some ideas for LWC and IWC retrievals in mixtures of rain and graupel/hail using Z DR and K DP were proposed by Carlin et al. [105] based on simulations of deep convective storms using a one-dimensional model of melting hail (Ryzhkov et al. [22]) and the two-dimensional Hebrew University cloud model (HUCM) coupled to a PRFO (Ryzhkov et al. [15]).

Illustration of the Polarimetric Microphysical Retrievals for Hurricane Harvey
Polarimetric microphysical retrievals in rain and ice have been used to investigate the microphysical structure of Hurricane Harvey (2017) as it made landfall on the Texas coast on 25 August 2017. The eyewalls of hurricanes and their outer rain bands have very different microphysical characteristics and the potential for flash flooding. Hurricane Harvey exemplifies this difference well. The composite vertical profiles of LWC/IWC, D m , and N t retrieved from columnar vertical profiles (CVP) (Murphy et al. [106]) of polarimetric radar variables have been generated in the eyewall area of the hurricane and its outer rain bands (Hu et al. [107]). The centers of the CVP columns are shown in the PPIs of Z H in Figure 15. The eyewall region was better sampled by the KCRP WSR-88D radar in Corpus Christi (left panel in Figure 15), whereas the outer rain band was better visible from the perspective of the KHGX WSR-88D radar in Houston (right panel in Figure 15). The corresponding CVPs in a height versus time format are shown in Figures 16 and 17.
A comparison of the two CVPs in Figures 16 and 17 indicates that the eyewall region and its periphery marked by the white box in the left panel of Figure 15 is characterized by a high concentration of very small ice particles aloft with Z H generally less than 20 dBZ above the melting layer. Vertical gradients of Z, LWC, and D m are obviously higher in rain below the melting layer in the eyewall area than in the outer rain bands, which tells that the bulk of precipitation in the eyewall is formed as "warm" rain, and ice aloft makes a relatively small contribution to the rainfall at the surface, although the IWC aloft is high. It is the warm rain mechanism that is primarily responsible for precipitation in the eyewall.  [106]) of polarimetric radar variables have been generated in the eyewall area of the hurricane and its outer rain bands (Hu et al. [107]). The centers of the CVP columns are shown in the PPIs of ZH in Figure 15. The eyewall region was better sampled by the KCRP WSR-88D radar in Corpus Christi (left panel in Figure 15), whereas the outer rain band was better visible from the perspective of the KHGX WSR-88D radar in Houston (right panel in Figure  15). The corresponding CVPs in a height versus time format are shown in Figures 16 and 17.
A comparison of the two CVPs in Figures 16 and 17 indicates that the eyewall region and its periphery marked by the white box in the left panel of Figure 15 is characterized by a high concentration of very small ice particles aloft with ZH generally less than 20 dBZ above the melting layer. Vertical gradients of Z, LWC, and Dm are obviously higher in rain below the melting layer in the eyewall area than in the outer rain bands, which tells that the bulk of precipitation in the eyewall is formed as "warm" rain, and ice aloft makes a relatively small contribution to the rainfall at the surface, although the IWC aloft is high. It is the warm rain mechanism that is primarily responsible for precipitation in the eyewall. The outer rain band is characterized by noticeably higher ZH in ice above the melting layer, larger ice particles, and lower Nt. The relative contribution of "warm" rain to the surface precipitation is lower, because vertical gradients of the polarimetric variables, IWC, and Dm are smaller compared to the eyewall region. Figure 17 shows that the LWC even decreases from the The combination of warm rain generated below the melting layer and rain resulting from the melting of ice originating in stronger convective updrafts makes an outer rain band of tropical cyclones more efficient at producing heavy rain compared to the eyewall. Extremely heavy rain from the stalled external rain bands was responsible for most of the flooding during Harvey.   The combination of warm rain generated below the melting layer and rain resulting from the melting of ice originating in stronger convective updrafts makes an outer rain band of tropical cyclones more efficient at producing heavy rain compared to the eyewall. Extremely heavy rain from the stalled external rain bands was responsible for most of the flooding during Harvey.   The outer rain band is characterized by noticeably higher Z H in ice above the melting layer, larger ice particles, and lower N t . The relative contribution of "warm" rain to the surface precipitation is lower, because vertical gradients of the polarimetric variables, IWC, and D m are smaller compared to the eyewall region. Figure 17 shows that the LWC even decreases from the melting layer down to the surface. This means that a good portion of precipitation at the surface is from the melting of ice. The coexistence of graupel-size ice, ice crystals, and supercooled water lofted in stronger updrafts facilitates electric charge separation and significant lightning discharges in the outer rain band. In contrast, no lightning was reported in the eyewall region of Hurricane Harvey.

Thermodynamic Retrievals
The combination of warm rain generated below the melting layer and rain resulting from the melting of ice originating in stronger convective updrafts makes an outer rain band of tropical cyclones more efficient at producing heavy rain compared to the eyewall. Extremely heavy rain from the stalled external rain bands was responsible for most of the flooding during Harvey.

Thermodynamic Retrievals
Hydrometeors undergoing a phase transition exchange latent heat with the environment, and these diabatic processes are a fundamental driver of atmospheric motion across a range of spatial and temporal scales. Condensation of water vapor in convective updrafts releases latent heat and drives convection through increases in buoyancy. Riming/accretion and refreezing represent another source of latent heat release that warms the environment. Evaporation and sublimation cool and moisten the environment, whereas the depositional growth of ice results in warming and drying. Melting of hydrometeors also causes cooling of the environment and, combined with evaporation and precipitation loading, may produce strong downdrafts associated with "cold pools" and microbursts. Polarimetric radar variables are very sensitive to the phase transitions of hydrometeors and, thus, can be very useful for thermodynamic retrievals in addition to the microphysical retrievals described in the previous section.
Carlin et al. [108] and Carlin and Ryzhkov [87] suggested a new paradigm of using polarimetric radar data for thermodynamic retrievals of warming and cooling rates associated with latent heat release and absorption. The rate of change of the environmental temperature due to latent heating can be described by (Grim et al. [109]): where N(D) is the size distribution of mixed-phase particles; dm i,subl /dt and dm i,melt /dt are the rates of ice loss due to sublimation and melting; dm w,evap /dt is the rate of water mass loss due to evaporation; and L s , L v , and L f are latent heat of sublimation, vaporization, and fusion. A similar equation can be written for the rate of change of the environmental water vapor mixing ratio q v : In Equations (17)- (18), T is the air temperature, T p is the temperature of the particle, c p is the specific heat of the air, ρ a is the air density, ρ a,d is the dry air density, k a is the thermal diffusivity of air, and dQ is the heat available for melting. The rates dm/dt are determined by environmental thermodynamic variables, whereas the size distributions of pure rain and pure ice can be retrieved from polarimetric radar measurements using Equations (6)- (14) and assuming µ = 0. These can be used directly to estimate the warming/cooling rates dT/dt and moistening/drying rates dq v /dt using Equations (17) and (18). Conceptually, if reliable SD retrievals from radar measurements are not possible through the full depth of the cloud (e.g., due to a muted polarimetric signature in regions of intense aggregation), SDs obtained in regions of reliable retrievals can instead be used to initialize 1D Lagrangian cloud models to evolve the SD in less-reliable regions. An example of such an approach is the analytical model of depositional growth and the aggregation of snowflakes of Passarelli [110] and Mitchell [111], which requires an initial intercept and slope parameter of the particle size distributions to evolve downward that could be derived from radar. In practice, quantifying the uncertainty of the SD profile obtained this way will require estimates of the uncertainty associated with both the initial SD retrieved from radar and the model used to evolve the SD downward, which should be a focus of future studies.
Diabatic warming due to the condensation of water vapor in convective updrafts is a primary driving force of deep convection. Unfortunately, it is hard to capture the condensation process directly with a radar because of the very small size of cloud water droplets, which are usually invisible to the radar at commonly used precipitation radar wavelengths. However, warming due to condensation in deep moist convections can be better quantified with polarimetric radar markers, such as Z DR column height, compared to more commonly used metrics such as Z H or the rain rate (also estimated from Z H ). Such an approach is exemplified by (though not limited to) the so-called "SLH" latent heating retrieval method described in Shige et al. [112,113]. In these approaches, convection-resolving models are coupled to radar operators to calculate reflectivity and run over a long period of time to generate a lookup table of mean vertical profiles of latent heating indexed by a radar observable such as a 10-dBZ height (e.g., see Figure 5 of [112]). These tables are then used to derive the vertical profiles of latent heating from observed reflectivity data. Such lookup table approaches have proved popular and continue to be developed and refined (e.g., Ahmed et al. [114]).
Our simulations of deep convective storms using the HUCM coupled with a polarimetric forward radar operator show that using the height of the 1-dB Z DR surface above the freezing level is a much better index of the vertical warming/cooling profiles than any Z-related parameter. This is illustrated in Figure 18, where these profiles are displayed as functions of the height of the 10-dBZ surface of Z H and 1-dB surface of Z DR in a deep hail-bearing convective storm (top panels a and b). It is obvious that the vertical cross-section of latent heating retrieved using the Z DR -based index (panel e) is in much better agreement with the "true" modeled distribution of latent heating (panel c) than the one retrieved using the Z H -based index (panel d). While the exact values of heating simulated in the lookup table shown in Figure 18b are expected to be somewhat sensitive to the model parameters used in the simulation (e.g., the environment, CCN, etc.), the correspondence between the latent-heating rate and Z DR column height overall appears robust and is in agreement with the results presented in Figure 6, as the latent-heating rate and w are inherently related.
Carlin et al. [108] demonstrated that the use of Z DR columns as proxies for convective updrafts has apparent advantages for radar data assimilation compared to the traditional utilization of Z H . One common technique of reflectivity data assimilation is using a cloud analysis, which inserts temperature and moisture increments, as well as hydrometeors deduced from Z H via empirical relations, to induce and sustain updraft circulations. In the study of Carlin et al. [108], the advanced regional prediction system's (ARPS) cloud analysis was modified from its original Z H -based formulation to provide moisture and latent heat adjustments based on Z DR columns. In this approach, observed Z DR columns are used to isolate where positive temperature and moisture perturbations are applied, and modest amounts of moisture are removed outside of these locations if/at/or near saturation.
In the top panel of Figure 19, 15-dBZ contours of Z H (in grey) and color-shaded (Z DR ) column depths between 2000 and 2300 UTC are displayed in 10-min intervals for the 19 May 2013 tornadic supercells in Oklahoma. In the two bottom panels of Figure 19, the ARPS model runs after assimilation of the Z H data only (single-pol run, left bottom panel), and the Z DR column information (dual-pol run, right bottom panel) are compared. The analyzed updraft tracks in the dual-pol run are more coherent and consolidated, with less spurious convection than in the single-pol run. The short-term forecast according to the dual-pol run shows a reduction in forward speed and northward position bias of the main updrafts, an issue encountered in many storm-scale modeling experiments. These results suggest that polarimetric radar data can provide more targeted increments of heating with positive forecast impacts.
In the melting layer and below, cooling due to the melting and evaporation of hydrometeors usually occurs. These cooling rates can be quantified using simplified 1D Lagrangian models with spectral bin microphysics coupled with a PRFO. In our analyses, we utilized 1D Lagrangian models for melting hail (Ryzhkov et al. [22]) and melting snow (Carlin and Ryzhkov [87]) and the PRFO of Ryzhkov et al. [15].
A typical vertical profile of the simulated cooling rate from our 1D model [22] below the freezing level for a hailstorm is shown in Figure 20 (right panel). This model includes melting, the soaking of meltwater into graupel and hailstones with densities less than that of solid ice, the spontaneous breakup of raindrops greater than 8 mm, and the shedding of meltwater into raindrops according to Rasmussen and Heymsfield [115]. In addition to the aforementioned processes described in [22], the sublimation of dry hail and evaporation of meltwater and shed water are included. The PRFO utilized T-matrix scattering computations and followed all assumptions for melting hail and graupel described in [15], including varying σ α as a linear function of f w and the aspect ratio during melting according to [115]. Since this spectral model explicitly and separately treats the melting and evaporation of graupel, hail, and rain, their relative contributions to the cooling rate can be estimated (left panel in Figure 20). As expected, the sublimation and melting of graupel are the dominant contributors to the cooling rate within the first kilometer beneath the melting level, whereas the melting of hail and evaporation of raindrops dominate the cooling rate below, with the evaporation of raindrops originating from the complete melting of graupel and hail particles, making the largest contribution near the surface. Combining these results with the radar operator [15], it can be shown that the vertical gradient of K DP is well-correlated with the intensity of diabatic cooling in the first two km below the freezing level.
(dual-pol run, right bottom panel) are compared. The analyzed updraft tracks in the dual-pol run are more coherent and consolidated, with less spurious convection than in the single-pol run. The short-term forecast according to the dual-pol run shows a reduction in forward speed and northward position bias of the main updrafts, an issue encountered in many storm-scale modeling experiments. These results suggest that polarimetric radar data can provide more targeted increments of heating with positive forecast impacts. The methodology for estimating cooling rates due to rain evaporation using a 1D spectral bin model and the data collected by dual-polarization radar and micro-rain radar has been developed and demonstrated in Xie et al. [116]. It was shown that polarimetric radar observations can be used to estimate the evaporative cooling rates from the retrieved DSD aloft and the evaporation model. Carlin and Ryzhkov [87] showed through simulations that the maximal cooling rates in the melting layer (ML) of stratiform precipitation can be determined from the maximal value of K DP within the ML. This is well-illustrated in Figure 21. The scatterplots of maximal cooling rates versus maximal Z H show tremendous variability caused by the diversity of size distributions of ice aloft and changes of vertical profiles of temperature and humidity within the ML (top panels in Figure 21), whereas the almost perfect linear dependencies of the max(dT/dt) on max(K DP ) (bottom panels in Figure 21) show great promise. The reason for the strong correlation between the cooling rates and K DP is that both parameters are dominated by contributions from smaller particles in the size spectrum within the melting layer, whereas Z is heavily weighted by a few large melting snowflakes. Carlin and Ryzhkov [87] show that the cooling rate due to melting is proportional to the moment of the order 1.3 of the size distribution, whereas K DP is proportional to the moment of the order 1.1.
assumptions for melting hail and graupel described in [15], including varying σα as a linear function of fw and the aspect ratio during melting according to [115]. Since this spectral model explicitly and separately treats the melting and evaporation of graupel, hail, and rain, their relative contributions to the cooling rate can be estimated (left panel in Figure 20). As expected, the sublimation and melting of graupel are the dominant contributors to the cooling rate within the first kilometer beneath the melting level, whereas the melting of hail and evaporation of raindrops dominate the cooling rate below, with the evaporation of raindrops originating from the complete melting of graupel and hail particles, making the largest contribution near the surface. Combining these results with the radar operator [15], it can be shown that the vertical gradient of KDP is wellcorrelated with the intensity of diabatic cooling in the first two km below the freezing level. (Top panels) Scatterplots of the maximal simulated cooling rates versus maximal Z within the melting layers for various size distributions of snow above the melting layers at the S, C, and X bands. (Bottom panels) Scatterplots of the maximal simulated cooling rate versus the maximal K DP within the melting layers for various size distributions of snow above the melting layers at the S, C, and X bands. Adapted from Carlin and Ryzhkov [87].

Conclusions
This paper provides an overview of several benefits that polarimetric weather radar may provide to the cloud-modeling community in an effort to improve the performance of numerical weather prediction (NWP) models and the microphysical parameterizations that they use. Polarimetric weather radar data can be leveraged by NWP models with the use of microphysical and thermodynamic retrievals or a polarimetric radar forward operator (PRFO). The retrieval concept essentially converts radar information into quantities that the model prognoses, such as liquid or ice water content, particle mean volume diameter, and number concentration, as well as the warming/cooling and moistening/drying rates. The PRFO path converts the model output into polarimetric radar fields for direct comparison with the radar observations.
A brief review of existing PRFOs is provided, with a list of necessary formulas for computing polarimetric radar variables from the model outputs presented in the Appendix A. The use of PRFOs with models employing spectral bin microphysics is relatively straightforward, whereas coupling the PRFOs with models that use bulk microphysical parameterization requires certain adaptations, such as the generation of artificial mixed-phase hydrometeor categories that are typically not produced in bulk schemes.
The techniques for microphysical polarimetric retrievals in rain are relatively well-established, as opposed to the retrievals in snow, which were introduced only recently. The retrieval equations for ice/snow have been derived assuming an inverse dependence of particle density on the equivolume diameter. Therefore, they cannot be used in areas dominated by graupel or hail. The ice microphysical retrievals heavily rely on the reliable estimation of the specific differential phase (K DP ), which may require substantial spatial averaging.
Thermodynamic radar retrievals are a new frontier of research and show good promise, as the areas of the phase transition of water, such as vapor condensation or snow melting, are usually associated with pronounced polarimetric radar signatures. The way to quantify the rates of diabatic warming/cooling and moistening/drying using polarimetric radar data is outlined in the article and illustrated by a couple of examples.
This study describes several modifications in the HUCM, such as enabling the model to reproduce realistic Z DR columns in the simulations, that result from the partnership between cloud modelers at the Hebrew University of Jerusalem and radar meteorologists at the University of Oklahoma and the National Severe Storms Laboratory. Such partnerships can be considered as a model to facilitate the use of polarimetric weather radars for optimizing the performance of NWP models in a future.
Author Contributions: A.V.R. and A.K. provided conceptualization, project administration, and supervision; A.V.R. wrote the original draft, and J.S. and J.T.C. did the review and editing; J.S. and M.P. developed the software for the forward operator and coupled it with the HUCM cloud model; and J.T.C. conducted research on the radar thermodynamic retrievals. All authors have read and agreed to the published version of the manuscript. .
where angular brackets indicate integration over particle size distribution N(D), e.g., The scattering amplitudes s (π,0) a,b are computed either using the Rayleigh formulas for uniformly filled and two-layer spheroids or numerical methods such as the T-matrix or others.
Attenuation effects are taken into account by using equations where a tilde overbar denotes the radar variables affected by attenuation/differential attenuation.