Assimilation of Polarimetric Radar Data in Simulation of a Supercell Storm with a Variational Approach and the WRF Model

: Polarimetric radar data (PRD) have potential to be used in numerical weather prediction (NWP) models to improve convective-scale weather forecasts. However, thus far only a few studies have been undertaken in this research direction. To assimilate PRD in NWP models, a forward operator, also called a PRD simulator, is needed to establish the relation between model physics parameters and polarimetric radar variables. Such a forward operator needs to be accurate enough to make quantitative comparisons between radar observations and model output feasible, and to be computationally efﬁcient so that these observations can be easily incorporated into a data assimilation (DA) scheme. To address this concern, a set of parameterized PRD simulators for the horizontal reﬂectivity, differential reﬂectivity, speciﬁc differential phase, and cross-correlation coefﬁcient were developed. In this study, we have tested the performance of these new operators in a variational DA system. Firstly, the tangent linear and adjoint (TL/AD) models for these PRD simulators have been developed and checked for the validity. Then, both the forward operator and its adjoint model have been built into the three-dimensional variational (3DVAR) system. Finally, some preliminary DA experiments have been performed with an idealized supercell storm. It is found that the assimilation of PRD, including differential reﬂectivity and speciﬁc differential phase, in addition to radar radial velocity and horizontal reﬂectivity, can enhance the accuracy of both initial conditions for model hydrometer state variables and ensuing model forecasts. The usefulness of the cross-correlation coefﬁcient is very limited in terms of improving convective-scale data analysis and NWP.


Introduction
The accuracy of the initial conditions of numerical weather prediction (NWP) models has a great impact on their forecast results, especially for high-resolution NWP. To build up the initial hydrometeor information which is the key to successful simulating and forecasting of convective storms, the inevitable process known as "spin-up" is needed for NWP models [1]. In order to reduce the negative influence of "spin-up" on the forecast skill of short-lived convective storms, Doppler weather radar observations, including the horizontal reflectivity (Z H ) and radial velocity (V r ), are commonly used to improve the initial condition with advanced data assimilation (DA) methods, such as the three-verification of ensemble probabilistic forecasts in a real case. Subsequently, by establishing look-up tables to increase the calculation efficiency of J10 during the DA procedure, Putnam et al. [17] first implemented the direct assimilation of Z DR through the EnKF method for a real supercell case in the literature. Through a series of observing system simulation experiments (OSSEs) of an idealized supercell storm, Zhu et al. [18] investigated the impact of assimilating Z DR within an EnKF framework by using J10 coupled with a DM MP scheme, and also discussed the potential influence of updating hydrometeor number concentrations for DA effect. Their conclusions show that the assimilation of Z DR can improve the accuracy of analyzed hydrometeor fields in terms of pattern and intensity, and that updating the number concentrations with mixing ratios is important for deciding whether the benefit of assimilating Z DR can be preserved. However, they are all limited to the use of partial PRD for assimilation, and significant errors in the magnitude of analyzed PRD do exist, indicating that the operational application of PRD assimilation remains challenging.
In addition, some polarimetric signatures found in convective storms provide the opportunity for indirect assimilation of PRD. Using the Advanced Regional Prediction System (ARPS) cloud analysis module [5,49,50], Carlin et al. [51] managed to assimilate PRD through the detection of Z DR column signatures [33,52], which often occur in the strong updraft region within intense thunderstorms, to modify the temperature and moisture of the ARPS model. Encouraging results were obtained with improved analyses (more coherent updraft) and forecasts (more realistic reflectivity structures with better skill scores) compared to the original cloud analysis scheme without polarimetric information included. There are limits to the widespread use of this approach, however, due to the absence of Z DR columns in relatively weak precipitation systems and its high dependence on empirical relationships between model state variables and PRD.
Recently, Zhang et al. [53] (hereafter Z21) developed a new set of simplified and parameterized PRD observation operators which can link NWP model state variables and PRD for DA use, and verified its validity and applicability by applying an ideal case and a real case respectively. They found that the parameterized operators yield results consistent with that of rigorous calculation in J10. However, Z21 increased computational efficiency significantly, resulting in costs less than one percent of the computing time of J10 to complete the same task, which makes it more suitable for DA applications. Thus, in this study, the applicability of Z21 is investigated under a variational DA framework, and the different DA effects of each polarimetric variable are studied by conducting various OSSEs.
The rest of this paper is organized as follows. The new parameterized PRD operators are briefly reviewed in Section 2 and the overview of the 3DVAR DA system is provided in the following section. In Section 4, the experimental design is given, and some preliminary experiment results are discussed in Section 5. A summary and conclusions for this study are presented in the last section.

Microphysics Models and Parametrization
SM and DM MP schemes are commonly used in NWP models [54][55][56][57][58]. In SM simulations of NWP models, water mixing ratios (q), which is directly relevant to water content (W = ρ a q, where ρ a is the air density), is the only prognostic physics variable for hydrometeor physics. For DM MP schemes, the number concentration (N t ) and the water mixing ratio (q) are both predicted variables. For a given two-parameter drop/particle size distribution (DSD/PSD) model and a specified hydrometeor x, both N tx and q x can be converted to mass-weighted diameter (D mx = 4 ρ a q x πρ x N tx 1/3 ) and water content (W x = ρ x q x ); then, many other related quantities such as DSD/PSD parameter, rainfall rate and reflectivity are calculated by these two state variables. The detailed description of the derivation can be found in Z21.

Parameterized PRD Operators
Polarimetric radar variables (Z H , Z DR , K DP , ρ hv ) are expressed by the integrals of DSD/PSD weighted by the scattering amplitudes [59,60]. The scattering amplitudes are normally calculated by the T-matrix method, and the integrals are calculated numerically for accurate results. This was done in J10. However, the complex integral form and resulting expensive computational cost make it not convenient for variational DA use, which usually needs the operators to be efficient and easy calculation of the derivative.
Mahale et al. [61] represented polarimetric variables as simple functions of water mixing ratio q r and mass/volume-weighted diameter D m in the polynomial form for rain. In the case of ice phases and mixtures such as dry/melting snow, hail, and graupel, the calculation and parameterization are more complicated than those of rain induced by the increased variability in density during the melting stage and irregular shape, as well as the orientation of the particles. For simplification, a simple function of the percentage of melting (γ x ) as following is used in Z21 to calculate the particle density ρ x for the melting hydrometeor species x: where γ x = q r /(q r + q x ) (the same as J08), ρ dx stands for the dry/pure state of ice hydrometeor corresponding to the melting state, and ρ w is the density of water. Except for the mean axis ratio and standard deviation of the canting angles, Z21 makes some adjustments to allow a larger dynamic range of Z DR and ρ hv . The shape and orientation of hydrometer particles follow the modeling and representation documented in J08. Similar to rain, for a given ice or mixed phase species x, including snow, hail, graupel, and their melting parts, PRD variables are calculated for a set of D m at a given γ x , and then parameterized as functions of D m as follows: where the reflectivity factor is Z x = 11.25 × 10 −3 ρ a q x πρ x D 3 m (D m , unit: mm). The fitting coefficients associated with γ x are provided in Tables 1 and 2 of Z21, and more detailed descriptions for calculations, assumptions, and result analysis can also be found there.

The 3DVAR DA System
To test the assimilation performance of Z21 in a variational DA scheme, in this study, a 3DVAR DA system designed especially for radar data assimilation at the convective scale is employed. This 3DVAR system was initially developed at the Center for Analysis and Prediction of Storms (CAPS) and subsequently refined at National Severe Storms Laboratory (NSSL). The standard 3DVAR cost function can be written as follows [3]: The cost function defined above contains three terms on the right-hand side. The first term is the background term, which defines the departure of the analysis vector x from the background vector x b weighted by the inverse of the NWP model background error covariance matrix B, where the analysis vector x contains the following variables: the three wind components (u, v, and w) and the mixing ratios of hydrometeors (rain water, snow, hail, graupel, cloud water, cloud ice). The background error covariance matrix is modeled by the product of a diagonal matrix of the standard deviation of the background error and a spatial recursive filter [3]. The standard deviations for the model variables are derived from the statistics of the Rapid Update Cycle (RUC) model's 3 h forecasts over several years from 2001 to 2004 [62]. The second term is the observation term, which defines the departure of the analysis vector x from the observation vector y o , where H(x) is the forward operator which builds a "bridge" between the model state variables and observation data. In this study, the new PRD operators (Z21), which can use both the mixing ratio and the number concentration from a DM MP scheme, will be discussed later. In the second term, R is called the observation error covariance matrix which includes both instruments and representativeness errors. The observation errors' standard deviation for PRD variables (V r , Z H , Z DR , K DP , and ρ hv ) is specified empirically as 2 m/s, 3 dBZ, 0.5 dB, 0.5 • /km, and 0.1 respectively. Finally, as the third term on the right-side, J c (x) represents the dynamic or equation constraints. This term has proved to be very important, particularly for this convective-scale 3DVAR system. For example, the mass continuity equation imposed as a weak constraint is helpful in producing more suitable wind analysis [2,3].
Radar radial velocity V r is regarded as a crucial variable in the DA procedure [3,4,10] because it provides important wind field information and dynamic characteristics of the precipitation event. In this 3DVAR DA system, the forward observation operator for radial velocity considering the effects of the Earth's curvature is expressed as: where h and s is radar beam height and distance, respectively, over and along the curving Earth's surface, r is the radial slant range, and ϕ is the radar azimuth angle. The 4/3effective Earth radius mode is assumed to follow for the propagation of radar beams [59]. The effect of hydrometeor fall velocity is neglected here, and will be considered in the follow-up real data case study. For PRD, after the calculations for each species x available in specified MP scheme through Equations (2)-(5), the final variables for the pixel are obtained by summing the contributions from all species as follows: For the 3DVAR method used here, considering the variations of observational spacing existing among different observations, multiple analysis passes are usually used to analyze different observation types with different filter scales. Indeed, applying multiple-pass with a recursive filter is proven superior to the common single-pass method theoretically [6,43]. Hence, in this study, each PRD variable is analyzed in a separate pass.
The new operators have been proved efficiently in PRD simulation and therefore particularly suitable for DA purposes (Z21), including both variational and ensemble DA. In such applications, the TL/AD operators, coupled with the forward operators, are required for the variational assimilation, which are not necessarily found in the EnKF [63]. As a consequence, complex operators such as those used in J10 are more applicable to be employed through the EnKF method for DA use. Despite this difficulty, due to the simple polynomial form of Z21, the programming of corresponding TL/AD operators is implemented smoothly according to the rules of adjoint code construction [64,65].
The AD operator is the transpose of the TL operator. Therefore, how well TL operator approximates the nonlinear forward operator determines how valid the gradient of the cost function is. Following [66], to verify the correctness of the TL operator, one can define a function R(a), given by where f (x) and f (x) are forward, tangent linear operator respectively. The x 0 represents a random picked base state and dx represents a perturbation state. Here, R(a) is calculated over all points in the domain. According to the theory, if R(a) in Equation (12) is close to 1 for any small values of a (see Table 1), meaning that the TL operator approximates the nonlinear forward operator well. The test results of the TL operator are listed in Table 1, as a keeps decreasing from 1 to 1E -14 , the R(a) values for all polarimetric variables (R_Z H , R_Z DR , R_K DP and R_ρ hv ) are approaching 1 uniformly and steadily in the case of double precision. However, when a goes down to an extra small value, the R(a) starts to increase, which is caused by the rounding error of the calculator floating-point arithmetic operation. The above outcomes enable us to conclude that the assimilation system with the Z21 integrated passes the validity test and that the subsequent experimental research can be carried out reasonably.

Experimental Design
As discussed previously, PRD shows potential for enhancing storm simulations in terms of intensity and structure and ensuing short-range forecasts [16][17][18]40,[42][43][44]51]. However the forward operators used thus far are either too simple [43,44] or too complicated [39,41]. The recently proposed parameterized PRD operators in Z21 can be used in different microphysics schemes and are relatively efficient.
In this study, a fully compressible and non-hydrostatic Advanced Research Weather Research and Forecasting (WRF-ARW) model (version 3.7.1) is used in a three-dimensional space for a supercell storm simulation. The horizontal resolution is 1 km with 121 grid points in the east-west and north-south directions. Simultaneously, 51 stretched levels up to 20 km (50 hPa, approximately) above ground level (AGL) are chosen. Open boundary conditions for the lateral and Rayleigh damping for the top boundary are used in this idealized case. The whole integration of the WRF-ARW model for simulation is 3 hours.
A modified sounding from a classic supercell event that occurred on 20 May 1977 in Del City, Oklahoma is used to provide the environmental field of the storm. A thermal bubble with a potential temperature perturbation of 3 K is added in the truth field to initiate convection [14,[67][68][69]. This warm bubble is centered at the location of (x = 90 km, y = 40 km, z = 1.5 km) and has a 10 km horizontal radius and 1.5 km vertical radius inside the model domain. The Milbrandt and Yau scheme [55], a DM MP scheme consisting of the predicted mixing ratios of six-class hydrometeors, such as cloud water (q c ), cloud ice (q i ), rain water (q r ), snow (q s ), hail (q h ), and graupel (q g ), and their corresponding number concentrations (N c , N i , N r , N s , N h , and N g , respectively), is applied in description of microphysical processes in the idealized case. The time step for the integration of the prognostic equations is set to be 3 seconds. The standard 1.5-order TKE closure scheme is chosen for the turbulence parameterization, and other configurations are the same as that of Ge et al. [67].
During the 3 h truth simulation, the cloud forms around 10 min, rain water appears at 15 min, ice hydrometeors are generated at 20 min, and a single convective storm develops in the first 30 min (not shown). The storm reaches its mature stage at 45 min, then starts to split, evolves and slightly weakens. At 2 hours into the model integration, the right-splitting cell tends to dominate, as indicated by a clear hook echo and strong updraft. Figure 1 shows the cycled DA and forecast scheme employed in OSSEs. Based on the truth simulation with WRF, all polarimetric radar observations (V r , Z H , Z DR , K DP , and ρ hv ) can be simulated and used in the iterations of 3DVAR analysis to update the model state variables (u a , v a , w a , and q a x ). After the analysis step, a 5 min forecast procedure is performed and the result is the background field (u b , v b , w b , and q b x ) as the input of the DA analysis step in the next cycle. Considering the mixing ratios of precipitation are too small at the early stage of storm initiation, the idealized data are generated after the 30 min of model integration. Two assumed S-band radars are located in the southeast and southwest corner, respectively, and are performed on the VCP11 scan mode of WSR-88D radar (14 tilts from 0.5 • to 19.5 • ). Idealized PRD is derived from the truth simulation output using the forward operators in Equations (7)-(11). All OSSEs presented in this study are listed in Table 2. The control run assimilating V r and Z H data, labeled as ExpVrZh, is considered as a reference for comparison purposes. Sensitivity experiments ExpVrZhZdr, ExpVrZhKdp, and ExpVrZhRhv are the same as the control run but with the additional assimilation of Z DR , K DP , and ρ hv , respectively. ExpVrZhPol is a special run with all PRD assimilated. These sensitivity experiments are performed to help understand the impact of PRD on the analysis of hydrometer variables. Therein, the number concentrations of hydrometeors will not be updated during the DA cycles on account of the complexity introduced by their wide dynamic range (from 0 to above 1E 12 m -3 ) and the problem of a significant increase of the degrees of freedom. The explanation of all abbreviations used in this paper can be found in Appendix A. The Root Mean Square Errors (RMSEs) are commonly used as a criterion for the accuracy judgment of DA [14,18,40,68,69]. To be consistent with the previous requirement of Z H assimilation, the RMSEs of analysis and forecast displayed in Figure 2 are only calculated over precipitation grids with Z H exceeding 10 dBZ.  , (b) v, (c) w, mixing ratios of (d) cloud water q c , (e) cloud ice q i , (f) rain water q r , (g) snow q s , (h) hail q h , and graupel q g , and the same calculating procedure as truth simulations for polarimetric radar variables (i) Z H , (j) Z DR , (k) K DP , and (l) ρ hv via Z21 forward operators.
At a first glance of Figure 2, for all model state variables and simulated radar observations, the RMSEs of analysis and forecast have a generally decreasing tendency over time for each OSSE. Since number concentrations are not updated with mixing ratios during the DA process, there may be a potential imbalance between them, leading to fluctuations in the early DA stages. Compared to ExpVrZh (black line), ExpVrZhZdr (red line) produces slightly better analyses and forecasts for most variables, where more pronounced improvements occur in the last few cycles, especially in ice hydrometeor species q i (Figure 2e), q h , and q g (Figure 2h). Since K DP is most sensitive to the presence of liquid water, the RMSE of q r (Figure 2f) of ExpVrZhKdp (blue line) has the most obvious difference. The better convergence through the whole DA cycle suggests that the satisfied rain water field much closer to the truth can be analyzed with additional K DP assimilation. In fact, as with ExpVrZhZdr, it also plays a positive role in obtaining a smaller RMSE in the last few cycles for q i , q h , and q g . For ExpVrZhRhv (cyan line), the relatively distinct variations only occur in the first or last few cycles, but as a whole, its RMSEs almost follow that of ExpVrZh with no significant discrepancies, indicating the limited effect of ρ hv assimilation.
In terms of observed radar variables, the RMSEs of ExpVrZhZdr converge better and faster, indicating its simulated PRD is superior to that of ExpVrZh, especially for Z DR (Figure 2j) and Z H (Figure 2i). The above conclusions agree well with those of Jung et al. [40], except for benefits shown earlier in DA cycles here. This is largely the same as the result of Zhu et al. [18]. Apparently, ExpVrZhKdp produces consistently better K DP (Figure 2k) in both analysis and forecast cycles, as expected. Moreover, there are some subtle advancements of Z H in the last few cycles, but no contributions to the simulation of Z DR and ρ hv (Figure 2l). As confirmed by the variation of the RMSEs, the assimilation of ρ hv cannot bring benefits in observation space as well. Although both Z H and ρ hv can converge reasonably, the influence of ρ hv assimilation is too small to be noticed, which may be caused by the relatively small variation range of itself (on the order of 0.01). Noted that the reverse changing problem (the RMSE decreases in the forecast cycle, but increases in the corresponding analysis cycle) occurring in Z DR and K DP can be resolved by assimilating the corresponding PRD variable.
In the aspect of RMSEs, different abilities and effects in improving the accuracy of the model state variables and simulation of PRD have been presented for each polarimetric variable. Superior and more balanced results are expected to be achieved by jointly assimilating all PRD. As shown in purple line in Figure 2, although some fluctuations arise in w (Figure 2c), q s (Figure 2g), q h and q g at certain cycles, the best convergent trends are obtained in ExpVrZhPol for almost all kinds of hydrometeors, especially in later DA cycles. For simulated PRD, ExpVrZhPol also acquires better RMSE convergence compared to other experiments, especially for Z DR and ρ hv . Actually, with respect to most variables, the sharpest downward variation range occurring in the first cycle of ExpVrZhPol indirectly reveals the greatest potential capability of revising the model analysis field.

Evaluation of PRD Assimilation
The impacts of PRD assimilation on the horizontal distribution of intensity and polarimetric characteristics of the simulated storm are further analyzed in Figures 3-6. Compared with the true Z H (Figure 3a), all OSSEs represent the hook echo signature and the strong Z H core at the leading edge of the storm. Among them, ExpVrZh (Figure 3b) and Ex-pVrZhRhv (Figure 3e) have very similar results with obviously high Z H (greater than 15 dBZ overestimated) in the rear stratiform precipitation area (hereafter SPA). However, the Z H intensity described by ExpVrZhZdr (Figure 3c) is superior to them, both in the front convective precipitation area (CPA) and the rear SPA. However, the weak echo (<20 dBZ) on the northern edge is over-adjusted to deviate from the truth. On the contrary, although the weak echo is well analyzed, a much smaller Z H exists in the SPA of ExpVrZhKdp (Figure 3d). In terms of intensity and structure, the overall result of ExpVrZhPol (Figure 3f) is much closer to the truth, particularly in the echo edge and rear SPA, except for a certain overestimation in the CPA. There is a small Z DR patch (<1 dB) in the southeast corner of the main storm in ExpVrZh (Figure 4b), and also the same in ExpVrZhRhv (Figure 4e). Beyond that, the Z DR in the SPA is significantly smaller than the truth (Figure 4a), with the largest difference over 1.5 dB. ExpVrZhKdp (Figure 4d) has even worse simulation with much smaller Z DR . However, the intensity of Z DR in ExpVrZhZdr (Figure 4c) is much enhanced in the fore-mentioned two underestimated areas, particularly in the SPA. However, to some extent, overestimation persists in the weak echo at the edge of storm. ExpVrZhPol (Figure 4f) equally gives the best analysis among all OSSEs, except for the underestima-tion that occurred in the south-east of the storm. Although the dominant distribution of K DP is well exhibited in ExpVrZh (Figure 5b), ExpVrZhZdr ( Figure 5c) and ExpVrZhRhv (Figure 5e) without K DP assimilation, the central strength of the leading edge and the hook echo are clearly weaker. The assimilation of K DP is beneficial that the simulated K DP shown in Figure 5d is perfectly matched with the truth (Figure 5a) in both intensity and structure. Undoubtedly, as a result of additional assimilation of K DP , ExpVrZhPol (Figure 5f) also has a similar good performance to ExpVrZhKdp. Based on the simulated result of ExpVrZh (Figure 6b) which is pretty close to the true ρ hv (Figure 6a), ExpVrZhRhv (Figure 6e) mainly brings improvement in details, whereas ExpVrZhZdr ( Figure 6c) and ExpVrZhKdp (Figure 6d) make further overestimation of ρ hv in certain areas. Affected by Z DR and K DP , ExpVrZhPol (Figure 6f) does not exhibit better ρ hv distribution, but instead has an underestimated area (<0.9).   Figure 3, but for the differential reflectivity Z DR (unit: dB) of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively.  Figure 3, but for the specific differential phase K DP (unit: • /km) of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively.  (Figure 7e). Additionally, Z H in these OSSEs are overestimated markedly in the middle level of the forward flank (5 km horizontally, the same hereafter), which is reduced properly in ExpVrZhKdp (Figure 7d). Additionally, the overestimation in the low level of the rear SPA (30~40 km) has been alleviated by assimilating Z DR , such that the intensity is much closer to the truth and the echo goes from ungrounded to grounded. ExpVrZhRhv shows no enhancement compared with ExpVrZh, or even worse results (overestimation in some areas, Figure 7e vs. Figures 7a and 7b). More than that, the ahead-developing convective cell (~6 km) and inner strong echo of the main storm are both better described in ExpVrZhKdp and ExpVrZhPol (Figure 7f). In comparison with the truth (Figure 8a), the signature of the Z DR column (~5 km) related to the updraft within the storm is more evidently discerned in ExpVrZdr (Figure 8c) and ExpVrKdp (Figure 8d), whereas the sagged high Z DR area (~25 km) attributed to the existence of a small number of large raindrops is well retained in ExpVrZh (Figure 8b), ExpVrZdr, and ExpVrRhv (Figure 8e). Modifications of amplitude and shape of Z DR can be noticed in the tail (~35 km) in ExpVrZhZdr and ExpVrZhKdp. With respect to above features, the best analysis result in terms of intensity and structure is given in ExpVrZhPol (Figure 8f). ExpVrZhKdp reveals an extremely high degree of similarity with the truth (Figure 9d vs. Figure 9a), including the K DP core associated with high liquid water content (LWC) within the storm, an overshooting structure linked to the strong updraft rushing through the melting layer, and even the small K DP area hanging ahead. The effect of K DP assimilation is maintained in ExpVrZhPol (Figure 9f), showing almost identical distribution characteristics to that of ExpVrZhKdp. Comparatively, other OSSEs (Figure 9b,c,e) nearly fail to reproduce these mentioned signatures. The melting layer associated with the melting hydrometeors is precisely described by a ρ hv drop in these OSSEs, but with a wider melting band and a few underestimated areas above. Other than a small positive contribution from the ρ hv assimilation (Figure 10e) in the low-value area located in the middle level of the forward edge (~5 km) of the storm, there are no notable improvements in the distribution of analyzed ρ hv for different OSSEs.   Figure 7, but for the differential reflectivity Z DR of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively. Figure 9. The same as Figure 7, but for the specific differential phase K DP of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively. Figure 10. The same as Figure 7, but for the cross-correlation coefficient ρ hv of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively. Figures 11-14 show the mixing ratios of cloud ice q i , rain water q r , hail q h , and graupel q g in the truth and analysis fields of OSSEs, respectively, through identical vertical slices as in Figure 7. These OSSEs with PRD assimilation partially (mainly in the right part) and slightly (within 0.1 g/kg) strengthen q i more than that in ExpVrZh (Figure 11b), which is still much weaker compared to the truth (Figure 11a). Besides the basic distribution of q r , well depicted in OSSEs, the outstretched structure corresponding to the convection initiation ahead (~5 km) is better analyzed in ExpVrZhZdr (Figure 12c) and ExpVrZhKdp (Figure 12d). The latter better portrays the rear (>30 km) relatively small q r (0.1~0.2 g/kg), while the former is slightly over-adjusted (>0.2 g/kg). The overshooting pattern (~10 km) resulting from the strong updraft within the storm and suggesting the presence of super-cooled water is clearly displayed in ExpVrZhKdp as well, completely consistent with the truth (Figure 12a). However, the more obvious q h (Figure 13c) core (>1.8 g/kg) and enhanced q g (Figure 14c) in the middle-(>1.5 g/kg) and upper level (>4.0 g/kg) can be seen in ExpVrZhZdr. As expected, these mentioned benefits induced by Z DR or K DP assimilation are also retained in ExpVrZhPol. In addition, there are many other improvements such as the preferable descriptions of relatively small q i (<0.2 g/kg), q h (<0.6 g/kg), and q g (<4.5 g/kg) in the upper left corner f, Figures 13f and 14f), the superior q r analysis much closer to the truth (Figure 12f vs. Figure 12a), and so on. However, apart from minor differences, all analyzed mixing ratios in ExpVrZhRhv are largely the same as the corresponding ones in ExpVrZh (Figures 11e, 12e, 13e and 14e vs. Figures 11b, 12b, 13b  and 14b, respectively). Additionally, analyses of cloud water q c and snow q s yield similar results as other hydrometeor variables, that is, q c is similar to q r , and q s is similar to other ice species (q i , q g , and q h ), which are not detailed here.

Evaluation of Hydrometeor Analysis
In general, the above results further demonstrate the potential positive effects of Z DR assimilation on the analysis of solid hydrometeor particles, more benefits of K DP assimilation in improving the simulation of rain water rather than ice hydrometeors, and the limited usefulness of ρ hv in terms of enhancing model hydrometeor analysis. Attributing to the combined positive impacts from PRD assimilations, more comprehensive hydrometeor analysis results in both intensity and pattern can be obtained. The conclusions are also in accordance with the findings of RMSE analysis. Figure 11. The same as Figure 7, but for the mixing ratio of cloud ice q i (unit: g/kg) of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively.  Figure 7, but for the mixing ratio of rain water q r (unit: g/kg) of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively. Figure 13. The same as Figure 7, but for the mixing ratio of hail q h (unit: g/kg) of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively. Figure 14. The same as Figure 7, but for the mixing ratio of graupel q g (unit: g/kg) of (a) truth, and analysis fields of (b) ExpVrZh, (c) ExpVrZhZdr, (d) ExpVrZhKdp, (e) ExpVrZhRhv, and (f) ExpVrZhPol, respectively.
Different PRD variables have different information contents and error effects corresponding to each hydrometeor. Z DR increases as the raindrop size increases or ice particles melt. Therefore, assimilation of Z DR yields improved analysis of all hydrometeors. K DP depends on both the raindrop size and the number concentration, hence, its assimilation improves the analysis of q r . Since ρ hv reduction occurs only in the melting/mixing regions which have a very small percentage of the storm and it has a small dynamic range with relatively large errors, the improvement with assimilation of ρ hv is minimal.

Evaluation of Forecast
To investigate the improvement effect of modified initial conditions from different DA schemes on model predictions, the 1 h forecasts beginning at the 120 min of a storm are made at 15 min intervals for these OSSEs. In addition, the mean RMSEs of forecasts are calculated as seen in Figure 2 and displayed in Figure 15. Unsurprisingly, the increasing trend of RMSE with time is exhibited in all model state variables and simulated PRD, but the variation characteristics for different OSSEs are still existing diversely.
Taking all parameters into account, ExpVrZhPol (purple line) generally outperforms other OSSEs. Although the differences in the horizontal wind of u and v (Figure 15a,b) are almost indistinguishable, it can be seen upon careful inspection that ExpVrZhPol has a relatively low RMSE on the whole. Furthermore, except for the better performance of ExpVrZhZdr (red line) in the intermediate stage (30-45 min) of prediction, ExpVrZhPol also has the best overall prediction effect on vertical velocity w (Figure 15c). An analogous situation appears in the perturbation potential temperature θ (Figure 15d) and perturbation pressure p (Figure 15e) as well. For hydrometeors, the remarkable advancements resulting from combined assimilation of PRD can be found in cloud water q c (Figure 15g), rain water q r (Figure 15i), snow q s (Figure 15j), and graupel q g (Figure 15l) almost throughout the entire prediction period, and mainly in both ends of the period for cloud ice q i (Figure 15h) and hail q h (Figure 15k). While taken individually, with respect to RMSEs, ExpVrZhZdr has a good performance in ice hydrometeors (q i , q h , and q g ) forecasts, suggesting the potential advantages of Z DR assimilation in increasing the prediction accuracy of ice species. Further, due to the high sensitivity of K DP to liquid water, ExpVrZhKdp (blue line) illustrates more superiorities in the forecasting of water vapor q v (Figure 15f), q c , and q r . Comparatively speaking, ExpVrZhRhv (cyan line) has fully comparable behaviors with ExpVrZh (black line), demonstrating the limited benefits of ρ hv assimilation. In observation space, ExpVrZhPol also manifests superior overall performance, thoroughly in Z H (Figure 15m), and mainly in the late forecast stage (after 30 min) of Z DR (Figure 15n), K DP (Figure 15o) and ρ hv (Figure 15p), which is of great significance for quantitative applications of PRD (precipitation estimation, hydrometeor classification, and so on). By contrast, ExpVrZhKdp is the second choice, having a suboptimal effect comprehensively. It is noteworthy that the RMSEs of ExpVrZhRhv nearly coincide with that of ExpVrZh for all PRD variables, further indicating the minor positive impact of ρ hv assimilation on the forecasts.

Summary and Conclusions
In this paper, based on the new simplified and parameterized polarimetric radar observation operators developed by Z21, the corresponding TL/AD operators necessarily required by the variational DA procedure are constructed. The gradient check of the cost function in the minimization step of analysis cycles proves the operators work reasonably. In fact, as mentioned previously, the updating of number concentrations which affect the DA effect [18] is not considered here and is worthy of further research.
To evaluate the influences of assimilating the differential reflectivity Z DR , specific differential phase K DP , cross-correlation coefficient ρ hv , as well as horizontal reflectivity Z H and radial velocity V r , cycled 3DVAR DA and forecast schemes are performed. The truth is generated with a 3 h simulation of an idealized supercell storm using the WRF model. The PRD is constructed from the idealized truth simulations via Z21 forward operators. Then, these observations are assimilated into the WRF model at a 5 min interval with the initial DA beginning at 30 min of model integration.
Using the experiment of assimilating both Z H and V r as the benchmark run, OSSEs with additional assimilation of Z DR , K DP , and ρ hv , respectively, are performed to examine the impact of each polarimetric variable on DA results. Additionally, an extra OSSE which assimilates all PRD is conducted to see whether a more balanced result can be achieved. The results suggest that the Z DR assimilation adjusts all model state variables to varying degrees, and helps reduce RMSEs at almost each DA cycle. Moreover, the assimilation process makes some polarimetric signatures, such as hook echo, Z DR column and melting layer, more obvious, and improves the distributions of hydrometeors and observations much closer to the truth as well. High sensitivity to liquid water of K DP is believed to be responsible for the more improved analysis of rainwater, and the intensity and pattern of PRD are also enhanced with additional assimilation of K DP , especially for Z H and K DP itself. The small discrepancy between OSSEs with and without ρ hv assimilation indicates its limited ability in adjusting model state variables and improving the analysis accuracy. In contrast, due to the accumulation of favorable DA effects of different polarimetric variables, assimilating all PRD obtains more balanced analysis results, as expected. The above polarimetric signatures are also well presented in observation space, and many other improvements have been made, including the optimization of rain water analysis which is perfectly matched with the reference truth, and better descriptions of some small value areas in mixing ratios of cloud water, snow, and graupel. These different assimilation results are attributed to the different information contents and error effects from PRD variables. It is noted that these characteristics of assimilation PRD can change depending on the type of storm and model microphysics as well as observation errors, which needs further study. A series of prediction experiments, which initiate at the 120 min of different OSSEs and run over a 1 h period at a 15 min interval, are conducted to further study the impact of improved initial conditions on ensuing forecasts. The RMSEs of forecasts are calculated and the results show that the Z DR assimilation can bring obvious benefits for forecasts of ice hydrometeor species including cloud ice, hail, and graupel and even for the vertical wind field and perturbation pressure in the late stage of prediction. Assimilating K DP is found to improve the prediction of water vapor, cloud water, and rainwater, which is consistent with its performance in the liquid water analysis. There is no visible promoting effect on the prediction of most variables, suggesting that ρ hv assimilation is not very helpful for the improvement of NWP. Similar to analysis results, the OSSE assimilating all PRD shows a more comprehensive performance in terms of RMSEs, having the most accurate (smallest RMSE) results of horizontal wind field, perturbation potential temperature, cloud water, rain water, snow, hail, and graupel in the end of forecast period. In addition, the forecasts in observation space are also better, with higher accuracy in the second half of period.
Although these findings illustrate the potential usefulness of PRD assimilation in improving convective-scale NWP, additional important research questions remain. A subsequent study will apply the PRD assimilation to real data cases to assess its benefits to real-world forecasts. Furthermore, the feasibility and capability of this new parameterized polarimetric radar operators under different MP schemes require investigation as well.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
All abbreviations and acronyms used in this paper are summarized briefly in Table A1.