A Polarimetric Radar Operator and Application for Convective Storm Simulation

: In this study, a polarimetric radar forward model operator was developed for the Weather Research and Forecasting (WRF) model that was based on a scattering algorithm using the T-matrix methodology. Three microphysics schemes—Thompson, Morrison 2-moment, and Milbrandt-Yau supported in the operator. This radar forward operator used the microphysics, thermodynamic, and wind ﬁelds from WRF model forecasts to compute horizontal reﬂectivity, radial velocity, and polarimetric variables including differential reﬂectivity ( Z DR ) and speciﬁc differential phase ( K DP ) for S-band radar. A case study with severe convective storms was used to examine the accuracy of the radar operator. Output from the radar operator was compared to real radar observations from the Weather Surveillance Radar–1988 Doppler (WSR-88D) radar. The results showed that the radar forward operator generated realistic polarimetric signatures. The distribution of polarimetric variables agreed well with the hydrometer properties produced by different microphysics schemes. Similar to the observed polarimetric signatures, radar operator output showed Z DR and K DP columns from low-to-mid troposphere, reﬂecting the large amount of rain within strong updrafts. The Thompson scheme produced a better simulation for the hail storm with a Z DR hole to indicate the existence of graupel in the low troposphere.


Introduction
Ground-based Doppler weather radars have been routinely used by the national weather agencies of different countries to monitor hazardous weather and related precipitation [1]. With advances in radar technology during the last few decades, it has been shown that radar-derived estimates of precipitation can be greatly improved using polarimetric techniques. Polarimetric Doppler radars transmit and receive both horizontal and vertical polarized signals. From the two different reflected power returns, more information on the convective storm dynamics, kinematics, and precipitation within the developing storms can be obtained. In addition to variables from the regular Doppler radar (horizontal reflectivity Z H and radial velocity V R ), variables derived from polarimetric radars typically include some of the following: differential reflectivity (Z DR ), differential phase (Φ DP ), specific differential phase (K DP ), correlation coefficient (ρ HV ), and linear depolarization ratio (LDR). These variables contain additional information that can be used to estimate the type, shape, size, density, and orientation of hydrometeor precipitation particles within the storms [2][3][4]. Of these observables, Z DR and K DP are two major polarimetric radar variables. Z DR represents the ratio of the backscattered horizontal-to-vertical power return. Large Z DR values normally indicate the presence of oblate particles or horizontally oriented particles. Φ DP is the difference in the phase shift between horizontal and vertical signals. K DP is one-half the range derivative of Φ DP. High values of K DP indicate significant amounts of liquid water content.
Ground-based S-, C-, and X-band Doppler weather radars (single-or dual-polarization) usually provide measurements every few minutes, and rapid-scan radars with a single or limited set of elevation angles can provide measurements every few seconds to help detect tornadoes, hail, and storm-scale processes [38][39][40][41][42][43][44][45][46][47][48][49]. However, these radar observations are usually not sufficient to fully describe the convective and precipitation processes and their variations within the developing clouds. Radar observations also lack information about thermodynamic conditions, and the gaps in spatial coverage and the limited domain often constrain the value of the radar observation in studying the key processes during the development of convective storms.
Numerical simulations from high-resolution model systems are important components used to forecast and analyze severe storms as well as other atmospheric phenomena, providing continuous temporal and spatial coverage of the atmospheric conditions and their variation on a regular grid. In addition, numerical simulations contain more variables that can be used to describe the dynamics, physics, microphysics, and precipitation processes as well as cross-scale processes when compared to observational data. By comparing the numerical simulations with real observations, a deeper understanding of the phenomena and the mechanisms of and factors affecting them can be obtained. On the other hand, real observations play an essential role in validating and assessing forecast accuracy, hence helping to optimize the numerical models [50].
For many instruments (e.g., radar, lidar, satellite, as well as some unconventional instruments such as telescopes) that are currently used in the field to monitor atmospheric phenomena from the large-scale systems to turbulence, numerical models often do not provide the variables that the instruments observe. Therefore, a simulator or observational operator is needed to achieve the comparison and validation [51][52][53][54][55][56][57][58]. Specifically, in most numerical models, radar quantities, especially the polarimetric variables, are not considered as standard model output variables. In order to use such model data for verification or to understand physical processes, a radar forward operator is needed to convert the model variables into simulated radar measurements. The polarimetric radar forward operator is a vibrant research area. Pfeifer et al. [51] created a simulator using the Tmatrix method to provide synthetic polarimetric parameters Z H , Z DR , and LDR for the bulk microphysics scheme with the Consortium for Small-scale Modeling (COSMO-DE) model. Ryzhkov et al. [52] presented a radar emulator to estimate polarimetric radar variables for the Hebrew University Cloud Model with spectral microphysics. The polarimetric variables were calculated with the scattering amplitudes for hydrometeors using T-matrix codes for resonance-sized particles and Rayleigh formulas for smaller particles. This emulator has been implemented in numerous studies such as Snyder et al. [59], Ilotoviz et al. [60], and Shpund et al. [61] to examine the characteristic of polarimetric signatures for a better understanding and forecast of updraft and graupel or hail. Jung et al. [53] developed a polarimetric radar forward operator to compute Z H , vertical reflectivity Z V and Z DR , reflectivity difference Z DP , and K DP for the single-moment Milbrandt and Yau microphysics scheme for the Advanced Regional Prediction System (ARPS) model. In a later study by Jung et al. [62], the radar forward operator was updated and the doublemoment Milbrandt and Yau microphysics scheme was adopted to form an improved simulation of polarimetric signatures. The application of this radar operator includes a number of studies focused on assimilation of polarimetric radar observations for deep convective systems such as supercell [63][64][65][66]. In Snyder et al. [67,68], a three-moment microphysics scheme was built upon Jung et al. [62] to examine differences in simulated polarimetric signatures for radars at different wavelengths (X-band vs. S-band). A recent study by Oue et al. [69] demonstrated their Cloud-resolving model Radar SIMulator (CR-SIM) that could be used to compute polarimetric radar variables for multiple instruments (multi-wavelength radars and lidars) with input from multiple cloud and convective storm resolving models. Although all of the above-mentioned radar operators used the T-matrix algorithm to compute radar variables, there were differences among the operators due to the focuses of implementation and parameters used in the operator as well as the microphysics schemes. It is desired to develop different radar forward operators in the research community to provide techniques for validation of different numerical models, evaluate the accuracy and sensitivity of different assumptions in microphysics schemes for various weather systems, assess the radar retrieval methods, and provide more tools that can be used in data assimilation. This served the motivation of the current study.
The community Weather Research and Forecasting (WRF) model [70] is commonly used to simulate and analyze weather systems of different regimes and scales [71][72][73][74][75][76][77]. The above-mentioned studies demonstrate the capability of the WRF model to successfully simulate convective storms. At the same time, it was also shown that the model simulations are quite sensitive to the physics parameterizations. In this study, we developed a standalone polarimetric radar forward operator for three 2-moment microphysics schemes by Morrison et al. [78], Milbrandt and Yau [79,80], and Thompson et al. [81] that are commonly adopted in simulations of convective storms. This operator read in the WRF output files and computed polarimetric radar variables on the WRF model grid. The main purpose of this paper was to introduce of the polarimetric radar forward operator and demonstrate its capability. Another component of this study involved the analysis of observations from the National Weather Service (NWS) Weather Surveillance 1988 Doppler (WSR-88D) radar system for a severe convective storm. By comparing the simulated radar variables with the real observations, this analysis enabled us to understand more about the differences and limitations of the WRF model and the microphysics schemes. This paper is organized as follows: Section 2 introduces the radar operator. In addition, a brief discussion is given on the observational data used, the model configuration, and the experiment design. In Section 3, the output from the polarimetric radar operator is discussed. A severe thunderstorm event was used as a case study to demonstrate the radar forward operator, and the output from the radar operator was compared with the real WSR-88D radar observations. Section 4 provides a summary and addresses some unresolved questions.

Polarimetric Radar Forward Operator
In this research study, we developed a portable polarimetric radar forward model operator coded in Fortran-90. This radar operator used output from WRF model simulations to compute idealized polarimetric radar observables for a S-band radar. The operator contained the following steps, as shown in Figure 1. In the first step, the hydrometeor mixing ratios, number concentrations, and atmosphere conditions including temperature, moisture, pressure, 3-D wind (speed and direction), air density, and model grid information were collected from the WRF output files. The radar location (longitude, latitude, and height above sea level) could be defined in the parameter file. The second step started with a module to compute the radar scan coordinates (elevation angle, range, and azimuth angle) for each model Cartesian grid (x, y, and z). The T-matrix algorithm was used to calculate scattering properties. The T-matrix scattering was calculated using look-up tables by numerically integrating forward and backward scattering amplitudes with hydrometeor sizes ranging from zero to the maximum sizes for each hydrometeor. Look-up tables of forward and backward scattering amplitudes were created based on Mishchenko's T-matrix code [82] for different hydrometeors with various sizes and shapes. The radar variables were thus computed for each individual type of hydrometeor. Finally, the radar variables for all types of hydrometeors were added together and the final variables were derived for each model grid box. The hydrometeors in the radar operator included cloud, rain, ice, snow, and graupel and/or hail. The output variables from the radar operator included radial velocity (V R ), horizontal reflectivity (Z H ), vertical reflectivity (Z V ), differential reflectivity (Z DR ), and specific differential phase (K DP ). This operator could be used for radars at different frequencies (S-, C-, X-, etc.). For the particular example shown in this paper, we used S-band radar with a wavelength of 10 cm and a frequency of 3 GHz. The radar beamwidth was set as 0.91 • and radar range resolution was 250 m.
Radar reflectivity factors were computed as the total of different hydrometeor species with an integration over the particle size distribution using the following equations as given by Ryzhkov et al. [52]: Here, Z hh , Z vv , and Z vh are expressed in mm 6 m −3 , and K DP is in degrees per kilometer. The corresponding logarithmic scale values of Z H in dBZ and Z DR in dB were obtained by: In the above equations, λ is radar wavelength, |K w | 2 = 0.92 is the dielectric factor of water, N i (D) gives the number of particles within the diameter interval between D and D + dD for the specific hydrometeor type i, M is the number of different species in the volume, and D is the equivalent volume diameter for the hydrometeor size bin. Re represents the real part of the complex number, and | | is its magnitude. The coefficients A 1 -A 7 represent the angular moments of orientation distribution. f a,b (π) and f a,b (0) represent the complex backward and forward scattering amplitudes along the major (a) and minor (b) axes, respectively. The subscript * represents conjugation. The scattering amplitudes of each hydrometeor with different diameters were pre-computed and stored in the look-up tables using Mishchenko's T-matrix code (available online on http://www.giss.nasa.gov/ staff/mmishchenko/t_matrix.html, accessed on 15 April 2022). The look-up tables were built for each hydrometeor type as a function of particle size, density, aspect ratio, and air temperature. It was assumed that all hydrometeor particle types have a two-dimensional axisymmetric Gaussian distribution of particle orientations. As in Ryzhkov et al. [52], it was assumed that the mean canting angle is 0 • and the standard deviation σ is 10 • for cloud, rain, and ice, and 40 • for snow and graupel or hail. Cloud particles were assumed to have a spherical shape, whereas the other hydrometeor types were assumed to be oblate spheroids. The aspect ratio for rain particles was a function of drop size following Brandes et al. [14], and the aspect ratio for graupel or hail followed the equations given by Ryzhkov et al. [52]. The scattering amplitudes were calculated for particles with equivolume diameters of 0.1-9 mm with bin widths of 10 µm for rain, 5 µm-50 mm with bin width of 10 µm for graupel, 5 µm-70 mm with bin width of 10 µm for hail, and 0.1-50 mm with bin width of 10 µm for snow.
The operator supported three 2-moment microphysics schemes: Milbrandt and Yau, Thompson, and Morrison. The related parameters of the hydrometeor particles were defined based on the specific microphysics scheme. For the Morrison scheme, a Gamma distribution with a fixed shape factor was assumed for particle size: where N 0 , µ, and λ are the intercept, shape, and slope parameters of the size distribution, respectively. D is the particle diameter. For ice, rain, snow, and graupel, µ = 0. The number concentration of cloud droplet was set as a fixed value of 250 cm −1 . The slope was given by: where Γ is the gamma function, Q is the hydrometeor mass mixing ratio, and N T is the total number concentration of the specific hydrometeor. Parameters c and d were given by the power law mass-diameter relationship m = cD d assuming fixed densities (500 kg m −3 for ice, 100 kg m −3 for snow, and 400 kg m −3 for graupel) for the hydrometeors.
The intercept was given by: The Milbrandt and Yau scheme contains six hydrometeors: cloud, rain, ice, snow, graupel, and hail. The Gamma distribution was assumed for the particle size: where N T is the total number concentration, ν and α are dispersion parameters, and λ is the slope. For cloud droplets, α = 1 and ν = 3. For other hydrometeors, ν = 1, so a more commonly used Gamma distribution was applied: where the intercept parameter N 0 was defined as: The slope was given by: where Q is hydrometeor mass mixing ratio. Parameters c and d were given by the power law mass-diameter relationship m = cD d . In the Milbrandt and Yau scheme, hydrometeor density is defined as 100 kg m −3 for spherical snow, 400 kg m −3 for graupel, and 900 kg m −3 for hail. The Thompson scheme contains five hydrometeors: cloud, rain, ice, snow, graupel, and hail. It predicts the mixing ratio for each hydrometeor species but the number concentrations for rain and ice only. The scheme also assumes a generalized Gamma distribution similar to Equation (8) for cloud droplets, cloud ice, and rain. The slope and intercept for these species were defined as in Equations (9) and (10).
For snow, the number density distribution was: where M n is the nth moment of the snow size distribution. κ 0 = 490.6, κ 1 = 17.46, Λ 0 = 20.78, Λ 1 = 3.29, and µ s = 0.6357. For graupel, a generalized gamma distribution was assumed with the slope and intercept as: where ρ d is the air density, and the particle density ρ h is defined as 500 kg m −3 for graupel. Further details of the microphysics schemes are contained in Morrison et al. [78], Milbrandt and Yau [79,80], and Thompson et al. [81]. In these microphysics schemes, the particles are either liquid or solid. Therefore, melting or mixed-phased hydrometeors were not considered when developing the radar simulator package.

Numerical Experiments
A severe convective storm event in the afternoon of 29 May 2018 was used for a case study as a means of comparing simulated to real polarimetric radar fields. Simulations for this event were conducted with the WRF model version 3.8.1. Three nested domains were centered on the location of the event (Figure 2). The model domains had horizontal resolutions of 6.25, 1.25, and 0.25-km with 46 vertical levels, in order to adequately resolve the detailed storm scale structures with an affordable computational cost. The following physical schemes were used: the rapid radiative transfer model (RRTM; Mlawer et al. [83]) long-wave radiation, Dudhia short-wave radiation [84], Betts-Miller-Janjic cumulus parameterization [85], Mellor-Yamada_Janjic (MYJ) planetary boundary layer (PBL) schemes [85], and the Unified Noah land-surface model [86]. The WRF model outputs were saved every 5 min and were used as inputs for the radar forward operator. In this paper, the result focused on the simulation for the innermost domain. The radar operator was applied for all horizontal grids and for vertical levels from the lowest to the highest levels. Polarimetric variables were simulated for the S-band WSR-88D radar KDDC (the black star in Figure 2, 99.969 • W in longitude, 37.761 • N in latitude, and 813 m in altitude) in Dodge City, KS.

Results
To demonstrate the capability of the radar forward operator, it was applied to the WRF model simulation for the case study, as described above. In this section, we discuss the results from the radar operator, and compare the result with the real observation from the WSR-88D radar to evaluate the performance of the microphysics schemes used in the WRF model simulations.

Case Event and WSR-88D Radar Observation
On the afternoon of 29 May 2018, several multicell clusters and small-line-type convective storms formed in association with a dry-line boundary that extended from Colorado The Level II radar base data from the WSR-88D KDDC was obtained from the NOAA National Center for Environmental Information (NCEI; https://www.ncdc.noaa.gov/ nexradinv/choosesite.jsp, accessed on 15 April 2022) and processed and analyzed for the severe storms simulated by the WRF model. In order to make comparisons with the model simulations, the Level 2 radar data were interpolated into the Cartesian coordinates (latitude, longitude, and height) with 250 m horizontal grid spacing and 500 m vertical spacing using the software Py-ART [87]. Figure 3 displays radar horizontal reflectivity at 1 km (where most particles were liquid water), 3 km (where mixed-phase and melting occurred), and 8 km (where most particles were solid) height at 2346 UTC 29 May when the storm reached its mature stage. As in Figure 3, a line of deep convective cells was developing in Oklahoma to Kansas with maximum reflectivity above 55 dBZ throughout the lower to upper troposphere. Within 10 min of this time, large hail with diameters up to 2.5 inches was reported at more than 10 places within this domain.  Figure 4 shows the WRF model output for the mixing ratio and number concentration of rain and graupel or hail at 1 km altitude when the simulated storm reached the mature stage at 2345 UTC 29 May. Since both graupel and hail are predicted by the Milbrandt and Yau scheme, the sum of the two species was plotted in Figure 4e,k. The Thompson scheme does not provide number concentration for graupel, so it was not plotted here. The three microphysics schemes performed quite differently at this time. WRF_Morr produced a broader area of rainwater above 0.1 g kg −1 when compared to the other schemes. In WRF_MY, the storm area indicated by rainwater was narrower than WRF_Morr, but had comparable or even larger values of rain number concentration. WRF_Th produced the largest value for the maximum rainwater mixing ratio (7.2 g kg −1 vs. 3.5 in WRF_MY, and 3.2 in WRF_Morr) but the storm area indicated by rainwater was much narrower. In terms of rain number concentration, WRF_Th produced much higher values than the other two experiments (10 6 vs. 10 4 kg −1 ). WRF_MY produced graupel or hail in the lower troposphere with a maximum mixing ratio of 1.3 g kg −1 and WRF_Th produced graupel in the lower troposphere with a larger maximum amount of 2.8 g kg −1 . In WRF_Morr, the graupel amount was quite small (10 −4 g kg −1 ) during the entire simulation period.  Figure 5 shows Z H as produced by the radar operator corresponding to the WRF model output shown in Figure 4. Consistent with Figure 4a-c, Z H in WRF_Morr displayed a broader storm area than the other two experiments. WRF_Th predicted higher reflectivity convective cores than the other two with the maximum Z H of 60 dBZ. WRF_MY produced a narrower storm area but with the maximum Z H greater than WRF_Morr and a larger area of reflectivity >50 dBZ. For differential reflectivity, WRF_MY generated much smaller Z DR values for regions with Z H lower than 45 dBZ, indicating smaller rain particles. In WRF_Morr, Z DR was typically greater than 1.5 dB at most rainy areas with Z H above 30 dBZ, indicating the existence of large rain particles. A few areas of Z DR greater than 3 dB with the maximum value of 4.4 dB were also found to the east of strong convective regions. For example, the area near 36.4 • N, 99.4 • W was located to the east of the strong convection with Z H > 45 dBZ. This indicates that this location contained a small number of large drops of rain and lacked small drops. WRF_Th behaved differently when compared to WRF_Morr and WRF_MY. Z DR was generally 2-3 dB for regions with Z H over 30 dBZ except for the locations where graupel existed (Figure 4f). For example, Z DR was lower than 0.5 dB near 37.2 • N, 99.0 • W where Z H was above 40 dBZ, which is examined further in Section 3.3. Specific differential phase is sensitive to the amount of liquid water. Therefore, the locations of large K DP values (>0.4) agreed well with high Z H and high rainwater mixing ratio.

Radar Operator Result
The mean mass diameter (MMD) is a parameter that can be used to easily evaluate the particle size distribution. Since the Z DR field is largely proportional to the median diameter of the rain particles within the detected volume, the MMD of rain could help us understand the behavior of Z DR in the three simulations. The MMD in the unit volume for a Gamma distribution was defined as: where ρ d is the air density and ρ r is density of liquid water, q r is rainwater mixing ratio, and N Tr is the total number concentration of rain droplets. As shown in Figure 6, the overall pattern of MMD r at 1 km height generally agreed with Z DR shown in Figure 5d (Figure 5f). Figure 7 shows the mixing ratios for rainwater, snow, and graupel and/or hail in the mid-troposphere (3 km) where liquid water and ice coexist. At 3 km height, liquid water still dominated the precipitation. Similar to Figure 4, rainwater at mid-level in WRF_Morr expanded over a much broader area than the rest of the experiments. WRF_Th produced the largest maximum value (7-8 g kg −1 ), but it was concentrated over small areas. WRF_MY produced the smallest amount of rainwater with maximum value of 4.5 g kg −1 . At 3 km, the snow mixing ratio (<0.5 g kg −1 ) was much smaller than rain and graupel or hail (both >1 g kg −1 ). Among all three experiments, WRF_Th produced the least amount of snow. WRF_Morr produced more snow for the northern part of the storm, while WRF_MY produced more snow for the southern part of the storm. For graupel and hail, WRF_MY produced a larger amount than WRF_Morr within the convective centers. WRF_Th generated the largest maximum value of 7.8 g kg −1 in graupel in the convective core near 37.2 • N, 98.97 • W.    Figure 8 shows the output from the radar operator for the WRF model output displayed in Figure 7. In the mid-troposphere where both liquid and frozen water exist, Z H over the convective centers became stronger than at lower levels due to the existence of graupel or hail. This could be seen in all three experiments. Z DR values were smaller here than at lower levels due to the same reason and the existence of snow. K DP values were generally smaller than Figure 5 except in the convective cores in WRF_Morr where larger amount of rainwater was predicted at mid-level than low level. Among all three experiments, WRF_MY produced the smallest values in Z DR , which was consistent with the patterns in the lower troposphere (Figure 5e) and the smaller amount of rainwater (Figure 7b) than the other experiments. In WRF_Morr, Z DR was generally greater than 0.5 dB at locations where there was more rain than graupel and below 0.5 dB at locations where graupel dominated, which was also the case for WRF_Th. In areas where the rainwater amount was large (Figure 7c), Z DR was around 2 dB and K DP values were above 0.4 deg km −1 . For regions with a large amount of graupel (Figure 7i), Z DR was below 0.5 dB and K DP was less than 0.1 deg km −1 . One of the major benefits of this forward operator package was to provide polarimetric radar variables for each hydrometeor. Therefore, we had the opportunity to check the contribution from each individual hydrometeor species to the radar properties. Figure 9 shows Z H and Z DR at 3 km height for rain, snow, and graupel from WRF_Morr. It was seen that strong reflectivity (>40 dBZ) and large Z DR (>1 dB) appeared in convective storms from rain. Reflectivity from snow was quite small (<25 dBZ in most areas) and Z DR was near zero. For graupel, the horizontal reflectivity exceeded 40 dBZ over the northeastern part of the storm where graupel mixing ratio was larger than 1 g kg −1 . Z DR from graupel was small (below 0.2 dB).  Figure 10 shows the mixing ratio and number concentration of snow and graupel and/or hail at 8 km height (near 350 hPa). Since the number concentration of snow and graupel are not predicted in Thompson scheme, these were not included here. The total of graupel and hail were plotted for Milbrandt and Yau scheme in Figure 10e,j. With a large area of snow >2 g kg −1 and a smaller area of graupel >1 g kg −1 , WRF_Th overall predicted the most snow and the least graupel at high-level troposphere among the three experiments. WRF_MY produced comparable amounts of snow and graupel and hail as WRF_Morr but with larger maximum values (9.4 vs. 8.2 g kg −1 ) in graupel and hail over the convective centers. The overall number concentration of snow in WRF_MY was about two times higher than WRF_Morr. However, the number concentration of graupel and/or hail in WRF_Morr was generally larger than in WRF_MY. Figure 11 shows the output from the radar operator for the WRF model fields shown in Figure 10. In the upper troposphere, Z H became smaller than mid-level due to the existence of large amounts of snow and ice. Among all three experiments, WRF_Th produced the strongest echo in convective centers with the maximum Z H exceeding 60 dBZ due to the large amount of graupel (Figure 8f). Z DR in all three experiments was generally less than 0.2 dB due to the dominance of frozen particles. K DP values in high tropospheric levels were generally below 0.1 deg km −1 . Note that there were two areas with K DP above 0.2 deg km −1 in WRF_Morr, indicating the existence of supercooled water (red contour lines shown in Figure 11g) brought upward by strong updrafts. Figure 10. Mixing ratios of snow (a-c) and graupel and/or hail (d-f) for WRF_Morr, WRF_MY, and WRF_Th, respectively. Total number concentrations of snow (g-h) and graupel and/or hail (i,j) for WRF_Morr and WRF_MY, respectively. All plots are for the fields at 2345 UTC 29 May 2018 at 8 km height. Figure 11. Z H (a-c), Z DR (d-f), and K DP (g-i) for WRF_Morr, WRF_MY, and WRF_Th, respectively. Red contour lines in K DP plot show the rain water mixing ratio of 0.5 g kg −1 . All plots are for the fields at 2345 UTC 29 May 2018 at 8 km height. The star signs in the plots show the location of the KDDC radar. Figure 12 plots the NDDC radar measured reflectivity and differential reflectivity of a hail storm within the region of 36.5-37.3 • N and 98.6-99.6 • W at 2321 UTC (see the red box in Figure 3). The Level II radar data were interpolated into Cartesian coordinates by Py-ART using a Cressman function with the radius of influence increasing with distance from the radar. Storm reports indicated that hail with diameter up to 6 cm at five different places were produced by this storm from 2250 to 2330 UTC. From Figure 12a, two high Z H (>55 dBZ) areas were produced. One was in the convective core near 36.  (Figure 12d), the Z DR hole could be seen from the surface up to the mid-troposphere. On both sides of the Z DR hole, Z DR columns (large Z DR up to 3.0 dB from surface to~4 km height) were found, indicating large particles of rainwater in the lower atmosphere.  Figure 13 shows the west-east vertical cross-section (the black lines shown in Figure 5a-c) of the simulated storm from different experiments. For the storm in WRF_Morr, three strong updrafts were created along the latitude line with Z H up to 55 dBZ from the low-to mid-level troposphere. In general, large values of Z DR (>2 dB) concentrated at the surface to 2 km height with extremely large Z DR values (>3 dB) appearing below 1 km height. At the updraft centers, Z DR columns (Z DR values above 2 dB from surface to 3-4 km height) were found, which indicated large rain particles brought upward by the updrafts. The location of 0.2 dB agreed well with the 0 • C isotherm, indicating the separation of frozen and liquid precipitation particles. K DP corresponded well with the Z H field. K DP columns (large K DP up to 3 deg km −1 from surface to mid-level) were found in the convective centers with Z H > 50 dBZ. Supercooled water was seen from 9 to 12 km which caused K DP exceeding 0.2 deg km −1 (red contour line in Figure 13g). For WRF_MY, Figure 13b captures two major updrafts with Z H near 60 dBZ. The sizes of the convective centers in both horizontal and vertical directions were larger in WRF_MY than WRF_Morr. However, the Z DR columns in WRF_MY were generally smaller than in WRF_Morr with 2-3 dB confined below 2.5 km height. Again, K DP columns were also found in WRF_MY but with smaller maximum values. In Figure 13h, small areas of total K DP up to 0.4 deg km −1 were found around 10 km height. This was largely due to the supercooled water drops (contribution indicated by the red contour line of 0.2 deg km −1 ) carried up by the strong updrafts. WRF_Th produced a storm that had graupel at the surface and low troposphere. There was only one strong updraft in the center of the storm with Z H above 60 dBZ from the mid-to upper-troposphere. Beneath the strong updraft, downdraft and Z H of 55 dBZ were found from the surface to 2 km height. The Z DR plot in Figure 13f shows a Z DR hole (small Z DR less than 0.2 dB surrounded by larger values of Z DR ) in the low levels. The small Z DR with large Z H near 98.9 • W indicated the existence of graupel from the surface to the upper troposphere (confirmed by Figure 5i, Figure 7i, and Figure 10f). Comparing with the simulated storms shown in Figure 13, the observed hail storm shown in Figure 12 had a smaller horizontal span with strong echoes (55-65 dBZ) existing mostly below 6 km height. The strong echoes in the simulated storms (Figure 13b,c) reached 9 km height and above, indicating much stronger updrafts and a larger amount of graupel or hail in the upper atmosphere of the simulated storms. The observed storm showed Z DR columns (Z DR > 2 dB) in the low-to mid-troposphere in regions with Z H exceeding 35 dBZ, which was also found in all simulated storms reflecting the large amount of large rain drops. The above results showed that the radar forward operator was developed properly and demonstrated value by the ability to produce representative polarimetric radar fields. To note, specific differential phase was not Level II base data from NWS WSR-88D radar; therefore, no comparison between observations and predictions was made for K DP .

Observation and Simulated Variables for a Hail Storm
In order to better understand the existence of the Z DR column and K DP column, scatterplots with q r were generated for the three experiments. In Figure 14a-c, large Z DR values (≥1.5 dB for grid points below freezing level) and their corresponding q r values were plotted. In WRF_Morr, large Z DR ranged from 1.5 to 3.8 dB with q r up to 8.5 g kg −1 . For q r above 4 g kg −1 , a large portion of the grid points had Z DR of 2.1-2.3 dB. For grid point with Z DR > 3 dB, q r was generally less than 1 g kg −1 . For WRF_MY, Z DR was between 1.5 and 2.8 dB with q r almost evenly distributed within 0-4 g kg −1 . In addition, most of the grid points with Z DR > 2.3 dB had q r lower than 2.0 g kg −1 . For WRF_Th, Z DR values were between 1.5 and 2.5 dB with q r up to 8.2 g kg −1 . For larger q r (>6 g kg −1 ), Z DR values were typically lower than 1.8 dB. Figure 14d-f plots large K DP values (≥1.0 deg km −1 for grid points below freezing level) and their corresponding q r . In WRF_Morr and WRF_MY, K DP increased with q r . Therefore, the largest K DP agreed with the largest value in q r . In WRF_Th, such increasing trend of K DP was not apparent. For high q r values (>6 g kg −1 ), K DP spread out between 1.0 and 3.8 deg km −1 . Figure 14. Scatterplots of large Z DR vs. q r in (a-c) and K DP vs. q r in (d-f) for WRF_Morr, WRF_MY, and WRF_Th, respectively, when temperature was equal to or above 0 • C.

Discussion
In this paper, we introduced a new polarimetric radar forward operator for three different microphysics schemes: 2-moment Morrison, 2-moment Milbrandt and Yau, and Thompson. The radar forward operator was developed based on scattering calculations using the T-matrix algorithm. The radar forward operator used thermodynamic, microphysics, and wind information from WRF model simulations as input and produced polarimetric radar variables including V R , Z H , Z DR , and K DP . A case study for a severe thunderstorm on 29-30 May 2018 was examined to demonstrate the capability of the radar forward operator.
High-resolution WRF simulations for the severe thunderstorms showed that the microphysics schemes performed significantly different in terms of the type, mass, and number concentration of different hydrometeors. The Morrison scheme produced a broader area with a large amount of rain in the low-to mid-troposphere. The Thompson scheme produced a much larger rain number concentration than both the Milbrandt and Yau scheme and Morrison scheme. In addition, the significant difference in the Thompson scheme from the other schemes was the presence of a much larger amount of graupel from the surface into the lower troposphere over the convective core. In the upper troposphere, the Thompson scheme produced the largest amount of snow and the least amount of graupel among all experiments. The snow number concentration from the Milbrandt and Yau scheme was about two times larger than that in the Morrison scheme, while the graupel and hail number concentration from the Milbrandt and Yao scheme was smaller than the Morrison scheme.
The output from the radar forward operator followed the distribution of the microphysics properties. From Z H , it was shown that the Morrison scheme produced a broader storm area at low-to mid-levels than the other two experiments. At upper levels, the Milbrandt and Yau showed stronger Z H at convective centers than the Morrison scheme, but the Thompson scheme generated the strongest radar echo among all three schemes. Regarding Z DR , the Milbrandt and Yau scheme produced generally smaller Z H values, which could be explained by the much smaller mean mass diameter of rain drops. Since each hydrometeor was processed separately in the radar forward operator, it was convenient to examine the contributions from different hydrometeors. This was especially useful in mid-levels where liquid and solid or ice particles coexist. This was seen in the polarimetric variables at 3 km height for the Morrison scheme, which showed that a major part of the high reflectivity and high Z DR values were contributed by rain. For K DP , all three experiments indicated a good correlation between large K DP values and the large amount of rain. K DP could also indicate the supercooled water at high altitudes. Z DR columns and K DP columns were found in all three model simulations corresponding to the large amount of rain within the strong updrafts in the low-to mid-troposphere. The Z DR hole was found to spatially closely correspond with the existence of graupel in the lower troposphere as produced by the Thompson scheme.
The results confirmed that the radar operator was properly built for the three microphysics schemes. Since these schemes do not allow mixed-phase particles, melting was not considered in the current operator. However, melting is an important process for snow and graupel or hail in the real atmosphere. The absence of melting processes may cause errors in polarimetric variables. For future developments, the melting process is a key component to add. In addition, two other polarimetric variables, the correlation coefficient ρ HV and the linear depolarization ratio (LDR), are still under development. ρ HV is the correlation coefficient between the horizontal and vertical polarized echoes. LDR is the ratio of the vertical power return from a horizontal pulse or the horizontal power return from a vertical pulse, which is sensitive to the existence of mixed types of precipitation particles. Both variables can bring important information on the characteristics of mixed-phase precipitation particles. We must also compare the output from our operator with other publicly available radar operators (e.g., pyDualPol software available from http://arps.ou.edu/downloadpyDualPol.html, accessed on 15 April 2022). The comparison would be toward an assessment of scheme accuracy relative to radar-observed storms to gain confidence in how well our method simulates polarimetric variables. Comparing results from the pyDualPol operator to the operator developed here would help the research community evaluate the overall value of such operators.
For future work, it is our plan to produce simulated polarimetric observables and compare them with the real radar observations to validate and evaluate the microphysics schemes. Currently, our case studies (besides the event shown in the paper, we are also investigating two other heavy precipitation cases in 2018) only focus on S-band radars. Since many C-or X-band radars are also available for both research and operational communities, it is valuable to expand the present examination and include these higher frequency radars. The current forward operator was developed for WRF model output, yet with modifications in the interface modules, this operator can be adapted for other regional forecast models as well. The comparisons shown here between KDDC WSR-88D radar and WRF model simulations are only from one case study, which may not represent the general behavior of different microphysics schemes. More case studies for different precipitation regimes (light, moderate, and heavy) need to be conducted for a thorough understanding of the microphysics properties and polarimetric signatures produced by the different schemes.
In addition, as an idealized research tool, the errors that occur in the model forecast is passed to the radar operator. The impact of these errors is two-fold: (1) the errors provide a great opportunity to explore the uncertainties in the WRF model microphysics schemes through sensitivity studies, which may be an important component of our future research; (2) they allow for the validation of the WRF-model-simulated radar variables against real observations since attenuations, errors, and noise are expected in real observations. In future studies, the community should explore techniques to enhance the simulator toward creating data closer to the real radar observations.