Potential Impacts of Lightning Flash Observations on Numerical Weather Prediction With Explicit Lightning Processes

Although various lightning observation platforms exist, lightning flash observations have not been fully utilized for numerical weather prediction (NWP), mainly due to complications related to modeling lightning processes. Herein, we investigated the potential impacts of lightning observations using an ensemble Kalman filter (EnKF) with an NWP model, including explicit lightning processes. The EnKF included flow‐dependent error covariances between lightning flashes and the atmospheric state. Observing system simulation experiments were performed for an idealized convective case. Ensemble‐based correlations with various observations indicated that lightning flashes provide valuable information below cloud tops where satellite infrared radiances cannot observe well due to the presence of anvil clouds. Moreover, a single‐step data assimilation experiment showed that lightning flashes positively impact the analysis and forecast of wind, hydrometeors, and electric charge. The results revealed that lightning flashes are potentially useful for an EnKF with an NWP model including explicit lightning processes.

Lightning observations are a potentially valuable tool for further improving NWP, as they are closely associated with solid particles within convective clouds (Takahashi, 1978). Lightning flashes are observed by a ground-based sensor network, such as the National Lightning Detection Network (NLDN, Cummins & Murphy, 2009) in the US, the Broadband Observation network for Lightning and Thunderstorm (BOLT, Yoshida et al., 2014) in Japan, and the World Wide Lightning Location Network (WWLLN, Lay et al., 2004), as well as visible sensors such as the Geostationary Lightning Mapper (GLM, Goodman et al., 2013) onboard . In particular, the GLM covers a wide area of land and ocean at an 8-km horizontal resolution and provides continuous two-dimensional (2D) locations of lightning flashes. As noted by Mansell (2014), lightning observations are not as useful as ground-based radar observations over land. However, they can be used to obtain crucial information over the ocean, where ground-based radars are unavailable.
Several studies have examined the assimilation of lightning flashes. Fierro et al. (2016Fierro et al. ( , 2019 and Hu et al. (2020) assimilated pseudo-water vapor observations derived from the GLM using a 3D variational (3DVAR) system and demonstrated beneficial impacts on short-term forecasts. Wang et al. (2017) developed a nudging-based assimilation method for retrieved graupel and associated latent heat releases and demonstrated that their method successfully improves short-term precipitation forecasts. Allen et al. (2016) used an ensemble Kalman filter (EnKF, Evensen, 1994;Houtekamer & Zhang, 2016) with statistical relationships between model variables and lightning flashes derived from offline simulations. They used an NWP model with the electrification method by Mansell et al. (2005) to derive statistical relationships; however, they turned off the electrification method for their data assimilation (DA) experiments. Similarly, Wang et al. (2018) successfully assimilated pseudo-observations of graupel with an EnKF system. More recently, using Allen et al. (2016)'s statistical relationships, Kong et al. (2020) successfully assimilated real-world GLM observations with an EnKF system. With these approaches, observation operators that convert model variables such as graupel volume into a lightning observation would need to be optimized for each event. Alternatively, the use of an NWP model, including the explicit electric processes, would allow the EnKF to automatically construct covariance among the model variables and lightning flashes. Mansell (2014) tested this approach and showed that lightning flashes improve the distributions of analyzed updraft and radar reflectivity in an isolated convective cloud. He also showed that spurious convection was successfully suppressed by assimilating zero flashes. However, he did not investigate the details of the covariances that can inform what information is obtained by assimilating lightning flashes.
Although previous studies have demonstrated positive results regarding lightning flash assimilations, our knowledge on the topic remains limited. In particular, it is unclear what information can be obtained from lightning observations on top of satellite IR radiances (e.g., . Previous studies have shown beneficial impacts of all-sky IR radiance DA on NWP, yet lightning flash DA might further improve NWP. Here, we aimed to investigate the potential of lightning observations with an EnKF system by performing observing system simulation experiments (OSSEs) for an idealized supercell case. This study used an NWP model with explicit lightning processes and investigated the ensemble-based correlations between model variables and lightning flashes, compared with those between model variables and satellite IR radiances. To the best of the authors' knowledge, such a comparison has not been shown thus far. As previous studies have pointed out (e.g., Mansell, 2014), lightning observations are not expected to outperform weather radar observations. Therefore, we focused on the comparison between lightning and IR radiance observations. First, we performed a nature run and simulated the observations. Furthermore, we performed OSSEs using an EnKF system with simulated radar observations and ran an ensemble forecast. We then compared the ensemble-based correlations between the model variables and simulated observations, including lightning flashes and satellite IR radiances. Finally, to prove the concept, we assimilated GLM observations, corresponding to the flash number in a model at a single step without cycling, and demonstrated their impacts on NWP for a single case. In contrast to Allen et al. (2016), we turned on the electrification process in all simulations. As this study aimed to present an essential first step to assess the potential impacts of lightning flash observations on top of all-sky satellite IR radiances, we chose to use perfect-model OSSEs instead of real GLM data. By doing so, correlation analyses could be designed to be straightforward comparisons between lightning flash and IR radiance observations because the model and observation biases for the real GLM data have not been assessed. Furthermore, to effectively assimilate lightning flashes in cycling experiments with an EnKF system, a method to initiate lightning flashes must be implemented, even if all ensemble members have zero flashes. The development of such a method is beyond the scope of this initial study and will be addressed in future research.
The rest of this study is organized as follows. Section 2 provides the experimental settings. Section 3 describes the nature run and shows the assimilation impacts of radar observations on electric variables. Section 4 presents an ensemble-based correlation analysis with various observation types. The impact of assimilating flash observations is presented in Section 5. Finally, Section 6 provides a summary.

SCALE-LETKF
In this study, we used the SCALE-LETKF system (Lien et al., 2017), consisting of an atmospheric Regional Model from the Scalable Computing for Advanced Library and Environment (SCALE-RM, Nishizawa et al., 2015;Sato et al., 2015) and the local ensemble transform Kalman filter (LETKF, Hunt et al., 2007;Miyoshi & Yamane, 2007). The SCALE-RM solves the set of fully compressible nonhydrostatic equations with prognostic variables of 3D momentums, total density, mass-weighted potential temperature, and mass concentration of tracers such as water vapor and hydrometeors. The SCALE-RM was set up with a horizontally 2-km and vertically 500-m mesh domain that covers 352 × 352 × 22.25 km. We used the microphysics scheme of Seiki and Nakajima (2014), with 10 prognostic hydrometeor variables comprising mixing ratios and number densities for five hydrometeor categories of cloud water, rain, cloud ice, snow, and graupel.
The SCALE-RM was recently extended by Sato et al. (2019) to include electric processes such as noninductive charge separation (Takahashi, 1978) and neutralization based on Fierro et al. (2013)'s parameterization. The prognostic variables for the electric processes are the charge density of each hydrometeor category, and the electric potential and electric field are diagnosed from charge density by solving the Poisson equation. Following Fierro et al. (2013), the discharge process occurs in a vertical cylinder (6-km radius) that centers the grid point in which the electric field is stronger than a prescribed threshold of 120 kV m −1 . If multiple grid points at different horizontal locations exceed the threshold of the electric field, the center of the vertical cylinder is randomly determined among the multiple grid points (Sato et al., 2019). The height of the vertical cylinder is the same as the domain height. Each discharge process reduces the charge amount in each 3D grid point inside the cylinder if the charge density exceeds 0.1 nC m −3 (Fierro et al., 2013). This discharge process is sequentially repeated for each vertical cylinder until the electric field in each grid point becomes weaker than the threshold. The number of discharge processes initiated in each grid point is a 3D output variable from the parameterization of Fierro et al. (2013), corresponding to the number of lightning flashes in each grid point. Fierro et al. (2013) reported that their parameterization gives reasonable lightning flashes that are qualitatively comparable to observations, and Sato, Hayashi, et al. (2021) and Sato, Miyamoto, et al. (2021) reported that the SCALE-RM reasonably simulates lightning flashes in an idealized tropical cyclone and heavy precipitation events in Japan.
The LETKF, a variant of the EnKF, assimilates observations using the ensemble-based forecast error covariances with an observation operator that converts the model state variables into the observed variables. The LETKF provides the best estimates of the state variables, including zonal and vertical winds, pressure, temperature, the water vapor mixing ratio in addition to microphysics-related variables such as the mixing ratios, number concentrations, and charge densities of five hydrometeors, and the electric potential. Here, the LETKF was extended to assimilate vertically accumulated lightning flash observations at an 8-km mesh or simply "GLM observations," which correspond to flashes observed by GLM rather than events or groups (Goodman et al., 2013). The horizontal resolution of 8 km corresponds to that of GLM at nadir (Goodman et al., 2013).

Observing System Simulation Experiments (OSSEs)
We performed OSSEs in which synthetic observations were generated from a nature run simulating a supercell thunderstorm using a 3-K elliptical warm bubble superimposed on a horizontally homogeneous convectively unstable base state. The bubble had a horizontal radius of 10 km and a vertical radius of 1.5 km and was placed at the center of the domain. The initial thermodynamic and wind profiles were based on Weisman and Klemp (1982) and Weisman and Rotunno (2000), respectively.
The initial ensemble perturbations were generated using the method of Aksoy et al. (2009), who perturbed the warm-bubble location, its amplitude, and wind profiles. The warm-bubble location was randomly perturbed by a Gaussian distribution with a mean of zero and a standard deviation of 20 km. The ensemble-mean location was the domain center, that is, the same as the nature run. The warm-bubble amplitude was also randomly drawn from a Gaussian distribution with a mean of 3 K and standard deviation of 0.5 K. If the random number was less than 1.5 or larger than 4.5, it was replaced by 1.5 or 4.5, respectively. Following Coffer et al. (2017), the horizontal-wind perturbations were randomly drawn from a Gaussian distribution with a mean of zero and a standard deviation of 2 m s −1 without correlations. The resulting perturbed wind profiles were used as the initial and boundary conditions for each ensemble member. Deterministic forecasts from the ensemble-mean used the unperturbed perfect wind profile as the boundary conditions.
Synthetic observations were assimilated from ground-based X-band radar and a GLM sensor. The X-band radar mimics the PAWR and was placed near the domain center. The PAWR observations included radar reflectivity and radial Doppler velocity and were simulated by the observation operator originally developed by Xue et al. (2009) and modified by Amemiya et al. (2020). For simplicity, no attenuation was considered. To effectively suppress spurious convection, this study used the method by Aksoy et al. (2009) in which radar reflectivity observations lower than 15 dBZ were replaced with 5 dBZ. Similar to Aksoy et al. (2009), to reduce the number of the radar observations, both reflectivity and radial velocity observations were thinned at three-by-three horizontal grid points. No thinning was applied in the vertical direction. Therefore, the PAWR provided high-vertical resolution observations every 5 min. We assumed that the PAWR observations were available in the entire domain, except above a 12-km height. To prepare synthetic GLM observations, we first vertically accumulated 3D lightning flashes outputted from the Fierro et al. (2013) parameterization. The resulting data, at a horizontally 2-km mesh, were further accumulated into an 8-km mesh. Following the success by Mansell (2014), this study assimilated both nonzero and zero flashes.
Observation errors were drawn from Gaussian distributions. Following Maejima et al. (2017), the standard deviations were selected as 3 m s −1 for the radial velocity and 10% of the observed values with a minimum threshold value of 2 dBZ for reflectivity. For the GLM observations, which are vertically accumulated flash numbers, the standard deviation was set to 1. However, no errors were added to the observations with zero flashes. Figure 1 shows the workflow of the DA experiments. Before assimilating the observations, we performed 65-min spin-up ensemble forecasts. The ensemble size was 320 and the assimilation interval was 5 min. Overall, we performed 12 assimilation cycles from t = 65 min to t = 120 min. Covariance inflation was achieved by the relaxation to prior perturbation method with a coefficient of 0.5 (F. Zhang et al., 2004). The standard deviations of the horizontal localization scales were selected as 2 km for the radar data (both reflectivity and Doppler velocity) and 30 km for the GLM data based on the results of several additional sensitivity experiments. These localization scales for the radar and GLM data corresponded to the radius of influence (where the observation impact becomes zero) of approximately 7.3 and 109.5 km, respectively. The assimilation of the GLM observations was still beneficial, even with a shorter horizontal localization scale of 10 km. However, the forecast errors for the vertical velocity slightly increased. The vertical localization scale for the PAWR observations was selected as 1 km, but no vertical localization was applied to the GLM observations because the GLM observations were 2D. The observation error standard deviations in the LETKF were selected as 5 dBZ for radar reflectivity, 3 m s −1 for radial velocity, and 1.0 for the GLM.
Four experiments, named NODA, RDA, GLMDA, and NO GLMDA, were conducted. NODA was a model-free run without DA. RDA assimilated every-5-min radar observations. GLMDA and NO GLMDA were initiated by the analysis ensemble of RDA at t = 90 min. GLMDA only assimilated GLM observations at t = 120 min, whereas NO GLMDA assimilated no observations at t = 120 min. In GLMDA and NO GLMDA, we also ran 30-min deterministic forecasts from the analysis ensemble-mean at t = 120 min. No IR radiance observations were assimilated in this study.
In the ensemble-based correlation analysis, we investigated three types of observations: PAWR, IR radiance, and GLM. Following  and Honda, Kotsuki, et al. (2018), IR radiances were calculated using a radiative transfer model, RTTOV 12.3 (Saunders et al., 2018). Previous studies showed that moisture-sensitive bands such as the 6.2 and 6.9 μm are closely correlated with the atmospheric state (e.g., F. Zhang et al., 2016;Y. Zhang et al., 2021). We focused on a moisture-sensitive band of 6.9 μm, which has been successfully assimilated by several previous studies (Honda, Kotsuki, et al., 2018;Honda et al., 2019). We performed 30-min ensemble forecasts from the analysis ensemble at t = 90 min in the RDA, and these observations were then simulated.

Nature Run
Figures 2 and 3 exhibit the typical structure of a supercell thunderstorm in the nature run, with a hook echo pattern near the updraft core (Figures 2a, 3a and 3b). The supercell comprised a widespread anvil cloud near the tropopause (Figures 2a and 3b). The satellite IR radiances displayed a broad cloudy area and could not identify the convective core ( Figure 2c). Indeed, in a real-world scenario, GOES-16 IR radiances observed a widespread anvil cloud in the mature stage of a convective system (Y. Zhang et al., 2018).
The GLM observation peaked at approximately 20 km away from the updraft core (Figure 2d), which is consistent with the observations of Calhoun et al. (2013). The charge structure of the storm was a typical tripole, and numerous lightning flashes occurred in the midtroposphere and upper troposphere (Figures 2b, 3c and 3d).
Among the various observation systems compared here, PAWR is the most useful one for capturing the 3D structure of the supercell because radar reflectivity directly relates to the mixing ratios of hydrometeors. In addition, the GLM plays an essential role in detecting the convective core in the matured convective cloud, especially where radars are not available due to the terrain or over the ocean, and satellite IR radiances cannot provide the information below the anvil clouds.

OSSEs
Here, we describe the results of OSSEs with the assimilation of the PAWR observations. Figure 4 compares the analyses from NODA and RDA with the nature run. In RDA, the horizontal distribution of radar reflectivity is very similar to that in the nature run, indicating that the assimilation of the PAWR observations significantly improved the quality of the hydrometeor field. This is consistent with previous studies on precipitation events (Miyoshi, Kunii, et al., 2016;. Although no electric observations were directly assimilated, the analyzed charge structure in RDA was close to the nature run because the charge amount is associated with the mixing ratio of hydrometeors. In NODA, the supercell structure is unclear owing to the large spread of the supercell location. Assimilating the PAWR observations improves analyses and forecasts of the supercell in terms of the distributions of hydrometeors and charges. Figure 5 compares the forecasts of NODA and RDA with the nature run. The forecast in RDA predicted two distant areas of high radar reflectivity and widespread lightning flashpoints, consistent with the nature run (Figures 5a-5c). Moreover, the charge distribution in the RDA forecast was significantly closer to that in the nature run than the NODA forecast (Figures 5d-5f).

Ensemble-Based Correlations With Various Observation Types
To explore the potential impacts of assimilating various observation types, we investigated ensemble-based correlations between each observation type and the atmospheric state. High correlations, either positive or negative, corresponded to large analysis increments created by the observation. As shown in Figure 1, to analyze ensemble-based correlations, we performed 30-min ensemble forecasts from the analysis ensemble at t = 90 min in RDA and simulated various observations, including the PAWR, GLM, and IR radiances. Here, the GLM observations corresponded to the accumulated number of lightning flashes in the previous 30 min. Figure 6 shows the ensemble-based correlations between each observation near the center of the storm and the total mixing ratio of hydrometeors. As expected, the radar reflectivity observation exhibited very high correlations with the hydrometeors (Figure 6a). The vast majority of the ensemble had anvil clouds, so that the IR radiance observation near the storm center exhibited no correlation with the amount of hydrometeors ( Figure 6b).
Alternatively, lightning flash observations showed high correlations below the cloud top, where IR radiances displayed no correlations. Below the cloud top, the charge structure varied among the ensemble members, so that the GLM observations exhibited higher correlations than the IR radiances. This suggests that assimilating GLM observations could potentially improve the hydrometeor distribution below the cloud top. The GLM observations exhibited clear relationships with not only the hydrometeors but also other variables such as wind velocity and the amount of water vapor. Figure 7 shows the ensemble-based correlations with these variables. As expected, the GLM observations showed a correlation with the vertical velocity, yet the correlation with the vertical velocity was lower than that with water vapor (Figures 6c and 7a). This is consistent with previous observations such as those of Basarab et al. (2015), who showed that the echo volume is a better indicator of lightning activity compared to updraft intensity. Moreover, the correlations between the GLM observations and the vertical velocity were marginally shifted from the GLM observation location (Figure 7a), indicating that too strong an updraft reduces lightning flashes (e.g., Calhoun et al., 2013).
The GLM observations exhibited clear correlations with the amount of water vapor and temperature. In particular, the correlations with water vapor were comparable to those with hydrometeors shown in Figure 6c, consistent with the effectiveness of the simple pseudo-water vapor assimilation method in Fierro et al. (2016Fierro et al. ( , 2019 and Hu et al. (2020). The GLM observations also displayed clear positive (negative) correlations with temperature in the upper level (lower level), reflecting convection and diabatic heating. Similarly, Kong et al. (2020) showed that the flash extent density derived from the graupel volume or mass highly correlates with both thermodynamic and dynamical model variables. In contrast, the IR radiance observations displayed no clear relationships with these three variables, especially below the cloud top. The IR radiance observations were associated with the amount of water vapor and the temperature above the cloud top. However, these results are case dependent. F. Zhang et al. (2016) showed that IR radiance observations have clear relationships with the wind field in a tropical cyclone. Therefore, further investigation of other phenomena such as tropical cyclones is required. Figure 8 shows the averaged ensemble-based correlations in the cloudy regions. The GLM observations have stronger correlations with the hydrometeors, vertical velocity, and temperature below the cloud top, compared to the IR radiance observations. The GLM observations were closely associated with the amount of water vapor below the midtroposphere, whereas the IR radiances were correlated with the water vapor and temperature only near the cloud top. As shown by Y. Zhang et al. (2021), the correlation characteristics of the IR radiances in the cloudy regions were independent of the choice of water vapor bands ( Figure S1). These results indicate that the GLM observations provide valuable information below the cloud top, especially for hydrometeors and moisture, unlike the IR radiance observations. The GLM observations were correlated with various hydrometeors. Figure 9 illustrates that the GLM observations have clear relationships with graupel, snow, and cloud ice in the upper level and rain and cloud water in the lower level. Among the five hydrometeors, graupel showed the highest correlation with the GLM observations, yet the other hydrometeors were also highly correlated, probably because developed convection with much graupel also contains a large amount of the other hydrometeors. In addition, a recent observational study by Hayashi et al. (2021) reported that the lightning activity showed a higher correlation with the ice-related volume that includes ice, snow, graupel, and hail compared to that with the graupel and hail volume. In contrast, the IR radiances have no clear correlations with any hydrometeors below the cloud top.
Finally, we investigated the effects of accumulating the GLM observations in time. Figure 10 shows the ensemble-based correlations of the GLM observations accumulated for 5, 10, and 30 min. The results indicate that the longer the accumulation period, the higher the correlations with hydrometeors. This would be related to the random process for the origin of lightning flashes in the present model parameterization of   Fierro et al. (2013), so that instantaneous lightning flashes are not as highly correlated with the atmospheric state. A 30-min accumulation for the GLM observations is applied in Section 5.

Assimilation Impacts of GLM Lightning Flash Observations
Here, we aimed to demonstrate the potential assimilation impacts of the GLM lightning flash observations and thus performed two additional experiments (GLMDA and NO GLMDA) with and without the assimilation of the GLM observations. In GLMDA, the total amount of hydrometeors in the vertical cross sections was increased compared to that of NO GLMDA and was more similar to the nature run (Figures 11a-11c). In addition, in GLMDA, the charge distribution was improved by the GLM assimilation (Figures 11d-11f). These results indicate that the assimilation of the GLM observations has the potential to improve the analyses of hydrometeor and charge distributions.
Assimilating GLM observations improved the forecast. Figure 12 compares the distributions of hydrometeors and the root mean square errors (RMSEs) relative to the nature run between forecasts from the ensemble-mean analyses in GLMDA and NO GLMDA. In GLMDA, a convective core in the northwest was weaker than that in NO GLMDA, which is consistent with the nature run (Figures 12a-12c). We also observed that the forecast improved for hydrometeors and charge density in most vertical layers (Figures 12d and 12f), except for the amount of hydrometeors initially. This was attributed to the poor quality of the first guess above the 11-km height, where no radar observations were assimilated. Vertical velocity was slightly improved (Figure 12e). Recent observational studies indicated that lightning flashes are not always associated with vertical velocity (Basarab et al., 2015;Carey et al., 2019). In summary, the single-cycle OSSEs demonstrated Figure 8. Radius-height cross sections of averaged ensemble-based correlations of (a-d) GLM and (e-h) IR radiances with (a, e) The mixing ratio of total hydrometeors, (b, f) Vertical velocity, (c, g) Meridional velocity, and (d, h) Temperature. The 3D ensemble-based correlations were computed using the 30-min ensemble forecast from the analysis at t = 90 min in data assimilation (DA) only where GLM flash >0 in at least 32 of the 320 members, and the column maximum mixing ratio of total hydrometeors q h > 1 g kg −1 was then averaged into the radius-height cross sections based on the horizontal distance from each observation.  that the assimilation of the GLM observations could potentially improve NWP, as revealed by the ensemble-based correlation analyses in Section 4.

Conclusions
To improve convective-scale NWP, it is essential to assimilate various observations and obtain improved initial conditions. Despite their close relationships with cloud particles within convective clouds, lightning flashes have not been systematically utilized for NWP. By conducting OSSEs for an idealized convective case with the LETKF coupled with the SCALE-RM, including the explicit electric processes, this study aimed to investigate the potential of assimilating GLM lightning flash observations to improve NWP.
To estimate the potential impacts of lightning observations, ensemble-based correlations between various observation types and atmospheric variables were compared. The GLM observations exhibited high correlations with the amount of hydrometeors, water vapor, and temperature in the clouds, even when the IR radiances had unclear correlations due to widespread anvil clouds. The correlations were larger when the GLM observations were accumulated for a longer period. Therefore, GLM observations provide information on the atmospheric state within the cloud, where the IR radiances cannot do so owing to the presence of anvil clouds.
The impacts of assimilating the GLM observations were further verified by additional OSSEs. The results indicated that the assimilation of GLM observations significantly improved the analyzed distributions of the hydrometeors and charge density. Moreover, the GLM assimilations contributed to improving the forecast skills in terms of vertical velocity as well as hydrometeors and their charge amount. This study highlighted the potential impacts of GLM lightning flash observations, which are expected to contain more direct information below the cloud top than IR radiance observations. Moreover, this study investigated a single idealized convective case without cycling DA. However, this study was limited to one-cycle DA and perfect-model OSSEs as an essential first step to assimilate real lightning flash data effectively. An important next step is to perform cycling DA for an extended period and include different cloud types. In addition, further research should focus on determining an improved method of assimilating lightning observations that takes non-Gaussian errors and biases from imperfect microphysics parameterizations into account. Another important future direction will be to assimilate lightning observations in the real world together with all-sky satellite IR radiances. Furthermore, it would be necessary to use a more sophisticated lightning parameterization such as Mansell et al. (2002) instead of the Fierro et al. (2013) scheme for future research using high-resolution NWP models.

Data Availability Statement
Data and source code required for simulations are stored at http://doi.org/10.5281/zenodo.5506403. , and (f) Total charge density (nC m −3 ) relative to the nature run between GLMDA and NO GLMDA forecasts. Cool (warm) colors correspond to improvement (degradation) by assimilating the GLM observations at t = 120 min.