# Evaluation of Lightning Prediction by an Electrification and Discharge Model in Long-Term Forecasting Experiments

Authors: Liangtao Xu; Shuang Chen; Wen Yao
Journal: Advances in Meteorology
 (2022)

---

## Abstract

Over nearly three rainy seasons of lightning activity in North China, numerical prediction experiments were carried out using the Weather Research and Forecasting model coupled with electrification and discharge schemes (WRF-Electric). The numerical forecast results were evaluated using the neighborhood-based equitable threat score (ETS) and fraction skill score (FSS) verification methods based on nationwide observational lightning data. An algorithm was used to generate the coverage of the total flash (intracloud and cloud-to-ground flashes) by fitting the cloud-to-ground flash data. The numerical results showed that the region of lightning activity could be well predicted by the mesoscale WRF-Electric model, particularly during a 6–12-hour forecasting period. The average ETS score of the 6–12-hour forecasting period was 0.34 for a 20 km neighborhood radius. The predictive skill of the model varied not only monthly but also diurnally. The model showed better forecasting skills during the main rainy season (June–July–August) and at 14 : 00–20 : 00 local time. The predictability of the model was enhanced with increasing thunderstorm scale. On the other hand, the coverage of predicted lightning activity was relatively concentrated, and the lightning flash density was higher than the observations. The main discrepancies in the model prediction were related to the design of the discharge parameterization. Thus, in discharge parameterization, the initial threshold for lightning should be modified according to the model resolution, while the magnitude of the neutralization charge in a single discharge should be referenced to the observational results.


---

## Body

## 1. Introduction

It is important to predict lightning activity better and in a timely fashion, because it plays a critical role in the aerospace, energy, power, and travel industries. For the 0–2-hour nowcasting of lightning activity, an extrapolation method based on multiple storm tracking algorithms or the machine-learning method can be applied using radar, satellite, and lightning detection datasets. Such an approach can successfully predict lightning activity within this timeframe [1, 2]. However, the accuracy of extrapolation-based nowcasting decreases rapidly after 2 hours [3]. In this case, a numerical weather prediction method is usually adopted instead or in combination to predict lightning activity over longer lead times [4, 5].Two kinds of the method are employed to predict lightning activity using a numerical weather model. One method is based on the predicted dynamic and microphysical fields in association with the lightning activity. This method is called “diagnosed parameterization of lightning” [6, 7] and involves establishing a lightning prediction proxy. A statistical model is run according to the characteristics of the lightning activity and large-scale atmospheric dynamic and thermal parameters (e.g., CAPE, CIN, K-index, or lifted index), and then the lightning activity can be forecasted using the prediction field obtained from the numerical simulation [8–13].Previous researchers have also built several diagnostic variables concerning lightning activity according to the relationship between lightning and cloud-scale properties [14–16]. McCaul et al. [17] used the model-derived graupel flux in convective clouds and the total ice content to obtain a statistical relationship between these two parameters and the total lightning flash density. Yair et al. [18] proposed a lightning potential index (LPI) to predict the potential trend of lightning occurrence. They indicated that the LPI could predict lightning activity better than the traditional thermodynamic indices. A machine-learning framework can also be employed to diagnose the occurrence of lightning based on numerical prediction fields [19, 20].The second method is referred to as “physical parameterization of lightning.” Predictive equations for the charge density of various hydrometeors are introduced into the numerical model, which predicts the evolution of the total charge density of the cloud. The electric field in the cloud is calculated using the predicted charge density of the thunderstorm. Meanwhile, parameterization of lightning is developed based on the physical process of lightning. This method conforms to the physical process of lightning occurrence. As a result, with the progress of numerical models, the physical parameterization method has a greater potential to predict lightning activity than the diagnosed parameterization method.A mesoscale numerical model coupled with electrical processes is generally referred to as mesoscale electrification and discharge model. Such models have been used to predict lightning activity in thunderstorms [21–24]. For example, Wang [25] introduced electrical processes into the GRAPES (Global/Regional Assimilation Prediction System) model and performed predictive experiments on previously observed regional lightning activity. Xu et al. [26, 27] coupled multiple electrification and discharge schemes to the Advanced Research version of the Weather Research and Forecasting (WRF-ARW) model, referred to as the WRF-Electric model. They performed predictive experiments on lightning activity using bulk lightning parameterization schemes. Fierro et al. [28] introduced a more detailed lightning parameterization scheme into the WRF model (known as the E-WRF model) and examined the results from predictive tests of various kinds of storms. The E-WRF model was evaluated in 10 high-impact weather cases, and the simulated flash origin density fields gave good agreement with the observations [29]. Lynn [30] and Wang et al. [31] indicated that lightning assimilation could significantly improve the ability of models to predict lightning activity. Lightning prediction experiments have also been carried out using RAMS (Regional Atmospheric Modeling System) coupled with electrical processes [32].A month-long lightning prediction experiment using the diagnosed parameterization in a forecasting model was conducted by Wilkinson [33]. However, lightning prediction experiments using mesoscale electrical models have mainly been applied in case studies. Such an approach is very useful in exploring lightning prediction methods, but it is not easy to evaluate the predictive skill of the model using only case studies. Therefore, based on the WRF-Electric model, a regional (North China) lightning forecasting platform was established at the beginning of June 2015, and lightning prediction experiments were carried out four times daily. In the present study, using the lightning prediction results during the rainy seasons (May to September) of 2015–2017, the predictive skill of the WRF-Electric model was objectively evaluated, providing a basis for further improvements in the electrical part of the model.The rest of this paper is organized as follows. The description of the model and setup is illustrated in Section2. The data and verification method are described in Section 3. The verification results and discussion are provided in Section 4. A summary and conclusions are given in Section 5.


## 2. Model and Setup

### 2.1. Model Description

The WRF-ARW model (version 3.4.1) [34] coupled with electrification and discharge schemes [26, 27, 35] was used for the simulation. Figure 1 shows a schematic illustration of the model. The discharge initiation threshold (defined as Ebreak) in the bulk lightning parameterization is the height-varying electric field, calculated as follows [36]:(1)Ebreak=±167ρairzβ,where ρairz is the (nondimensional) air density, which varies with height. The units of Ebreak are kV m−1. β is a limit coefficient. In the high-resolution (hundreds of meters) cloud model, β is usually set to 1.0. When the resolution of the model is several kilometers, the high-value area of the electric field may be smoothed, so the value of β in this study is set to 0.5.Figure 1
Schematic illustration of the WRF-Electric model.The discharge process is triggered when the vertical electric field in one grid point exceeds the height-varying initiation threshold, and then the model will reduce the charge density at the grid point where the absolute electric field exceeds 30 kV m−1 by 20%. The discharge scheme does not distinguish between lightning types (intracloud or cloud-to-ground flash) or between possible multiple stroke processes in a single lightning flash process, so we link one discharge process in the model to a complete lightning flash process in the observation.The two electrical models developed by Fierro et al. [28] and Xu et al. [26, 27] are established by introducing the physical processes of electrification and discharge into the WRF-ARW model. The adopted electrification processes are mainly noninductive charging schemes based on classic laboratory studies. The major difference between the two models are the coupled microphysical scheme and the discharge scheme. The NSSL two-moment scheme and a relatively complex discharge scheme were used by Fierro et al. [28] whereas the Milbrandt two-moment scheme and a bulk discharge scheme were used by Xu et al. [26, 27].


### 2.2. Settings for the Continuous Prediction Experiments

An automatic forecasting system was built based on the WRF-Electric model. Forecasting results for the 24-hour regional lightning activity are provided every 6 hours by the system. The forecast result is mainly the accumulated grid point flash number (FN, for total flash) per hour. The forecasting system has been run continuously since 1 June 2015, and nearly three rainy seasons of forecasting data have been accumulated. The verification work was conducted based on these forecasting data.The National Centers for Environmental Prediction Global Forecast System data at a resolution of 1° and 3 h intervals [37] were used for the initial and lateral boundary conditions. Two-way nesting domains were used in the experiment, with mesh sizes of 124 × 124 and 160 × 160 and horizontal grid resolutions of 12 km and 4 km, respectively. The center of the domains was (40°N, 116.2°E), and both domains used standard static nesting without a prescribed set of moves. Figure 2 shows the terrain height and locations of the two domains. In the vertical regime, 28 uneven levels were used for all the meshes. The model top was defined as 50 hPa. Only two layers of nested domains were used here, and the vertical resolution is relatively low. It is expected that the computational cost can be reduced under relatively high horizontal resolution.Figure 2
Terrain height and locations of the two domains. Red points denote the locations of the three SAFIR-3000 substations. The abbreviation of BJ, HB, and TJ represent Beijing, Hebei, and Tianjin.The physical options included the Milbrandt two-moment microphysical parameterization [38, 39] with the noninductive charging scheme adapted from research by Gardiner et al. [40], Ziegler et al. [41], and Tan et al. [42] (TGZ scheme), the Rapid Radiative Transfer Model (RRTM) longwave radiation scheme [43], the Dudhia shortwave radiation scheme [44], the Kain–Fritsch cumulus scheme [45], and the Yonsei University (YSU) planetary boundary layer scheme [46]. No cumulus parameterization was used in the 4 km mesh. The two meshes were integrated for 24 h with time steps of 15 s and 5 s (see Table 1 for more details). Only the high-resolution forecast results (4 km) were used in this study.Table 1
Model design.
ParameterDomain 1Domain 2Dimensions (x, y)124 × 124160 × 160Grid size (km)124Time step (s)4515Integration hours2424Boundary layer schemeYSUYSUMicrophysics schemeMilbrandt two-momMilbrandt two-momCharging schemeTGZTGZCumulus schemeKain–FritschNoLateral boundary conditionGFS 1° × 1°NestedInitial conditionGFS 1° × 1°NestedLongwave radiationRRTMRRTMShortwave radiationDudhiaDudhia


## 2.1. Model Description

The WRF-ARW model (version 3.4.1) [34] coupled with electrification and discharge schemes [26, 27, 35] was used for the simulation. Figure 1 shows a schematic illustration of the model. The discharge initiation threshold (defined as Ebreak) in the bulk lightning parameterization is the height-varying electric field, calculated as follows [36]:(1)Ebreak=±167ρairzβ,where ρairz is the (nondimensional) air density, which varies with height. The units of Ebreak are kV m−1. β is a limit coefficient. In the high-resolution (hundreds of meters) cloud model, β is usually set to 1.0. When the resolution of the model is several kilometers, the high-value area of the electric field may be smoothed, so the value of β in this study is set to 0.5.Figure 1
Schematic illustration of the WRF-Electric model.The discharge process is triggered when the vertical electric field in one grid point exceeds the height-varying initiation threshold, and then the model will reduce the charge density at the grid point where the absolute electric field exceeds 30 kV m−1 by 20%. The discharge scheme does not distinguish between lightning types (intracloud or cloud-to-ground flash) or between possible multiple stroke processes in a single lightning flash process, so we link one discharge process in the model to a complete lightning flash process in the observation.The two electrical models developed by Fierro et al. [28] and Xu et al. [26, 27] are established by introducing the physical processes of electrification and discharge into the WRF-ARW model. The adopted electrification processes are mainly noninductive charging schemes based on classic laboratory studies. The major difference between the two models are the coupled microphysical scheme and the discharge scheme. The NSSL two-moment scheme and a relatively complex discharge scheme were used by Fierro et al. [28] whereas the Milbrandt two-moment scheme and a bulk discharge scheme were used by Xu et al. [26, 27].


## 2.2. Settings for the Continuous Prediction Experiments

An automatic forecasting system was built based on the WRF-Electric model. Forecasting results for the 24-hour regional lightning activity are provided every 6 hours by the system. The forecast result is mainly the accumulated grid point flash number (FN, for total flash) per hour. The forecasting system has been run continuously since 1 June 2015, and nearly three rainy seasons of forecasting data have been accumulated. The verification work was conducted based on these forecasting data.The National Centers for Environmental Prediction Global Forecast System data at a resolution of 1° and 3 h intervals [37] were used for the initial and lateral boundary conditions. Two-way nesting domains were used in the experiment, with mesh sizes of 124 × 124 and 160 × 160 and horizontal grid resolutions of 12 km and 4 km, respectively. The center of the domains was (40°N, 116.2°E), and both domains used standard static nesting without a prescribed set of moves. Figure 2 shows the terrain height and locations of the two domains. In the vertical regime, 28 uneven levels were used for all the meshes. The model top was defined as 50 hPa. Only two layers of nested domains were used here, and the vertical resolution is relatively low. It is expected that the computational cost can be reduced under relatively high horizontal resolution.Figure 2
Terrain height and locations of the two domains. Red points denote the locations of the three SAFIR-3000 substations. The abbreviation of BJ, HB, and TJ represent Beijing, Hebei, and Tianjin.The physical options included the Milbrandt two-moment microphysical parameterization [38, 39] with the noninductive charging scheme adapted from research by Gardiner et al. [40], Ziegler et al. [41], and Tan et al. [42] (TGZ scheme), the Rapid Radiative Transfer Model (RRTM) longwave radiation scheme [43], the Dudhia shortwave radiation scheme [44], the Kain–Fritsch cumulus scheme [45], and the Yonsei University (YSU) planetary boundary layer scheme [46]. No cumulus parameterization was used in the 4 km mesh. The two meshes were integrated for 24 h with time steps of 15 s and 5 s (see Table 1 for more details). Only the high-resolution forecast results (4 km) were used in this study.Table 1
Model design.
ParameterDomain 1Domain 2Dimensions (x, y)124 × 124160 × 160Grid size (km)124Time step (s)4515Integration hours2424Boundary layer schemeYSUYSUMicrophysics schemeMilbrandt two-momMilbrandt two-momCharging schemeTGZTGZCumulus schemeKain–FritschNoLateral boundary conditionGFS 1° × 1°NestedInitial conditionGFS 1° × 1°NestedLongwave radiationRRTMRRTMShortwave radiationDudhiaDudhia


## 3. Data and Methodology

### 3.1. ADTD Dataset

At present, there is no total lightning detection network covering the whole prediction period or region (inner domain of Figure2). Therefore, the nationwide cloud-to-ground (CG) lightning observation data from the Advanced Time of Arrival and Direction (ADTD) system network of China were used to compare with the numerical model results. The ADTD system was initiated by the Chinese Academy of Sciences and the China Meteorological Administration in 2003 and is based on the very-low-frequency/low-frequency time difference direction hybrid positioning technique. Li et al. [47] have given a detailed introduction to the development of the ADTD system. The ADTD system has been operated since 2003 and includes 357 stations covering most of China. Many studies of CG flash characteristics in China are based on ADTD data [48–50].The ADTD system is generally assumed to have a detection efficiency of 90% and a positioning accuracy of 500 m over a range of 150 km [48], although there is a lack of systematic research into these features. Also, the system provides measurements of the time, strike location, polarity, and peak current of CG flashes. The Chinese Academy of Meteorological Sciences has conducted artificially triggered lightning experiments in South China. According to the data from a triggered lightning test (six return strokes), the average positioning accuracy of the ADTD system for the return stroke is 3768 m.The prediction experiments were run at a horizontal grid size of 4 km. To better compare the prediction results, the observed lightning location data were also processed on a grid of 4 km, and the number of lightning flashes in every 4 km × 4 km grid was counted to be used as the observed value for comparison with the predicted results. The flash density is defined as the number of lightning flashes in each 4 km × 4 km grid in this study, for which the unit is flashes/16 km2.


### 3.2. SAFIR Dataset

The SAFIR-3000 system is a total flash positioning system based on very high frequency (VHF) and low frequency [51]. It had three substations in the Beijing–Tianjin area (as shown in Figure 2). The system was well maintained during the period 2005–2008. The theoretical detection efficiency is 95%, and the positioning accuracy is 500 m. It is generally accepted that the positioning accuracy in the core detection zone is 2 km and the detection efficiency is up to 90% [52].The data period used in this study was from June to September in 2007–2008. The analysis period was 0600–1200 UTC each day (1400–2000 BTC, local Beijing time). This is the period during which thunderstorms frequently occur. The analysis region covered was (39.3–40.3°N, 116–118°E), which is located in the core detection zone. First, the VHF radiation sources and CG return strokes were separately grouped into flashes [51]. The SAFIR_CG data represent the CG flashes, and the SAFIR_Tot data represent the total flashes.


### 3.3. A Fitting Algorithm

We aimed to design a simple algorithm to obtain the coverage of the total flash by fitting the CG flash data. The fitted lightning flash occurrence region is a better representation of the actual total lightning occurrence region. Although a few of the thunderstorms contained only intracloud (IC) flashes or only CG flashes, most of the thunderstorms contained both IC flashes and CG flashes, and the occurrence of IC flashes and CG flashes is closely linked. Lightning will occur when there is an accumulation of charge in clouds and IC flashes often develop before CG flashes [53]. Whether the lightning develops into an IC flash or a CG flash depends on the charge structure in the clouds, although it may involve some stochastic processes. The proportion of CG flashes to total flashes in thunderstorms is about 20–50% with regional variations [52, 54, 55]. Total flashes occur over larger areas than CG flashes.When a CG flash occurs at a statistical point (e.g., the red point in Figure3) on the timescale of a thunderstorm (cumulative over several hours), we assume that some IC flashes may occur (blue points in Figure 3(a)) around the CG flash point. The blue points are defined as preliminary guesses for IC flashes. When a grid point is guessed at least twice, then it is assumed that an IC flash occurs in this grid (purple points in Figure 3(b)). A new dataset can be generated by processing (fitting) the CG flash data. The fitting algorithm should first be evaluated with a total flash locating system.Figure 3
Sketch map of the simple algorithm for fitting total flash.
(a)(b)A new fitting flash dataset referred to as SAFIR_CG_fitted can be obtained by processing the SAFIR_CG data. No specific fitting analysis is performed when the number of grids with lightning flash occurrences is ≤2 during a 6-hour period.Figure4 shows the fitting results for a 6-hour period (06 : 00–12 : 00 UTC) on 30 June 2008. The area of occurrence of the SAFIR_CG data is smaller than that of SAFIR_Tot, although the occurrence area of the SAFIR_CG data is often located in the kernel areas of the SAFIR_Tot data. The algorithm expands the CG flash to match the relationship between the CG flash and the total area of flash occurrence. Compared with the SAFIR_CG data, the area of occurrence of the SAFIR_CG_fitted data is more consistent with that of the SAFIR_Tot data. Therefore, from qualitative analysis, the fitting algorithm is effective without introducing too many false signals. The relative detection efficiency (RDE) was defined to measure the performance of the CG datasets.(2)RDE=NA∩NBNB,where NA represents the flash area covered by dataset A (CG or fitted CG) and NB represents the flash area covered by dataset B (total flash). The RDE can give the fraction of all flashes detected by B that were also detected by A. Some false lightning signals will be introduced using CG flash data to fit the total flash data. Thus, the concept of the false signal ratio (FSR) was introduced:(3)FSR=NA−NA∩NBNB.Figure 4
The SAFIR_Tot, SAFIR_CG, and SAFIR_CG_fitted data accumulated from 06 : 00 to 12 : 00 UTC on 30 June 2008.
(a)(b)(c)The fitting algorithm was evaluated using the SAFIR dataset from June to September in 2007–2008 with a 6 h interval. Compared with the original SAFIR_CG data, the SAFIR_CG_fitted data showed a significant improvement in the RDE. The average RDE of the SAFIR_CG data concerning SAFIR_Tot was 35%. The fitting algorithm can increase the RDE by an average of 20% after introducing a small FSR (about 10%). When the number of grids with lightning occurrence exceeded 25, the RDE increased by about 26–30% with a smaller FSR. The fitting algorithm can effectively bridge the spatial gap between the CG flash and the total flash.Because the ADTD system provides only CG flash data, the fitting algorithm generates the total flash from the CG flash. In the next part of forecast verification, CG_original represents the CG flash detected by ADTD, and CG_fitted represents data generated from the CG_original by the fitting algorithm.


### 3.4. Verification Method

It is difficult to assess the actual quality of the high-resolution model using traditional metrics (point-to-point method) because slight displacement errors often result in “double penalties” (i.e., errors observed but not forecast, or forecast but not observed). To avoid these double penalties, the neighborhood-based equitable threat score (ETS) reported by Clark et al. [56] was adopted in this study. This is an approach to validating forecasts on a high-resolution numerical grid. The verification is based on the gridded observed lightning flash data and model forecast flash data.The traditional method of computing the ETS [57] uses a 2 × 2 contingency table of possible forecast outcomes at individual grid points where the table elements are hits (the correct forecast of an event), misses (an observed but not forecast event), false alarms (a forecast but not observed event), and correct negatives (correct forecast of nonevent). Using these elements, the ETS is expressed as(4)ETS=hits−chancehits+misses+false alarms−chance,where(5)chance=hits+misseshits+false alarmshits+misses+correctnegative+false alarms.The average ETSs were computed by summing (i.e., aggregating) contingency table elements over all cases. The miss rate (MR) and false alarm ratio (FAR) were also calculated as follows:(6)MR=misseshits+misses,FAR=false alarmshits+false alarms.To compute a neighborhood-based ETS, the criteria for hits are relaxed by considering adjacent grid points within a specified radius of each grid point. Following Clark et al. [56], the ETS was calculated for lightning forecasts at radii of 0, 5, 10, 20, 30, 40, 60, and 90 km for thresholds referring to 1 flash/16 km2 (lightning flash density), which means that only the region in which lightning occurs is quantitatively evaluated (a yes or no evaluation). When the neighborhood radius is set to 0, the neighborhood-based approach reduces to the traditional form of ETS.Fraction skill scores (FSSs) [58] were also used in this study. FSS is a neighborhood-based verification method in which the fractional coverage of predicted and observed grid area to the total area is compared. The FSS can be more tolerant of small displacement errors than the ETS, making it more suitable for evaluating a simulation with fine-scale grids [31]. In the actual calculation, the threshold q needs to be set first to preprocess the data, and q is used to convert the observation field (Or) and forecast field (Mr) into binary fields Io and IM. All grid squares exceeding the threshold have a value of 1 and all others have a value of 0:(7)Io=1Or≥q0Or<q,IM=1Mr≥q0Mr<q.Io and IM are the observed and predicted fields processed in the space window, respectively. For a certain neighborhood radius, Io and IM can be calculated by selecting q, and FSS can be further calculated. In this study, FSS was calculated for lightning forecasts at radii of 4, 12, 20, 28, 36, 44, 52, 60, 68, 76, 84, 92, and 100 km for thresholds referring to 1 flash/16 km2 (flash density). The CG_fitted flash data of ADTD are used in the verification. Because the observed CG flashes (CG_original) have been counted in the 4 km × 4 km grid, according to the number of grids in which CG flashes have been observed (N) in the 6-hour period, all the inspection periods are divided into seven categories (N < 50, 50 ≤ N < 100, 100 ≤ N < 500, 500 ≤ N < 1000, 1000 ≤ N < 2000, 2000 ≤ N < 3000, and N ≥ 3000). All periods are evaluated for classification by the FSS method.Our study focused on the forecasting ability of the numerical model in a 6-hour period. Based on the differences in forecasting time, the 6–12-hour, 12–18-hour, and 18–24-hour lightning forecasts were examined. The initial 0–6-hour forecasts were discarded from most analyses and the model was spinning up during this period, so the lightning forecasts may be inaccurate. Meanwhile, the model is initialized at four different times a day so that there are four forecasting results for the same 6 h observation time with a 6 h interval.


## 3.1. ADTD Dataset

At present, there is no total lightning detection network covering the whole prediction period or region (inner domain of Figure2). Therefore, the nationwide cloud-to-ground (CG) lightning observation data from the Advanced Time of Arrival and Direction (ADTD) system network of China were used to compare with the numerical model results. The ADTD system was initiated by the Chinese Academy of Sciences and the China Meteorological Administration in 2003 and is based on the very-low-frequency/low-frequency time difference direction hybrid positioning technique. Li et al. [47] have given a detailed introduction to the development of the ADTD system. The ADTD system has been operated since 2003 and includes 357 stations covering most of China. Many studies of CG flash characteristics in China are based on ADTD data [48–50].The ADTD system is generally assumed to have a detection efficiency of 90% and a positioning accuracy of 500 m over a range of 150 km [48], although there is a lack of systematic research into these features. Also, the system provides measurements of the time, strike location, polarity, and peak current of CG flashes. The Chinese Academy of Meteorological Sciences has conducted artificially triggered lightning experiments in South China. According to the data from a triggered lightning test (six return strokes), the average positioning accuracy of the ADTD system for the return stroke is 3768 m.The prediction experiments were run at a horizontal grid size of 4 km. To better compare the prediction results, the observed lightning location data were also processed on a grid of 4 km, and the number of lightning flashes in every 4 km × 4 km grid was counted to be used as the observed value for comparison with the predicted results. The flash density is defined as the number of lightning flashes in each 4 km × 4 km grid in this study, for which the unit is flashes/16 km2.


## 3.2. SAFIR Dataset

The SAFIR-3000 system is a total flash positioning system based on very high frequency (VHF) and low frequency [51]. It had three substations in the Beijing–Tianjin area (as shown in Figure 2). The system was well maintained during the period 2005–2008. The theoretical detection efficiency is 95%, and the positioning accuracy is 500 m. It is generally accepted that the positioning accuracy in the core detection zone is 2 km and the detection efficiency is up to 90% [52].The data period used in this study was from June to September in 2007–2008. The analysis period was 0600–1200 UTC each day (1400–2000 BTC, local Beijing time). This is the period during which thunderstorms frequently occur. The analysis region covered was (39.3–40.3°N, 116–118°E), which is located in the core detection zone. First, the VHF radiation sources and CG return strokes were separately grouped into flashes [51]. The SAFIR_CG data represent the CG flashes, and the SAFIR_Tot data represent the total flashes.


## 3.3. A Fitting Algorithm

We aimed to design a simple algorithm to obtain the coverage of the total flash by fitting the CG flash data. The fitted lightning flash occurrence region is a better representation of the actual total lightning occurrence region. Although a few of the thunderstorms contained only intracloud (IC) flashes or only CG flashes, most of the thunderstorms contained both IC flashes and CG flashes, and the occurrence of IC flashes and CG flashes is closely linked. Lightning will occur when there is an accumulation of charge in clouds and IC flashes often develop before CG flashes [53]. Whether the lightning develops into an IC flash or a CG flash depends on the charge structure in the clouds, although it may involve some stochastic processes. The proportion of CG flashes to total flashes in thunderstorms is about 20–50% with regional variations [52, 54, 55]. Total flashes occur over larger areas than CG flashes.When a CG flash occurs at a statistical point (e.g., the red point in Figure3) on the timescale of a thunderstorm (cumulative over several hours), we assume that some IC flashes may occur (blue points in Figure 3(a)) around the CG flash point. The blue points are defined as preliminary guesses for IC flashes. When a grid point is guessed at least twice, then it is assumed that an IC flash occurs in this grid (purple points in Figure 3(b)). A new dataset can be generated by processing (fitting) the CG flash data. The fitting algorithm should first be evaluated with a total flash locating system.Figure 3
Sketch map of the simple algorithm for fitting total flash.
(a)(b)A new fitting flash dataset referred to as SAFIR_CG_fitted can be obtained by processing the SAFIR_CG data. No specific fitting analysis is performed when the number of grids with lightning flash occurrences is ≤2 during a 6-hour period.Figure4 shows the fitting results for a 6-hour period (06 : 00–12 : 00 UTC) on 30 June 2008. The area of occurrence of the SAFIR_CG data is smaller than that of SAFIR_Tot, although the occurrence area of the SAFIR_CG data is often located in the kernel areas of the SAFIR_Tot data. The algorithm expands the CG flash to match the relationship between the CG flash and the total area of flash occurrence. Compared with the SAFIR_CG data, the area of occurrence of the SAFIR_CG_fitted data is more consistent with that of the SAFIR_Tot data. Therefore, from qualitative analysis, the fitting algorithm is effective without introducing too many false signals. The relative detection efficiency (RDE) was defined to measure the performance of the CG datasets.(2)RDE=NA∩NBNB,where NA represents the flash area covered by dataset A (CG or fitted CG) and NB represents the flash area covered by dataset B (total flash). The RDE can give the fraction of all flashes detected by B that were also detected by A. Some false lightning signals will be introduced using CG flash data to fit the total flash data. Thus, the concept of the false signal ratio (FSR) was introduced:(3)FSR=NA−NA∩NBNB.Figure 4
The SAFIR_Tot, SAFIR_CG, and SAFIR_CG_fitted data accumulated from 06 : 00 to 12 : 00 UTC on 30 June 2008.
(a)(b)(c)The fitting algorithm was evaluated using the SAFIR dataset from June to September in 2007–2008 with a 6 h interval. Compared with the original SAFIR_CG data, the SAFIR_CG_fitted data showed a significant improvement in the RDE. The average RDE of the SAFIR_CG data concerning SAFIR_Tot was 35%. The fitting algorithm can increase the RDE by an average of 20% after introducing a small FSR (about 10%). When the number of grids with lightning occurrence exceeded 25, the RDE increased by about 26–30% with a smaller FSR. The fitting algorithm can effectively bridge the spatial gap between the CG flash and the total flash.Because the ADTD system provides only CG flash data, the fitting algorithm generates the total flash from the CG flash. In the next part of forecast verification, CG_original represents the CG flash detected by ADTD, and CG_fitted represents data generated from the CG_original by the fitting algorithm.


## 3.4. Verification Method

It is difficult to assess the actual quality of the high-resolution model using traditional metrics (point-to-point method) because slight displacement errors often result in “double penalties” (i.e., errors observed but not forecast, or forecast but not observed). To avoid these double penalties, the neighborhood-based equitable threat score (ETS) reported by Clark et al. [56] was adopted in this study. This is an approach to validating forecasts on a high-resolution numerical grid. The verification is based on the gridded observed lightning flash data and model forecast flash data.The traditional method of computing the ETS [57] uses a 2 × 2 contingency table of possible forecast outcomes at individual grid points where the table elements are hits (the correct forecast of an event), misses (an observed but not forecast event), false alarms (a forecast but not observed event), and correct negatives (correct forecast of nonevent). Using these elements, the ETS is expressed as(4)ETS=hits−chancehits+misses+false alarms−chance,where(5)chance=hits+misseshits+false alarmshits+misses+correctnegative+false alarms.The average ETSs were computed by summing (i.e., aggregating) contingency table elements over all cases. The miss rate (MR) and false alarm ratio (FAR) were also calculated as follows:(6)MR=misseshits+misses,FAR=false alarmshits+false alarms.To compute a neighborhood-based ETS, the criteria for hits are relaxed by considering adjacent grid points within a specified radius of each grid point. Following Clark et al. [56], the ETS was calculated for lightning forecasts at radii of 0, 5, 10, 20, 30, 40, 60, and 90 km for thresholds referring to 1 flash/16 km2 (lightning flash density), which means that only the region in which lightning occurs is quantitatively evaluated (a yes or no evaluation). When the neighborhood radius is set to 0, the neighborhood-based approach reduces to the traditional form of ETS.Fraction skill scores (FSSs) [58] were also used in this study. FSS is a neighborhood-based verification method in which the fractional coverage of predicted and observed grid area to the total area is compared. The FSS can be more tolerant of small displacement errors than the ETS, making it more suitable for evaluating a simulation with fine-scale grids [31]. In the actual calculation, the threshold q needs to be set first to preprocess the data, and q is used to convert the observation field (Or) and forecast field (Mr) into binary fields Io and IM. All grid squares exceeding the threshold have a value of 1 and all others have a value of 0:(7)Io=1Or≥q0Or<q,IM=1Mr≥q0Mr<q.Io and IM are the observed and predicted fields processed in the space window, respectively. For a certain neighborhood radius, Io and IM can be calculated by selecting q, and FSS can be further calculated. In this study, FSS was calculated for lightning forecasts at radii of 4, 12, 20, 28, 36, 44, 52, 60, 68, 76, 84, 92, and 100 km for thresholds referring to 1 flash/16 km2 (flash density). The CG_fitted flash data of ADTD are used in the verification. Because the observed CG flashes (CG_original) have been counted in the 4 km × 4 km grid, according to the number of grids in which CG flashes have been observed (N) in the 6-hour period, all the inspection periods are divided into seven categories (N < 50, 50 ≤ N < 100, 100 ≤ N < 500, 500 ≤ N < 1000, 1000 ≤ N < 2000, 2000 ≤ N < 3000, and N ≥ 3000). All periods are evaluated for classification by the FSS method.Our study focused on the forecasting ability of the numerical model in a 6-hour period. Based on the differences in forecasting time, the 6–12-hour, 12–18-hour, and 18–24-hour lightning forecasts were examined. The initial 0–6-hour forecasts were discarded from most analyses and the model was spinning up during this period, so the lightning forecasts may be inaccurate. Meanwhile, the model is initialized at four different times a day so that there are four forecasting results for the same 6 h observation time with a 6 h interval.


## 4. Results and Discussion

### 4.1. Quantitative Verification of the Forecast Results

To investigate the performance of the WRF-Electric model in forecasting lightning, the results of continuous prediction experiments for three rainy seasons from May to September in 2015–2017 (only May in 2015 was missing) were evaluated using the neighborhood-based ETS and FSS approaches. The ETS, FSS, MR, and FAR for lightning were calculated for the 6–12-hour (PRE12), 12–18-hour (PRE18), and 18–24-hour (PRE24) forecast periods for different neighborhood radii. The forecasts were compared with the 6-hourly observed values of the CG_fitted and CG_original data of the ADTD system.Figure5 shows that the ETS increased rapidly as the radius increased from 5 to 30 km and that the MR and FAR also decreased sharply over the same range. The variation got smaller when the neighborhood radius was >30 km. For almost all neighborhood radii, the ETS score was the highest in PRE12 and the lowest in PRE24. The corresponding MR and FAR of PRE12 were the lowest among the three analysis periods. PRE24 had the highest MR. There was no apparent difference in the FAR between PRE18 and PRE24. Based on a comparison of the skill scores of the three periods, the model predictions were the best during PRE12.Figure 5
The ETS, MR, and FAR values for different model forecasting periods (PRE12, 6–12-hour period; PRE18, 12–18-hour period; and PRE24, 18–24-hour period) for different neighborhood radii. The solid lines represent the results compared with the CG_fitted data, and the dashed lines represent the results compared with the CG_original data.
(a)(b)(c)The ETS value of PRE12 for a neighborhood radius of 20 km was 0.34, with corresponding MR and FAR values of 0.38 and 0.55, respectively. The most significant problem in the forecast is the high MR values. For a lightning flash density forecast skill score at a threshold of 1 flash/16 km2, the high MR means that the predicted region in which lightning occurred was smaller than the observed region. We also compared the 6-hourly observed values of the original cloud-to-ground lightning data (CG_original) with PRE12, and the results are shown as dashed lines (PRE12_original) in Figure 5. The MR of PRE12 was higher than that of PRE12_original because the coverage area of the CG_original lightning was smaller than that of CG_fitted.The high MR is probably the overly high electric field threshold of lightning initiation (Ebreak) in the discharge scheme, which means that discharge cannot occur at the region where the electrification process is weak. A total of 28 uneven vertical levels and 4 km horizontal grids were used for the inner mesh in the model, which is a relatively low spatial resolution for the simulation of lightning discharge. The strong electric field in clouds could be smoothed due to the relatively coarse spatial resolution, with the result that the model cannot reflect the discharge activity of a weak electrified area of a thunderstorm. Although the β coefficient in (1) was set to 0.5 based on the model resolution, Ebreak was still high. The simulated thunderstorms are also sensitive to horizontal resolution [59]. Therefore, according to the model resolution, Ebreak should be modified to obtain the correct coverage of model prediction for lightning activity. The discharge scheme needs to be improved to predict more lightning in a weak electrification area of a thunderstorm, thereby reducing the MR.In North China, the predictability of lightning could be dependent on typical weather systems. The lightning forecasts from May to September for different forecast periods (PRE12, PRE18, and PRE24) were compared based on the neighborhood ETS approach to investigate the monthly variations in the model performance. The results showed that the forecasting skill of PRE12 was the best among the three forecast periods for all months (not shown). The monthly ETS scores were therefore compared using the PRE12 forecast results.Figure6 shows the ETS values of PRE12 for each month. In June, July, and August, the ETS scores are higher than those in May and September. From May to September, the average number of grids in which CG flashes were observed in the 6 hours was 81, 265, 278, 298, and 233, respectively. The systematic reduction of ETS in May and September may be due to the variation of thunderstorm scale. When the scale of thunderstorm becomes smaller, the predictability becomes worse, which is also reflected in FSS verification. The main rainy season in North China runs from June to August, and thunderstorms are relatively weak and less frequent in May and September. If Ebreak is higher than in the real situation, the model could lead to high MR values for weak thunderstorms, explaining the low forecasting skill in May and September. The predictability of lightning may be affected by convective intensity.Figure 6
ETS values for different neighborhood radii for different months in 6–12-hour forecasting periods (PRE12).Thunderstorms are characterized not only by monthly variations but also by diurnal variations. Therefore, the forecast results were compared to determine the local time period in which the model showed the best forecasting skill. Because every forecast contains 24-hour forecasting results, there are four forecast results for the 6-hour local standard time (LST = UTC + 08 : 00) periods selected (02 : 00–08 : 00, 08 : 00–14 : 00, 14 : 00–20 : 00, and 20 : 00–02 : 00). Figure7 shows the ETS values in a neighborhood radius of 20 km for those periods. Regardless of the forecasting time, the model shows the best forecasting skill for local time 14 : 00–20 : 00. The ETS value of PRE12 for local time 14 : 00–20 : 00 was the highest (>0.4).Figure 7
ETS values for a neighborhood radius of 20 km for different local times in different forecasting periods (PRE6, 0–6-hour period; PRE12, 6–12-hour period; PRE18, 12–18-hour period; PRE24, 18–24-hour period).In fact, except for the difference in monthly climatic background, the characteristics in terms of thunderstorm type will have an impact on the predictability of thunderstorms. Figure8 shows the FSS for neighborhood radii and numbers of grids in which lightning is observed, classified into seven categories in PRE12. The prediction FSS value is lower for a period with a smaller number of grids in which there were flashes (N < 500), and it can reach 0.2 for a neighborhood radius of 100 km. However, for a period with more grids in which there were flashes (2000 ≤ N < 3000 and N ≥ 3000), the FSS value is high, and the score can reach 0.5 at the neighborhood radii of 28 km and 36 km. The FSS value increases with neighborhood radii and the number of grids in which lightning is observed.Figure 8
The fraction skill scores (FSSs) for different neighborhood radii for different categories in 6–12-hour forecasting periods (PRE12). The forecasting periods are divided into seven categories by using the number of grids in which lightning is observed (N). Curves of different colors correspond to different categories (N < 50, 50 ≤ N < 100, 100 ≤ N < 500, 500 ≤ N < 1000, 1000 ≤ N < 2000, 2000 ≤ N < 3000, and N ≥ 3000).Although the number of grids in which lightning is observed is not counted based on individual thunderstorms, it can still reflect the thunderstorm scale to some extent. When the number of grids in which lightning is observed is small, the thunderstorm tends to be more localized within the 6 hours, and when the number is large, the thunderstorm tends to be more widespread. The verification results show that the predictability of the model also increases with the thunderstorm scale.


### 4.2. Qualitative Analysis of the Forecast Results

One goal of this section is to show how the lightning forecasts vary for forecast periods qualitatively. The model’s ability to forecast the density of lightning flashes was analyzed. For three typical cases, Figure9 shows the three cases of observed and predicted lightning flash density in 2015: 17 May (Figures 9(a)–9(e)), 10 June (Figures 9(f)–9(j)), and 15 July (Figures 9(k)–9(o)). The synoptic-scale circulations (figure omitted) of the three episodes were analyzed at 850 hPa, 500 hPa, and 200 hPa.Figure 9
Observations and experiments at different prediction times: ((a)–(e)) case on 17 May, ((f)–(j)) case on 10 June, and ((k)–(o)) case on 15 July. The first and second columns show the observed cloud-to-ground lightning (CG_original) and fitted lightning (CG_fitted) data ((a)–(b) 06 : 00–12 : 00 on 17 May, (f)–(g) 06 : 00–12 : 00 on 10 June, and (k)–(l) 18 : 00–24 : 00 on 15 July 2015). The third to fifth columns represent the predicted results at different forecasting periods: the 6–12-hour, 12–18-hour, and 18–24-hour forecast results.The first episode was associated with the low-level vortex and midlevel short-wave trough eastward. Three bands appeared in the lightning distribution and the choice of this period better reflects the ability of the model to predict the location of the occurrence of lightning. The second episode was induced by the deep vortex system at low and middle levels. The coverage of the lightning activity was widespread and the observed lightning flash density was relatively high. The last episode was caused by the vortex at middle and high levels consistent with easterly flow in low level. The third episode was selected because the region where lightning occurred was relatively localized and the lightning activity was weak.The coverage of the CG_fitted data was more continuous than that of the CG_original data; meanwhile, the coverage of the observed lightning activity (CG_fitted data) was larger than that of the predicted results. Besides the main distribution area of lightning activity, other instances of observed lightning were more scattered. The predicted lightning activity was relatively concentrated and was only representative of an area of high flash density. The model successfully captured the main distribution regions of lightning activity in these three cases (Figures9(c), 9(h), and 9(n)), while the scattered lightning flashes were hardly predicted. Under these conditions, the model was prone to misprediction in areas where lightning activity was weak.In all three cases, the model’s best predictive results were in PRE12 or PRE18, although the results showed some differences for forecasting periods. On 17 May and 10 June 2015, the distribution of lightning activity in PRE12 was consistent with the observations (Figures9(c) and 9(h)). For 17 May 2015, the model predicted the same three belts of lightning activity as the observed actions. On 15 July 2015, the best forecast result was during PRE18 (Figure 9(n)). In PRE12 (Figure 9(m)), the model did not predict the basic pattern of regions with observed lightning.The predicted flash density was higher than the observations, particularly on 10 June 2015 (Figures9(g)–9(h)), when the lightning activity was stronger and broader. The observed flash density was mainly <10 flashes/16 km2, whereas the predicted density reached 100 flashes/16 km2. The IC and CG lightning flashes were not distinguished in the discharge scheme. The ratio of CG lightning to total lightning is about 20% in North China [52], and therefore, the predicted flash density was unreasonably high. In the model discharge scheme, the magnitude of charge density is proportionally reduced arbitrarily, which may distort the amount of neutralization charge in a single lightning flash. The current discharge scheme is likely to have a problem in that the amount of neutralization charge in a single lightning strike is too small, resulting in a higher density of lightning flashes in the forecast. More observations of lightning discharge should be used in the discharge scheme to maintain a reasonable amount of neutral lightning charge in the simulation.The model successfully predicted the region of lightning activity. However, the lightning parameterization scheme needs to be improved to predict more lightning flashes in the weak electrification areas of a thunderstorm and reduce the density of lightning flashes in the strong electrification areas.In the previous studies [25, 29, 32], the lightning forecasting experiments with an electrical model were mainly carried out in one or several thunderstorm cases. This study performs long-term (three rainy seasons) lightning forecasting experiments with an electrical numerical prediction model and aim to verify it with lightning observations, adjusted to include contributions from intracloud lightning. However, the forecast skill score of this study is relatively lower.The bulk discharge scheme in this study or E-WRF [28] does not distinguish between lightning types (i.e., CG or IC lightning), which is a shortcoming of the electrical WRF model. More consideration should be given to how to divide the types of lightning. When predicting lightning using a mesoscale numerical model, it is difficult to determine the type of lightning in a bulk lightning scheme. Lynn et al. [60] proposed a dynamic algorithm to produce forecasting maps for cloud-to-ground and intracloud lightning. The algorithm used the dynamic and microphysics fields from the WRF cloud-resolving model to calculate the electrical potential energy for different types of lightning. Tan et al. [42] also indicated that the electric potential at the initial discharge time affects the type of lightning, which will help discriminate the type of lightning in a bulk parameterization discharge. Furthermore, positive CG and negative CG must be distinguished in the discharge scheme [61], even in a bulk scheme.


## 4.1. Quantitative Verification of the Forecast Results

To investigate the performance of the WRF-Electric model in forecasting lightning, the results of continuous prediction experiments for three rainy seasons from May to September in 2015–2017 (only May in 2015 was missing) were evaluated using the neighborhood-based ETS and FSS approaches. The ETS, FSS, MR, and FAR for lightning were calculated for the 6–12-hour (PRE12), 12–18-hour (PRE18), and 18–24-hour (PRE24) forecast periods for different neighborhood radii. The forecasts were compared with the 6-hourly observed values of the CG_fitted and CG_original data of the ADTD system.Figure5 shows that the ETS increased rapidly as the radius increased from 5 to 30 km and that the MR and FAR also decreased sharply over the same range. The variation got smaller when the neighborhood radius was >30 km. For almost all neighborhood radii, the ETS score was the highest in PRE12 and the lowest in PRE24. The corresponding MR and FAR of PRE12 were the lowest among the three analysis periods. PRE24 had the highest MR. There was no apparent difference in the FAR between PRE18 and PRE24. Based on a comparison of the skill scores of the three periods, the model predictions were the best during PRE12.Figure 5
The ETS, MR, and FAR values for different model forecasting periods (PRE12, 6–12-hour period; PRE18, 12–18-hour period; and PRE24, 18–24-hour period) for different neighborhood radii. The solid lines represent the results compared with the CG_fitted data, and the dashed lines represent the results compared with the CG_original data.
(a)(b)(c)The ETS value of PRE12 for a neighborhood radius of 20 km was 0.34, with corresponding MR and FAR values of 0.38 and 0.55, respectively. The most significant problem in the forecast is the high MR values. For a lightning flash density forecast skill score at a threshold of 1 flash/16 km2, the high MR means that the predicted region in which lightning occurred was smaller than the observed region. We also compared the 6-hourly observed values of the original cloud-to-ground lightning data (CG_original) with PRE12, and the results are shown as dashed lines (PRE12_original) in Figure 5. The MR of PRE12 was higher than that of PRE12_original because the coverage area of the CG_original lightning was smaller than that of CG_fitted.The high MR is probably the overly high electric field threshold of lightning initiation (Ebreak) in the discharge scheme, which means that discharge cannot occur at the region where the electrification process is weak. A total of 28 uneven vertical levels and 4 km horizontal grids were used for the inner mesh in the model, which is a relatively low spatial resolution for the simulation of lightning discharge. The strong electric field in clouds could be smoothed due to the relatively coarse spatial resolution, with the result that the model cannot reflect the discharge activity of a weak electrified area of a thunderstorm. Although the β coefficient in (1) was set to 0.5 based on the model resolution, Ebreak was still high. The simulated thunderstorms are also sensitive to horizontal resolution [59]. Therefore, according to the model resolution, Ebreak should be modified to obtain the correct coverage of model prediction for lightning activity. The discharge scheme needs to be improved to predict more lightning in a weak electrification area of a thunderstorm, thereby reducing the MR.In North China, the predictability of lightning could be dependent on typical weather systems. The lightning forecasts from May to September for different forecast periods (PRE12, PRE18, and PRE24) were compared based on the neighborhood ETS approach to investigate the monthly variations in the model performance. The results showed that the forecasting skill of PRE12 was the best among the three forecast periods for all months (not shown). The monthly ETS scores were therefore compared using the PRE12 forecast results.Figure6 shows the ETS values of PRE12 for each month. In June, July, and August, the ETS scores are higher than those in May and September. From May to September, the average number of grids in which CG flashes were observed in the 6 hours was 81, 265, 278, 298, and 233, respectively. The systematic reduction of ETS in May and September may be due to the variation of thunderstorm scale. When the scale of thunderstorm becomes smaller, the predictability becomes worse, which is also reflected in FSS verification. The main rainy season in North China runs from June to August, and thunderstorms are relatively weak and less frequent in May and September. If Ebreak is higher than in the real situation, the model could lead to high MR values for weak thunderstorms, explaining the low forecasting skill in May and September. The predictability of lightning may be affected by convective intensity.Figure 6
ETS values for different neighborhood radii for different months in 6–12-hour forecasting periods (PRE12).Thunderstorms are characterized not only by monthly variations but also by diurnal variations. Therefore, the forecast results were compared to determine the local time period in which the model showed the best forecasting skill. Because every forecast contains 24-hour forecasting results, there are four forecast results for the 6-hour local standard time (LST = UTC + 08 : 00) periods selected (02 : 00–08 : 00, 08 : 00–14 : 00, 14 : 00–20 : 00, and 20 : 00–02 : 00). Figure7 shows the ETS values in a neighborhood radius of 20 km for those periods. Regardless of the forecasting time, the model shows the best forecasting skill for local time 14 : 00–20 : 00. The ETS value of PRE12 for local time 14 : 00–20 : 00 was the highest (>0.4).Figure 7
ETS values for a neighborhood radius of 20 km for different local times in different forecasting periods (PRE6, 0–6-hour period; PRE12, 6–12-hour period; PRE18, 12–18-hour period; PRE24, 18–24-hour period).In fact, except for the difference in monthly climatic background, the characteristics in terms of thunderstorm type will have an impact on the predictability of thunderstorms. Figure8 shows the FSS for neighborhood radii and numbers of grids in which lightning is observed, classified into seven categories in PRE12. The prediction FSS value is lower for a period with a smaller number of grids in which there were flashes (N < 500), and it can reach 0.2 for a neighborhood radius of 100 km. However, for a period with more grids in which there were flashes (2000 ≤ N < 3000 and N ≥ 3000), the FSS value is high, and the score can reach 0.5 at the neighborhood radii of 28 km and 36 km. The FSS value increases with neighborhood radii and the number of grids in which lightning is observed.Figure 8
The fraction skill scores (FSSs) for different neighborhood radii for different categories in 6–12-hour forecasting periods (PRE12). The forecasting periods are divided into seven categories by using the number of grids in which lightning is observed (N). Curves of different colors correspond to different categories (N < 50, 50 ≤ N < 100, 100 ≤ N < 500, 500 ≤ N < 1000, 1000 ≤ N < 2000, 2000 ≤ N < 3000, and N ≥ 3000).Although the number of grids in which lightning is observed is not counted based on individual thunderstorms, it can still reflect the thunderstorm scale to some extent. When the number of grids in which lightning is observed is small, the thunderstorm tends to be more localized within the 6 hours, and when the number is large, the thunderstorm tends to be more widespread. The verification results show that the predictability of the model also increases with the thunderstorm scale.


## 4.2. Qualitative Analysis of the Forecast Results

One goal of this section is to show how the lightning forecasts vary for forecast periods qualitatively. The model’s ability to forecast the density of lightning flashes was analyzed. For three typical cases, Figure9 shows the three cases of observed and predicted lightning flash density in 2015: 17 May (Figures 9(a)–9(e)), 10 June (Figures 9(f)–9(j)), and 15 July (Figures 9(k)–9(o)). The synoptic-scale circulations (figure omitted) of the three episodes were analyzed at 850 hPa, 500 hPa, and 200 hPa.Figure 9
Observations and experiments at different prediction times: ((a)–(e)) case on 17 May, ((f)–(j)) case on 10 June, and ((k)–(o)) case on 15 July. The first and second columns show the observed cloud-to-ground lightning (CG_original) and fitted lightning (CG_fitted) data ((a)–(b) 06 : 00–12 : 00 on 17 May, (f)–(g) 06 : 00–12 : 00 on 10 June, and (k)–(l) 18 : 00–24 : 00 on 15 July 2015). The third to fifth columns represent the predicted results at different forecasting periods: the 6–12-hour, 12–18-hour, and 18–24-hour forecast results.The first episode was associated with the low-level vortex and midlevel short-wave trough eastward. Three bands appeared in the lightning distribution and the choice of this period better reflects the ability of the model to predict the location of the occurrence of lightning. The second episode was induced by the deep vortex system at low and middle levels. The coverage of the lightning activity was widespread and the observed lightning flash density was relatively high. The last episode was caused by the vortex at middle and high levels consistent with easterly flow in low level. The third episode was selected because the region where lightning occurred was relatively localized and the lightning activity was weak.The coverage of the CG_fitted data was more continuous than that of the CG_original data; meanwhile, the coverage of the observed lightning activity (CG_fitted data) was larger than that of the predicted results. Besides the main distribution area of lightning activity, other instances of observed lightning were more scattered. The predicted lightning activity was relatively concentrated and was only representative of an area of high flash density. The model successfully captured the main distribution regions of lightning activity in these three cases (Figures9(c), 9(h), and 9(n)), while the scattered lightning flashes were hardly predicted. Under these conditions, the model was prone to misprediction in areas where lightning activity was weak.In all three cases, the model’s best predictive results were in PRE12 or PRE18, although the results showed some differences for forecasting periods. On 17 May and 10 June 2015, the distribution of lightning activity in PRE12 was consistent with the observations (Figures9(c) and 9(h)). For 17 May 2015, the model predicted the same three belts of lightning activity as the observed actions. On 15 July 2015, the best forecast result was during PRE18 (Figure 9(n)). In PRE12 (Figure 9(m)), the model did not predict the basic pattern of regions with observed lightning.The predicted flash density was higher than the observations, particularly on 10 June 2015 (Figures9(g)–9(h)), when the lightning activity was stronger and broader. The observed flash density was mainly <10 flashes/16 km2, whereas the predicted density reached 100 flashes/16 km2. The IC and CG lightning flashes were not distinguished in the discharge scheme. The ratio of CG lightning to total lightning is about 20% in North China [52], and therefore, the predicted flash density was unreasonably high. In the model discharge scheme, the magnitude of charge density is proportionally reduced arbitrarily, which may distort the amount of neutralization charge in a single lightning flash. The current discharge scheme is likely to have a problem in that the amount of neutralization charge in a single lightning strike is too small, resulting in a higher density of lightning flashes in the forecast. More observations of lightning discharge should be used in the discharge scheme to maintain a reasonable amount of neutral lightning charge in the simulation.The model successfully predicted the region of lightning activity. However, the lightning parameterization scheme needs to be improved to predict more lightning flashes in the weak electrification areas of a thunderstorm and reduce the density of lightning flashes in the strong electrification areas.In the previous studies [25, 29, 32], the lightning forecasting experiments with an electrical model were mainly carried out in one or several thunderstorm cases. This study performs long-term (three rainy seasons) lightning forecasting experiments with an electrical numerical prediction model and aim to verify it with lightning observations, adjusted to include contributions from intracloud lightning. However, the forecast skill score of this study is relatively lower.The bulk discharge scheme in this study or E-WRF [28] does not distinguish between lightning types (i.e., CG or IC lightning), which is a shortcoming of the electrical WRF model. More consideration should be given to how to divide the types of lightning. When predicting lightning using a mesoscale numerical model, it is difficult to determine the type of lightning in a bulk lightning scheme. Lynn et al. [60] proposed a dynamic algorithm to produce forecasting maps for cloud-to-ground and intracloud lightning. The algorithm used the dynamic and microphysics fields from the WRF cloud-resolving model to calculate the electrical potential energy for different types of lightning. Tan et al. [42] also indicated that the electric potential at the initial discharge time affects the type of lightning, which will help discriminate the type of lightning in a bulk parameterization discharge. Furthermore, positive CG and negative CG must be distinguished in the discharge scheme [61], even in a bulk scheme.


## 5. Summary and Conclusions

Experiments to predict the regional lightning activity from 2015 to 2017 using the WRF-Electric model coupled with electrification and discharge schemes were carried out in this study. By employing the neighborhood-based ETS and FSS verification method, the ability to predict regional lightning activity in the numerical model was evaluated objectively. The model’s main problems were identified through the verification, which provides a basis for its further improvement. The main conclusions are as follows:(1)
The major region of lightning activity can be predicted well by the mesoscale electrification and discharge model. The quantitative verification results over North China showed that the predicted results were best in the 6–12-hour forecasting period. In the neighborhood-based ETS verification, the average ETS score in the 6–12-hour forecast period was 0.34 for a 20 km neighborhood radius. The prediction skill of the model varies not only monthly, but also diurnally. The model shows better forecasting skill in the main rainy season (June–July–August) and at 14 : 00–20 : 00 local time. The model results are more valuable for predicting lightning flashes for the 6–18-hour forecast period. The FSS verification results show that the predictability of the model increases with an increase in the thunderstorm scale. For the periods with a greater number of grids in which there were flashes (2000 ≤N < 3000 and N ≥ 3000), the FSS score can reach 0.5 at neighborhood radii of 28 km and 36 km. The forecast results do not perform well based on the FSS verification when the thunderstorm scale is small. The predictability of lightning may be also affected by convective intensity.(2)
The coverage of the lightning activity predicted by the model was small and relatively concentrated compared with observation. Some areas with scattered lightning activity were missed. The threshold of the initial lightning should be modified according to the model resolution. For this study, the lightning parameterization scheme needs to be improved to predict more lightning in thunderstorms with weak electrification to reduce the MR. The density of lightning predicted by the model was higher than the observed actions, which may be caused by the unreasonable amount of charge neutralized by a single simulated lightning flash. To improve the prediction of the flash density, the magnitude of the neutralization charge of single lightning in the discharge scheme should be determined from the observational results.Currently, the ability to quantitatively predict the lightning flash density by mesoscale electrification and discharge models lags far behind. Although the distribution of predicted lightning activity also depended on the different noninductive charging schemes, the main discrepancies in the current models are in the design of discharge parameterization. Parameterizing lightning reasonably in a mesoscale model is still an unresolved and highly challenging problem [62].The neighborhood-based ETS and FSS methods are used to evaluate lightning activity forecasts. These verification methods are widely used. However, they were originally developed for precipitation forecasting. Due to the unique characteristics of lightning activity, a recent study [33] has proposed a new method for verifying lightning forecasts in terms of three properties: coverage of the domain, distance to the observations, and intensity of the lightning forecast. The new method should be given more attention in future researches.Compared with the evaluation conducted by Dafis et al. [29] and Gharaylou et al. [63], the resolution of the initial field and simulation used in this study is low, but the current results are of reference value for improving the discharge scheme. If the initial field is improved and the resolution is increased, it is believed that the ability to predict the lightning in WRF-Electric will be significantly improved.



---

*Source: 4583030-2022-02-15.xml*