Evaluating the Performance of a Convection-Permitting Model by Using Dual-Polarimetric Radar Parameters: Case Study of SoWMEX IOP8

: In this study, a dual-polarimetric radar observation operator is established and modiﬁed for the Taiwan area for the purpose of model veriﬁcation. A severe squall line case during the Southwest Monsoon Experiment Intensive Observing Period 8 (SoWMEX IOP#8) on 14 June 2008, is selected and examined. Because the operator is adopted from the use of the midlatitude region, sensitivity tests are performed to obtain the optimal setting of the operator in the subtropical region. To accurately capture the dynamic structure of the squall lines, the ensemble-based data assimilation system, which assimilates both radial wind and reﬂectivity data, is used to obtain the optimal analysis ﬁeld on the mesoscale for evaluating the performance of model simulation. The characteristics of two microphysics schemes are investigated, and the results obtained using the schemes are compared with the S-band dual-polarimetric radar observations. The horizontal and vertical cross-sections show that the analyses resemble the observations. Both schemes can replicate the polarimetric parameter signature such as Z DR and K DP columns. When comparing model simulation with polarimetric parameters through the drawing of contour frequency by altitude diagrams (CFADs), the results reveal that the single moment microphysics scheme performs better than the double moment scheme in this case. However, the reﬂectivity ﬁeld in the stratiform area is more accurately captured when using the double moment scheme. Furthermore, validation with polarimetric variables ( Z H , Z DR and K DP ) histograms shows underestimation of the K DP ﬁeld in both schemes. Overall, this study indicates the beneﬁt of assimilating radial wind and reﬂectivity data for the analyses of severe precipitation systems and the necessity of assimilating polarimetric parameters for the accuracy of microphysical processes, especially complex microphysics schemes in subtropical region. percentile. results indicate that if double moment scheme applied in the numerical model simulation, dual-polarimetric as and K DP , to be assimilated for obtaining better hydrometeor variables.


Introduction
Meteorological radar can make observations of high spatial and temporal resolution, which aid weather surveillance and nowcasting. In the traditional Doppler weather radar, reflectivity (Z H ) and radial wind (V R ) observations provide crucial information for analyzing weather systems, enabling comprehension of the formation and structure of severe weather [1]. Furthermore, radar data also improve forecasts when applied to echo extrapolation nowcasting [2,3] and data assimilation (DA) techniques [4][5][6]. Dual-polarimetric weather radar offers more information than Doppler radar, such as the differential reflectivity (Z DR ), specific differential phase (K DP ), and co-polar correlation coefficient (ρ HV ). These polarimetric radar parameters provide information such as the shape of hydrometeors and liquid water content distribution in the atmosphere; they also enable the distinguishing of meteorological radar data simulator (PRDS), which can convert the hydrometeor mixing ratio and total number concentration into polarimetric parameters. A supercell cell case was tested in observing system simulation experiments (OSSEs). The results revealed that the location of the hail mixing ratio maximum overlapped occurrence of the minimum Z DR , and the position of the rain mixing ratio maximum and K DP maximum were the same, indicating that the PRDS successfully converted the model output variables into polarimetric parameters. In the same year, scholars [27] used a system called synthetic polarimetric radar (SynPolRad), in which the T-matrix method was used to simulate the propagation of radar electromagnetic waves in the atmosphere. SynPolRad can convert model outputs into Z DR , the linear depolarization ratio (L DR ), and reflectivity. In comparison with observations, the performance of liquid-phase hydrometeors more closely matched the characteristics of weather systems. However, the distribution of ice-phase hydrometeors contained clear position errors. The study [28] suggested the application of a forward operator to a spectral (bin) microphysics scheme and tested it using a real case. The results demonstrated that the simulation conducted using the bin model of polarimetric parameters was superior to the bulk models. However, the bin model microphysics scheme was too complex and thus the computational cost was immense. In another study [29], a moment-based warm rain observation operator was designed. The relationship between the moments and polarimetric parameters was established on the basis of more than 2 million disdrometer data collected globally. The meaningful statistical data and low-uncertainty data showed that the combination of the sixth moment and ninth moment outperformed the combination of the zeroth moment and third moment in bulk microphysics schemes. These forward operators provided a connection between the model variables and dual parameters. With the dual-polarimetric forward operators in particular, the model simulation results could be converted into polarimetric variables. Consequently, comparison of microphysics schemes has become a new research spotlight. Several studies [30][31][32] have indicated that double-moment schemes can simulate polarimetric signatures such as the Z DR arc and Z DR column in supercells. The fixed intercept parameter in single moment schemes inhibits DSD flexibility; therefore, the true observational DSD cannot be mimicked. Additionally, some scholars [33] indicated that microphysics schemes underestimate the liquid water content compared with the observed K DP . Therefore, by comparing with polarimetric variables, microphysics schemes have been discovered to still have uncertainties that must be improved. By assimilating polarimetric parameters, several studies [34][35][36] have discovered that polarimetric signatures such as the Z DR arc and Z DR column and DSD-related variables are more consistent with observations.
In Taiwan, the weather radar network will soon be completely upgraded to dual-polarimetric radars. Various studies regarding the application of Doppler radars have reported the benefit to rainfall evaluations of rainfall estimation [37] and DA techniques [38,39]. The application of polarimetric radar data has mainly focused on QPE [40], the PID technique [41], and DSD retrieval [42,43] for microphysical studies. Few studies have applied these data to the validation of model forecasts and DA when using dual-polarimetric parameters. The study [44] directly employed a PRDS to assimilate polarimetric radar data in a typhoon case in Taiwan. The present study selects the case of a squall line severe weather system. Sensitivity tests first illustrate the necessity of a setup threshold for the PRDS when examining model outputs to observations over the Taiwan area. Then, radial wind and reflectivity observations are assimilated to ensure the precipitation system is close to reality, and the optimal setting of the PRDS is employed to evaluate the performance of a convection-permitting model when using two microphysics schemes and dual-polarimetric variables. The article is organized as follows. The model setting and DA strategy are detailed in Section 2. The PRDS and sensitivity test are presented in Section 3. The results of simulations are examined in Section 4. Section 5 is the discussion. Finally, the conclusions are given in Section 6.

Case Overview
The case of an MCS during Southwest Monsoon Experiment intensive observing period 8 (SoWMEX-IOP8) in 2008 is investigated. SoWMEX began on 15 May and ended on 30 June 2008.
Around May and June, the moist and warm southwesterly frequently causes heavy rainfall events in Taiwan. Consequently, the purpose of SoWMEX is to determine (1) the development mechanism of MCSs caused by the southwesterly; (2) the interaction between the southwesterly and complex terrain in Taiwan; (3) the dynamic, thermodynamic characteristics and microphysics processes during the southwesterly; and (4) development of the radar DA technique and improvement of the quantitative precipitation forecast (QPF) in Taiwan. Therefore, abundant observations are available for the experimental period, such as data from the National Center of Atmospheric Research (NCAR) S-band dual-polarimetric radar (S-POL) deployed on the southwestern coast of Taiwan. IOP8 began at 0000 UTC on 14 June and ended at 0000 UTC on 17 June. The stationary front was located along the southeastern coast line of China (Figure 1a), and southwesterly flow notably influenced Taiwan (Figure 1b). The skew-T diagram obtained at 0600 UTC during IOP8 (not shown) reveals a strong southwesterly flow (30-40 knots) with abundant water vapor content (approximately 15 g/kg), and thus Taiwan was situated in an unstable area (convective available potential energy (CAPE) > 1600 J/kg). In addition, the 500-hpa chart (Figure 1c) revealed that a trough was located near western Taiwan, fueling the development of severe systems. The 200-hpa chart showed that the difluence area was located in the upper air of Taiwan ( Figure 1d). In summary, all conditions at the synoptic scale indicated a favorable environment for the development of severe weather. The case of an MCS during Southwest Monsoon Experiment intensive observing period 8 (SoWMEX-IOP8) in 2008 is investigated. SoWMEX began on 15 May and ended on 30 June 2008. Around May and June, the moist and warm southwesterly frequently causes heavy rainfall events in Taiwan. Consequently, the purpose of SoWMEX is to determine (1) the development mechanism of MCSs caused by the southwesterly; (2) the interaction between the southwesterly and complex terrain in Taiwan; (3) the dynamic, thermodynamic characteristics and microphysics processes during the southwesterly; and (4) development of the radar DA technique and improvement of the quantitative precipitation forecast (QPF) in Taiwan. Therefore, abundant observations are available for the experimental period, such as data from the National Center of Atmospheric Research (NCAR) S-band dual-polarimetric radar (S-POL) deployed on the southwestern coast of Taiwan. IOP8 began at 0000 UTC on 14 June and ended at 0000 UTC on 17 June. The stationary front was located along the southeastern coast line of China (Figure 1a), and southwesterly flow notably influenced Taiwan (Figure 1b). The skew-T diagram obtained at 0600 UTC during IOP8 (not shown) reveals a strong southwesterly flow (30-40 knots) with abundant water vapor content (approximately 15 g/kg), and thus Taiwan was situated in an unstable area (convective available potential energy (CAPE) > 1600 J/kg). In addition, the 500-hpa chart (Figure 1c) revealed that a trough was located near western Taiwan, fueling the development of severe systems. The 200-hpa chart showed that the difluence area was located in the upper air of Taiwan (Figure 1d). In summary, all conditions at the synoptic scale indicated a favorable environment for the development of severe weather.   Figure 2a) were captured by the NCAR S-POL. Squall line A was made landing and causing rainfall for southern Taiwan. Squall lines B and C were developing and moving toward Taiwan. These squall-line systems formed and affected southern Taiwan, causing extremely heavy rainfall that lasts 6 h. The accumulated rainfall was 150 mm per day, and some mountainous areas received more than 200 mm per day on 14 June ( [41] Figure 2). Studies have focused on the mechanism of the severe weather system [45], analysis of microphysics properties through observations [41], application of QPE [46], and QPF [47,48]. The present study is the first to simulate polarimetric variables and validate model behavior using these variables.  Figure 2). Studies have focused on the mechanism of the severe weather system [45], analysis of microphysics properties through observations [41], application of QPE [46], and QPF [47,48]. The present study is the first to simulate polarimetric variables and validate model behavior using these variables.

Model Configuration
In this study, a non-hydrostatic and compressible model, Weather Research and Forecast (WRF, [49,50] Figure 3a). In total, 52 layers with terrain-following coordinates are given in the vertical direction, and the model top is set at 10 hPa. The following physical parameterizations are used in this study: the Rapid Radiative Transfer Model (RRTM) longwave scheme [51], Dudhia shortwave scheme [52], Yonsei planetary boundary layer scheme [53], and Grell-Freitas cumulus parameterization scheme [54]. In addition, two microphysics schemes are selected to examine the performance: the Goddard cumulus ensemble (GCE) scheme [16] and Morrison scheme [13]. Both the GCE and Morrison schemes consider five classes of hydrometeors (cloud water, rain, cloud ice, snow, and graupel/hail). In both microphysics schemes, graupel is selected as the rimed category species for the default setting in WRF. The GCE scheme is a single-moment microphysics scheme that predicts the hydrometeor mixing ratio only, whereas the Morrison scheme is a complete double-moment microphysics scheme in which the mixing ratio and total number concentration are prognostic variables. Initial and lateral boundary conditions (every 6 h) are provided by National Centers for Environmental Prediction-Global Forecast System (NCEP-GFS) 1 • × 1 • reanalysis data obtained from 0000 UTC on 14 June to 0000 UTC on 15 June 2008.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 25 210), respectively ( Figure 3a). In total, 52 layers with terrain-following coordinates are given in the vertical direction, and the model top is set at 10 hPa. The following physical parameterizations are used in this study: the Rapid Radiative Transfer Model (RRTM) longwave scheme [51], Dudhia shortwave scheme [52], Yonsei planetary boundary layer scheme [53], and Grell-Freitas cumulus parameterization scheme [54]. In addition, two microphysics schemes are selected to examine the performance: the Goddard cumulus ensemble (GCE) scheme [16] and Morrison scheme [13]. Both the GCE and Morrison schemes consider five classes of hydrometeors (cloud water, rain, cloud ice, snow, and graupel/hail). In both microphysics schemes, graupel is selected as the rimed category species for the default setting in WRF. The GCE scheme is a single-moment microphysics scheme that predicts the hydrometeor mixing ratio only, whereas the Morrison scheme is a complete double-moment microphysics scheme in which the mixing ratio and total number concentration are prognostic variables. Initial and lateral boundary conditions (every 6 h) are provided by National Centers for Environmental Prediction-Global Forecast System (NCEP-GFS) 1° × 1° reanalysis data obtained from 0000 UTC on 14 June to 0000 UTC on 15 June 2008.

WRF Local Ensemble Transform Kalman Filter Radar Assimilation System (WLRAS)
An ensemble DA system, called the WRF Local ensemble transform Kalman filter Radar Assimilation System (WLRAS) [38], is used. The WLRAS is based on the concept of the local ensemble transform Kalman filter (LETKF) [55]. WRF is used to conduct model forecasting during the forecast step in DA cycles. The WLRAS enables multiple radar data to be assimilated and adjusts the model's prognostic variables-such as horizontal and vertical wind, geopotential height, potential temperature, mixing ratio, and number concentration of hydrometeors-through cross-variable background error correlations. Furthermore, different localization radii are applied to model prognostic variables in the WLRAS. Horizontal and vertical localization radii are established to select particular observations within a certain distance to ensure more precise correction of variables and increase the computational efficiency. The horizontal and vertical localization radii are given, referring to a previous study [38], in this research ( Table 1). The observation errors for reflectivity and radial wind are assumed to be 5 dBZ and 3 m/s, respectively. Before assimilation, model variables are converted into observation variables and interpolated into observational space using the observation operators (described in Section 3.4).

WRF Local Ensemble Transform Kalman Filter Radar Assimilation System (WLRAS)
An ensemble DA system, called the WRF Local ensemble transform Kalman filter Radar Assimilation System (WLRAS) [38], is used. The WLRAS is based on the concept of the local ensemble transform Kalman filter (LETKF) [55]. WRF is used to conduct model forecasting during the forecast step in DA cycles. The WLRAS enables multiple radar data to be assimilated and adjusts the model's prognostic variables-such as horizontal and vertical wind, geopotential height, potential temperature, mixing ratio, and number concentration of hydrometeors-through cross-variable background error correlations. Furthermore, different localization radii are applied to model prognostic variables in the WLRAS. Horizontal and vertical localization radii are established to select particular observations within a certain distance to ensure more precise correction of variables and increase the computational efficiency. The horizontal and vertical localization radii are given, referring to a previous study [38], in this research ( Table 1). The observation errors for reflectivity and radial wind are assumed to be 5 dBZ and 3 m/s, respectively. Before assimilation, model variables are converted into observation variables and interpolated into observational space using the observation operators (described in Section 3.4).

Radar Data QC and Process
In this study, five radars-RCWF, RCCG, RCKT, S-POL, and NCU-are employed, the locations of which are shown in Figure 3b. RCWF, RCCG, and RCKT are S-band ground-based Doppler radars belonging to the Central Weather Bureau (CWB). NCAR S-POL is a dual-polarization S-band radar. Finally, NCU is a C-band dual-polarimetric radar operated by National Central University. For each radar, reflectivity and radial wind data are subjected to QC. For data from the CWB radars, ground cluster, sea cluster, and noise signals have been removed [56]. For the S-POL and NCU data, the value of ρ HV over 0.85 (NCU) and 0.92 (S-POL) is applied to filter out non-meteorological signals. In addition, the radial wind fields of all radars are unfolded by the Rakit (Radar Kit), which is developed by Radar Laboratory of National Central University.
After QC, the radar data are processed in the form of super-observations [38,57,58], which comprise the following steps: (1) Divide the PPI coverage into several areas. The coverage is partitioned every 5 km along the radial direction and every 5 • of the azimuth angle. (2) Split into small areas. The values of variables in a particular area are calculated by using an inverse-distance weighting method with respect to the center of the area. The super-observation data are distributed more densely near the center of the radar and become sparse away from the radar center. Nevertheless, the super-observations retain the characteristic of weather systems and reduce computational costs.

Observation Operator
As mentioned in Section 3.2, converting the model variables into observation variables in observational space is necessary for the purpose of assimilation. Spatial interpolation is performed by finding eight grid points that enclose the observation point with inverse-distance weights. The model variables are then converted into observation variables. The following introduces the radial wind and reflectivity operator.
The formula for converting wind components in the model into radial wind is presented in Equation (3); x, y, and z are the distances from the center of radar in the three dimensions, and u, v, and w are the zonal, meridional and vertical wind. Whereas v t is terminal velocity of rain, which is calculated as in Equation (4). In (4), p 0 is the pressure at the surface, p is mean pressure at that height, ρ a is the air density, and q r is the rain mixing ratio.
The observation operator of dual-polarimetric parameters used in this research was first proposed by [26]. This operator converts model variables, such as the mixing ratio and total number concentration, of hydrometers into radar variables (reflectivity, differential reflectivity, and specific differential phase). Several characteristics are considered to match the observations. Here, we briefly introduce some features of the operator, but the details can be found elsewhere [26,44]. (i) Considering the oblate shape of raindrops, a raindrop axis ratio relationship that changes with drop diameter is applied in the T-matrix method which can simulate the propagation of radar electromagnetic waves in the atmosphere. Besides, the raindrop axis ratio is obtained from the observation of continental US. Consequently, the operator simulates the horizontal and vertical reflectivity as dual polarimetric radar observations. In addition, the Rayleigh scattering approximation is employed to simulate the reflectivity of ice-phase hydrometeors in the horizontal and vertical directions, and a constant axis ratio of 0.75 is set for snow and hail particles. After retrieving the backscattering amplitudes, a power law function is fitted for calculation purposes. The backscattering amplitude of hydrometeors in horizontal and vertical directions can be expressed by Equations (5) and (6).
The f a,x and f b,x represent for the horizontal and vertical backscattering amplitudes, x represents for the type of hydrometeor (rain, snow graupel/hail, rain-snow and rain-graupel/hail), and α and β stand for the coefficients of fitting power law function. The value of α and β please refer to section 3.C of [26]. (ii) The melting effect is considered in the operator. As the ice-phase particles fall into the region with temperature higher than 0 • C, they begin to melt. Therefore, a thin layer of water covers the surface of the ice-phase particles. Because of the differing ability to reflect electromagnetic waves, radars may misidentify the melting particles as large rain drops, resulting in enlarged reflectivity signals. This is the so called "brisectionght band" effect, which is relatively ubiquitous in radar observations. The operator checks a certain grid point and examine whether rain and ice-phase particles coexisted. Afterward, the conversion fraction according to the amounts of rain and snow/hail particles is calculated, and the mass of both types of particle merge together, so the mixture "rain-snow" or "rain-hail" is decided. (iii) The canting angle effect is considered in the operator. The particles themselves tumble and tilt when they fall to the ground. This effect has an impact on the horizontal and vertical backscattering energy received by radars. Considering the aforementioned effects, the model variables are converted into observation variables using the following formulas: Z v,r = 4λ 4 α rb 2 N 0r Equations (7) and (8) are the formulas of the rain species for horizontal and vertical reflectivity while Equations (9) and (10) represent for the snow, graupel/hail, rain-snow and rain-graupel/hail species. Equation (11) shows the formula for specific differential phase. λ is the radar wavelength (10.7 cm in this research); K w is the dielectric coefficient of water; A, B and C are the coefficients for canting angle effect; α k and β k are the nondimensional coefficients. Λ x and N 0x are the slope parameter and intercept parameter in the DSD, respectively, and are calculated using Equations (12)- (14): Remote Sens. 2020, 12, 3004 where q x and N tx are the model outputs, representing the hydrometeor mixing ratio and total number concentration, respectively; µ is the shape parameter; and Γ is the gamma function. For single moment schemes with a fixed intercept parameter, Equation (12) can be used to obtain the slope parameter.
In the final step, the contributions from rain, snow, hail, rain-snow, and rain-hail are summed. To determine the reflectivity and differential reflectivity, the contributions are summed in a linear scale and then converted into a logarithmic scale. The specific differential phase is simply summed to yield the total contribution.

Experimental Design
where and are the model outputs, representing the hydrometeor mixing ratio and total number concentration, respectively; is the shape parameter; and Γ is the gamma function. For single moment schemes with a fixed intercept parameter, Equation (12) can be used to obtain the slope parameter. In double moment schemes, Equations (13) and (14) must be employed.
In the final step, the contributions from rain, snow, hail, rain-snow, and rain-hail are summed. To determine the reflectivity and differential reflectivity, the contributions are summed in a linear scale and then converted into a logarithmic scale. The specific differential phase is simply summed to yield the total contribution.

Results
The observation operator designed in a previous study [26] is based on climatological information (e.g., temperature and relationship of raindrop axis ratio) in the mid-latitude region. Because microphysics processes and characteristics could be different in mid-latitude and subtropical areas [42], examining the setup of polarimetric radar observation operators when applying them over Taiwan area is necessary. Therefore, some parameters and thresholds are

Results
The observation operator designed in a previous study [26] is based on climatological information (e.g., temperature and relationship of raindrop axis ratio) in the mid-latitude region. Because microphysics processes and characteristics could be different in mid-latitude and subtropical areas [42], examining the setup of polarimetric radar observation operators when applying them over Taiwan area is necessary. Therefore, some parameters and thresholds are investigated through sensitivity tests, and the optimal setting is obtained for the subsequent experiments.
After modification of the operator, the reflectivity and radial wind are assimilated to enable correction and adjustment of the dynamic and thermodynamic fields so that the analysis field is close to the observation field. The differences in microphysics schemes and performance can then be investigated in simulations of the polarimetric radar parameters.

Sensitivity Tests of the Operator
The deterministic forecast mentioned in Section 3.5 is used to demonstrate the result of sensitivity tests. The deterministic forecast at 1100 UTC on 14 June shows that the location of the squall-line system was the center of the Taiwan Strait (not shown). Compared with radar observation (Figure 2d), the model simulation reveals both the temporal and spatial phase errors. Therefore, a strategy as [62,63], which ignored the position error and selected the main precipitation area between model and observation, is applied to compare the model simulation and observation. As a result, without assimilating radar data to correct the position error, we select a domain with area the same as the red box in Figure 1d and that covers the main convection and stratiform zone for the sensitivity tests. Contoured frequency by altitude diagrams (CFADs) [64] of reflectivity are employed to examine the vertical structure of the weather systems. The data selection and grouping threshold used when manufacturing the CFADs are shown in Table 2. The domain selected in observation fields is illustrated in Figure 2d (S-POL), indicated by the red box. The size of the selected domain in deterministic forecasts is the same as the selected area in S-POL, and the domain covers both the convection and stratiform areas of squall lines. S-POL data are interpolated to 3-and 0.25-km horizontal and vertical resolution, respectively. The analysis height is from 1 to 15 km above sea level (ASL). By investigating the CFAD, the signal of the bright band structure in reflectivity fields is revealed, and whether the conversion ratio of mixed-phase hydrometeors in the operator is more suitable can be evaluated. The default value of the transfer ratio in this study [26] is 0.5 and 0.4 for rain-snow and rain-hail, respectively. Table 2. Data selection range and interval for reflectivity, differential reflectivity, and specific differential phase.  Figure 5 shows the CFAD of observed and simulated reflectivity at lower than 8-km height at 1100 UTC on 14 June 2008. Following the method in [65], the 25% (solid line), 50% (long-dashed line), and 75% (short-dashed line) accumulated frequency percentages are presented to enable statistical evaluation of data dispersion. Clearly, the observed reflectivity indicates a slight increase in the percentage at approximately 5 km, which indicates the height of 0 • C. When applying the default setting of the observation operator [26], the simulated reflectivity in Figure 5b displays the bright band structure around the same height (~5 km). In this article, the bright band structure obtained using the Morrison scheme is clearer than that obtained using the GCE scheme; thus, the Morrison scheme is selected for the subsequent sensitivity tests. However, a distinct discontinuity distribution is also revealed. The enhanced reflectivity signal is result from the melting snow and graupel. The mixing phase hydrometers are often regarded as large raindrop due to the surface of which shedding a thin layer of water. Since the fraction of mixing phase is tunable in the operator, the maximum conversion fraction of the mixing phase is modified and examined if the structure between 3 and 5 km can be improved to be closer to the observation. Figure 5c,d show the CFADs of modification obtained by halving the maximum conversion fraction of both mixing-phase hydrometeors (reducing from 0.5 to 0.25 for rain-snow and from 0.4 to 0.2 to rain-hail species) and considering no mixing effect. Figure 5c shows that the peak values of the 25%, 50%, and 75% lines are reduced by approximately 3-5 dBZ. However, a discontinuous distribution remains within the 3-5-km heights. Moreover, no bright band structure exists in Figure 5d, which was obtained by considering the absence of melting particles. The result of the sensitivity tests of the maximum conversion fraction suggests that the fraction influences the reflectivity near the melting level.

Radar Variables
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 25 3-5-km heights. Moreover, no bright band structure exists in Figure 5d, which was obtained by considering the absence of melting particles. The result of the sensitivity tests of the maximum conversion fraction suggests that the fraction influences the reflectivity near the melting level. In double moment schemes, total number concentration is another factor affecting reflectivity. In the situation of the same mixing ratio, a larger total number concentration results in a smaller mean mass diameter of particles. Here, the minimal total number concentration of mixed-phase hydrometeors is investigated. Figure 6 shows the results of a sensitivity test for the adjustment of the total number concentration; Figure 6a-c represent the minimum total number concentration as 10/m 3 , 100/m 3 , and 1000/m 3 , respectively. As long as the mixed-phase hydrometeors are considered, a portion of the total number concentration must also be converted. The conversion proportion is the same as the mixing ratio conversion fraction calculated using coexisting hydrometeors. However, the total number concentration of mixed-phase hydrometeors is only converted from the total number concentration of snow or hail species. The reason is that when ice-phase hydrometeors melt, the liquid water simply attaches to the surface of the ice particles without liquid particle numbers decreasing. In addition, a specified threshold of the total number concentration is essential for preventing an unrealistic total number concentration or excessively simulated reflectivity. If the total number concentration does not exceed the threshold, the mixed-phase hydrometeors are not considered at that grid point. Figure 6a shows that the bright band structure is less discontinuous than that presented in Figure 5b, proving that some of the extreme reflectivity values have been In double moment schemes, total number concentration is another factor affecting reflectivity. In the situation of the same mixing ratio, a larger total number concentration results in a smaller mean mass diameter of particles. Here, the minimal total number concentration of mixed-phase hydrometeors is investigated. Figure 6 shows the results of a sensitivity test for the adjustment of the total number concentration; Figure 6a-c represent the minimum total number concentration as 10/m 3 , 100/m 3 , and 1000/m 3 , respectively. As long as the mixed-phase hydrometeors are considered, a portion of the total number concentration must also be converted. The conversion proportion is the same as the mixing ratio conversion fraction calculated using coexisting hydrometeors. However, the total number concentration of mixed-phase hydrometeors is only converted from the total number concentration of snow or hail species. The reason is that when ice-phase hydrometeors melt, the liquid water simply attaches to the surface of the ice particles without liquid particle numbers decreasing. In addition, a specified threshold of the total number concentration is essential for preventing an unrealistic total number concentration or excessively simulated reflectivity. If the total number concentration does not exceed the threshold, the mixed-phase hydrometeors are not considered at that grid point. Figure 6a shows that the bright band structure is less discontinuous than that presented in Figure 5b, proving that some of the extreme reflectivity values have been corrected. Figure 6c shows that the threshold is so high that the bright band is almost absent. In Figure 6b, we infer that the bright band structure is more suitable for this case because the CFAD is less discontinuous than that in Figure 6a. When performing DA in this experiment, the setting of the threshold is 100/m 3 .
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 25 corrected. Figure 6c shows that the threshold is so high that the bright band is almost absent. In Figure 6b, we infer that the bright band structure is more suitable for this case because the CFAD is less discontinuous than that in Figure 6a. When performing DA in this experiment, the setting of the threshold is 100/m 3 .
(a) (b) (c)    (Figure 2d) show two squall lines approaching Taiwan. Squall line B landed in southwestern Taiwan, and strong reflectivity (>45 dBz) appears as a Y-shaped feature, indicating a strong convection area along the southwestern coast. Squall line C is located offshore of Taiwan behind squall line A. Both squall lines have strong reflectivity, which is up to 45-50 dBZ, and broad stratiform areas are distributed to the west of the main convection area. Because the reflectivity field is assimilated into the convection-permitting model, the simulated reflectivity field at the analysis time (1100 UTC) is expected to resemble the observations; the simulated reflectivity obtained using the GCE (Figure 7a) and Morrison (Figure 7b) schemes favorably captures the Y-shaped strong convection structure of squall line B near the coast. In addition, simulation using the GCE scheme obtains clearer squall line structure than that using the Morrison scheme. In the analysis using the Morrison scheme, the Y-shaped structure is slightly different from the observation, and the southern part of the precipitation system is excessively weak.

Performance of the Analysis Against the S-POL
(a) (b)   Taiwan, and strong reflectivity (>45 dBz) appears as a Y-shaped feature, indicating a strong convection area along the southwestern coast. Squall line C is located offshore of Taiwan behind squall line A. Both squall lines have strong reflectivity, which is up to 45-50 dBZ, and broad stratiform areas are distributed to the west of the main convection area. Because the reflectivity field is assimilated into the convection-permitting model, the simulated reflectivity field at the analysis time (1100 UTC) is expected to resemble the observations; the simulated reflectivity obtained using the GCE (Figure 7a) and Morrison (Figure 7b) schemes favorably captures the Y-shaped strong convection structure of squall line B near the coast. In addition, simulation using the GCE scheme obtains clearer squall line structure than that using the Morrison scheme. In the analysis using the Morrison scheme, the Y-shaped structure is slightly different from the observation, and the southern part of the precipitation system is excessively weak.

Performance of the Analysis Against the S-POL
Next, the vertical structure of the squall line is evaluated. Because the simulated reflectivity reflects the precipitation system favorably when using either the GCE or Morrison scheme, Z DR and K DP are further examined. Figure 8 shows the average vertical cross-section (indicated by a rectangle enclosed by dashed black boxes in Figures 2d and 7) of reflectivity Z H (Figure 8a,d,g), differential reflectivity Z DR (Figure 8b,e,i), and the specific differential phase K DP (Figure 8c analysis time (1100 UTC) is expected to resemble the observations; the simulated reflectivity obtained using the GCE (Figure 7a) and Morrison (Figure 7b) schemes favorably captures the Y-shaped strong convection structure of squall line B near the coast. In addition, simulation using the GCE scheme obtains clearer squall line structure than that using the Morrison scheme. In the analysis using the Morrison scheme, the Y-shaped structure is slightly different from the observation, and the southern part of the precipitation system is excessively weak. Near the main convection area in S-POL (Figure 8a), the distributions of ( Figure 8b) and (Figure 8c) form a vertical column-like structure that extends 6 km. We suspect that this structure is owing to strong vertical velocity that brings large raindrops to upper levels. A sufficiently strong upward motion keeps large raindrops in the air, resulting in high and . In both analyses, the high-area (Figure 8e,f) corresponds to the main updraft area, proving that the column signature is reproduced in both analyses. Additionally, the strong column structure is simulated in both schemes (Figure 8h,i) and occurs in the main updraft area. The column is near the strong convection area, which is reasonable because of the high liquid water content near that area. The results illustrate that the polarimetric signature can be replicated in the model to some extent.

Statistical Verification of System Structure Using the CFAD
To evaluate the performance of the model statistically, CFADs are plotted that illustrate the  (Figure 8g), the characteristics of the squall-line system are captured after radar DA. The bright band structure in the Morrison scheme result is much stronger than that in the GCE scheme result. Both schemes can handle the position of strong convection areas, and the system can develop up to 9-km height.
In the S-POL observations, the stronger Z DR (Figure 8b) and K DP (Figure 8c) signals overlap near the strong convection area. The peak of Z DR is approximately 1.2 dB, revealing that the raindrops are mainly spherical. The Z DR distributions obtained using the GCE (Figure 8e) scheme resemble the Z H distribution obtained using the GCE scheme (Figure 8d), because these radar parameters are converted using the hydrometeor mixing ratio. In the Morrison scheme (Figure 8f), the polarimetric parameters are composed of the mixing ratio and total number concentration, and thus the Z DR pattern does not resemble the Z H pattern (Figure 8c). The benefit of the double moment is that the size sorting effect is considered. The slightly deviated position of the high Z DR value from Z H can be seen at the main core of the convection.
The K DP in the S-POL observation (Figure 8c) can be up to 1 • /km, indicating that the region has high liquid water content. Above 6-km height, K DP is enhanced, implying that this phenomenon is caused by the irregular shape of rimming crystals and large raindrops. The K DP distributions obtained using both schemes (Figure 8h,i) are almost the same and located in the core convective area. However, both GCE ( Figure 8h) and Morrison (Figure 8i) simulated K DP above 6-km ASL could not duplicate similar patterns as the S-POL observations (Figure 8c). This is due to the current operator considered the shape of ice-phase hydrometeors as elliptic as previous study [26], leading difficulties to presenting polarimetric variables in cold-rain processes.
Near the main convection area in S-POL (Figure 8a), the distributions of Z DR (Figure 8b) and K DP (Figure 8c) form a vertical column-like structure that extends 6 km. We suspect that this structure is owing to strong vertical velocity that brings large raindrops to upper levels. A sufficiently strong upward motion keeps large raindrops in the air, resulting in high Z DR and K DP . In both analyses, the high-Z DR area (Figure 8e,f) corresponds to the main updraft area, proving that the Z DR column signature is reproduced in both analyses. Additionally, the strong K DP column structure is simulated in both schemes (Figure 8h,i) and occurs in the main updraft area. The K DP column is near the strong convection area, which is reasonable because of the high liquid water content near that area. The results illustrate that the polarimetric signature can be replicated in the model to some extent.

Statistical Verification of System Structure Using the CFAD
To evaluate the performance of the model statistically, CFADs are plotted that illustrate the characteristics of the model simulation and observations. The coverage of selected data for CFAD is shown by the red boxes in Figure 2d (S-POL observation) and Figure 7 (Analysis field). Besides, the polarimetric data range and interval are shown in Table 2. Figure 9 shows the S-POL CFAD of Z H (Figure 9a), Z DR (Figure 9b), and K DP (Figure 9c) at 1100 UTC on 14 June 2008. Figure 10 displays the radar variables at 1100 UTC on 14 June 2008, of the CFAD for the DA experiments when using different microphysics schemes. When examining the variables Z DR and K DP , the CFAD is only demonstrated below 6-km ASL. This is because complex cold rain processes result in the observation operator having difficulty simulating the higher level Z DR and K DP . Therefore, we follow previous studies that have focused on the warm rain region [32,33].  (a) (b) (c) Figure 9. CFADs of S-POL observations at 1100 UTC on 14 June 2008: (a) reflectivity, (b) differential reflectivity, and (c) specific differential phase. Figure 10. CFAD of the analysis field (DA) at 1100 UTC on 14 June 2008: (a-c) reflectivity, differential reflectivity, and specific differential phase, respectively, in the GCE scheme. (d-f) same as (a-c) but for the Morrison scheme.

Inspection of the Differences in Microphysics
In this section, we discuss the different behavior during DA. In the GCE scheme, only the hydrometeor mixing ratios are predicted. Therefore, during the analysis step of DA cycling, the DA system modifies the mixing ratio directly. Additionally, the polarimetric radar observation operator only utilizes the mixing ratio for the radar parameter conversion in the GCE scheme. Consequently, the modification of radar parameters is clear and direct. Figure 11 presents the scatter plots (with correlation coefficient) between rain mixing ratio increment and radar parameter increment in the GCE scheme (Figure 11a-c represent , , and , respectively). The -and relationships are similar; positive (negative) increment corresponds to increases (decreases) in Figure 10. CFAD of the analysis field (DA) at 1100 UTC on 14 June 2008: (a-c) reflectivity, differential reflectivity, and specific differential phase, respectively, in the GCE scheme. (d-f) same as (a-c) but for the Morrison scheme.
In Figure 9a, focusing on the warm rain area first, the 50th percentile below 5-km ASL is approximately 28 dBZ. The reflectivity is concentrated between 21 and 37 dBZ. The maximum reflectivity is approximately 48 dBZ. At 5-km ASL, the probability of high reflectivity is greater, indicating that the bright band structure is located near 5-km ASL. In the upper layer, the echo is maintained within 10-20 dBZ. For the Z H CFADs after DA (Figure 10a,d), the 25th, 50th, and 75th percentiles are close to the observation, although the 50th and 75th percentiles are slightly underestimated in both schemes. Furthermore, the bright band structure is enhanced at 5-km ASL in the model, with increasing Z H probability indicating the existence of mixed-phase hydrometeors. In the mid-level (5-9 km), Z H is overestimated by approximately 5 dBZ in the GCE scheme ( Figure 10a) and 10 dBZ in the Morrison scheme (Figure 10d). In higher air (above 9 km), the analyses show that the estimated probabilities of Z H are near approximately 10 dBZ, underestimated when compared with the observation. Nonetheless, the analysis fields in both schemes are relatively reliable in this research, indicating that the thermodynamic fields match the observation after DA.
The Z DR CFAD in S-POL (Figure 9b) shows the 50th percentile at approximately 0.5 dB in the warm cloud region, indicating that the raindrops tend to be spherical. The maximum Z DR is approximately 1.5 dB. Therefore, judging from the Z H and Z DR values, the DSD in this event may indicate small but highly concentrated raindrops. In the Z DR CFADs (Figure 10b,e), both schemes overestimate Z DR after DA. The GCE scheme (Figure 10b) maintains the Z DR value at 1.2 dB, whereas the Z DR in the Morrison scheme (Figure 10e) increases with decreasing height ASL. The behavior of Z DR can be explained as follows. In the GCE scheme (Figure 10b), when calculating the polarimetric variables, only the hydrometeor mixing ratio is considered because of the single moment setting. Consequently, the mixing ratio increases during DA cycles. Because of the fixed intercept parameter and increasing mixing ratio, the overall raindrop size increases, as reflected by an increasing mean droplet diameter. Therefore, the outcome leads to increased Z DR . In the Morrison scheme (Figure 10e), the increased mixing ratio causes increased Z DR , whereas the extra parameter, the total number concentration, also affects Z DR . Apart from the mixing ratio, the default setting of the PRDS causes the overestimation of Z DR . As the operator is being constructed, the raindrop axis ratio relationship is applied in the T-matrix method to retrieve the backscattering amplitude functions. We mentioned in Section 3 that the default axis ratio relationship is derived from the observations made in the contiguous United States. The raindrops are elliptical in shape and large in diameter. Considering the observations made in Taiwan, such as the S-POL Z DR CFAD (Figure 9b), the differences in cloud microphysics processes and DSD features between the United States and Taiwan can be obtained.
Regarding the K DP CFAD in S-POL (Figure 9c), half of the data are limited to 0.1 • /km, and the 75th percentile is at approximately 0.2 • /km. K DP can be treated as a proxy of liquid water content; however, it is not sensitive in the low rainfall rate region for S-band radar, such as the stratiform area [66]. Therefore, we simply distinguish the convective and stratiform areas in accordance with another study [67]. The result shows that the K DP of the convective area is between 0.3 • /km and 0.5 • /km (not shown). Although the low K DP may be affected by non-meteorological signals, we still assume that the K DP data can represent the true atmospheric situation, because the radar data are processed with sufficient QC. In the K DP CFAD analyses (Figure 10c,f), the frequency of high K DP is captured by the DA system although the 75th percentile is constrained below 0.05 • /km. In this section, polarimetric variables such as the differential reflectivity and specific differential phase are investigated. The general pattern in the simulated Z H may be the same, but further inspection of the other polarimetric variables such as Z DR and K DP can uncover the details in each microphysics scheme.

Inspection of the Differences in Microphysics
In this section, we discuss the different behavior during DA. In the GCE scheme, only the hydrometeor mixing ratios are predicted. Therefore, during the analysis step of DA cycling, the DA system modifies the mixing ratio directly. Additionally, the polarimetric radar observation operator only utilizes the mixing ratio for the radar parameter conversion in the GCE scheme. Consequently, the modification of radar parameters is clear and direct. Figure 11 presents the scatter plots (with correlation coefficient) between rain mixing ratio increment and radar parameter increment in the GCE scheme (Figure 11a-c represent Z H , Z DR , and K DP , respectively). The q r -Z H and q r -Z DR relationships are similar; positive (negative) increment corresponds to increases (decreases) in Z H and Z DR . The outcome is intuitive because all the raindrops size are growing (shrinking) as mixing ratio increasing (decreasing). The q r -K DP increment relationship is a strong correlation. K DP is proportional to the liquid water content, which is similar to the rain mixing ratio in the model. This is one reason the analysis field of the GCE scheme is closer to S-POL observations than that of the Morrison scheme. Figure 12 presents the Morrison scatter plot between the increment of mixing ratio or total number concentration and the radar parameters. Figure 12a-c show the q r -Z H , q r -Z DR , and q r -K DP increment relationships, whereas Figure 12d-f show the N t -Z H , N t -Z DR , and N t -K DP counterparts. The relationships between q r and the polarimetric variables resemble those in the GCE scheme. However, the Morrison scheme is a double-moment scheme, meaning the use of two model variables during conversion to the radar parameters. The q r -Z H (Figure 12a) and q r -Z DR (Figure 12b) relationships are thus less strongly correlated compared with the relationships in the GCE scheme. At first, we suspect that the double moment scheme would outperform the single moment scheme owing to the flexibility of DSDs. However, the analysis did not meet our expectation. The total number concentration and radar variables are nearly uncorrelated, especially N t -Z H (Figure 12d) and N t -Z DR (Figure 12e). The result illustrates that as reflectivity is assimilated, the impact upon N t is relatively weak. However, the N t -K DP (Figure 12f) shows stronger correlation compare to Z H and Z DR , indicating that the modification of N t could improve the K DP value slightly. When we investigated the CFAD in the previous section, the Z DR CFAD in the Morrison scheme becomes worse after DA. We infer that when the double moment scheme is applied, more radar observations such as Z DR and K DP are required to provide extra information. That is, if a sophisticated double moment scheme is used, assimilating the reflectivity only is insufficient. correlation compare to and , indicating that the modification of could improve the value slightly. When we investigated the CFAD in the previous section, the CFAD in the Morrison scheme becomes worse after DA. We infer that when the double moment scheme is applied, more radar observations such as and are required to provide extra information. That is, if a sophisticated double moment scheme is used, assimilating the reflectivity only is insufficient. correlation compare to and , indicating that the modification of could improve the value slightly. When we investigated the CFAD in the previous section, the CFAD in the Morrison scheme becomes worse after DA. We infer that when the double moment scheme is applied, more radar observations such as and are required to provide extra information. That is, if a sophisticated double moment scheme is used, assimilating the reflectivity only is insufficient. Although the analysis field in the Morrison scheme is not optimal, we can obtain the characteristic of microphysics schemes by simulating the polarimetric radar parameters. In Section 4.2, the comparison in the horizontal direction shows a similar pattern to the reflectivity. However, in reality, the same reflectivity may be obtained for fewer large drops or numerous small drops. That may represent a different microphysics feature if we investigate the polarimetric parameters further. Figure 13 demonstrates the vertical distribution of the area-mean hydrometeor variables (mixing ratio, total number concentration, and mean mass diameter which is calculated by Eqaution 15) at cycle 5, 1100 UTC June 14. The area is marked as the red rectangle in Figure 6. The total number concentration in the GCE scheme is calculated from the mixing ratio by using Equations (12) and (14). On average, the rain mixing ratio in the GCE scheme (Figure 13a) is twice that in the Morrison scheme (Figure 13d), whereas the total number concentration is also larger in Although the analysis field in the Morrison scheme is not optimal, we can obtain the characteristic of microphysics schemes by simulating the polarimetric radar parameters. In Section 4.2, the Z H comparison in the horizontal direction shows a similar pattern to the reflectivity. However, in reality, the same reflectivity may be obtained for fewer large drops or numerous small drops. That may represent a different microphysics feature if we investigate the polarimetric parameters further. Figure 13 demonstrates the vertical distribution of the area-mean hydrometeor variables (mixing ratio, total number concentration, and mean mass diameter which is calculated by Eqaution 15) at cycle 5, 1100 UTC June 14. The area is marked as the red rectangle in Figure 6. The total number concentration in the GCE scheme is calculated from the mixing ratio by using Equations (12) and (14). On average, the rain mixing ratio in the GCE scheme (Figure 13a) is twice that in the Morrison scheme (Figure 13d), whereas the total number concentration is also larger in the GCE scheme (Figure 13b), approximately triple that in the Morrison scheme (Figure 13e). The total number concentration is constant in the GCE scheme ( Figure 13b) but decreases as the height decreases in the Morrison scheme (Figure 13e). Using the mixing ratio and total number concentration, the mean mass diameter can be calculated as follows: The mean mass diameter combines the DSD parameters and is more intuitive to evaluation of the difference between microphysics schemes. The results indicate that the raindrops in the Morrison scheme (Figure 13f) are larger than those in the GCE scheme (Figure 13c). Furthermore, the mean mass diameter increases as ASL decreasing in the Morrison scheme (Figure 13f), indicating that the collision-coalescence process is more obvious in the Morrison scheme. By contrast, the mean mass diameter is held as constant in the GCE scheme ( Figure 13c). The outcome is relatively similar to the Z DR CFAD: a larger D mr corresponds to larger Z DR . Similar situations are discovered near the melting layer to 8-km ASL. The larger size of snow and graupel results in higher reflectivity near the melting level to 8-km ASL in the Morrison scheme (Figure 10d). In addition, the K DP CFAD in the GCE scheme ( Figure 10c) is larger than that in the Morrison scheme (Figure 10f) owing to the higher mean rain mixing ratio in the GCE scheme (Figure 13a).  Figure 13. Vertical distribution of the mean mixing ratio (g/kg), total number concentration (/m 3 ), and mean mass diameter (mm) with respect to rain, snow, and graupel hydrometeors at 1100 UTC on 14 June 2008. The area is indicated by the red box in Figure 6. The sky blue, red, and green lines represent rain, snow, and graupel, respectively. The upper and lower rows present the results of the GCE and Morrison schemes, respectively.

Histograms
CFAD analyses provide a general pattern of validating model performance. The simulated polarimetric parameters are evaluated by histograms which provide a qualitative comparison between observations and model simulations. Besides, the precipitation system is divided to convective and stratiform area for further evaluation. The method used herein refers to that in a Figure 13. Vertical distribution of the mean mixing ratio (g/kg), total number concentration (/m 3 ), and mean mass diameter (mm) with respect to rain, snow, and graupel hydrometeors at 1100 UTC on 14 June 2008. The area is indicated by the red box in Figure 6. The sky blue, red, and green lines represent rain, snow, and graupel, respectively. The upper and lower rows present the results of the GCE and Morrison schemes, respectively.

Histograms
CFAD analyses provide a general pattern of validating model performance. The simulated polarimetric parameters are evaluated by histograms which provide a qualitative comparison between observations and model simulations. Besides, the precipitation system is divided to convective and stratiform area for further evaluation. The method used herein refers to that in a previous study [33]. The method is described as follows. (1) The observation data are rearranged from low to high values and the total number of observational data counted. Then, the 0th, 10th, up to the 100th percentile data are selected for grouping criterion. (2) The model data are sorted into 10 groups by the selected criterion in the previous step (i.e., the first group contains the model data whose value is between 0th and 10th, and so on). (3) The model data number in each group is counted and summed to obtain the total model grid point data number. (4) This total model number is divided by 10, and a number for the reference number in each group is obtained. The number in each group should be close to the model reference grid point number. If a certain group number exceeds (is smaller than) the reference number, the variables are overestimated (underestimated). Finally, we determine the analysis performance in each group and whether the analysis has under-or overestimated in the group. For simplicity, we normalize the group number by dividing with reference number. Consequently, we can check much more easily and quickly in each group. Furthermore, by dividing the total observation number with the total model grid point number, the ratio can indicate whether the entire model makes an underor overestimation. If the ratio is close to 1, the analysis field approaches the observation distribution in whole aspect. If particular groups are selected, for example, data over 80% define a strong or convective area. We can employ this method to understand the characteristics of different weather system classification regions. Figure 14 shows the histogram validation of the GCE and Morrison analysis fields in combination with three heights ASL: 1.5, 2.5, and 3.5 km. Three heights are selected for the reason of increasing the sample, and the warm rain area structures are investigated. The performance of Z H in both scheme analysis fields is almost identical, and the distributions in each interval are overestimated slightly. In the Morrison scheme, former 80% intervals behave much better than in the GCE scheme. Furthermore, the first 10% interval indicates that the Morrison scheme has the flexibility to simulate the stratiform area DSD, which has been reported by many studies [13,31]. However, both schemes underestimate the convection area, specifically the reflectivity echo at higher than 40 dBZ, showing that the analysis fields of both schemes still require improvement. In the Z DR field, both schemes overestimate the value, and the last interval has 7.68 and 8.47 times the reference grid points. The reason for the overestimation of Z DR is mentioned in Section 4.3. In the K DP field, the two schemes underestimate the amount, indicating that the liquid water content is much lower compared with in the observations. Underestimation of liquid water content in model simulation was also reported in another study [33]. When checking the last 10% of data with values higher than 0.42 • /km, representing the convection area, the GCE scheme performs better than Morrison scheme. This result is mainly due to the rain mixing ratio and K DP being strongly related to the liquid water content. The DA system can provide important distributions, even though we only assimilated radial wind and reflectivity. much lower compared with in the observations. Underestimation of liquid water content in model simulation was also reported in another study [33]. When checking the last 10% of data with values higher than 0.42°/km, representing the convection area, the GCE scheme performs better than Morrison scheme. This result is mainly due to the rain mixing ratio and being strongly related to the liquid water content. The DA system can provide important distributions, even though we only assimilated radial wind and reflectivity.

Discussion
In this case study, sensitivity tests illustrated that for evaluating the performance of model simulation objectively, it is necessary to setup the threshold of parameters in the dual-polarimetric

Discussion
In this case study, sensitivity tests illustrated that for evaluating the performance of model simulation objectively, it is necessary to setup the threshold of parameters in the dual-polarimetric observation operator. The studies [42,68] indicate that the DSD characteristic in Taiwan is closer to maritime precipitation, which featuring high number concentration and medium drop size. After excluding total number concentration <100/m 3 of model grid point in mixing-phase area, the simulated bright band level is smoother and closer to S-POL observation.
Examining the performance of numerical simulations, the analysis fields of GCE and Morrison schemes show similar structure and intensity of squall lines to radar reflectivity after radar data is assimilated. Results show that simulated Z H of GCE (single moment) could be as good as the Morrison (double moment) scheme (Figures 7a and 8a,d,g), and simulated Z DR of GCE has less bias than that of the Morrison scheme. In addition, the vertical cross-sections of Z DR and K DP illustrates that the polarimetric signatures in the main convection area could be replicated in both schemes. However, compared to the simulated radar parameters with the S-POL CFAD, both schemes overestimated Z DR , especially the Morrison scheme. This is due to the hydrometeor variables of mixing ratio and number concentration is dependent in the single moment GCE scheme. This is confirmed when the increment of model variables and simulated dual-polarimetric parameters are examined (Figures 11 and 12). The simulation with Morrison scheme, whose total number concentration is prognostic, shows an unrealistic results comparing to retrieved total number concentration field (not shown). Besides, the distribution of the K DP CFAD in GCE more resemble to S-POL observation than that in Morrison scheme when both scheme underestimate the 75th percentile. Overall, results indicate that if double moment scheme is desired to be applied in the numerical model simulation, dual-polarimetric parameters, such as Z DR and K DP , are essential to be assimilated for obtaining better hydrometeor variables.
Further investigation on the vertical distribution of mixing ratio, total number concentration and mean mass diameter, the results show the characteristics in microphysics schemes. The comparison of the mean mass diameter revealed that the mean drop size was greater in the Morrison scheme ( Figure 13). These results explain the larger bias of Z DR CFAD in the analysis with Morrison scheme (Figure 9 versus Figure 10). In addition, the mean mass diameter in Morrison is larger than that in GCE while the amount of rain mixing ratio in GCE is larger than that in Morrison ( Figure 13). Therefore, K DP CFAD in the analysis with GCE scheme has greater value than the Morrison scheme (Figure 10c versus Figure 10f). By separating the squall line system into convective and stratiform areas through histograms (Figure 14), different characteristics of physical processes in the two regions could be identified. The histogram results revealed that the Morrison scheme can simulate the stratiform more precisely than the GCE scheme. Conversely, the GCE scheme exhibited superior performance in convective areas compared with the Morrison scheme, because the data distribution was closer to the observed distribution. One should notice that whereas the Z DR and K DP value are compared with S-Pol CFAD, we mainly focused on the performance of warm-rain in the analysis. The reason is that the assumption of polarimetric operator on the shape of snow and graupel/hail particles resulted in the difficulties in simulating the Z DR and K DP value correctly.

Conclusions
In this study, a sophisticated dual-polarimetric observation operator for radar data was applied to Taiwan. A squall line case during the SoWMEX in 2008 was investigated, and the results of the convection-permitting model were verified using NCAR S-POL data. To accurately capture the precipitation system, radar reflectivity and radial wind were assimilated in the WRF model at 3-km resolution. In addition, two microphysics schemes were applied to evaluate the performance of simulations using radar observations.
In this application of the polarimetric radar observation operator in Taiwan, sensitivity tests were conducted to obtain the proper settings for this weather system in Taiwan. When examining the Morrison scheme, the maximum conversion fraction of mixing hydrometers was discovered to be the minor reason for discontinuity in the CFAD reflectivity, and adjusting the total number concentration played a major role in improving the bright band structure. The sensitivity test shows that while calculating the total number concentration of mixing-phase hydrometeors, a particular threshold can prevent unrealistic and discontinuous bright band structure happening. The proper threshold (i.e., 100/m 3 ) was given and used in the remainder of the experiments for model verification.
Based on the proper setup of the operator in this event, more evaluations were performed in the two microphysics schemes when comparing the model simulation with dual-polarimetric variables. The structures of reflectivity were similar between the GCE and Morrison schemes, and both schemes resembled the S-POL CFAD. However, both schemes overestimated Z DR compared to the observations, especially for the Morrison scheme. The investigation of the increment of total number concentration further revealed that assimilating radial wind and reflectivity data could only modify number concentration slightly. Moreover, the comparison of the mean mass diameter illustrated that the mean drop size was greater in the Morrison scheme. Therefore, the Z DR CFAD had larger bias in the Morrison scheme. The histogram of polarimetric variables provided a quantitative comparison and also described the difference of GCE and Morrison scheme by separating convective and stratiform area. The Morrison scheme matched the stratiform area of observation more precisely than the GCE scheme. Whereas, the GCE scheme performed better in convective areas compared with the Morrison scheme, owing to K DP distribution was closer to the observation.
In one case study, it was discovered that assimilating radial wind and reflectivity data may not be sufficient for enabling the optimal mesoscale analysis, especially for a double moment microphysics scheme. In future work, dual-polarimetric parameters will be assimilated, and the performance of numerical models both for analysis and short-term forecasting will be examined. In addition,