Evaluation of Forward Operators for Polarimetric Radars Aiming for Data Assimilation

In the preparation for polarimetric radar data assimilation, it is essential to examine the accuracy of forward operators based on different formulations. For this purpose, four forward operators that focus on warm rain conditions are compared both with each other and actual observations with respect to their performance for C-band dual polarimetric radars. These operators mutually consider radar beam broadening and climatological beam bending. The first operator derives polarimetric parameters assuming an exponential raindrop size distribution obtained by the models and is based on fitting functions against scattering amplitudes. The other three converters estimate the mixing ratio of rainwater from the measured polarimetric parameters. The second converter uses both the horizontal reflectivity (ZH) and the differential reflectivity (ZDR), the third uses the specific differential phase (KDP), and the fourth uses both KDP and ZDR, respectively. Comparisons with modeled measurements show that the accuracy of the third converter is superior to the other two. Another evaluation with actual observations shows that the first converter has slightly higher fractions skill scores than the other three. Considering the attenuation effect, the fitting function and the operator only with KDP are found to be the most suitable for data assimilation at C-band.


Introduction
Assimilation of radar reflectivity data provides significant improvements in quantitative precipitation forecasts (QPFs) as shown, for instance, by Sun and Crook (1997), Kawabata et al. (2011), Wang et al. (2013), and Schwitalla and Wulfmeyer (2014).They directly assimilated the radar reflectivity factor at horizontal polarization Z H by using the so-called "Z-Qr" relationship (a power law function of Z H ) as an observational operator.This type of operator is used to estimate the rainwater content Q rain (g m −3 ), which is predicted by numerical models, from the measured Z H .
Because their measurements were mostly made with conventional single-polarization radars, which have been used in operational radar networks so far, the operators were formulated with a single input parameter (Z H ).
Recently, dual-polarization radars have been deployed for operational purposes.For instance, the operational radar network in the United States was upgraded with polarimetric S-band radars (e.g., Seo et al. 2015).In France, polarimetric radars in the S-, C-, and X-bands were installed as parts of the operational network (Figueras i Ventura et al. 2012).The German Weather Service (DWD) designed a new radar network composed of polarimetric C-band radars (Helmert et al. 2014).Moreover, an operational X-band polarimetric radar network has been deployed in Japan since 2010 (Maesaka et al. 2011).Dual-polarization radars have advantages over conventional ones on quantitative precipitation estimates (QPEs) because of their superior reliability in estimating the rainfall rate (Bringi and Chandrasekar 2001).
Polarimetric parameters measured with dual-polarization radars include the differential reflectivity Z DR .The utilization of Z DR in combination with the radar reflectivity factor Z H helps to estimate rainwater content (e.g., Jameson 1991;Anagnostou et al. 2008), especially for S-band radars.However, the algorithms that rely on Z H and Z DR are less efficient at both the Cand X-bands because Z DR can be strongly biased by differential attenuation, especially in heavy rainfall (Ryzhkov et al. 2014).A method based on the specific differential phase K DP has also been used extensively in many studies (e.g., Jameson and Caylor 1994;Ryzhkov and Zrnić 1995;Kim et al. 2010) because it is immune to attenuation and radar miscalibration (Sachidananda and Zrnić 1986), even for radars operating at either the C-or X-bands.However, K DP is less reliable at lower rainfall rates because it becomes noisy and increasingly susceptible to the variability of the raindrop size distribution (DSD; Chandrasekar et al. 1990).Most of the aforementioned methods assume a priori power law form relations between polarimetric parameters and rainwater content with coefficients determined by fitting radar observables to rain parameter observations on the ground.In contrast, Adachi et al. (2015) proposed a method to retrieve rain parameters, including the DSD and rainwater content, from polarimetric radar measurements based on scattering theory without empirical fitting functions.
On the other hand, many operators for numerical models, which calculate the scattering amplitudes directly from the DSD retrieved by the model, have been proposed using Mie scattering theory, the Rayleigh approximation (e.g., Ryzhkov et al. 2011), and the T-matrix method (e.g., Pfeifer et al. 2008;Jung et al. 2010;Augros et al. 2015).Because the T-matrix approach (e.g., Vivekanandan et al. 1991;Mishchenko et al. 1996) is capable of calculating the scattering matrix for spheroidal particles such as raindrops, this method is superior to Mie theory, which is valid only for spherical particles, although the computational cost of the T-matrix is much larger.
Some of these operators have been applied to QPF with data assimilation (DA).Jung et al. (2008a, b) adopted fitting functions with the T-matrix method in warm rain conditions (Zhang et al. 2001) and the Rayleigh approximation for ice particles to calculate polarimetric parameters at the S-band.They implemented this method in the Advanced Regional Prediction System (ARPS) (Xue et al. 2000) and applied it to an ensemble Kalman filter.Their observation-system-simulation experiment showed good agreement with a perfect model experiment.Li andMecikalski (2012, 2013) implemented power laws with Q rain by using Z H , K DP , and Z DR for the three-dimensional variational system (3D-Var) of the Weather Research and Forecasting Model (WRF: Skamarock et al. 2008;WRF 3D-Var: Barker et al. 2012).Their studies focused on warm rain at the C-band and showed improvements in short-range forecasts in real cases.In addition, Kawabata et al. (2011) and Sun et al. (2013) showed positive impact on short-range forecasts using 4D-and 3D-Var systems with warm rain schemes by assimilating only Z H .These studies were done without any considerations of rainfall types like convective or stratiform rain.
As mentioned above, many convertors with polarimetric parameters have been proposed.These convertors are classified into two types of operators: the first calculates polarimetric parameters from rain parameters, including the DSD in numerical models, while the second derives rain parameters from observed polarimetric parameters.It is therefore necessary to evaluate their performance so that the best choice can be made for implementation in a DA system.In this study, we implemented four methods of the two types (a fitting function and three power law methods).Moreover, because weather radar networks are operated at the S-, C-, and X-bands, the coefficients used for the three power laws in our methods consider all three frequencies (Bringi and Chandrasekar 2001).Furthermore, a fitting function originally proposed by Zhang et al. (2001) for the S-band was applied and expanded to the C-band in the present study.This means that we focus on C-band radars in the present study, because of available observational data to evaluate the proposed methods.However, our methods are also applicable for S-and X-bands by replacing the coefficients in formulations implemented in the present study.
We first evaluated the accuracy of the four operators by comparing our numerical model results with radar observations at the C-band.This evaluation is the first step to utilize the best observation operator for polarimetric radars in DA systems.The final goal of this study includes an implementation of the operator into the WRF Var system, which has been operated as a rapid-update cycle system to evaluate precipitation forecasts over Germany by one of the present authors (Schwitalla et al. 2011).Because the WRF Var system is limited to assimilating observations relating to rainwater and is unable to consider ice particles, we focused on warm rain.There have been many studies that have shown significant improvements of rainfall forecasts with direct assimilation of warm rain, as mentioned above.In addition to these studies, we use more accurate measurements of rain with polarimetric radars.
Other motivations for this study include improving the quality of both QPE and QPF by using numerical models.High-quality QPE data are essential to understand the whole water cycle in a certain area, as studied in the Catchments As Organized System (CAOS) project (Zehe et al. 2014;Bauer et al. 2015; http:// www.caos-project.de/).Moreover, QPF with high reliability is indispensable for understanding the mechanisms that organize, develop, and decay mesoscale convective systems.Therefore, the present study can also contribute to research on small-scale heavy convective rainfalls, such as the Tokyo Metropolitan Area Convection Study for Extreme Weather Resilient Cities (TOMACS) project (Nakatani et al. 2015).This direction should lead us toward more effective disaster prevention related to convective heavy rainfall events.As a first step toward achieving this goal, we developed forward operators for polarimetric radars and evaluated them in this study.
This paper is organized as follows: In Section 2, forward operators for polarimetric radars are described.Simulations using the WRF model and case descriptions are given in Section 3. Evaluations are presented in Section 4. Our discussions and conclusions are given in Sections 5 and 6, respectively.The theoretical background of polarimetric parameters is shown in the appendix.

Space interpolator
Like many studies of DA (e.g., Seko et al. 2004;Zhang et al. 2009), radar observations were averaged and then interpolated to the model grid points along the radar beam in radial directions, because their resolution was finer than that of the model and highfrequency signals were eliminated.This interpolated observation, so called "super observations", remained original azimuthal and vertical resolutions.While corresponding modeled values were interpolated to the same points under considering Gaussian broadening of the radar beam and the vertical gradient of the refractive index calculated with the International Standard Atmosphere (Doviak and Zrnić 1993).Furthermore, the modeled values were derived with valuable convertors described in Section 2.2.In all evaluations in this study, the super observations and interpolated modeled values were used.This means that the observations and simulations were compared on the plan position indicators (PPIs), but the radial resolutions were modified to that of the model grid spacing.In addition, beam blocking areas were defined using highresolution topography data.

Variable converters
A variable converter that estimates polarimetric parameters such as Z H (mm 6 m −3 ), vertical reflectivity Z V , Z DR (dB), K DP (° km −1 ), specific attenuation A h (dB km −1 ), and specific differential attenuation A DP (dB km −1 ) from the DSD in the model (see appendix for details of the polarimetric parameters) was developed.Simply speaking, this converter uses fitting functions against the scattering amplitude calculated by the T-matrix (Zhang et al. 2001;hereafter FIT).Note that the radar reflectivity factor Z H is converted to conventional reflectivity Z h (dBZ) by Eq. (B2).Details are given in the next subsection.The parameters A h and A DP estimated with the model were used to simulate the attenuation effects in Z h and Z DR in the model, respectively.We also utilized three other converters that estimate rainwater content from observed polarimetric parameters as described below.
a. Calculation of polarimetric parameters with the modeled rain parameters Although cloud microphysics schemes assume that raindrops in the models are spherical, we first introduce an axis ratio of a raindrop r given by Brandes et al. (2002Brandes et al. ( , 2005) ) to the variable converters as follows: where D is the equivolume diameter (mm).
FIT is based on fitting functions.Following Zhang et al. (2001), the backscattering amplitude is represented by a power law function as follows: where subscript h and v represent horizontal and vertical polarization, and α h, v and β h, v are coefficients that are determined by fitting the values of equivolume diameters D to those of the scattering amplitudes, respectively.| S h, v | is calculated with the T-matrix method (Fig. 1).The difference between the forward scattering amplitudes is represented as where f n (D) and f v (D) represent forward scattering amplitudes, and α k and β k are determined by fittings.
We made these three fittings using the T-matrix and Eqs.
(1) -( 3) and determined the coefficients shown in Table 1.Furthermore, based on the assumption of Eq.
(2), Zhang et al. (2001) derived the following equations for Z H, V and K DP from N 0 and the slope parameter Λ as follows: and where the subscripts H and V represent horizontal and vertical polarization, and K w is a constant defined as K w = (ε -1)/(ε + 2) in which ε is the complex dielec- tric constant of water.
The parameters A H and A DP are calculated from K DP as where α H , β H , α d , and β d are coefficients that follow Bringi and Chandrasekar (2001), which are also listed in Table 1.Finally, Z H and Z DR are affected by the attenuation effect A H and A DP , respectively (Eqs.B7, 8).
The fitted curvature for S v (D) is larger than that for S h (D) in case D is less than 1 mm (Fig. 1), although this rarely happens for warm rain in nature.In fact, values of S v obtained with the T-matrix are almost identical to those of S h in that range in the figure.Therefore, we simply replace α v and β v with α h and β h of Eq. (2) in the range less than or equal to FIT calculates the polarimetric parameters by using the exponential DSD (Eq.A1 at μ = 0).It was used in Jung et al. (2008a, b) for DA with S-band radars, as mentioned in the introduction.Although Jung et al. (2008a, b) used an axis ratio proposed by Green (1975) and applied it to S-band radar data, we applied the method to C-band with another axis ratio by Brandes et al. (2002Brandes et al. ( , 2005) ) and the WRF model under warm rain conditions (no ice hydrometers).As for the axis ratio, the model proposed by Brandes et al. (2002,   10 -6      Green (1975), as shown in Thurai and Bringi (2005).In addition, when new fittings at Sand X-bands are calculated, FIT can be used for these frequencies, as well.
b. Calculation of rainwater from observed polarimetric parameters In addition to FIT, described above, which calculates polarimetric parameters from the DSD in the model, we also used three other existing converters that estimate the rainwater content Q rain (g m −3 ) from the observed polarimetric parameters (Bringi and Chandrasekar 2001).These convertors are the socalled "Z-Qr" relationships: the first one estimates Q rain from Z H and Z DR , represented by and the second estimates it from K DP as where f is the frequency of the radars in GHz.The third estimates Q rain from K DP and Z DR as a 1, 3 , b 1, 2, 3 , and c 1, 2, 3 are the empirical constants described in Bringi and Chandrasekar (2001; Table 2 in this study).We used coefficients for C-band in this study, however, since Bringi and Chandrasekar (2001) proposed coefficients for S-and X-bands as well, Eqs. ( 8) -( 10) are applicable for other frequencies.This is the most important for our choice among a lot of estimators proposed in other studies, because it allows for assimilating radar data at all three frequencies with the same operator only by replacing the coefficients.
As far as we know, these equations are the one and only candidates for this purpose at present.Equations ( 8) and (10) were utilized in Li andMecikalski (2012, 2013), and Eq. ( 9) was applied in Yokota et al. (2016) for DA experiments with actual observations.As described above, FIT emulates polarimetric radars and Q(Z H , Z DR ), Q(K DP ), and Q(K DP , Z DR ) are estimators of the actual rainwater content.However, in this paper, we simply call all four together the "forward operator", because this term is used for "observation operator" in the cost function of variational methods, which contrasts with the "adjoint (backward) operator", which will be developed in a future study.Meanwhile, Q rain in the model is calculated from the mixing ratio of rainwater (kg kg −1 ) and air density (kg m −3 ).

Case 1: 14 August 2014
On 14 August 2014, many convective systems that had been initiated over Germany moved eastward and were distributed sporadically during daytime.Figure 2a shows the composite reflectivity map using the PPI data at the lowest elevation of each C-band radar data, which may include data above the melting layer.Because a large number of convective systems were observed, this case is useful for statistical evaluation.Numerical simulations were conducted with the WRF model in this study.In the WRF model, we used two-moment (Morrison et al. 2009; see Appendix A) instead of one-moment microphysics schemes for better representations of clouds and precipitation because the former predicts the number density and mixing ratio of water and ice particles explicitly, reducing the number of assumptions in the calculation of the parameters for the gamma distribution of raindrop size.
First, a simulation with 5-km grid spacing was done over Western Europe with initial conditions provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) analysis with 16-km grid spacing starting at 0600 UTC on 14 August 2014.Then, the six hourly analysis data were used to drive lateral boundary conditions until 1800 UTC to reduce forecast errors throughout the entire forecast period.In addition to the Morrison scheme, the Kain-Fritsch convective cumulus parameterization (Kain 2004) was applied in this simulation.After that, a downscaling simulation with 1.67-km grid spacing and 57 layers in vertical up to 50 hPa was conducted without any cumulus parameterizations from 0900 to 1800 UTC.The domains for the simulation are shown in Fig. 3a.Other important physical schemes used in these simulations are the NOAH-MP (Niu et al. 2011) as a land-surface model and the YSU planetary boundary scheme by Hong et al. (2006).
The simulation results with 1.67-km grid spacing were compared with the DWD operational radar observations at Neuheilenbach, Offenthal, Feldberg, and Tuerkheim (circles in Fig. 3a) from 1000 to 1800 UTC.The wavelength was 5.3 cm (C-band), the antenna gain was around 45 dB, and the maximum range was approximately 180 km.The original radar observations were made with a range resolution of 125 m and a beam width of 1°.However, data that were averaged over eight gates (1000 m) along the radial were used in the experiments.PPIs at elevation angles of 0.5°, 1.5°, 2.5°, and 3.5° within the maximum ranges were used in this study.We obtained the observational data of Z h , Z v , and ϕ DP without attenuation correction after operational quality controls such as eliminating ground clutters by DWD (personal communication).K DP was calculated from ϕ DP as follows: (1) ϕ DP with the value of Z h more than 1 dBZ were selected; (2) negative ϕ DP were rejected; (3) running means of ϕ DP with five bins were calculated, and then removed outliers by comparing with each observation and the mean; (4) afterwards, the running means were re-applied, and gradients were calculated; (5) finally, the gradients with positive value were stored as K DP data.
In this evaluation, data were limited to altitudes lower than those of the freezing level.These freezing levels were defined on the PPIs and determined with the model prediction data varying in time and space.It might be problematic that observations just below the freezing level may include melting ice particles.The observed melting layer and any errors in the altitude of the freezing level by the model should be considered as observational errors in the DA.
As further quality controls for observations, data were confined to the ranges from 0 to 70 dBZ for reflectivity, from 0 to 4 dB for Z DR , and from 0 to 5° km −1 for K DP .Regarding K DP , two additional quality controls were made to avoid contamination from low-quality observations; (1) difference between modelled and observed K DP should be less than 3.0° km −1 ; (2) difference of rain water contents calculated with K DP using Eq. ( 9) and with Z H and Z DR using Eq. ( 8) should be less than 0.3 g kg −1 .These thresholds were determined after an initial investigation using the data set and methods described in this study.

Case 2: 30 June 2012
Because it is expected that the characteristics of rainwater (e.g., the drop size distribution) in small and medium-size cloud systems are different than those in strong thunderstorms, and there was no strong thunderstorm in Case 1 and even in the entire summer, 2014, we conducted another simulation with 9-km horizontal grid spacing using the ECMWF analysis with 16-km grid spacing from 0000 to 1800 UTC on 30 June 2012.We then downscaled it to 3-km grid spacing and further to 1-km grid spacing (Fig. 3b) in the case of a supercell (see the conventional radar observation in Fig. 2b).However, because there were no polarimetric observations available in this case, we decided to make model-to-model comparisons.
For this purpose, a virtual C-band radar was put near the thunderstorm at 48.9°N, 9°E.The configuration of this radar was the same as for the DWD radars, except for its location, which was determined under consideration of intense rain regions in the supercell.Then simulated radar data was created with FIT using the WRF prediction described above.Details of this comparison are given in the next section.Settings of both above-mentioned simulations were determined under consideration of possible DA systems in our future plans.

Evaluations
We evaluated the four operators described above with actual observations, simulations, and each other.First, simulated rainwater in the model was converted to polarimetric parameters with FIT and compared with the actual observations (OBS) in Case 1 for evaluation (see Section 4.1).This was done to see errors of FIT, including representative errors of the model simulation, i.e., errors of locations, cloud microphysics, ice phase in convective systems, and the inaccuracy of the forward operator.Observations of all four radars were used in this evaluation.
Second, for "Z-Qr relationships", we used FIT as a radar simulator to provide polarimetric parameters for the "Z-Qr" relationships (Section 4.2).Schematically speaking, input data (hereafter A) are converted to output data (hereafter C) by a convertor (hereafter B).When we evaluate the accuracy of the convertor, the output data C are compared with other data (D).In our evaluation, first we converted the modeled rainwater content (A) by FIT (B) to the polarimetric parameters (C).Then C data are re-converted by the "Z-Qr" relationships (B¢) to the rainwater content (A¢).This convert-inverse process provides some differences between A and A¢.Therefore, we are able to estimate the inaccuracies of these "Z-Qr" relationships, although they include the common error of FIT.This method enables us to compare the converted rainwater with the simulated rainwater (TRUE) directly, without any observations.Simulation data both in Case 1 and 2 were used in this evaluation.
Finally, FIT and the "Z-Qr relationships" were compared with actual observations in Case 1 using fractions skill scores (FSSs; Section 4.3).In general, it is difficult to make a direct comparison between different variables: polarimetric parameters and rainwater, which are derived with two types of converters (model variables to observations and observations to model variables).For this comparison, we applied the FSS, which is used for verification of rainfall forecasts with certain criteria and can help us to understand the difference in accuracy between convertors.Although the FSS evaluates only the accuracy of the predicted rainfall area against observations, it nevertheless illustrates the accuracy for each criterion when we carefully choose the criteria for the scores (bins of precipitation amounts).This is based on the idea, which, in cases where the numbers of cases at the certain criterion for each operator are the same or similar, the probability density of the events for nature (all cases) are also the same or similar.Therefore, the FSS is able to measure skills of different variables using the criterion.Furthermore, "grid-to-grid" scores (e.g., threat score, standard deviation) only capture forecast skills when the exact locations of the observations and forecasts coincide.However, when the forecasts suffer from displacement error, the "grid-to-grid" scores are penalized by the so-called "double penalty" problem.To avoid this problem, the FSS was suggested by Roberts and Lean (2008).First, they defined a certain size of verification domains over a forecast field, and then counted the number of grids in the domain that exceeded the criteria for verification.Dividing that number by the total number of the grid mesh in the domain, the percentiles exceeding the criteria were obtained (P fcst and P obs in Eq. 11).The FSS is defined as where n represents the total number of grids in the verification domains, and P fcst and P obs are the forecast and observed percentiles in the domain, respectively (see Fig. 8a as a reference).The advantages of the FSS include identifying spatial fluctuations in forecasts and evaluating valid forecasts that depend on a horizontal scale by varying the size of the verification domains.For instance, Duc et al. (2013) extended the FSS to the time dimension and clarified that forecast skills against small-scale precipitation systems were related to short time scales.

Evaluation of FIT with radar observations in
Case 1 The simulation captured important features of the convective systems, for instance, two major systems at almost the same locations and intensities as the observations: a rain band in the north-south direction existed to the east of the radar site, a convective system in the northwestern region of the radar obser-vation, and sporadic precipitation embedded in these systems (upper-left in Fig. 4).However, the simulation (FIT) produced more widely distributed Z h than OBS.This means that the simulation missed small-scale structures that appeared in the observations.One of the causes of these distributions being so widespread is the model grid spacing: numerical models can explicitly represent phenomena only at resolutions coarser than the size of at least three times that of the model grid cells.Furthermore, this can be affected by low-quality initial conditions due to poor resolution and fewer observations, such as radar reflectivity.
Z DR in FIT shows wider areas but more moderate values than in OBS.K DP in FIT captured smaller areas and also more moderate values than in OBS (middle and right columns of Fig. 4).
Next, we focus on the statistical analysis of 1000 -1800 UTC in Case 1, which included 30,000 -100,000 points.Quality controls for each parameter determine the number of the points.The averages and standard deviations of Z h calculated by FIT and OBS are similar, and the difference between FIT and OBS is quite small (Table 3).Moreover, the standard deviation of the difference between FIT and OBS is 15.65 dBZ, which is similar to other studies (e.g., 15 dBZ by Kawabata et al. 2011).Regarding Z DR and K DP , there are negative biases and small standard deviations, which cause relatively large standard deviations in FIT minus OBS.A similar impression may be obtained from Fig. 4, but the difference in K DP is notably larger.
Frequencies of occurrence with value classes for each parameter for both OBS and FIT in Case 1 are shown in Fig. 5.The lowest and highest boundaries, and the differences of each bin range of these graphs, are chosen to represent the statistical features of the data.Z h by FIT and OBS shows similar distributions over all the classes, although the lowest class by OBS shows a much larger value, and the classes from 11 -41 dBZ by FIT show slightly larger values than OBS.The distributions of Z DR by FIT and OBS show similar numbers over all the classes, but the value of the lowest class in OBS is smaller than that in FIT.In contrast, the distribution of K DP by FIT concentrates in the lowest class, but that by OBS spreads from 0 to 3.3° km −1 .This results in a negative bias and small standard deviation in FIT.It is not clear whether this difference is caused by errors in FIT or OBS, because K DP is a differential value from the differential phase ϕ DP and is therefore well known to be difficult to cal- culate.

Comparison of
Z DR ) As shown in Fig. 6 with respect to the Offenthal radar, Q(Z H_FIT , Z DR_FIT ) and Q(K DP_FIT , Z DR_FIT ) appear to be in fairly good agreement compared to TRUE in terms of distributions, but slightly stronger Q rain values than TRUE are seen in intense rainfall regions.Q(K DP_FIT ) shows good agreement with TRUE in terms of both distributions and intensities.Note that, in this section, we only discuss simulation data as described in Section 3.3.The subscript "_FIT" represents a parameter provided by the FIT method.Since the locations of each rain region are mostly the same in Figs.6a -d, differences in intensity are important in this evaluation, because sensitivities for each polarimetric parameter appear in different locations, even in a convective system.
Statistical analysis with all four radar data in Case 1 and the virtual radar in Case 2 also shows that the averages and standard deviations of TRUE and Q(K DP_FIT ) are similar, while Q(Z H_FIT , Z DR_FIT ) and Q(K DP_FIT , Z DR_FIT ) are larger by a factor of from 1.5 to 2.0 (Table 4).Furthermore, the difference between TRUE and Q(K DP_FIT ) is much smaller than that of TRUE against Q(Z H_FIT , Z DR_FIT ) and Q(K DP_FIT , Z DR_FIT ).Thus, it can be concluded that Q(K DP_FIT ) has a smaller error than those of the others in the formulations.
Given the formulations of the polarimetric parameters (Eqs.B7, 8), Q(Z H_FIT , Z DR_FIT ) and Q(K DP_FIT , Z DR_ FIT ) are affected by attenuation, but Q(K DP_FIT ) is not.Thus, in Case 2, it is clear that more attenuation occurred in the northeastern area of the supercell, which is behind strong intensities along the radar beams (Figs.7a, c), compared with Q(K DP_FIT ) and TRUE.The attenuation could be problematic, although it was not confirmed by the statistical analysis because it happened only behind regions of wide and/or intense rain, and the regions were relatively narrow.0.1-0.5 0.5-0.90.9-1.31.3-1.7 1.7-2.1 2.1-2.5 2.5-2.9 2.9-3.(a) Fig. 6.Same as Fig. 4 but for rainwater (g m -3 ).

Fractions skill score
In this section, we compared observed polarimetric parameters with FIT, and rainwater converted by "Z-Qr relationships" from observed polarimetric parameters with simulated rainwater for Case 1. Observations of all four radars were used in this evaluation.FIT (Fig. 8b) shows a clear dependency on horizontal scale and intensity: the scores of the lowest intensity at 2 dBZ decrease from 0.49 to 0.35 as the size of the verification domain decreases; likewise, scores for the 5-km verification domain decrease from 0.35 to 0.14 as intensity increases.This is quite reasonable when we consider the relationship between horizontal scale and intensity of convective systems.For instance, Wang et al. (2013) showed that the FSSs by WRF without DA were 0.0 -0.3 for 5 mm h −1 and 25 km radius, while our FSS with downscaled initial conditions is 0.19 for 25 dBZ and a verification size of 15 km.It can be said that our simulation has similar FSSs with that.Namely, as an observational fact, strong echo regions of a thunderstorm appear only in narrow areas, while weak echoes spread to wide areas.This verification illustrates that the numerical simulation describes natural convective activity well.However, it is also clear that the highest score of 0.49 is not high, because, as pointed out by Bryan et al. (2003), the 1.67 -2.0 km resolutions is still insufficient for representing each convective system.In addition, we can point out that the quality of the initial conditions is quite low because they do not include any information on hydrometers at high resolutions.
The FSSs in Q(Z H , Z DR ), Q(K DP ), and Q(K DP , Z DR ) are close to that of FIT in the weakest intensity region, while they decrease rapidly with increasing intensity.However, we can still recognize the scale dependency in the results.Comparing these three, there is no large difference.It may be noted that Q(K DP ) is worst in the regions over 0.3 g m −3 , but the values are too small to be meaningful.In conclusion, FIT is better than the "Z-Qr" relationships in moderate and strong intensity regions.

Robustness of fitting coefficients
To confirm the robustness of the fitting coefficients in FIT (Eqs.4, 5), we varied α h, v, k and β h, v, k (Table 1) with errors between −15 % and 15 %.Equation (2) with varied α h and β h are shown in Fig. 9.It is clear that Eq. ( 2) is robust to varied α h but sensitive to varied β h ; moreover, we obtained similar results for Eq. ( 2) with varied α v and β v , and Eq. ( 3) with varied α k and β k (not shown).Using these varied coefficients, we calculated Z h with varied α h and β h , Z DR with varied α v and β v , and K DP with varied α k and β k (Fig. 10).Note that to avoid vertical reflectivity larger than horizontal reflectivity, we calculated Z DR with negative errors in α v and β v .
Similar to Fig. 9, varied αh affected both the stan- dard deviation and bias of Z h more weakly than did varied β h (Fig. 10a).However, because the errors in the standard deviation are between −3 % to 3 % as compared with −15 % to 15 % errors of the varied αh and β h and relatively large biases are seen, it can be said that Eq. 2 is robust with regard to the fitting accuracy but slightly sensitive to the bias.On the other hand, errors in Z DR reach 70 % when β v has −15 % error.For K DP (Eq.5), the biases and the standard deviation error for varied α k did not vary, but the standard deviation error for β k showed slightly larger values, from −8 % to 8 %.

Discussion
For the three "Z-Qr" relationships, Q(K DP ) showed better accuracy than Q(Z H , Z DR ) and Q(K DP , Z DR ) (Table  4).In addition, correlations of TRUE against Z h , K DP , and Z DR were 0.76, 0.87, and 0.38, respectively.These polarimetric parameters were provided by FIT.This also implies that K DP or a combination of Z h and K DP are superior to other combinations.Meanwhile, Fig. 7 reminds one of the necessity to consider the attenuation effect in the model, even though it does not appear in the statistical analysis.The FSSs showed that these three converters have similar skills over all intensities and horizontal scales.When our methods  are applied to the X-band, since this frequency is quite easily attenuated, it would be better to use only K DP in the "Z-Qr" relationships for DA.
The results of the FSS suggest that it is better for DA to adopt converters from models to observations.However, it is noteworthy that Gaussianity and linearity are also important for DA.The Gaussianity of both methods (FIT and the "Z-Qr" relationships) are quite low (not shown), which is often seen for prognostic variables related to water substances (e.g., Fig. 3 in Li and Mecikalski 2010 and Fig. 1 in Kawabata et al. 2011).Concerning linearity, FIT includes nonlinearity to the model prognostic variables in its formulations.For instance, Eq. ( 4) includes a power law with Λ, which is derived from the number density and mixing ratio of rainwater in the model.The "Z-Qr" relationships are completely linear (Eqs.8 -10) because radar observations are converted to the model prognostic variable of rainwater, and there is no sensitivity to the model variable of Q r .Moreover, linearity is related to cloud microphysics: the Morrison two-moment scheme handles the number density and mixing ratio of rainwater.FIT coincides with both variables, while the "Z-Qr" relationships are only connected to the mixing ratio.Increasing the number of prognostic variables in cloud microphysics may provide larger nonlinearity, while larger impacts on DA are expected.WRF-VAR considers perturbations to the mixing ratio of rain and cloud water but not to the number density.This means that assimilating information related to the number density is treated in the cost function but not in its gradient.This indirect utilization can lead to a smaller impact of the FIT assimilation.To understand which of the Q(K DP ) or FIT methods is better, it is necessary in the future to investigate impacts on DA results.
In Fig. 5c, the modelled K DP derived with FIT only appears in the lowest bin, while the observed K DP does in all bins.This is the same feature as in Bringi and Chandrasekar (2001;Fig. 7.85); They showed that observed and simulated K DP in connection with reflectivity and the observed K DP is widely distributed, while the simulated K DP appeared in very narrow ranges.This would confirm that the observed K DP represents rich diversity in nature, such as existence of various size of raindrops in an observation bin, while the simulated one does less.Therefore, assimilating widely distributed K DP would result in a rich diversity in the modeled rainwater.This is one of the reasons why we have to try assimilating polarimetric parameters.

Summary and conclusion
To utilize polarimetric radar observations for DA, four forward operators for polarimetric radars were implemented, and their accuracies for C-band dual polarimetric radars were examined.First, a space interpolator was developed that considers the beam broadening and the climatological bending effect.Then, two types of variable converters were implemented: one converts rainwater contents Q rain and number density in numerical models to polarimetric parameters (reflectivity Z h , differential reflectivity Z DR , and specific differential phase K DP ), and the others convert observed polarimetric parameters to Q rain , the so-called "Z-Qr" relationships.The former converter was fitting functions with the T-matrix method (FIT), and the latter three, Z h and Z DR (Q(Z H , Z DR )), K DP (Q(K DP )), and K DP and Z DR (Q(K DP , Z DR )), were connected to Q rain .
Two numerical simulations with the WRF model were conducted: a sporadic rainfall case on 14 August 2014 (Case 1) and a supercell case on 30 June 2012 (Case 2).Actual polarimetric observations were available only for Case 1.
FIT showed good agreement with observations for Z h and Z DR , but not for K DP .This disagreement should be investigated for both fitting and observational accuracies.We used FIT as a radar simulator to provide polarimetric parameters from the model rainwater to the "Z-Qr" relationships for evaluation.
A statistical comparison between the "Z-Qr" relationships showed that Q(K DP ) had the highest accuracy.The reason may be that correlations of the modeled rainwater against Z h and K DP were higher than that against Z DR .Although an examination with FSS showed no difference among the three applied operators, because Q(K DP ) was not affected by the attenuation, especially in the case of a strong thunderstorm (Case 2), it is suggested that Q(K DP ) is the best among these "Z-Qr" relationships for DA of warm rain.Meanwhile, we installed sophisticated quality controls for calculating K DP (Section 3.1).However, the observed K DP was quite different from the modeled one (Figs. 4,5).Thus, it might be necessary to consider other K DP estimations for each of observations and FIT.
An investigation on robustness of the fitting coefficients was conducted.The results showed that Z h is robust to errors in the coefficients, K DP is slightly affected, and Z DR has quite large sensitivity.
A comparison of two types of converters demonstrated that FIT was slightly better than the "Z-Qr" relationships; however, because the difference was not very large, it is difficult to determine which one is better at this stage.In addition, when considering Gaussianity, linearity, and the prognostic variables in the model, further investigations using tangent linear models and assimilation results are necessary.Furthermore, we have to investigate the impacts on rainfall forecast with Q(K DP ) and FIT through DA experiments.In this next step, the adjoint code of FIT will be developed and be assimilated polarimetric parameters directly, while the rainwater converted by Q(K DP ) will be assimilated as a retrieval method.These are issues for the second stage of our study.The issues may also include the dependency of the proposed methods on radio wave frequencies; X-and S-bands, because sensitivities to rain intensity of each polarimetric parameter observed at different frequencies would differ from the ones at C-band discussed in the present study.The specific differential phase (Smyth and Illingworth 1998) K DP (° km −1 ), the specific attenuations for horizontal polarization, A H (dB km −1 ), and for vertical polarization, A V , are defined (Oguchi 1973) as where Z h (x) and Z DR (x) represent the unattenuated reflectivity and the differential reflectivity at a range x (km), respectively.

Fig. 1 .
Fig. 1.Scattering amplitude against rainwater size and the fitting functions at C-band.Open circles and triangles represent backscattering amplitudes calculated by the T-matrix, and the bold solid and bold dashed lines are fitted with power-law functions (Eq.2).Cross marks represent the difference between the horizontally and vertically polarized forward scattering amplitudes, with the fitting function (Eq. 3) shown by the dotted line.

Fig. 3 .
Fig. 3. Simulation domains for (a) Case 1 on 14 August 2014 and (b) Case 2 on 30 June 2012.The smallest boxes in (a) and (b) represent the domains of the finest resolution simulations that were used in this study.The four black circles in (a) show the observation ranges of the DWD radars in this study.

Fig. 4 .
Fig. 4. Z h (left), Z DR (middle), and K DP (right) on a PPI at 0.5° for (a) observed and (b) simulated data of the Offenthal radar at 1700 UTC on 14 August 2014.White areas show the shadow regions due to topography or that were out of observational range.

Fig. 8 .
Fig. 8. Fractions skill score (FSS).A schematic of the FSS calculation presenting displacement errors in precipitation systems is shown in (a).In "grid-to-grid" methods, the forecast has no skill in this case.While in the FSS with verification grids defined as 5×5, both percentiles of the radar observation and forecast are 6/25.The FSS in this case is 1.0.Calculated FSSs in FIT, Q(Z H , Z DR ), Q(K DP ), and Q(K DP , Z DR ) are shown in (b), (c), (d), and (e), respectively.Vertical and horizontal axes of (b)-(e) illustrate horizontal scales and criteria for the calculations, respectively.

Fig. 10 .
Fig. 10.Standard deviations of FIT minus OBSs (same with Table3).(a) Z h for various α h and β h with errors between −15 % and 15 %, (b) Z DR for various α v and β v with errors between −15 % and 0 %, (c) K DP for vari- ous α k and β k with errors between −15 % and 15 %.

A
, v (m) are the forward scattering amplitudes.In addition, the specific differential attenuation A DP (dB km −1 ) is defined as H and A DP are used to estimate the attenuated Z h att and Z DR att , respectively, from the model as

Table 3 .
Averages and standard deviations of Z h , Z DR , and K DP according to FIT and FIT minus OBS in Case 1.

Table 4 .
Averages and standard deviations of TRUE, Q(Z H_FIT , Z DR_FIT ), Q(K DP_FIT ), and Q(K DP_FIT , Z DP_FIT ), and their differences with TRUE.The values are in kg m −3 .Same as Fig.6but for Case 2. The virtual radar is located at 48.9°N, 9°E.