Sensitivities of Quantitative Precipitation Forecasts for Typhoon Soudelor (2015) near Landfall to Polarimetric Radar Data Assimilation

: Based on the preciousness and uniqueness of polarimetric radar observations collected near the landfall of Typhoon Soudelor (2015), this study investigates the sensitivities of very short-range quantitative precipitation forecasts (QPFs) for this typhoon to polarimetric radar data assimilation. A series of experiments assimilating various combinations of radar variables are carried out for the purpose of improving a 6 h deterministic forecast for the most intense period. The results of the control simulation expose three sources of the observation operator errors, including the raindrop shape-size relation, the limitations for ice-phase hydrometeors, and the melting ice model. Nevertheless, polarimetric radar data assimilation with the unadjusted observation operator can still improve the analyses, especially rainwater, and consequent QPFs for this typhoon case. The di ﬀ erent impacts of assimilating reﬂectivity, di ﬀ erential reﬂectivity, and speciﬁc di ﬀ erential phase are only distinguishable at the lower levels of convective precipitation areas where speciﬁc di ﬀ erential phase is found most helpful. The positive e ﬀ ect of radar data assimilation on QPFs can last three hours in this study, and further improvement can be expected by optimizing the observation operator in the future


Introduction
Enhancing the accuracy of quantitative precipitation forecasts (QPFs) has been a relentless goal and challenge for worldwide weather research and operational units. Accurate QPFs are crucial for various downstream hydrological applications, which can serve as yardsticks for governments to anticipate measures against potential water shortage or heavy rainfall hazards. In pursuit of this goal, the appropriate forecasting method can vary at different temporal and spatial scales. Taking tropical cyclones that approach the island of Taiwan for example, global numerical weather prediction (NWP) models are the first tool to predict their cyclogeneses and subsequent tracks within a medium range of 3-10 days [1,2]. After they intensify into typhoons and move closer to land, regional NWP models equipped with sophisticated physics can weigh in to provide convective-scale forecasts within a short range of 0.5-3 days [3,4]. Both global and regional NWP models require data assimilation components to optimize initial conditions by assimilating routine meteorological observations. For regional NWP models, weather radars are now widely used to provide the richest wind and hydrometeor information of precipitation systems. Thus, once typhoon rainbands impinge on the land of Taiwan, regional NWP models with radar data assimilation are the mainstream method for QPFs within a very short range of 2-12 h [5,6]. As for precipitation nowcasts within 0-2 h, radar echo extrapolation is currently a more Inspired by these recent successes, which focus on supercell and MCS cases whose scales span few to dozens of kilometers, this study explores for the first time the impact of polarimetric radar data assimilation for a tropical cyclone case that spans hundreds of kilometers. The goal is to examine how the polarimetric radar observation operator developed by the authors of Reference [45] behaves in the tropical cyclone case and whether the assimilation of polarimetric variables further benefits very short-range QPFs, which are crucial for early warning against heavy rainfall hazards. The studied case is Typhoon Soudelor (2015), which has been the most completely observed typhoon by the only ground-based S-band polarimetric radar in Taiwan up to the present. The antenna and radome were even damaged by its gusts and under repair for more than one year. Afterward, no other typhoons have impinged on northern Taiwan as severely as Typhoon Soudelor, and therefore, the polarimetric radar observations near the landfall of Typhoon Soudelor are precious and unique in this study. This study uses a DM MP scheme and incorporates the observation operator of Reference [45] into a WRF local ensemble transform Kalman filter (LETKF [64]) radar assimilation system to investigate the sensitivities of 6 h QPFs to assimilating polarimetric variables. The 6 h forecast length is the temporal scale at which the assimilation of a single S-band radar could possibly influence the forecast in the typhoon scenario according to a previous OSSE study [5]. In Section 2, the model and assimilation configurations of this system are described as well as the detailed procedures of the observation operator and data preprocessing. In Section 3, the case of Typhoon Soudelor is introduced, and the experimental design of polarimetric radar data assimilation is presented. Section 4 is devoted to the results for all the experiments. Section 5 provides in-depth discussions on the sources of the observation operator errors. Section 6 is devoted to the summary and future prospects.

Model and Assimilation Configurations
The WRF-LETKF radar assimilation system (WLRAS) was originally developed and evaluated through OSSEs for a typhoon scenario in Reference [5], which demonstrated the capability of this system to improve very short-range QPFs by assimilating V r and Z H . In this study, WLRAS is employed again but coupled with a later version (V4.0) of the Advanced Research WRF model and the observation operator of Reference [45] to address the challenge of real case polarimetric radar data assimilation. Referring to the configuration of an operational WRF run from Taiwan's Central Weather Bureau (CWB), this study uses two grids with two-way nesting, where the fine grid covers the whole observation area of Taiwan's weather radar network. Figure 1 shows the simulation domains of both grids, spanning 280 × 280 and 331 × 331 horizontal grid points at 15 km and 3 km resolutions, respectively. There are 45 vertical levels with a top of 30 hPa. The physics schemes in use include the WRF double-moment 6-class MP scheme (WDM6 [65]), rapid radiative transfer model (RRTM) model (RRTM) longwave radiation scheme [66], Goddard shortwave radiation scheme [67], fifth generation Penn State University and the National Center for Atmospheric Research mesoscale model (MM5) surface layer scheme [68], Noah land surface model [69], Yonsei University planetary boundary layer scheme [70], and Kain-Fritsch cumulus parameterization scheme ( [71]; applied to the coarse grid only). The WDM6 scheme, which contains three ice categories (cloud ice, snow, and graupel), is evolved from the WSM6 scheme by adding the prediction of number concentrations for cloud water, rainwater, and cloud condensation nuclei. Hence, prognostic state variables in this study contain the three-dimensional wind components (u, v, and w), perturbation potential temperature ( ), perturbation geopotential ( ), perturbation surface pressure of dry air ( ), mixing ratios of water vapor ( ), cloud water ( ), rainwater ( ), cloud ice ( ), snow ( ), and graupel ( ), and number concentrations of cloud water ( ), rainwater ( ), and cloud condensation nuclei ( ). The LETKF assimilation scheme is a variant of the square-root filter (also known as the deterministic EnKF), which does not require perturbing observations at analysis steps. Its distinguishing feature is that every grid point can be analyzed independently by assimilating surrounding observations within a covariance localization radius simultaneously. The complete Figure 1. The simulation domains of the two grids with two-way nesting, spanning 280 × 280 and 331 × 331 horizontal grid points at 15 km and 3 km resolutions, respectively. The green dot and circle are respectively the location and maximum unambiguous range of the Wufenshan (RCWF) radar. The blue dashed line is the best track of Typhoon Soudelor analyzed by CWB, where the red solid segment represents the simulation period (1200 UTC 7 August-0200 UTC 8 August) in this study.
The LETKF assimilation scheme is a variant of the square-root filter (also known as the deterministic EnKF), which does not require perturbing observations at analysis steps. Its distinguishing feature is that every grid point can be analyzed independently by assimilating surrounding observations within a covariance localization radius simultaneously. The complete algorithm can be found in Reference [64]. A noteworthy upgrade is made to the parallel computing framework of LETKF used in Reference [5] to improve assimilation efficiency for dense radar data. Suppose there are P processors available. In Reference [5], the grid to be analyzed was horizontally divided into P rectangular sub-grids of the same size, computed by separate processors. In this study, the grid to be analyzed is cut into y-z slices along the x direction, and the P processors take turns computing one slice at a time; that is, processor 1 computes slices 1, P + 1, 2P + 1, and so on. Both parallel computing frameworks give comparable efficiency if observations are evenly spread around the whole grid. However, if observations cluster in small areas as radar observations usually do, the latter framework gives higher efficiency for more balanced workload assigned to each processor. The configuration of LETKF in this study is largely inherited from Reference [5], which is optimized for the typhoon scenario. LETKF updates all prognostic state variables except µ , N tc , N tr , and N tn . µ is excluded because of the difficulty in updating a two-dimensional vertically integrated field. N tc , N tr , and N tn are left un-updated due to the concern about their highly non-Gaussian error statistics. Instead, these four variables are dynamically adjusted by the updated variables at the forecast steps of assimilation cycles. Applying the mixed localization method, which was first proposed in Reference [5], the horizontal covariance localization radii are 36 km for u and v and 12 km for the other updated variables to reflect the larger-scale error structure of horizontal winds in the typhoon scenario. The vertical covariance localization radius is uniformly 4 km. An ensemble size of 40 is used for limited computing resources, and a fixed multiplicative covariance inflation factor of 1.08 is used for the systematic underestimation of ensemble-based error covariances.

Observation Operator and Data Preprocessing
In Reference [5], the observation operator of WLRAS interpolated state variables from grid points to observation points and then converted them to V r and Z H . Following Reference [72], the conversion to V r summed up the projections of u, v, w, and rainwater terminal velocity (v t ) on the radar beam, and the conversion to Z H only considered the contribution from rainwater (q r ) with an assumed Marshall-Palmer DSD [73], as where x, y, and z are Cartesian coordinates with the origin at the radar site, and ρ a is the density of air. In this study, the observation operator of Reference [45] replaces the original Z H operator (Equation (2)) to additionally consider the contribution from ice-phase hydrometeors, the bright band signature of the melting layer, and conversions to Z DR and K DP . Firstly, a melting ice model is used to simulate the bright band signature resulting from mixed-phase hydrometeors that most bulk MP schemes do not deal with. Mixed-phase hydrometeors, such as rain-snow, rain-graupel, and rain-hail mixtures, are generated at a grid point by collecting a fraction of mass from rainwater and ice-phase hydrometeors if both coexist at that point, as where x denotes the ice-phase hydrometeor species, and F max is the maximum fraction, which empirically depends on the species. Secondly, backscattering amplitudes along horizontal and vertical axes are simulated as power-law functions of hydrometeor size that fit T-matrix scattering calculations [74] for rainwater and Rayleigh scattering calculations for ice-phase and mixed-phase hydrometeors. Reflectivity factors at horizontal and vertical polarizations (Z h and Z v ) contributed from each hydrometeor species are then calculated by integrating backscattering cross-sections over a gamma DSD, and finally expressed as where λ is the radar wavelength, K w is the dielectric constant of water, q and ρ are respectively the mixing ratio and density of the hydrometeor species, N 0 and µ are respectively the intercept and shape parameters of the gamma DSD, α and β are respectively the constant and exponent of the fitted power-law functions with subscripts a and b denoting horizontal and vertical axes, and A, B, and C are coefficients that redeem the canting behavior of falling hydrometeors. Thirdly, K DP contributed from each hydrometeor species is calculated by integrating the real part of the difference between forward scattering amplitudes along horizontal and vertical axes over the gamma DSD, and finally expressed as where C k is also a coefficient that redeems the canting behavior. It is noteworthy that q is prognostic in all SM, DM, and TM bulk MP schemes. N 0 is constant in SM schemes or diagnostic in DM and TM schemes as where N t is prognostic in DM and TM schemes, and µ is constant in DM schemes or diagnostic in TM schemes with prognostic reflectivity. Hence, Equations (4)-(7) are solvable and applicable to all bulk MP schemes. Lastly, the contributions from all the hydrometeor species observable by weather radars to Z h , Z v , and K DP are summed up to yield total values, in which total Z h and Z v yield Z H and Z DR . Taking the WDM6 scheme, for example, K DP = K DP,r + K DP,s + K DP,g + K DP,rs + K DP,rg .
Besides the refined observation operator, appropriate data preprocessing and quality control to reduce observation errors also play an important role in real-case data assimilation. The polarimetric radar data used in this study come from the RCWF S-band polarimetric radar of CWB, located on the top of Mount Wufen (121.7725 • E, 25.0728 • N, 770 m altitude) on the northern coast of Taiwan ( Figure 1). This radar is identical to the WSR-88D network in the U.S. with a maximum unambiguous range of 230 km, a volume scan interval of 7.5 min, and nine plan position indicator (PPI) elevation angles ranging from 0.5 • to 19.5 • . Quality control procedures, such as removing non-meteorological echoes, dealiasing V r data, and correcting signal attenuation, are fundamental to the preprocessing of raw radar data. Firstly, a high-resolution terrain file is used to filter out primary ground clutter, and then the remaining echoes are removed if ρ HV < 0.8, an empirical value for the RCWF radar to distinguish between meteorological and non-meteorological echoes. Secondly, V r data are dealiased automatically with an algorithm that checks temporal and spatial continuity. Failed dealiasing is then manually unfolded. Thirdly, Φ DP data are dealiased as well and then smoothed to estimate K DP values, which are discarded if negative. Lastly, attenuation correction is applied to Z H and Z DR data with a single-coefficient Φ DP -based method [75], although the correction is almost negligible for the S band. Figure 2 shows the Z H , Z DR , K DP , and ρ HV PPI observations, for the case of Typhoon Soudelor near landfall, at the lowest elevation angle posterior to the foregoing quality control procedures. It can be seen that most data-void areas and remaining Z DR noise are in the southwest quadrant where the higher Central Mountain Range is situated. The RCWF radar observations are thinned and interpolated to 19 constant altitude plan position indicator (CAPPI) altitudes ranging from 1 to 10 km at a 3 km horizontal resolution (corresponding to the fine grid) and 0.5 km vertical resolution prior to data assimilation. Observation errors are given as 1 m s −1 for V r , 2 dBZ for Z H , 0.2 dB for Z DR , and 0.5 • km −1 for K DP referring to Reference [47].
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 26 from 1 to 10 km at a 3 km horizontal resolution (corresponding to the fine grid) and 0.5 km vertical resolution prior to data assimilation. Observation errors are given as 1 m s −1 for , 2 dBZ for , 0.2 dB for , and 0.5° km −1 for referring to Reference [47].

Case of Typhoon Soudelor
Typhoon Soudelor was the most devastating tropical cyclone in the Northwest Pacific region in 2015, resulting in dozens of fatalities and a tremendous loss of property on its way through the Northern Mariana Islands, Taiwan, and eastern China. Its westbound best track analyzed by CWB is shown in Figure 1, where the red solid segment represents the simulation period (1200 UTC 7 August-0200 UTC 8 August) in this study. This tropical cyclone reached tropical storm intensity on July 30, achieved Category 5 intensity with a minimum central pressure of 900 hPa and 10 min maximum sustained winds of 59.7 m s −1 on August 3, and struck the island of Taiwan with

Case of Typhoon Soudelor
Typhoon Soudelor was the most devastating tropical cyclone in the Northwest Pacific region in 2015, resulting in dozens of fatalities and a tremendous loss of property on its way through the Northern Mariana Islands, Taiwan, and eastern China. Its westbound best track analyzed by CWB is shown in Figure 1, where the red solid segment represents the simulation period (1200 UTC 7 August-0200 UTC 8 August) in this study. This tropical cyclone reached tropical storm intensity on July 30, achieved Category 5 intensity with a minimum central pressure of 900 hPa and 10 min maximum sustained winds of 59.7 m s −1 on August 3, and struck the island of Taiwan with Category 3 intensity on August 7-8. During the hours around landfall on the eastern coast of Taiwan at 2040 UTC 7 August, strong gusts and rainbands north of the typhoon center severely impinged on northern Taiwan (Figure 2a) and damaged the RCWF radar by blowing its radome and antenna away. Nevertheless, Typhoon Soudelor was the first typhoon well observed by the RCWF radar after its polarimetric upgrade, and its valuable data obtained before the damage motivated this study.
At 2000 UTC 7 August, the eye of Typhoon Soudelor was located 110 km south of the RCWF radar ( Figure 2a). The northern semicircle of Typhoon Soudelor, the so-called dangerous semicircle, consisted of several huge rainbands with Z H reaching 50 dBZ. In contrast, fewer rainbands were evident and only close to the eye in the southern semicircle, and the most inner one featured a significant amount of large raindrops distinguished by high Z DR and K DP values reaching 2.7 dB and 1.8 • km −1 , respectively (Figure 2b,c). Z DR and K DP values inside the northern huge rainbands were lower over the ocean but increased over the land, especially K DP . This indicates that orographic lifting affected rainbands moving toward terrain and resulted in stronger convection and precipitation. ρ HV values were mostly higher than 0.98 within a range of 170 km but decrease beyond that (Figure 2d). This indicates that the melting layer was located at a 4 km altitude where mixed-phase hydrometeors existed. With respect to subsequent rainfall (Figure 3), three heavy rainfall areas (labeled A-C) were observed on the windward slopes of northern Taiwan with 3 h rainfall exceeding 150 mm during 2000-2300 UTC 7 August. Their locations correspond well to each rainband. During the next three hours, the rainfall decreased in northern Taiwan but increased in central Taiwan (area D) as the typhoon center gradually moved through central Taiwan. A 6 h maximum rainfall of 393 mm was reached in area A where massive disasters including floods, landslides, road closures, and the disruption of power and water supply occurred.

Experimental Design of Polarimetric Radar Data Assimilation
For the case of Typhoon Soudelor near landfall, a series of experiments are carried out to investigate the sensitivities of 6 h QPFs to polarimetric radar data assimilation. Figure 4 shows the schematic design of these experiments, aiming at improving a 6 h deterministic precipitation forecast for the most intense period since 2000 UTC 7 August. NoDA (no data assimilation) is a control experiment, which only performs a cold-start simulation initialized with the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) 0.25 degree forecast at 1200 UTC 7 August. The initial eight and final six hours of the simulation are regarded as a model spin-up period and benchmark forecast, respectively. DA represents all the experiments that assimilate RCWF observations with the LETKF scheme, using a 40-member ensemble perturbed around the NoDA initial condition via the random-cv (random control variables) facility of the WRF 3DVar system. The initial ensemble goes through a 6 h spin-up period and nine assimilation cycles at a 15 min interval, and then the analysis ensemble mean is used to generate the 6 h deterministic forecast. Radar data assimilation only takes place in the fine grid, whose analysis increments are fed back to the coarse grid via two-way interaction. The predicted rainfalls on land in all the experiments are verified against the observed rainfall, interpolated from more than 750 surface rain gauges to the fine grid points. The locations of all the rain gauges are shown in Figure 5 as well as terrain contours and the verification area, the land within the 230 km maximum unambiguous range of the RCWF radar. The temporal resolution of rain gauge data is 10 min, and therefore, hourly accumulations can be calculated.  Table 1). The green circle is the maximum unambiguous range of the RCWF radar. Uppercase A, B, C, and D represent individual heavy rainfall areas, among which area A reaches a 6 h maximum rainfall of 393 mm. The black rectangle highlights areas A-C for easier visual comparison.

Experimental Design of Polarimetric Radar Data Assimilation
For the case of Typhoon Soudelor near landfall, a series of experiments are carried out to investigate the sensitivities of 6 h QPFs to polarimetric radar data assimilation. Figure 4 shows the schematic design of these experiments, aiming at improving a 6 h deterministic precipitation  Table 1). The green circle is the maximum unambiguous range of the RCWF radar. Uppercase A-D represent individual heavy rainfall areas, among which area A reaches a 6 h maximum rainfall of 393 mm. The black rectangle highlights areas A-C for easier visual comparison.
back to the coarse grid via two-way interaction. The predicted rainfalls on land in all the experiments are verified against the observed rainfall, interpolated from more than 750 surface rain gauges to the fine grid points. The locations of all the rain gauges are shown in Figure 5 as well as terrain contours and the verification area, the land within the 230 km maximum unambiguous range of the RCWF radar. The temporal resolution of rain gauge data is 10 min, and therefore, hourly accumulations can be calculated.  This study focuses on the sensitivities of QPFs to the assimilation of different radar variables, and therefore various combinations of radar variables are constructed for the DA experiments, as listed in Table 1. Uppercase V, Z, D, and K in the experiment names stand for assimilating , , , and , respectively. Because is linearly correlated to three-dimensional wind components and consequently beneficial for their analyses, as proved in Reference [5], it is assimilated in all the DA experiments. For , , and , only positive values are assimilated. Negative-value observations are rejected by the data assimilation system to avoid adverse impact because negative is around the limit of detection and negative and are against the assumption of oblate raindrops in the observation operator. The results of all the experiments are shown next in three This study focuses on the sensitivities of QPFs to the assimilation of different radar variables, and therefore various combinations of radar variables are constructed for the DA experiments, as listed in Table 1. Uppercase V, Z, D, and K in the experiment names stand for assimilating V r , Z H , Z DR , and K DP , respectively. Because V r is linearly correlated to three-dimensional wind components and consequently beneficial for their analyses, as proved in Reference [5], it is assimilated in all the DA experiments. For Z H , Z DR , and K DP , only positive values are assimilated. Negative-value observations are rejected by the data assimilation system to avoid adverse impact because negative Z H is around the limit of detection and negative Z DR and K DP are against the assumption of oblate raindrops in the observation operator. The results of all the experiments are shown next in three parts: (1) the NoDA simulation and observation operator is assessed by examining the horizontal pattern and vertical distribution of simulated radar variables, (2) the analyzed radar variables and state variables in experiments VZ, VD, and VK are compared to distinguish the different impacts of assimilating Z H , Z DR , and K DP , and (3) the sensitivities of QPFs to different assimilated radar variables are revealed qualitatively and quantitatively.

Assessment of NoDA Simulation and Observation Operator
The NoDA simulation is first explored by comparing the forecast of Z H , Z DR , and K DP CAPPI fields with the RCWF radar observations at a 3.5 km altitude at 2000 UTC 7 August (Figure 6). A black rectangle is added to highlight the aforementioned areas A-C for easier visual comparison. Horizontally, NoDA exhibits significant differences of the rainband pattern and intensity from the observations. The simulated rainbands have position errors and look more discrete while the observed rainbands are connected with weaker echoes between. The intensity inside the simulated rainbands is generally overestimated for Z DR and K DP . Despite these inconsistencies, the observed orographic lifting effect that strengthens convection and precipitation in areas A-C is still visible in the simulation. For overall vertical distribution, the contoured frequency by altitude diagrams (CFADs) with the median curve for Z H , Z DR , and K DP in convective precipitation areas (CPAs; defined by Z H > 40 dBZ at a 3 km altitude) within the 230 km maximum unambiguous range at 2000 UTC 7 August are compared between the RCWF radar observations and NoDA (Figure 7). A blue rectangle is added to highlight the majority of the observed values at lower levels (below 5 km) for easier visual comparison, with ranges of 35-50 dBZ, 0.7-1.4 dB, and 0.2-0.7 • km −1 , respectively. The medians of the observed Z H , Z DR , and K DP all exhibit an increasing trend as the altitude decreases, which reflects the size-sorting mechanism that larger raindrops mostly exist at lower levels.
The CFADs in stratiform precipitation areas (SPAs; defined by 10 dBZ < Z H < 30 dBZ at a 3 km altitude) are also examined (Figure 8), although typhoon rainfall is dominated by the rainbands located in CPAs. The observed Z H , Z DR , and K DP values are similar to those in CPAs at higher levels but noticeably lower below, and the layer of the size-sorting mechanism is shallower (below 2.5 km). This result indicates that the absolute majority of raindrops in SPAs are smaller and more spherical. The NoDA simulation in SPAs has better consistency with the observations than in CPAs at lower levels, supporting Reference [76], that found the theoretical oblateness more accurate for smaller raindrops in typhoon cases. At higher levels in both CPAs and SPAs, the limitations of the observation operator for ice-phase hydrometeors are seen. It is noteworthy that the bright band signature, supposed to be clear near the melting layer of SPAs, is visible in the observations but excessive in the simulated Z H and Z DR , especially an unrealistic Z DR spike between 4 and 6 km altitudes. This excess is proved attributable to the melting ice model of the observation operator because it vanishes if the model is ignored (Figure 9). Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 26  The CFADs in stratiform precipitation areas (SPAs; defined by 10 dBZ < < 30 dBZ at a 3 km altitude) are also examined (Figure 8), although typhoon rainfall is dominated by the rainbands located in CPAs. The observed , , and values are similar to those in CPAs at higher levels  Figure 7, but in stratiform precipitation areas (SPAs). Figure 8. As in Figure 7, but in stratiform precipitation areas (SPAs).

Comparison of the VZ, VD, and VK Analyses
To distinguish the different impacts of assimilating Z H , Z DR , and K DP , the final analyzed radar variables in experiments VZ, VD, and VK are compared with the foregoing NoDA simulation (Figures 6-8). From the CAPPI fields at a 3.5 km altitude (Figure 6), all the three DA experiments correct the discrete rainband pattern in NoDA, but not all of them correct the overestimated intensity for Z DR and K DP . Taking a closer look at the black rectangle highlighting areas A-C, the overestimation of Z DR is unchanged in VZ and VD but slightly mitigated in VK. In contrast, the overestimation of K DP is better corrected in VZ and VK, the latter of which exhibits a very similar pattern to the observations. From the CFADs in CPAs (Figure 7), all three DA experiments rectify the excessive increasing trend for K DP at lower levels rather than for Z DR , and VK shows the most realistic dispersion for the whole spectra of Z H and K DP in comparison to the observations. From the CFADs in SPAs (Figure 8), the three DA experiments resemble one another closely and correct the unrealistic bright band signature for Z H and Z DR . These results about the analyzed radar variables can be generalized that, with the unadjusted observation operator, polarimetric radar data assimilation is still beneficial for analyses in this typhoon case. Besides, the different impacts of assimilating Z H , Z DR , and K DP are only distinguishable at the lower levels of CPAs where larger raindrops mostly exist, with a ranking order of K DP , Z H , and Z DR from best to worst. To further understand the variations of the analyzed state variables in the WRF model, Figure 10 shows the horizontal averages of w, N tr , q v , q r , q c , q i , q s , and q g by altitude in CPAs within the 230 km maximum unambiguous range at 2000 UTC 7 August for experiments NoDA, V, VZ, VD, and VK. The results reveal that the NoDA simulation has particularly high values for w, q r , q i , q s , and q g , while the variations among the other DA experiments are more obvious for w, q r , q c , and q g . As q r is the state variable most related to rainfall, its ranking order of VK, VZ, V, VD, and NoDA from lowest to highest is to be compared with the skills of QPFs discussed next.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 26 Figure 9. As the NoDA forecasts in Figures 7 and 8, but the melting ice model of the observation operator is neglected.

Comparison of the VZ, VD, and VK Analyses
To distinguish the different impacts of assimilating , , and , the final analyzed radar variables in experiments VZ, VD, and VK are compared with the foregoing NoDA simulation (Figures 6-8). From the CAPPI fields at a 3.5 km altitude (Figure 6), all the three DA experiments correct the discrete rainband pattern in NoDA, but not all of them correct the overestimated intensity for and . Taking a closer look at the black rectangle highlighting areas A-C, the overestimation of is unchanged in VZ and VD but slightly mitigated in VK. In contrast, the overestimation of is better corrected in VZ and VK, the latter of which exhibits a very similar pattern to the observations. From the CFADs in CPAs (Figure 7), all three DA experiments rectify the excessive increasing trend for at lower levels rather than for , and VK shows the most realistic dispersion for the whole spectra of and in comparison to the observations. From the CFADs in SPAs (Figure 8), the three DA experiments resemble one another closely and correct the unrealistic bright band signature for and . These results about the analyzed radar variables can be generalized that, with the unadjusted observation operator, polarimetric radar data assimilation is still beneficial for analyses in this typhoon case. Besides, the different impacts of

Sensitivities of QPFs to Different Assimilated Radar Variables
The predicted rainfalls during the first three, second three, and total six hours of deterministic forecasts in experiments NoDA, VZ, VD, and VK are compared with the aforementioned observed rainfall (Figure 3). The black rectangle highlights areas A-C for easier visual comparison. During the first three hours, the patterns of areas A-C in the three DA experiments resemble the observations more than in NoDA. Outside the rectangle, the observed 3 h rainfall is generally lower than 50 mm (in blue), while NoDA and VD have more areas higher than 50 mm (in green and yellow) than VZ and VK. During the second three hours, all the NoDA and DA experiments have similar underestimates in areas A-C and overestimates in area D. This result implies that the positive effect of radar data assimilation on QPFs can last three hours in this study. As for the total 6 h accumulations, the improved pattern of areas A-C is diluted but still distinguishable in the three DA experiments. For quantitative verification, the root-mean-square error (RMSE), spatial correlation coefficient (SCC), and equitable threat score (ETS) on a threshold of 15 mm for the predicted hourly rainfall on land within the 230 km maximum unambiguous range during 2000 UTC 7 August-0200 UTC 8 August are compared in two groups of experiments ( Figure 11). All the curves appear to be separate in the first three hours but tangled in the second three hours, and therefore reflect the 3 h positive impact of radar data assimilation again. Focusing on the first three hours, the ranking order from best (lowest RMSE, highest SCC, or highest ETS) to worst for the first group is VK, VZ, V, VD, and NoDA. This order is the same as the foregoing order from lowest to highest q r , which suggests that QPFs for this typhoon case highly depend on rainwater analyses. Moreover, the ranking order from best to worst for the second group is VZK, VZDK, VZD, VZ, and NoDA. This result suggests that, when V r and Z H are already assimilated, further QPF improvement is still achievable by the additional assimilation of polarimetric variables, especially K DP .

Sensitivities of QPFs to Different Assimilated Radar Variables
The predicted rainfalls during the first three, second three, and total six hours of deterministic forecasts in experiments NoDA, VZ, VD, and VK are compared with the aforementioned observed rainfall (Figure 3). The black rectangle highlights areas A-C for easier visual comparison. During the first three hours, the patterns of areas A-C in the three DA experiments resemble the observations more than in NoDA. Outside the rectangle, the observed 3 h rainfall is generally lower than 50 mm (in blue), while NoDA and VD have more areas higher than 50 mm (in green and yellow) than VZ and VK. During the second three hours, all the NoDA and DA experiments have similar underestimates in areas A-C and overestimates in area D. This result implies that the positive effect of radar data assimilation on QPFs can last three hours in this study. As for the total 6 h accumulations, the improved pattern of areas A-C is diluted but still distinguishable in the three DA experiments. For quantitative verification, the root-mean-square error (RMSE), spatial correlation coefficient (SCC), and equitable threat score (ETS) on a threshold of 15 mm for the predicted hourly rainfall on land within the 230 km maximum unambiguous range during 2000 UTC 7 August-0200 UTC 8 August are compared in two groups of experiments ( Figure 11). All the curves appear to be separate in the first three hours but tangled in the second three hours, and therefore reflect the 3 h positive impact of radar data assimilation again. Focusing on the first three hours, the ranking order from best (lowest RMSE, highest SCC, or highest ETS) to worst for the first group is VK, VZ, V, VD, and NoDA. This order is the same as the foregoing order from lowest to highest , which suggests that QPFs for this typhoon case highly depend on rainwater analyses. Moreover, the ranking order from best to worst for the second group is VZK, VZDK, VZD, VZ, and NoDA. This result suggests that, when and are already assimilated, further QPF improvement is still achievable by the additional assimilation of polarimetric variables, especially .

Discussion
When examining the vertical distribution by CFADs in Figure 7, the increasing trend in NoDA is excessively simulated for and , especially , in accord with the overestimation shown in CAPPI. This serious overestimation inferentially originates from an unsuitable raindrop shape-size relation applied by the observation operator rather than a model wet bias for the two following reasons. One is that the overestimation is maximal for and minimal for

Discussion
When examining the vertical distribution by CFADs in Figure 7, the increasing trend in NoDA is excessively simulated for Z DR and K DP , especially Z DR , in accord with the overestimation shown in CAPPI. This serious overestimation inferentially originates from an unsuitable raindrop shape-size relation applied by the observation operator rather than a model wet bias for the two following reasons. One is that the overestimation is maximal for Z DR and minimal for Z H respectively, having the most direct and indirect relations to the raindrop shape among the three variables. The other is that, according to Reference [76], examining in-situ 2D video disdrometer observations for 13 typhoons that struck Taiwan, the raindrops larger than 1.5 mm under high horizontal winds tend to be more spherical than the theoretical ones in References [37] and [77], referred to by Reference [45] and in this study.
At higher levels (above 5 km), the observed Z H , Z DR , and K DP values generally become lower, except for a few outliers. This result corresponds with the observations in Reference [78] for the case of Typhoon Nida (2016), asserting that lower values are mainly contributed from widespread snow, while a few higher values are contributed from supercooled rainwater and mixed-phase hydrometeors in updraft regions. However, the simulated Z DR and K DP are seriously underestimated with overwhelming near-zero values. This serious underestimation inferentially originates from the limitations of the observation operator for ice-phase hydrometeors in the three following aspects. Firstly, snow, graupel, and hail are all treated as ice ellipsoids with the same axis ratio during Rayleigh scattering calculations, while the true shapes are all different. Secondly, Mie scattering calculations may be more realistic for giant ice-phase hydrometeors like hailstones. Lastly, ice crystals, which are not considered in the observation operator, could be observable by weather radars when they are vertically aligned or canted by electric fields [79].
Hence, the NoDA simulation exposes three sources of the observation operator errors, including the raindrop shape-size relation, the limitations for ice-phase hydrometeors, and the melting ice model. These error sources all deserve further studies of optimization for typhoon cases.

Summary and Future Prospects
Typhoon Soudelor has been the most completely observed typhoon by the RCWF S-band polarimetric radar in Taiwan up to the present. No other typhoons have impinged on northern Taiwan as severely as Typhoon Soudelor did, and therefore, the polarimetric radar observations near the landfall of Typhoon Soudelor are precious and unique. This study uses the WDM6 microphysical scheme and incorporates the observation operator of Reference [45] into the WRF-LETKF radar assimilation system to investigate the sensitivities of 6 h QPFs for Typhoon Soudelor near landfall to the assimilation of the RCWF radar observations. The configuration of the WRF model refers to an operational run from Taiwan's Central Weather Bureau, and the configuration of LETKF is largely inherited from Reference [5]. After the quality control, thinning, and interpolation of the radar observations, a series of experiments assimilating various combinations of radar variables are carried out and compared with a control simulation without radar data assimilation for the purpose of improving a 6 h deterministic precipitation forecast for the most intense period since 2000 UTC 7 August. The results of the control simulation expose three sources of the observation operator errors, including the raindrop shape-size relation, the limitations for ice-phase hydrometeors, and the melting ice model. Nevertheless, this study inherits the default observation operator of Reference [45] and concludes that polarimetric radar data assimilation with the unadjusted observation operator can still improve the analyses of state variables, especially rainwater, and consequent QPFs for this typhoon case. The different impacts of assimilating Z H , Z DR , and K DP are only distinguishable at the lower levels of convective precipitation areas, with a ranking order of K DP , Z H , and Z DR from best to worst. Finally, the positive effect of radar data assimilation on QPFs can last three hours in this study, and further QPF improvement is achievable by the additional assimilation of polarimetric variables, especially K DP , on top of V r and Z H . For future prospects, further studies as stated above can be focused on suppressing the three sources of the observation operator errors exposed by the control simulation. The theoretical raindrop shape-size relation can be easily replaced with the in-situ one for typhoons in Reference [76] to re-simulate the scattering amplitudes via the T-matrix scattering calculations for rainwater. In contrast, the limitations for ice-phase hydrometeors are too difficult to deal with, and therefore ignoring the assimilation of polarimetric variables at higher levels can be a good strategy to adopt. As for the melting ice model, adjustment can be more meaningful than disuse because it is theoretically tenable to consider mixed-phase hydrometeors for the bright band signature. For example, a study joined by the authors recently found that, for an MCS case in Taiwan, the excessively simulated bright band signature can be mitigated when a minimum threshold of number concentrations for mixed-phase hydrometeors is applied to the melting ice model [80]. With a more accurate observation operator, the contribution from polarimetric variables can be more robust and also expectable for different types of precipitation systems. The results will be more representative of general performance if more cases are available and evaluated.