## 2 Methods
### WRF Model Configuration and Experimental Design in Data Assimilation  
In this study, the WRF model (version 3.8.1; Ref. [38]) and a data assimilation system based on WRF Variational Analysis (WRFVAR, version 3.9.1; Ref. [39]) were used for forecasting. The Doppler lidar wind profiler data were assimilated into the WRF-3 DVAR system for prediction of the AT event on 5 August 2018. The WRF model domain is centered on Taiwan, and has 451 \(\times\) 451 horizontal grid points with horizontal resolutions of 2 km \(\times\) 2 km, covering southeast China (19.5\({}^{\circ}\)-27.8\({}^{\circ}\) N, 116.4\({}^{\circ}\)-125.6\({}^{\circ}\) E). The domain has 52 vertical layers from the surface to a 20-hPa isobaric layer using a terrain following sigma coordinate. The physics schemes used in all domain include the Rapid Radiative Transfer Model longwave radiation and shortwave radiation [40], the Yonsei University planetary boundary layer scheme [41], the Unified Noah land surface scheme [42], and the Goddard microphysics scheme [43] instead of the cumulus scheme. Figure 2: Comparisons between wind observations of (**a**) the initial profile and (**b**) thinning profile at Pingtung Airport from 0100 to 0400 UTC (0900 to 1200 LST) on 5 August 2018. Abnormal winds are represented by black circles in (**a**). Radar data and surface observations including those from manually operated meteorological stations and AWSs provided the initial and boundary conditions of the WRF model, namely the Radar WRF operated by the Taiwan CWB (CWB RWRF). This WRF-3 DVAR system assimilating both radar and surface observations is an assimilation forecasting system for convective scale and its lateral boundary forcing is updated every 1 h with a 13-h length. The WRF-3 DVAR system can be categorized as a data assimilation system that provides an analysis \(x^{a}\) through minimizing the cost function J(x) [44] as follows. \[\mathrm{J}(\mathbf{x})=\frac{1}{2}\Big{(}x_{f}^{\prime}\Big{)}^{T}B^{-1}\Big{(} x_{f}^{\prime}\Big{)}+\frac{1}{2}(y-H[x])^{T}R^{-1}(y-H[x])=J_{b}+J_{o} \tag{1}\]  
\[x_{f}^{\prime}=x-x^{b} \tag{2}\]  
where J(x) is the cost function of the atmospheric state. The background term \(J_{b}\) plays a key role in the performance of data assimilation and it represents the analysis increment \(x_{f}^{\prime}\). The observation term \(J_{o}\) represents the distance between the vector of observation \(y^{o}\) and the analysis observation vector transformed by the observation operator \(H[x]\), which maps state space to observation space. The analysis increment \(x_{f}^{\prime}\) is defined as the distance between the analysis state \(x\) and the initial background field value \(x^{b}\). The terms \(B\) and \(R\) are the covariance matrix of background and observation errors, respectively. In practice, \(R\) is often assumed to be static and diagonal, and is provided prior to the assimilation [45]. \(B\) is calculated using the approach of the National Meteorological Center [46] with option 3 (CV3; Ref. [47]) to build the background error covariance in 3 DVAR. To evaluate the impact of the data assimilation methods of the lidar wind profiler on AT prediction near Pingtung Airport, three groups--C1, C2 and C3--were used in our experimental design. In the C1 group, the experiments did not include data assimilation (Table 1) and the initial and lateral boundary conditions with a horizontal resolution of 2 km were obtained from the Taiwan CWB RWRF. Subsequently, 12-h forecasts on a 2-km grid were launched at the starting times of 0100, 0200, 0300, and 0400 UTC on 5 August 2018, for the experiments NDA01, NDA02, NDA03, and NDA04, respectively (0100-1300, 0200-1400, 0300-1500, and 0400-1600 UTC, respectively; Table 1 and Figure 3). Figure 3: Schematic diagram of data assimilation-forecast strategy for the 2-km mesh for experiments of C1(NDA01, NDA02, NDA03 and NDA04), C2(DA01, DA02, DA03 and DA04) and C3 (CDA01, CDA02 and CDA03). The black box represents the major impact period of the thunderstorm event.