Direct Assimilation of Radar Reflectivity Data Using Ensemble Kalman Filter Based on A Two-Moment Microphysics Scheme for the Analysis and Forecast of Typhoon Lekima (2019)

: We investigate the impact of directly assimilating radar reflectivity data using an ensemble Kalman filter (EnKF) based on a double-moment (DM) microphysics parameterization (MP) scheme in GSI-EnKF data assimilation (DA) framework and WRF model for a landfall typhoon Lekima (2019). Observations from a single operational coastal Doppler are quality-controlled and assimilated. Compared with the baseline experiment initialized by GFS analysis, the reflectivity data assimilation experiment (Z-DA) resulted in an obvious improvement in both structural analysis and typhoon forecast skills in terms of intensity, precipitation, and track. Sensitivity experiments were conducted to evaluate the ability of EnKF to update certain state variables considering that the degree of freedom of analytical variables increased with a DM MP scheme. When either the total number concentration or other large-scale state variables that are not directly linked to reflectivity observations via the observation operator are not updated, the tendency of RMSIs and PS to be imbalanced is significantly increased during DA cycles compared to those of Z-DA with updating a full set of state variables, resulting in increased intensity and track forecast errors. These results indicate that the reliable ensemble covariance could handle the underconstraint issue associated with the DM scheme, and helps in obtaining more physically balanced analytical fields.


Introduction
Tropical cyclones (TC) can cause heavy losses of property and casualties through damaging winds and extreme precipitation.China has an average annual landfall of 7-9 TCs, rendering it one of the countries with the most severe TC disasters.Although the forecast accuracy of TC track has been improved in recent decades due to the continuous development of numerical weather prediction (NWP) models, the improvement in TC intensity and precipitation forecasting remains limited.
There are few observation platforms in the TC inner-core region.Radar is capable of observing TC inner-core structures and circulations at high spatial and temporal resolutions.Radar radial velocity and reflectivity could provide effective information on the wind structure and microphysical process at the TC vortex and convective scale [1][2][3].The assimilation of radar data would help in improving the initialization of the TCs.Many studies have attempted to define convective structures in TC initializations using a three- Academic Editor: Xiaolei Zou dimensional variational (3DVAR) method to assimilate radar radial velocity, leading to improved TC forecasts [4][5][6][7][8].However, using static background error covariance, 3DVAR is less able to describe the evolution of TCs and obtain correlations among cross-variables.Compared to 3DVAR, the ensemble Kalman filter (EnKF) has shown encouraging success in radar radial velocity with flow-dependent ensemble covariance.The EnKF method for assimilating radar velocity data is useful in TC vortex initialization, and contributes to improving TC forecast skill in intensity, structure, and rainfall prediction [3,[9][10][11].
Unlike radar radial velocity, a measure of wind field, reflectivity has knowledge about hydrometers that depends on scatter concentration, phase, size, and other properties.The assimilation of radar reflectivity is more complex and much more challenging than the assimilation of radial velocity [12,13].Cloud analysis is an early approach to assimilating reflectivity data into an NWP model, to adjust some of the hydrometer variables and in-cloud temperatures.Many studies showed the positive impact of using cloud analysis to initialize convective scale NWP and to mitigate the spin-up problem in shortrange forecasts [14][15][16][17].Cloud analysis is a promising tool for TC track and intensity forecasts [7].However, a number of empirical relationships between the hydrometer variables and the reflectivity are included in the cloud analysis method, introducing many additional uncertainties.
The VAR method to directly assimilate reflectivity has two common difficulties.One is associated with the tangent linear and adjoint of the strong nonlinear observation operator.Sun and Crook [18] found that, in their 4DVAR scheme, a large gradient of reflectivity cost function term associated with a small rainwater mixing ratio prevents effective convergence.To reduce this negative effect of nonlinearity, they retrieved a rainwater mixing ratio from reflectivity and then assimilated the rainwater mixing ratio into their 4DVAR scheme.Only warm-rain microphysics was used, and only rainwater mixing ratios were considered in their framework.Due to the neglect of many other hydrometers, especially ice-related hydrometers, this method would degrade the quality of analysis and forecast, especially in strong convective areas.Wang et al. [19] also encountered a similar problem, so they had to use reflectivity data of less than 55 dBZ.As other ways to mitigate the nonlinear problem, some kinds of the transformation of the hydrometer (rainwater) mixing ratio were introduced [20][21][22][23].Another difficulty of direct reflectivity assimilation in variational schemes is related to static background covariance, which is not only limited in describing correlations among cross-variables, but also likely to produce unphysical analyses of hydrometeors due to its static and homogeneous background variances [13].One potential solution is to employ a hybrid ensemble-variational (EnVAR) method [21] [24] to consider the model constraints in the form of ensemble covariance.Although the hybrid method is able to include flow-dependent background covariance, it still faces the nonlinearity challenge since it obtains the analysis in a variational framework.
In contrast with VAR or hybrid methods, the ensemble Kalman filter (EnKF) is a powerful approach for directly assimilating radar reflectivity [25][26][27].It is easy for EnKF to use a highly nonlinear observation operator (e.g., radar reflectivity) and models with complex ice microphysics because the adjoint model is no longer required.Moreover, it is straightforward for EnKF to utilize covariances between reflectivity and unobserved model variables.Several encouraging results for direct radar reflectivity assimilation using EnKF were reported [28].However, EnKF is vulnerable to model errors [29].The microphysics parameterization (MP) process is one of the most important sources of model errors.In recent years, studies indicated that prediction errors can be reduced by adopting the double-moment (DM) MP scheme [30][31][32].If a DM scheme is adopted in a model, it would be better to use the same MP scheme in the data assimilation stage.That is, the reflectivity observation operator is also based on the same DM MP scheme, and the number concentration of hydrometers is updated in the process of assimilation.Though using DM MP scheme in EnKF reflectivity DA is theoretically promising, it has a potential challenge since it increases the number of analytical variables to be estimated, which might further aggravate the underconstraint problem because only reflectivity observations are assimilated to update the number of analytical variables.The studies of EnKF to assimilate reflectivity based on the DM MP scheme are very limited [33][34][35].So far, reflectivity has not been assimilated by using EnKF based on DM MP scheme for TC simulation.Therefore, more studies are needed in order to better use EnKF to directly assimilate reflectivity to improve TC forecasts.
This study investigates for the first time the performance of EnKF reflectivity DA based on a partial DM MP scheme for a real TC case.We selected supertyphoon Lekima that occurred on August 2019 in the Northwest Pacific.The impact of the EnKF assimilation of radar reflectivity on analyses and TC forecasts is examined.Additional sensitivity experiments were conducted to evaluate the impact of updating the number concentration of hydrometers and other large-scale cross variables not directly related to the reflectivity observation operator.
The rest of this paper Is organized as follows: in Section 2, the TC Lekima case, the forecast model, radar observations, EnKF assimilation system, and experimental setup are described.Section 3 discusses the impact of reflectivity data on analyses, forecasts of intensity, track, and precipitation.Sensitivity to the impacts of the analytical variables is presented in Section 4. The conclusions are given in Section 5.

Case Overview
Supertyphoon Lekima (2019) was the most intense TC during the 2019 Northwest Pacific TC season.It made landfall in Wenling, Zhejiang province, reaching 16 force (52 m/s), at 1745 UTC, 9 August.After the first landfall in Zhejiang, Lekima passed over Jiangsu, Shandong and made its second landfall in Qingdao at 1250 UTC, 10 August, with maximal force 9 (23 m/s).Over its inland path, Lekima brought heavy precipitation and damage to 11 Chinese provinces, rendering it the second costliest TC in China since 2000, only after TC Fitow (2013).This study investigates the impact of directly assimilating radar reflectivity using EnKF.

Numerical Model
Predictions during EnKF cycles are performed using the Advanced Research Weather Research and Forecast ( WRF-ARW) model [36] version 3.9.1.Prognostic state variables include velocity components (u, v and w), perturbation potential temperature (θ), perturbation dry air mass in column (MU), and water vapor mixing ratio (qv), and microphysical state variables.In this study, the DM MP scheme of Thompson [37] is used.Microphysical variables were the mixing ratios of cloud water, rainwater, ice, snow, and graupel, and the total number concentrations of rain.The forecast domain had 548 × 452 grid points with a horizontal resolution of 3 km and 51 vertical levels from the surface to 20 hPa (Figure 1).Other physical parametrizations included the Yonsei University (YSU) planetary boundary layer scheme [38], the unified Noah land surface model [39], and the rapid radiative transfer model for global circulation model (RRTMG) shortwave and longwave schemes [40].

Assimilation Data and Algorithms
Data assimilated in this study were radar reflectivity factor (Z) from a single coastal WSR-88D radar at Wenzhou station in Zhejiang province.The radar data coverage is shown in Figure 1.For quality control, the original raw data were processed to remove ground clutter, speckles, and other artefacts, and converted into a mosaic of observations on a three-dimensional grid with a horizontal resolution of 3 km and 33 vertical levels.Only reflectivity data in precipitation regions (>=5 dBZ) were assimilated.
In this study, reflectivity observations were directly assimilated with a GSI EnKF system using 41 ensemble members.The EnKF scheme used here is commonly called an ensemble square root filter (EnSRF) [41].The length scale for horizontal covariance localization is 12 km, and the scale height (unit of −log (P/Pref) for vertical localization is 0.1.The relaxation to a prior spread (RTPS) multiplicative inflation scheme [42] was applied to account for the deficiency of the ensemble spread in the EnKF: where ′ is the posterior (analysis) perturbation for the i-th ensemble member; , are the prior (background) and posterior (analysis) standard deviation (spread) at each grid point; and is the inflation parameter, which is 0.95 in this study.For radar reflectivity assimilation, the gross check in the GSI EnKF system was applied.The observation error's standard deviation for reflectivity was specified as 5 dBZ.If the difference of observation and background reflectivity is greater than 7 times of the standard deviation of observation error, this observation fails the gross check and is not assimilated.
Unlike most previous research on direct reflectivity assimilation where reflectivity operators were based on a relatively simple formula that used Rayleigh scattering approximation and a single-moment (SM) MP scheme [13,21,43], this study adopted a more advanced and complex observation operator that was initially developed by Jung et al. [27], who used a fitted approximation to the T-matrix [44] for rain, and the Rayleigh approximation for ice species.This reflectivity operator takes into account nonspherical shapes for large precipitation species, and the non-Rayleigh scattering effect for raindrops.This operator was updated by Xue et al. [33] to add the total number of the concentrations of hydrometers as input variables, which is consistent with a DM MP scheme.
Using a fitted approximation to the T-matrix for the scattering amplitude calculation, the formula of rain reflectivity at the horizontal polarization is as follows: where is the radar wavelength, which is approximately 10.7 cm for the S-band radars.= 0.93 is the dielectric factor for water.Subscript r is the rain, = 4.28 × 10 , and = 3.04 are fitting values of the T-matrix calculation for rain.
is the intercept parameter for rain, which is consistent with Thompson DM MP.Slope parameter for rain is established from the intercept parameter and mixing ratio of rain.Γ is the complete gamma function of particle size distribution (PSD).
For other hydrometer species, the forms of radar reflectivity at the horizontal polarization are listed below: where where is the mean canting angle, and is the standard deviation of the canting angle.= 0 was assumed for all species, and = 20 °, 60 ° for snow and graupel, respectively.Subscript x can be dry snow (ds), dry graupel (dg), rain-snow mixture (rs), or raingraupel mixture (rg).
, are the fitting results for calculating scattering amplitudes from Rayleigh scattering approximation.

where
= / + / is a rain fraction calculated from the mixing ratio of rain, snow, or graupel.
For more details on the radar observation operators, the readers are referred to Jung et al. [27].

Experimental Design
Figure 2 illustrates the workflow of all the experiments.A control forecast (CNTL) serving as the baseline was initiated with NCEP GFS analysis at 0600 UTC, 09 August, and no radar data were assimilated.In all DA experiments, a 41-member ensemble was first generated at 0600 UTC.The initial fields for Member 41 were taken directly from GFS analysis, same with that in CNTL, and the initial fields for Members 1-40 were created by adding random perturbations generated from WRFDA 3DVAR with the cv3 background error covariance option [45] to the GFS analysis.Then, a 5 h ensemble forecast was launched for spin-up that would hopefully sample the mesoscale uncertainties of models.The lateral boundary conditions were generated from the GFS forecasts.To evaluate the impact of radar reflectivity data on the TC analyses and forecasts, in each DA experiment, radar reflectivity data were assimilated every 12 min for 5 cycles.The DA cycles started at UTC 1100 and ended at 1200, 9 August, about 5 h before Lekima's first landfall.The background fields for each cycle were provided by the 12 min ensemble forecast from the previous cycle.At the end of DA cycles, a 12 h deterministic forecast was launched from the analysis of Member 41 rather than from the ensemble mean in order to exclude the influence of ensemble forecasting to conduct a fair comparison with CNTL.To examine the impact of updating different variables, a set of sensitivity experiments were performed.Table 1 lists all experiments presented in this paper.Experiment Z-DA updated a full set of state variables, namely, u, v, w, MU, θ, qv, hydrometeor mixing ratios (qc, qi, qr, qs, qg), and the total number concentrations of rain water (Nr), consistent with the Thompson MP scheme in the WRF model.Experiment Z_NoNr excluded the update of Nr compared to experiment Z-DA.Increasing the analytical variables (i.e., Nr) in DM MP schemes may aggravate the underconstraint problem, since only reflectivity observations are assimilated, and the model error in a real case may lead to filter divergence with unreliable multivariate covariance.Experiment Z_NoNr was designed to examine the ability of EnKF to simultaneously estimate all possible state variables associated with the DM MP scheme.Two additional experiments, Z_NoUV and Z_NoCross (Table 1), were the same as experiment Z-DA but excluded the update of either horizontal wind components or all large-scale variables (u, v, w, θ, qv,), respectively.These two experiments were designed to investigate the impact of updated large-scale variables that were not directly linked to reflectivity observations via the observation operator.To verify the simulation results, the best track data from the China Meteorological Administration (CMA) [46,47] were used as the truth.The precipitation observations were the multisource merged hourly analysis product with a resolution of 0.1 °× 0.1 ° from CMA [48].

Impact on TC Analysis and Forecast
At the end of assimilation cycles (12,000 UTC), the forecasted and analyzed composite reflectivity of experiments CNTL and Z-DA are presented in Figure 3b,c.There was an obvious asymmetrical double eye-wall structure with the strongest reflectivity in the northwestern quadrant in the observed truth (Figure 3a).This feature was missing in the CNTL, but well-captured in experiment Z-DA.Furthermore, CNTL showed stronger and wider inner-core circulations, especially an overestimated precipitation area that occurred in the northern quadrant.In experiment Z-DA, a weaker and tighter inner core appeared compared to that in CNTL, the overestimated bias was corrected, and the rainband structures were much closer to the truth.The positive impact of radar assimilation on the analysis extended to the following forecasts (Figure 3d-l).For the 12 min forecast after DA (Figure 3d-f), CNTL still exhibited wilder and more asymmetric inner rain bands, and was missing a double eye-wall.In contrast, the forecast initiated from Z-DA analysis was closer to radar observations in bith the inner core and the region of outer spiral rain bands with weaker and more evenly distributed precipitation on the eastern coastal area in Zhejiang province.These improvements in precipitation lasted in the following forecast times (Figure 3g-i).
Figure 4 shows the horizontal wind speed and pressure of the analyses in the final cycle (1200 UTC, 9 August 2019) for experiments CNTL and Z-DA at a height of 1 km above the ground and in the vertical west-east cross-section through the individual vortex center of each experiment.Compared with CNTL, the wind speed was significantly greater, and the central pressure was lower in experiment Z-DA (Figure 4a,b).For experiment CNTL, the vortex was asymmetric with the strongest wind speed on the eastern side of the TC center (Figure 4c,d).With radar data assimilation, experiment Z-DA had a maximal wind of over 60 m/s, compared with 56 m/s in the CNTL.The height of the wind speed exceeding 44 m/s was over 3 km in the western side in experiment Z-DA, whereas it reached only 1 km in CNTL.Wind speeds exceeding 36 m/s extended to 6 km on both sides of the vortex in experiment Z-DA.These indicate that a deeper and more symmetrical vortex is obtained when reflectivity data are assimilated.

Intensity and Track Forecast
Figure 5 shows a comparison of the track and intensity between experiments CNTL and Z-DA for a 12 h forecast from 1200 UTC to 0000 UTC, 10 August, together with the observed best CMA track.Although the maximal wind speeds were underestimated, a clear improvement was found with reflectivity data assimilation compared to CNTL.Starting from slightly stronger maximal surface winds, experiment Z-DA had a persistent smaller intensity error throughout the forecast period (Figure 5a), indicating that radar reflectivity significantly impacted both the analysis and forecast.Regarding the track, CNTL showed an obvious westward bias compared to the best track.With radar data assimilation, experiment Z-DA reduced the westward displacement error and yielded better track results.Reflectivity data covered only the vortex circulation in the examined case, and the environmental field (a dominant factor for TC track) was actually not updated due to the relatively short localization length.With more radar data in the environmental field, the track forecast might be further improved.

Precipitation Forecast
Figure 6 presents the 3 h accumulated precipitation from 1800 UTC to 2100 UTC, 9 August, for both CNTL and Z-DA along with CMA precipitation observations.As verified against the observations, CNTL had a broader maximal rain band with a westward displacement around the inner core, and underestimated the precipitation north of the vortex, especially in the area between the city of Shanghai and Jiangsu province.The assimilation of reflectivity in Z-DA helped in decreasing both the overestimation in the inner core and underestimation in the north, adjusting the rainfall pattern closer to observations.

Sensitivity Experiments
The ability of EnKF to update certain state variables needs to be further examined, considering that the degree of freedom of analytical variables increased with a DM MP scheme.The analysis and forecast root-mean-square innovations (RMSIs) of sensitivity tests without updating certain state variables are displayed in Figure 7, along with the Z-DA and CNTL experiments as a benchmark or reference.All the experiments with radar data assimilation were obviously better than CNTL.Among them, the RMSIs of Z_No-Cross without updating all large-scale cross-variables, including horizontal and vertical wind components, potential temperature, and water vapor, were obviously the largest for all five assimilation cycles.If only excluding the horizontal wind, the result of Z_noUV was between those of Z-DA and Z_NoCross, indicating that our assimilation could better update all the large-scale variables through ensemble background error covariance.Furthermore, the RMSIs of Z_NoNr without updating the number concentrations of rain were also clearly larger than those of Z-DA at all analytical cycles, suggesting that an imbalance between mixing ratios and the corresponding total number concentrations may exist when only mixing ratios are updated.Considering that Z-DA had the lowest RMSIs, it is suggested to update both microphysical and large-scale variables in the current EnKF system.The impacts of not updating certain state variables on the imbalance introduced by analysis were further investigated.In Figure 8, the domain-averaged absolute surface pressure tendency, an important criterion for evaluating imbalance, was calculated during assimilation cycles and the first hour of forecast [49].The larger the pressure tendencies were, the less the balance was.Among all the sensitivity experiments without updating certain variables, pressure tendencies had higher values and more dramatic oscillation than those in Z-DA.These indicate that updated state variables missing in sensitivity tests could introduce more noise and provoke greater imbalance into analyses.Experiment Z-DA showed the smallest pressure tendencies, suggesting that updating all state variables leads to much more physically consistent analyses in the current reflectivity assimilation.We also examined the track and intensity forecast for all sensitivity experiments (Figure 9).Overall, experiment Z-DA with all variables updated had the relatively smallest errors for both intensity and track forecast, probably because the Z-DA experiment was not only less unbalanced, but was also able to properly update all state variables.When only hydrometer variables were updated, most of the benefits of assimilating reflectivity data were lost, and the intensity forecast was even worse than that in CNTL without assimilating radar reflectivity.

Conclusions and Discussion
Using the WRF model and GSI-EnKF system with radar DA capabilities, the impacts of direct radar reflectivity assimilation associated with a DM MP on analysis and forecast for typhoon Lekima (2019) were investigated.The same Thompson MP scheme was employed in both the WRF simulation and the observational operator for EnKF DA.As the first attempt to directly assimilate radar reflectivity by EnKF using a DM MP scheme for a real TC case, the reflectivity observations from one coastal WSR-88Ds radar were assimilated every 12 min over a 1 h DA window.Experiments with and without assimilating reflectivity data were conducted.Results show that the direct assimilation of reflectivity significantly reduced the RMSEIs of both analyses and forecasts.Compared to those in experiment CNTL without assimilation, analyses in the Z-DA experiment at the final cycle were much more accurate, with the circulations and structures much closer to the observations with double eye-wall structures and stronger wind speeds.With the better analyses, experiment Z-DA was able to improve Lekima's forecast skills in terms of track, intensity, and precipitation.
Additional sensitivity experiments were conducted to investigate the underconstraint issue due to the increased degree of freedom associated with the DM MP scheme.The results show that experiment Z-DA with all analytical variables updated produced the best analytical fit to the reflectivity observations and the smallest forecast error growth rate.The three sensitivity experiments (Z_NoCross, Z_NoUV, and Z_NoNr) without updating certain analytical variables produced worse results than those of experiment Z-DA for both analyses and forecasts, while differences among the three sensitivity experiments were relatively significant.The worst results in Z_NoCross were likely due to the lack of dynamic consistency in the analyzed fields.To understand the impact of updating different state variables, the mean absolute values of surface pressure tendency were examined.Experiment Z-DA, updating both hydrometers and other large-scale variables not directly involved in observation operator of reflectivity, had the smallest imbalance, resulting in the largest improvement in the forecast skill of intensity and track.This result indicates that our EnKF framework could handle the increased degree of freedom issue associated with the DM MP scheme, which is inconsistent with that in Dong and Xue [28], who used reflectivity to update only pressure and hydrometers due to the negative impact found when updating other large-scale variables.One possible reason is that more reliable ensemble covariance was obtained in the current simulation; further investigation and research are needed in the future.
As the first attempt to directly assimilate radar reflectivity based on a DM MP for TC application, this study is meaningful for exploring the potential of radar data in TC forecast.However, only one case of typhoon was studied, so more experiments with more TC cases are essential to obtain robust conclusions.Further improvement may also be possible by investigating many other DA factors, such as the observation errors associated with radar operator and adaptive vertical localization.

Figure 1 .
Figure 1.Domain of the WRF model at 3 km horizontal resolution with radar reflectivity observations (shaded; units: dBZ) at 1200 UTC, 9 August; the rectangle indicates the location of the Wenzhou radar station.

Figure 2 .
Figure 2. Workflow of the CNTL and the reflectivity DA experiments.

Figure 5 .
Figure 5. Predicted (a) maximal wind speed and (b) track for 12 h forecasts from 1200 UTC, 9 August, to 0000 UTC, 10 August 2019, for CNTL (blue lines), and Z-DA (red lines).The best CMA track is shown in black lines.

Figure 7 .
Figure 7. RMSIs of simulated reflectivity in experiments Z-DA (red line), Z_NoUV (purple line), Z_NoNr (orange line), and Z_NoCross (green line).The RMSI for experiment CNTL is shown in black line as a reference.

Figure 8 .
Figure 8. Mean absolute pressure tendency during the five DA cycles and the first 1 h forecast for the four DA sensitivity experiments.

Table 1 .
List of experiments.