Assimilation of X‐Band Phased‐Array Radar Data With EnKF for the Analysis and Warning Forecast of a Tornadic Storm

Abstract

The impact of assimilating China's operational X‐band Phased‐Array radar's (X‐PAR) data on the analysis and warning forecast of the vortex structure and intensity of the June 8, 2018 Foshan, Guangdong province, tornadic storm was investigated for the first time using an Ensemble Kalman Filter (EnKF) data assimilation system. Both radar radial velocity (Vr) and reflectivity (Z) from two S‐band operational radars and one X‐PAR were assimilated. Deterministic forecasts were launched every 6 min from 05:42 UTC (20 min before the tornado touched down) to 06:00 UTC from the EnKF mean analysis field. Five experiments were conducted to examine the added capability of Z assimilation of the EnKF system, and to investigate the impact of assimilating X‐PAR data on the analysis and prediction of the tornadic storm. Compared to the experiment without Z assimilation, the assimilation of Z reduced the analysis error and greatly reduced the forecast error of Z. The assimilation of X‐PAR data greatly improved the vortex structure of the tornadic storm at low levels, and improved the intensity of the rear inflow of the tornadic storm, especially with a higher assimilation frequency. Compared to the experiments without X‐PAR data assimilation, assimilating X‐PAR data improved the predictability of tornadic storm.

Introduction

The average annual number of tornadoes in China is not more than 100 from 1961 to 2010 (Fan & Yu, 2015). However, recent significant tornado events such as the EF4 Funing tornado, the EF4 Caihong tornado, and the EF3 Kaiyuan tornado have drawn great attention of social media in China (Xue et al., 2016). Despite of the advances in tornado analysis and large eddy simulation over the past decades (Mashiko, 2016; Mashiko et al., 2009; Sun et al., 2019), the prediction of tornadoes remains difficult. One of the reasons is the lack of instruments that can capture the internal structures of tornadic storms, especially the low-level atmosphere conditions (McLaughlin et al., 2009).

Doppler weather radar is the only available observational tool that can monitor the internal structure of tornadic storms. However, it takes the operational Weather Surveillance Radar-1988 Doppler (WSR-88D) radars about 5–6 min to perform a volume scan, which is too slow to effectively capture the rapid evolution of tornadic storms. In contrast, Phased-Array radars (PARs) need only about 2 min to scan a volume. PARs can detect short-lived storms more quickly via the “electronic scanning” approach with agile beam steering and multibeam scanning, making PARs more applicable to practical nowcasting and early warnings (Nai et al., 2017; Wurman et al., 2012). Some observational experiments have demonstrated the advantages of PARs compared with traditional weather radars when analyzing the structure and evolution of storm-scale weather disasters. Liu et al. (2018) observed the inflow and outflow associated with the formation and dissipation of the hook echo in a supercell event from a PAR in South China, which cannot be captured in the analysis based on traditional weather radars. Kuster et al. (2015) summarized the advantages of PARs in the depiction of storm-scale process, including: (a) variations of the low- and mid-level vortex intensity and associated changes in the extent of storm inflow; (b) fast strengthening of the low-level vortex; and (c) the abrupt motion changes of tornado circulations.

Despite the progress of radar observations, the study of PAR data assimilation is still limited. Supinie et al. (2016) showed that assimilating Mobile Weather Radar, 2005 X-band, Phased-Array data increases the mean low-level vorticity in a convective storm analysis. Supinie et al. (2017) performed a four-dimensional asynchronous implementation of the ensemble square root filter assimilation and conducted forecast experiments to examine the impact of high-frequency PAR observations on a tornadic supercell storm case compared with that from the traditional WSR-88D. They found that the rapid scanning capability of PARs can potentially extend severe weather warning lead time.

To produce dynamically consistent storms from PAR observations, a suitable assimilation method that can take advantage of the high temporal resolutions of PAR data is necessary. The Ensemble Kalman Filter (EnKF) technique, which initially developed by Evensen (1994), estimates the flow-dependent background error covariance through a group of forecasts. This technique has been widely examined with many weather systems and used in real-case situations. Some studies demonstrated the advantage of EnKF by comparing with other data assimilation techniques for convection-allowing forecasts. For example, Johnson et al. (2015) performed a case study of an upscale growing mesoscale convective system and showed that forecasts initialized using Gridpoint Statistical Interpolation (GSI)-based EnKF are more skillful than those using GSI-based three-dimensional variational data assimilation (3DVAR) with radar data assimilation because of the flow dependence and cross-variable correlation in the EnKF background error covariance. In addition, Caya et al. (2005) showed the EnKF typically produces more accurate analyses than four-dimensional variational data assimilation after several assimilation cycles through a simulated radar data assimilation experiments in a cloud-resolving model. Several studies also used the EnKF technique to capture storm development and understand its associated physical mechanism(s). Yokota et al. (2016) found that the low-level convergence forward of the storm and humidity in the rear of the storm had an important effect on the strength of a low-level mesocyclone when assimilating surface and radar data by local ensemble transform Kalman filter (LETKF). Markowski et al. (2012) showed tornadoes could be maintained only if partially embedded within or located near the convergence along rear-flank gust fronts by examining EnKF analyses. Zeng, Janjic, et al. (2021) compared the relative importance of reflectivity (Z) and radial velocity (Vr) data in convective-scale data assimilation with LETKF and showed that assimilating Z requires a shorter spin-up time to reconstruct convective systems but assimilating radial wind is more important to reconstruct the mesocyclone structure of supercells.

Guangzhou (GZ) is the core hub city of the Pearl River Delta region with a population of about 15 million. As one of the major precipitation centers in China, Guangzhou frequently suffers from severe storms which induce tornado, hail, local short-time heavy rainfall and so on. However, operational S-band weather radars in GZ, even though with a wide detection range, are unable to capture such small-scale weather systems and their evolution. Therefore, in order to effectively capture severe weather systems in GZ, China Meteorological Bureau, Guangzhou Meteorological Bureau and Zhuhai Naruida Technology Co., Ltd. have been conducting X-band dual-polarization PAR observation experiments since 2017 based on the statistical analysis of echo structure of severe convective storms in GZ during 2015–2017. At present, four X-band Phased-Array radars (X-PARs) have been deployed within or near GZ, covering the central urban area of GZ, and the number of X-PARs will increase to 42 in the future. Among these four radar sites, three of them (Huadu, Panyu, and Baiyun) are located in Guangzhou and the other one (Nanhai) is located in Foshan, which is close to GZ. Examining the influence of assimilating X-PAR data on the analysis and prediction of severe storms is one of the main purposes of this study.

Here, we adopted the Pennsylvania State University EnKF data assimilation system developed by Meng and Zhang (2008a, 2008b) to assimilate the X-PAR data. This study explores for the first time the impact of X-PARs data assimilation on tornadic storm analysis and prediction in China using a single case study. The advantage of X-PARs over operational S-band radars will be examined. In addition, the impact of assimilation frequencies of X-PAR data on the analysis of wind field structure and the prediction of a tornadic storm will be investigated.

The rest of the study is organized as follows. Section 2 provides the case overview and the experimental design. The analysis results are presented and discussed in Section 3 while the prediction results are shown in Section 4. A summary and conclusions are presented in Section 5.

Case Overview and Experimental Design Overview of the June 8 Foshan Tornado

On June 8, 2018, an EF1 tornado hit Foshan, Guangdong, which was embedded within a rainband of typhoon Ewiniar and was about 66.6 km away from the typhoon center (see Figure 1b for the geographical location). The tornado lasted for 30 s, leaving severe damage along its moving path of about 380 m (Figure 2). It was found that the high-level divergence and suction, a surface mesoscale convergence line, and uplift triggering of the small-scale vortex provided favorable conditions for severe weather development during the event (not shown). Local convection parameters observed by the Hong Kong (HK) sounding station indicated the tornadic storm happened in an environment containing moderate convective available potential energy (CAPE, 1694.4 J kg⁻¹), weak convective inhibition (CIN, 2.1 J kg⁻¹), low lifting condensation level (116 m), large 0–1 km storm-relative helicity (SRH, 237 m² s²) and strong 0–1 km vertical wind shear (13.1 × 10⁻³ s⁻¹). Based on the tropical cyclone (TC) tornado proximity sounding climatology summarized by Schneider and Sharp (2007) (CAPE > 500 J kg⁻¹, SRH > 100 m² s²), the environmental parameters for the Foshan tornado mainly fell into the high-tornado-threat category.

Figure 1. (a) The Weather Research and Forecasting simulation domains. The small solid black box (pointed by the arrow) indicates the inner domain D2, where radar observations were assimilated. The blue dashed line represents the best track of typhoon Ewiniar marked at 3 hr intervals from 06:00 UTC June 7, 2018 to 12:00 UTC June 8, 2018 (from south to north). The eye of typhoon Ewiniar at 06:00 UTC June 8, 2018 is marked by the red point. (b) Enlarged view of the dashed rectangle area in (a). The black dot at the center denotes the location of the tornado according to local tornado report. Location of the NH X-band Phased-Array radar (X-PAR) is marked by the red triangle. The red solid circle indicates the NH X-PAR's observation range with a diameter of 41 km. Locations of the S-band GZ radar and S-band ZQ radar are marked by the two black triangles. The black solid circles indicate the S-band radar's observation range with diameters of 460 km. The HK sounding station is marked by the dark cross.

Figure 2. The Foshan tornado track (blue dashed line) on satellite images with some damage photos along the tornado track. The A-E are the corresponding locations where we took these photos.

The tornadic storm was well captured by the Nanhai (NH) X-PAR. Figures 3 and 4 show the evolution of the convective storm which produced the tornado. From 05:42 UTC to 05:48 UTC, the storm intensified quickly, with maximum Z reaching 60 dBZ. At 05:42 UTC, the Vr map clearly shows the positive and negative velocity couplet associated with a mesovortex (Figure 4a). At 05:48 UTC, the height of the cyclonic vortex where the velocity couplet disappeared increased to 1.7 km (see the black line in Figure 4f). Corresponding to the velocity couplet, an echo overhang appeared (Figure 3f), which is an updraft indicator (Tanamachi et al., 2012). During 05:48 UTC–05:54 UTC, the storm continued to move northward and was about to enter Foshan (Figure 4c). The echo overhang was still obvious (Figure 3g), and the height of the vortex remained at 1.7 km AGL (Figure 4g), indicating the strength of the convective storm was maintained. From 05:54 UTC to 06:00 UTC, the tornadic storm moved into Foshan. The echo overhang almost disappeared at this time (Figure 3h), indicating the updraft might become weak at 06:00 UTC. The height of the cyclone decreased to about 1 km AGL (Figure 4h). The height of the maximum intensity of the vortex also decreased (see the height of the symbols “$\odot$” and “$\otimes$” in Figure 4h), perhaps creating a good chance for the vortex to touch the ground. Here, the position of the maximum intensity of the vortex is determined by the position of the maximum of azimuthal shear (Zhao et al., 2017). Two minutes later, Foshan reported the occurrence of the tornado on the ground.

Figure 3. The composite radar Z of GZ radar, ZQ radar, and NH X-band Phased-Array radar at 1 km AGL with 6 min interval (first row) and the corresponding vertical slice along the black line AB marked in panels a–d (second row), at 05:42 UTC (a, e), 05:48 UTC (b, f), 05:54 UTC (c, g), and 06:00 UTC (d, h). The rotation symbol in the first row represents the location of the tornadic storm. The marker “+” in (d) is the location of the reported tornado.

Figure 4. NH X-band Phased-Array radar (X-PAR) observed Vr. The location of the NH X-PAR is marked by the black triangle. The symbols “⊙” and “⊗” in the second row represent the direction away from and toward to the radar station, respectively, and the height of the symbols “⊙” and “⊗” represents the height of the maximum intensity of the vortex. The horizontal black lines in the second row represent the height of the cyclonic vortex.

Experimental Design

In this study, radar data from two operational S-band Doppler radars located at GZ and Zhaoqing (ZQ), and the NH X-PAR were used. The locations and observational ranges of the three radars are shown in Figure 1b. For the two S-band Doppler radars, the radar specifications are the same as that shown in Table 1 of Huang et al. (2018). In particular, their volume scans consist of nine elevations between 0.5° and 19.5° with a maximum Doppler range of 460 km. It takes 6 min for these two S-band Doppler radars to scan a volume and their spatial resolution is about 250 m. The NH X-PAR observes at 17 elevations between 0.9° and 29.7° with a maximum Doppler range of 41.6 km. It takes 2 min for the NH X-PAR to scan a volume and its spatial resolution is about 30 m. At the storm location, the vertical observation height ranges of GZ radar, ZQ radar, NH X-PAR are 1.23–8.52 km AGL, 0.67–20.96 km AGL, and 0.22–5.78 km AGL, respectively. For the X-PAR data, the attenuation in Z was corrected by parameterizing specific attenuation (A_H) and differential attenuation (A_DP) from specific differential phase (K_DP) linearly (Huang et al., 2016), which is given by [Image Omitted. See PDF] [Image Omitted. See PDF]where $\alpha$ and $\beta$ are 0.32 dB km⁻¹ and 0.059 dB km⁻¹, respectively. Both parameters were obtained from the scattering simulation using the raindrop size distribution data from Huang et al. (2018). The Doppler velocity dealiasing was achieved using the region-based algorithm given by the Atmospheric Radiation Measurement Climate Research facility Radar Toolkit (Helmus & Collis, 2016). For the other details of the quality control procedures of radar data, one can refer to Section 4 of Huang et al. (2016, 2018).

The numerical weather prediction model used in this work is the Weather Research and Forecasting (WRF) Model version 3.9.1. Two one-way nested domains were used. The outer domain D1 had 600° × 600° horizontal grid points with 2.5 km grid spacing, which covered the whole wind structure of typhoon Ewiniar. A much smaller inner domain D2 was used as the analysis domain for the tornadic storm for the purpose of saving computational cost, which had 71° × 71° grid points and 0.5 km grid spacing in the horizontal (Figure 1a). Both domains had 51 vertical levels with model top at 10 hPa. The boundary and initial conditions were provided by the 6-hourly and 1° × 1° final analyses of the Global Forecast System from the National Centers for Environmental Prediction. The physics options used include the Thompson microphysics scheme, the Rapid Radiative Transfer Model longwave radiation scheme, Dudhia shortwave radiation scheme, revised Monin-Obukhov surface layer scheme, unified Noah land-surface model, and Yonsei State University boundary-layer scheme. No cumulus scheme was used for D1 and D2.

The WRF-based EnKF system used in this work is based on that in L. Zhu et al. (2016) and Zhang et al. (2009), which was originally developed by Meng and Zhang (2008a, 2008b). Here, we added the capability to assimilate Z to this system. To simplify the design and also keep consistency with the WRF model, we use Z calculated by WRF. The innovation of Z is then calculated within the EnKF by adding Z to the list of analysis variables within the EnKF framework. The observation error for Vr and Z were set to 1 m s⁻¹ and 1 dBZ, respectively.

The observation operator for Vr is [Image Omitted. See PDF]where $u$, $v$, and $w$ are model-predicted zonal, meridional, and vertical wind speeds (m s⁻¹) interpolated to the observation locations, respectively; while $\phi$ and $\theta$ are the azimuth angle and elevation angle of the radar beam.

The observation operator for Z is from the Advanced Research version of WRF post algorithm based on Stoelinga (2005) and Thompson et al. (2008). The following are the formulas for radar reflectivity factors: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where the indices of a, r, s, g represent air, rain, snow and graupel, respectively, and $\rho$ is density, $q$ is mixing ratio, and $N_0$ is the intercept parameter which can be calculated as: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where $t_c$ is the temperature in Celsius. The equivalent reflectivity value associated with each hydrometeor mixing ratio is then calculated at model grid points, and a total equivalent reflectivity can be obtained by summing these values together: [Image Omitted. See PDF]

In this study, 45 initial ensemble members were created by adding perturbations generated by WRF 3DVAR random coefficient of variation option 3 (refer to Barker et al., 2004 for details) to the operational GFS analyses. The perturbed variables included potential temperature, horizontal wind components and mixing ratio for water vapor, with standard deviations of 1 K for temperature, 2 m s⁻¹ for wind, and 0.5 g kg⁻¹ for water vapor mixing ratio, which were similar to those applied by Meng and Zhang (2007). For inflating the background error, the relaxation-to-prior perturbation method (Zhang et al., 2004) was used with a relaxation coefficient of 0.5.

The following WRF variables were updated when assimilating Z and Vr: perturbation potential temperature (T), surface skin temperature (TSK), temperature at 2 m (T2), perturbation geopotential (PH), perturbation dry air mass in column (MU), surface pressure (PSFC), and perturbation pressure (P). Additionally, water vapor mixing ratio (QVAPOR), cloud water mixing ratio (QCLOUD), rain water mixing ratio (QRAIN), cloud ice mixing ratio (QICE), snow mixing ratio (QSNOW), graupel mixing ratio (QGRAUP), and rain number concentration (NRAIN) were updated only when assimilating Z, while vertical velocity (W), horizontal wind components (U and V), and 10-m horizontal wind (U10 and V10) were only updated when assimilating Vr. The purpose of these different choices of updated model variables was to reduce false correlation (Dowell et al., 2011).

The successive covariance localization method (Zhang et al., 2009) was used. The horizontal localization radius of influence (ROI) was set to 2 km for both Vr and Z. The vertical ROI was set to three vertical model grids for Vr and two vertical model grids for Z, using the fifth-order correlation function proposed by Gaspari and Cohn (1999). This method has been also used in several previous papers (Zhang & Weng, 2012, 2015; Zhang et al., 2009; K. Zhu et al., 2013). Due to the limitation of computational cost, the data thinning coefficient was set to 10, which means $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$10$}\right.$ of the total observations are randomly chosen and assimilated, corresponding to approximately 300-m and 2.5-km spatial resolution for X-band and S-band radar data, respectively.

Figure 5 shows the flowchart of the experiments. Ensemble forecasts were first run for 5 hr for spin-up. Starting from 05:00 UTC 8 June, radar data were then assimilated in D2 until 06:00 UTC 8 June with various assimilation frequencies (explained later), while no data assimilation was performed in D1.

Figure 5. (a) Data assimilation and forecast timelines of the experiments where deterministic forecasts were initialized from each experiment every 6 min from 05:42 UTC to 06:00 UTC and (b) configurations of the five data assimilation experiments where the upward pointing arrows indicate the times when radar data were assimilated. Note that no data were assimilated at 05:50 UTC because observation was not available at this time.

To examine the impact of assimilating operational S-band radar and X-PAR data, five experiments were carried out. The first experiment, Vr_S6, only assimilated Vr data from the two operational S-band radars (GZ and ZQ radar) every 6 min from 05:00 UTC to 06:00 UTC. The second experiment, VrZh_S6, was the same as Vr_S6 except that Z data from the two S-band radars were also assimilated. The purpose of this experiment is to examine the impact of the newly added observation Z. The third and fourth experiment (VrZh_S6_XZ6 and VrZh_S6X6) were conducted to investigate the impact of X-PAR on the analysis of the tornadic storm and also to investigate the different impact of NH X-PAR Vr and Z. Experiment VrZh_S6_XZ6 assimilated Z data from NH X-PAR at a 6-min interval starting from 05:00 UTC to 06:00 UTC in addition to those assimilated in the experiment VrZh_S6, while experiment VrZh_S6X6 assimilated both Vr and Z from NH X-PAR. In the fifth experiment, VrZh_S6X2, the assimilation frequency of X-PAR Vr and Z data was increased to every 2 min after 05:42 UTC (when observation was available) in order to show the impact of increased assimilation frequency. Deterministic forecasts from ensemble mean analysis of experiments were launched every 6 min since 05:42 UTC (20 min before the tornado touched down). All the forecasts ended at 06:02 UTC when the tornado touched down to the ground. The last forecast starting time of each experiment was 06:00 UTC, in which only 2-min of forecasts was generated.

The Polarimetric Radar Data (PRD) emulator software (Jung et al., 2008) developed by the Center for Analysis and Prediction of Storms (CAPS) of Oklahoma University in the US was used to emulate Z and Vr from WRF-model outputs, which is the same as the EnKF observation operator of Z and Vr mentioned above. The analyses and forecasts were then evaluated against GZ radar and NH X-PAR observations.

Analysis Results The Impact of Z From S-Band Radars on the Analysis of Storm

First, the impacts of the newly added Z assimilation on the analysis results of the tornadic storm were examined. Figures 6 and 7 show the analyses of storm-relative wind fields, vertical vorticity, and Z at 1 and 0.2 km height at the last assimilation time, while Figure 8 shows the cross sections of them along the lines in Figures 6a and 7b. The vertical vorticity was calculated using the model wind at 500 m resolution. Here, we focus on Figures 6a, 6b, 7a and 7b, while the rest will be discussed in the next subsection. For the dynamic structure, due to their limited resolution, Vr_S6 and VrZh_S6 could only resolve the large but weak vortex. There were clear vortex centers at the 1 km height in both Vr_S6 and VrZh_S6, but both of them failed to reproduce the vortex at the 200 m height. From the vertical cross section, it could be seen that the analysis vortexes were mainly existed above 200 m (see Figures 8a and 8b). This is not surprising since the S-band radar cannot observe the low-level wind structure.

Figure 6. The analyzed Z and storm-relative wind vectors at the height of 1 km AGL at 06:00 UTC, June 8, 2018 from experiments (a) Vr_S6, (b) VrZh_S6, (c) VrZh_S6_XZ6, (d) VrZh_S6X6 and (e) VrZh_S6X2. The positive (white solid contours) and negative (white dashed contours) vertical vorticity is shown at ±5, ±10, and ±15 × 10−3 s−1. The marker “+” is the location of the reported tornado. The black line CD in panels (a) (not shown in other panels) marks the cross section shown in Figure 8. The dashed box represents the area in the first line of Figures 3 and 4.

Figure 7. Same as Figure 6, but for the height of 200 m AGL. Remarkable vortexes (positive vertical vorticity >10 × 10−3 s−1) at 200 m height are emphasized by the black ellipse.

Figure 8. Cross sections of vertical wind (black contours, 2 m s−1 intervals, negative values dashed), vertical vorticity (white contours at ±5, ±10, and ±15 × 10−3 s−1), Z (color shaded), and wind (vectors, m s−1) at the analyze time 06:00 UTC from experiments (a) Vr_S6, (b) VrZh_S6, (c) VrZh_S6_XZ6, (d) VrZh_S6X6, and (e) VrZh_S6X2 along the lines in Figures 6a and 7a. The thick black dotted lines represent the corresponding heights of GZ radar 2.4° tilt (higher one) and NH X-band Phased-Array radar 0.9° tilt (lower one).

Figure 9 is the “sawtooth” plot of analysis and forecast errors through the assimilation cycles against GZ radar (2.4° tilt) and NH X-PAR (0.9° tilt) radar observations. Those two tilts were chosen as observations because the vortex of the tornado in the upper and lower levels could be well captured by them. The RMSEs were calculated at grid points where observed Z exceeded 10 dBZ (K. Zhu et al., 2020). In general, the RMSE differences of Vr between Vr_S6 (blue lines) and VrZh_S6 (black lines) are small with no more than 0.01 m s⁻¹ (Figures 9a and 9c), which is consistent with the vorticity structure analysis described above. This result indicates that the assimilation of Z from S-band radar has a limited effect on the upper-level dynamic field. VrZh_S6 got slightly smaller RMSEs than Vr_S6 at lower level for the last few analysis cycles (Figure 9c), which could be benefitted from the more balanced microphysics and dynamical structure of convection after the assimilation of Z. The assimilation of Z helps to reconstruct convective cells and suppresses spurious ones (Zeng, Janjic, et al., 2021) especially at the low levels (see Figures 8a and 8b), which in turn accelerates the adjustment of the wind structure and helps to reduce the spin-up time of model convection. That explains why experiment VrZh_S6 got consistently smaller RMSEs for both Vr and Z than that of Vr_S6 for the later analysis cycles.

Figure 9. (a) The RMSEs of model Vr against 2.4° tilt of GZ radar observation before and after each analysis from Vr_S6 (blue), VrZh_S6 (black), VrZh_S6_XZ6 (purple), VrZh_S6X6 (green), and VrZh_S6X2 (red). (b) Similar as (a) but for RMSEs of Z. (c) and (d) are the same as (a) and (b) but for RMSEs calculated at the 0.9° tilt of NH X-band Phased-Array radar observation.

In terms of Z structure, without the assimilation of S-band radar Z observation, Vr_S6 underestimated the Z intensity of the tornadic storm while overestimates the intensity outside the rainband (Figure 3d). The assimilation of Z in VrZh_S6 greatly increased the intensity in both the horizontal (Figures 6a, 6b, 7a and 7b) and vertical (Figures 8a and 8b) cross sections of the tornadic storm, while reduced the Z intensity in areas outside the rainband because of the assimilation of no reflectivity data (≤5 dBZ) (Zeng, Janjic, et al., 2021) in this region. The suppression of convection outside of typhoon rainband clearly improved the tornadic storm structure (Figures 6a and 6d). In terms of RMSEs, the average RMSE of Z in VrZh_S6 (black line) was 1.0 dBZ lower than that in Vr_S6 (blue line) at the tilt of 2.4° (Figure 9b). The RMSEs of Z in VrZh_S6 was also 1.0 dBZ lower than that of Vr_S6 (Figure 9d) at the tilt of 0.9° for the last four analysis cycles. To summarize, the additional S-band radar Z assimilation could help to reconstruct the convective system and reach its balance state faster than with only Vr assimilation. In terms of RMSEs of Vr and Z, experiment VrZh_S6 gets smaller value than that of Vr_S6 for the later analysis cycles, indicating VrZh_S6 fits better with the observations.

The Impact of X-PAR Data Assimilation

For the dynamic structure, the assimilation of X-PAR data greatly improved the low-level wind structure (Figures 6d and 7d). This is because X-PAR data well compensates the S-band radar observations at the low level. It can be seen that the vertical vorticities were enhanced in experiment VrZh_S6X6 (Figures 6d and 7d). At 1 km height, the maximum intensity of vortex reached more than 18 × 10⁻³ s⁻¹ in VrZh_S6X6, while no more than 15 × 10⁻³ s⁻¹ in VrZh_S6. Compared to VrZh_S6, the vortex center in VrZh_S6X6 was closer to the location of the observed tornado (see the marker “+” in Figure 6). The most significant difference can be found at 200 m height. The VrZh_S6X6 successfully reproduced a vortex, while no vortex at low levels was produced in VrZh_S6, which is consistent with the vertical slice in Figures 8b and 8d. Without assimilating X-PAR data, the vortex in VrZh_S6 did not reach the ground, while the vortex center extended to the ground with the X-PAR data assimilation. The significant difference was also found in the rear inflow region (see the lower left of Figures 8b and 8d). With X-PAR data assimilation, the rear inflow was greatly enhanced, especially at the low levels.

To see different impact of X-PAR Vr and Z on the analyzed dynamic structure of the storm, experiment VrZh_S6_XZ6 with only X-PAR Z assimilation was also plotted in Figures 6–8. The assimilation of X-PAR Z increased the intensity of the vortex. The maximum vertical vorticity increased from 14 × 10⁻³ s⁻¹ in experiment VrZh_S6 to 17 × 10⁻³ s⁻¹ in experiment VrZh_S6_XZ6 at 800 m AGL. However, the overall differences are very small and limited. Without NH X-PAR Vr assimilation, the dynamic structure of the tornadic storm at the near surface levels was hard to reproduce. It could be seen that the vortex from experiment VrZh_S6_XZ6 failed to extend to the ground (Figure 8c). The further assimilation of X-PAR Vr, as discussed above, greatly increased the intensity of vertical vorticity for all levels (Figure 8d), confirming the improvement of the storm's dynamic structure was mainly coming from the X-PAR Vr assimilation. This is consistent with the conclusion of Zeng, Janjic, et al. (2021), which also shows that assimilating Vr data is more important to reconstruct the dynamic structure of the storm.

The Z structure was also improved by the assimilation of X-PAR data. The most significant difference was in the rear of the tornadic storm where the convection was suppressed by inflow. It can be seen that at both heights, the Z intensity in the rear inflow region from experiment VrZh_S6X6 (Figure 6d, blue to green color) was clearly lower than that of VrZh_S6 (Figure 6b), and the Z structure of experiment VrZh_S6X6 was much closer to the observation (Figure 3d). The difference of analyzed vertical structure was also obvious. Without X-PAR data assimilation, the overall intensity of convection was overestimated (Figures 7b and 8b), while the assimilation of X-PAR data greatly reduced the Z intensity at both heights (Figures 7d and 8d). By comparing experiment VrZh_S6, VrZh_S6_XZ6 and VrZh_S6X6 (Figures 6b–6d and 7b–7d), it could be concluded that the reduction of Z intensity in this region can be mostly attributed to the assimilation of additional Z from NH X-PAR. The assimilation of additional X-PAR Vr further improved the convective system structure in the rear inflow region (see Figures 8c and 8d). The assimilation of X-PAR Vr improved the dynamic structure at low levels, resulting stronger rear inflow and stronger vertical vorticity. The spurious convection in the rear inflow region was further suppressed.

The RMSE results are consistent with the qualitative analysis (Figure 9). For the verification against Vr observation, it is shown that the assimilation of X-PAR data has little impact on the 2.4° tilt (Figure 9a), where the black and green lines overlay each other. Not surprisingly, the averaged RMSEs of Vr was reduced by 1.4 m s⁻¹ at the low levels after assimilation of X-PAR data (Figure 9c). As discussed above, this reduction of analysis error was mostly contributed by X-PAR Vr, as we could see the RMSEs of VrZh_S6_XZ6 are even higher than VrZh_S6 (Figure 9c). As to Z, the average RMSE of Z in VrZh_S6X6 was 0.5 dBZ lower than that of VrZh_S6 at 2.4° tilt, and 1.0 dBZ lower at 0.9° tilt (Figures 9b and 9d). For the RMSEs of Z, the performance of VrZh_S6_XZ6 and VrZh_S6X6 were similar, which further confirmed the improvement of Z structure was from the assimilation of Z. In all, the X-PAR data well compensated S-band radar observation at lower levels. The assimilation of X-PAR data greatly improved the low-level wind and storm structure.

The Impact of X-PAR Data Assimilation Frequency

The biggest advantage of X-PAR is its high observation frequency. In this study, as we mentioned above, X-PAR data are available every 2 min. Two experiments, VrZh_S6X6 and VrZh_S6X2, were conducted to examine whether the increased assimilation frequency helps to improve the analysis and forecast accuracy. Both experiments show obvious cyclonic vortices extending to the ground near the tornado report location within about a 4 km diameter (Figures 7d and 7e), demonstrating the assimilation of X-PAR data could improve the analysis of the tornadic vortex. The increased assimilation frequency helped to improve the vortex intensity, as we could see the maximum vortex intensity was stronger in experiment VrZh_S6X2 than in experiment VrZh_S6X6. The biggest difference was the wind structure in the rear inflow region. A downdraft could be clearly seen from 0.5 to 3 km AGL in the rear inflow region (Figure 8e), which did not appear in VrZh_S6X6 (Figure 8d). This increased descending flow in turn further suppressed the convection in that area. It can be seen that the Z of experiment VrZh_S6X2 was much weaker than that of VrZh_S6X6 in the rear inflow region (Figures 6d, 6e, 8d and 8e), which was more consistent with the observations (Figures 3d and 3h).

The RMSE results are consistent with the qualitative analysis. The increased assimilation frequency of X-PAR data improved the low-level analysis accuracy. At the 0.9° tilt, the average RMSE of Vr in experiment VrZh_S6X2 (red line) from 05:42 UTC to 06:00 UTC was 0.5 m s⁻¹ lower than that in VrZh_S6X6 (green line), and the average RMSE of Z in VrZh_S6X2 from 05:42 UTC to 06:00 UTC was 1.5 dBZ lower than that in VrZh_S6X6 (Figures 9c and 9d). As we explained above, the increased downdraft in the rear inflow region could suppress convection, which in turn improved Z structure. At the 2.4° tilt, the average RMSEs of Z and Vr in experiment VrZh_S6X2 were lower than that of experiment VrZh_S6X6, but the reduction of RMSEs were not as obvious as the result against NH X-PAR 0.9° tilt (Figure 9b). The increase of assimilation frequency of NH radar has limited impact on the upper level storm structure. Overall, the increased assimilation frequency could further improve the accuracy of low-level wind structure and the convective storm structure.

Forecasting Results

The tornado occurred at 06:02 UTC and lasted for 30 s. To see the impact of X-PAR data assimilation on the tornadic storm forecast, short deterministic forecasts from the ensemble mean analysis every 6 min since 05:42 UTC were produced. Here, we only examined the forecasting results of experiment Vr_S6, VrZh_S6, VrZh_S6X6 and VrZh_S6X2. Experiment VrZh_S6_XZ6 performed worse than VrZh_S6X6 on the analysis. Therefore, it is not shown. PRD emulator (Jung et al., 2008) developed by CAPS was used to simulate the Z and Vr from WRF forecasts. Here, the NH X-PAR radar data were used as observations.

Figures 10 and 11 show the forecast Vr and Z. All the forecasts were at 06:02 UTC when the tornado touched down. Here, we only plotted the forecast Vr initialized from 05:54 UTC (8-min forecast, first row of Figure 10) and 06:00 UTC (2-min forecast, second row of Figure 10). Other forecasts initialized from earlier analysis time did not show any sign of vortex structure and are not presented. For the Vr observation at 0.9° tilt, a strong circulation existed to the southeast of the radar site at 06:02 UTC (Figure 10i, see the positive and negative radial velocity couplet at the center). Without X-PAR data assimilation, both experiment Vr_S6 and VrZh_S6 failed to predict the vortex structure both in the 8-min forecast (Figures 10a and 10b) and in the 2-min forecast (Figures 10e and 10f). With X-PAR data assimilation, the 8-min forecast of VrZh_S6X6 still failed to predict the tornadic storm vortex (Figure 10c), while the 2-min forecasts produced a positive and negative radial velocity couplet similar to the observations (Figure 10g), albeit with relatively weaker intensity. This suggests that the low-level wind information was not well retained by the forecast model and was filtered out as soon as the model launched the forecast. With higher assimilation frequency, the predicted vortex structure of VrZh_S6X2 is clearly improved. For the 8-min forecast, a weaker radial velocity couplet was predicted at the place of tornado occurrence (Figure 10d). Similar results could also be found in Figure 12, which shows the 8-min forecasts at 200 m AGL. Among all the experiments, only VrZh_S6X2 exhibited a remarkable vortex at low levels (see black ellipse in Figure 12d), demonstrating the advantage of increasing assimilation frequency. For the 2-min forecast, the intensity of vortex was much closer to the observations (Figure 10h), making it the best among all the forecasts. Additionally, the predicted radial wind structure outside the vortex center was also improved, demonstrating the low-level wind information was well absorbed by the forecast model.

Figure 10. Forecasted NH X-band Phased-Array radar (X-PAR) Vr (0.9° tilt) at 06:02 UTC starting from 05:54 UTC from different experiments (a) Vr_S6, (b) VrZh_S6, (c) VrZh_S6X6, and (d) VrZh_S6X2. (e–h) are the same as (a–d) but for the forecasts starting from 06:00 UTC. (i) is the NH X-PAR observation at 06:02 UTC. The marker “+” near the domain center indicates the location of the reported tornado.

Figure 11. Observed Z at 06:02 UTC from 0.9° tilt of NH X-band Phased-Array radar (shaded), and the predicted Z (>40 dBZ, black contour) from different experiments (a) Vr_S6, (b) VrZh_S6, (c) VrZh_S6X6, and (d) VrZh_S6X2. All the forecasts for (a–d) were initialized at 05:54 UTC. (e–h) are the same as (a–d) but for the forecasts start from 06:00 UTC. The marker “+” near the domain center indicates the location of the reported tornado.

Figure 12. The forecast Z, and storm-relative wind vectors at 200 m AGL, valid at 06:02 UTC, from experiments (a) Vr_S6, (b) VrZh_S6, (c) VrZh_S6X6, and (d) VrZh_S6X2. All the forecasts started from 05:54 UTC. White contours, marker “+,' black ellipse and dashed box are the same as those in Figure 6.

In Figure 11, we also examined the forecast Z structure. The filled color contours are the observed Z at 06:02 UTC while the black line is the predicted Z larger than 40 dBZ. Here, similar as radial wind, the observations from the NH X-PAR were used for the forecast verification. By comparing the area enclosed by the black line and the red contour area, it can be seen that experiment VrZh_S6X2 is the best among all the experiments for both 8-min and 2-min forecast (Figures 11d and 11h), followed by VrZh_S6X6 (Figures 11c and 11g), while Vr_S6 is the worst (Figures 11a and 11e). Similar to the radial wind evaluations, the 2-min forecasts of radar reflectivity were better than that of the 8-min forecasts as indicated by a better overlap of the black line with the red shading area for the 2-min forecasts.

Figure 13 shows RMSEs of Z and Vr from forecasts starting from different analysis time against NH X-PAR data at 06:02 UTC. For example, one point at 05:42 represents the RMSE of 20-min forecast starting from 05:42 UTC, while one point at 05:48 represents the RMSE of 14-min forecast starting from 05:48 UTC. For the RMSEs of forecast Vr, Experiment VrZh_S6X2 is the best since 05:48 UTC, followed by VrZh_S6X6, while Vr_S6 is the worst. These results were consistent with previous subjective evaluations. Note that the data assimilation frequency was improved starting from 05:42 UTC in VrZh_S6X2. It was not surprising that experiments VrZh_S6X6 and VrZh_S6X2 starting from 05:42 UTC has gotten similar forecast errors. As is explained above, the assimilated low-level wind information of the tornadic storm was probably filtered out during the forecast if observations were assimilated only once. With continued assimilations, the model was able to maintain the correct low-level wind structure. The higher the assimilation frequency, the quicker the model adjusted to the observations. For the forecast starting from 05:48 UTC, the 14-min forecast RMSEs of Vr in VrZh_S6X2 was 0.03 m s⁻¹ lower than that of 2-min forecast in VrZh_S6, which means the assimilation of high frequency X-PAR data was able to predict better dynamic tornadic storm structure 12 min earlier than experiments with only S-band radar data assimilation. The observed low-level wind information by X-PAR well compensated the missing information of S-band radar and is thus the key for a successful tornadic storm prediction.

Figure 13. Forecast RMSEs at 06:02 UTC starting from different analysis time of (a) Vr and (b) Z. The RMSEs were calculated against 0.9°–4.5° tilt of the NH X-band Phased-Array radar observations. X-labels represents the different analysis time when we launched the deterministic forecasts. For example, the point at 05:42 means the 20-min forecast RMSEs while the point at 05:48 means 14-min forecast RMSEs.

The improvement in Z was not as obvious as Vr when X-PAR data were assimilated. It is until the last assimilation cycle that the forecast RMSEs of VrZh_S6X2 was smaller than that of VrZh_S6 (Figure 13b), which is probably due to the large uncertainties in the microphysics scheme or imbalance caused by rapid update (Bick et al., 2016; Zeng, de Lozar, et al., 2021). The microphysics process takes more time to spin up and adjust, but the direct assimilation of Z is still an effective way to improve the prediction of Z. It can be seen that the average RMSE of Z in VrZh_S6 is 1.2 dBZ lower than that in Vr_S6. As for the Z structure, the assimilation of 2-min X-PAR data still has 4–5 min time in advance in terms of producing a meaningful forecast when compared to the experiments without X-PAR data assimilation. In all, with the assimilation of high frequency X-PAR data, it is possible to improve the predictability of tornadic storms a few minutes in advance.

Summary

On June 8, 2018, a tornado embedded within a rainband of typhoon Ewiniar touched down near Foshan, Guangdong. The evolution and dissipation of the tornadic storm was well captured by one quasi-operational X-PAR and two operational S-band radars. In this study, the impact of X-PAR data assimilation on the analysis and forecast of the tornadic storm was investigated by using the WRF-based EnKF system developed by Meng and Zhang (2008a, 2008b). Furthermore, we added the capability of Z assimilation into this WRF-based EnKF system to investigate the impact of assimilating Z. Forty-five ensemble members were used, each of which had two one-way nested WRF model domains with grid spacing of 2.5 km for the outer domain and 500 m for the inner domain. The radar observations were only assimilated in the inner domain starting from 05:00 UTC.

Five experiments with the purpose of examining the impact of Z assimilation (Vr_S6, VrZh_S6) and the extra X-PAR data assimilation (VrZh_S6_XZ6, VrZh_S6X6, and VrZh_S6X2) on the analysis of the tornadic storm were conducted. The assimilation of S-band radar Z improved the overall Z structure of the tornadic storm. With both Vr and Z assimilation, the convective system could reach its balance state faster, resulting in smaller RMSEs than the experiments with only Vr assimilation for the later analysis cycles. The X-PAR observations well compensated the S-band observation at the low levels. Without X-PAR data assimilation, experiment VrZh_S6 failed to reproduce the tornadic vortex at 200 m height. With extra X-PAR data assimilation, the vortex signature was well reproduced, and the predicted vortex was greatly enhanced. Moreover, the vortex extended to the ground (a key signature for tornado touch down). Similar to the conclusion of S-band radar data assimilation, the assimilation of X-PAR Z improved the accuracy of the low-level Z structure, especially in the rear inflow weak echo region, while the assimilation of X-PAR Vr improved the accuracy of the low-level wind structure and further suppressed the convection in the same area. This was more obvious when the X-PAR data was assimilated at a higher frequency. A downdraft could be identified in the rear inflow region for experiment VrZh_S6X2. Verifications against radar observations indicated that experiments with X-PAR data assimilation had lower analyzed RMSEs, especially at the low levels. Experiment VrZh_S6X2, which assimilated X-PAR data every 2-min, had the lowest RMSEs.

Deterministic forecasts were also launched every 6 min from the cycling analysis and forecast system since 05:42 UTC to see if extra X-PAR data assimilation could improve the predictability of tornadic storms. All those forecasts ended at 06:02 UTC, when the tornado occurred. With a 2-min X-PAR data assimilation frequency, experiment VrZh_S6X2 was able to predict the positive and negative radial velocity couplet 8 min in advance, though the intensity was a little weaker than the observed one. Neither experiments with either lower assimilation frequency (VrZh_S6X6) nor experiments without X-PAR data (VrZh_S6) assimilation predicted this key signature of the tornadic storm. Base on the forecast RMSEs of Vr, the experiment with assimilation of high frequency X-PAR data could improve the predictability of tornado 12 min in advance than the experiment without X-PAR data assimilation. However, this advantage would be shortened to 4–5 min based on the forecast RMSEs of Z.

Overall, the assimilation of X-PAR is helpful for building a realistic tornadic storm structure, especially at the low levels where operational S-band radars do not have data. In this study, the tornadic storm forecasts were also improved based on the more realistic analysis fields.

Finally, we note that although the conclusions of this study are based on one single case, detailed examination of the results leads us to believe that our results have general meaning, at least for low-level small-scale weather disasters when a similar data assimilation technique such as the one described in this paper is used. We should also note that the tornado in this case is very weak, so its predictability might be limited. Therefore, even in the best experiment, the tornado was predicted only 8 min in advance. For stronger tornado cases, it is anticipated that assimilating PAR data with high assimilation frequency might be able to have longer warning lead time. Moreover, if we want to look into longer forecast lead time in the future, filter method that can reduce model imbalance, such as the integrated mass-flux adjustment filter (Zeng, de Lozar, et al., 2021), perhaps will need to be considered. In summary, this study demonstrates the potential ability and advantage for X-PAR in improving the prediction of small-scale weather disasters such as tornadoes.

Acknowledgments

This work was primarily supported by the National Natural Science Foundation of China (grants 42025501, 41875053, 61827901), the National Key Research and Development Program of China (grant number 2017YFC1501703), the 5th “333 High-level Personnel Training Project” of Jiangsu Province (BRA2019037) and the Open Research Program of the State Key Laboratory of Severe Weather (2020LASW-A01). We acknowledge Guangzhou Meteorological Bureau and Zhuhai Naruida Technology Co., Ltd. for collecting and archiving the radar data.

Data Availability Statement

The data which supports the analysis and conclusions of this paper are available at https://doi.org/10.5281/zenodo.4384032. For more detailed data, please contact us via e-mail.