Tomioka et al.
Progress in Earth and Planetary Science(2023) 10:60
https://doi.org/10.1186/s40645-023-00592-w

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Advantage of bulk lightning models
for predicting lightning frequency over Japan
Takumi Tomioka1,2, Yousuke Sato3,4* , Syugo Hayashi5 , Satoru Yoshida5   and Takeshi Iwashita6   

Abstract
This study examined the performance of an explicit bulk lightning model coupled with a meteorological model
for forecasting lightning by numerical weather prediction over Japan. The evaluation was conducted by comparing the lightning predicted by the explicit bulk lightning model, diagnosed empirically by the numerical model,
and observed by ground base measurements. From the results, the bulk lightning model performed better in termsof lightning frequency than did the diagnostic scheme, which overestimated the lightning frequency, although there
were no appreciable differences in the score of each method for the geographical distribution and time correlation compared with the observations. These results suggest that the explicit bulk lightning model is advantageous
for predicting lightning frequency. The sensitivity of the simulated lightning to the choice of cloud microphysicalmodel was also examined by using a two-moment and a one-moment bulk microphysical scheme. Sensitivity experiments on the choice of microphysical model indicated that the two-moment bulk scheme reproduced the observed
lightning well, while the one-moment bulk scheme overestimated the lightning frequency. Analyses suggested
that the overestimation of the lightning in the one-moment bulk scheme originated from active charge separation
by riming electrification, in which graupel was produced more frequently and was assumed to fall faster. These results
suggest that the explicit bulk lightning model with the two-moment bulk microphysical scheme offers an alternative
to conventional lightning prediction methods.
Keywords Lightning, Bulk lightning model, Cloud microphysical model, Numerical weather prediction
1 Introduction
Lightning is a threat to human life because it can injure
and even kill people. Every year, lightning strikes killseveral people in Japan and more than 1000 people*Correspondence:
Yousuke Satoyousuke.sato@sci.hokudai.ac.jp
1
Graduate School of Science, Hokkaido University, N10W8 Kita‑Ku,
Sapporo, Hokkaido 0600810, Japan
2
Present Address: Tokio Marine dR Co., Ltd., Tokyo, Japan
3
Faculty of Science, Hokkaido University, N10W8 Kita‑Ku, Sapporo,
Hokkaido 0600810, Japan
4
RIKEN Center for Computational Science, 7‑1‑26
Minatojima‑Minami‑Machi, Kobe, Hyogo 6500047, Japan
5
Meteorological Research Institute, 1‑1 Nagamine, Tsukuba, Ibaraki
3050052, Japan
6
Information Initiative Center, Hokkaido University, N11W5 Kita‑Ku,
Sapporo, Hokkaido 0600811, Japanworldwide (e.g., Fujibe 2017; Curran et al. 2000). Lightning is also a serious threat to today’s highly electrified
society by damaging electronics and other equipment.
Each year in Japan, the cost of lightning damage exceeds
100 billion yen (Yokoyama 2002), and the cost have been
increasing based on the insurance claim costs (https://​www.​c ity-​n et.​o r.​jp/​w p-​conte​nt/​u ploa​d s/​2015/​02/​0 0_​
kouky​ou_​raiga​itais​aku.​pdf, accessed on 2022/10/12).
However, such damage can be prevented by evacuationand equipment protection if the lightning is predicted inadvance, which is why accurate lightning prediction isrequired to reduce the damage from lightning.
Operational numerical prediction of lightning is not
yet feasible, but lightning can be diagnosed empirically based on observations. For example, the National
Severe Storms Laboratory and the European Centre for
Medium-Range Weather Forecasts diagnose lightning© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://creativecommons.org/licenses/by/4.0/.

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

based on hydrometeors in the clouds calculated bytheir forecasting models (Lopez 2018; https://​www.​nssl.​
noaa.​gov/​e duca​tion/​s vrwx​101/​light​ning/​forec​a sting/,
accessed on 2022/11/28). Meanwhile, the Japan Meteorological Agency (JMA) provides lightning guidance(Tsuchida 2018) by forecasting the probability of lightning using logistic regression based on explanatory variables considered to be related to lightning, such as theShowalter stability index, convective available potential
energy (CAPE), and precipitation amount, simulated byits model for operational numerical weather prediction(NWP). JMA also provides thunder nowcasting (Kasahara 2010a, b, 2011) by diagnosing lightning activity up
to 1 h ahead based on the observed lightning and precipitation. These diagnostic methods are accurate and havehelped greatly to reduce lightning damage, but they cannot predict the frequency of lightning. In addition, JMA’s
lightning guidance predicts the probability of lightning ineach 10-km or 20-km grid for every 3 h and so cannot be
used for forecasts finer than its grid size or time step, andthe lead time of JMA’s thunder nowcasting is between
one and a few hours, which is short compared to the
NWP lead time. McCaul et al. (2009) developed a scheme
to diagnose lightning frequency based on the microphysical properties of lightning clouds, such as upward flux ofgraupel in the mixed-phase region at − 15 °C and vertically integrated ice mass. This scheme can even forecast
the frequency of lightning with the NWP lead time, but it
uses empirical equations, which are not always applicablefor predicting lightning in storms that are not includedin the dataset used to fit the empirical equations, such as
storms that have never been experienced before, or in the
future climate.
Numerical simulation by an explicit bulk lightningmodel (BLM) (e.g., Barthe et al. 2005, 2012, 2016; Bovalo
et al. 2019; Hayashi 2006; Helsdon et al. 2001; MacGorman et al. 2001; Mansell et al. 2002, 2005; Mansell andZiegler 2013; Takahashi 1984) is a powerful tool for predicting lightning based on the predicted charge density ofhydrometeors. However, explicit BLMs are not suitablefor operational NWP because of their high computational cost, and they are yet to be used operationally. This
lack of experience with prediction by BLMs means that
their advantages are as yet unclear. The computationallyexpensive parts of BLMs are calculating the discharge
and solving the inverse matrix to calculate the electricfield from the charge density.
Fierro et al. (2013) developed a lightning discharge
scheme with low computational cost and that performed
well in previous studies (e.g., Fierro et al. 2013; Dafis et al.
2018; Chen et al. 2020; Lu et al. 2022; Sato et al. 2022).
The calculation for solving the inverse matrix should
also be accelerated. Most of those studies were focusedPage 2 of 22

on high-flash-rate events such as summer thunderstorms
and winter thunderstorms with strong convection caused
by low-pressure systems. However, in Japan, in addition
to the above types of lightning, some of the lightnings
associated with low cumulonimbus clouds along the
Sea of Japan in winter have special characteristics. This
type of lightning event occurs in other limited areas of
the world, such as around the Great Lakes in the UnitedStates along the west coast of both Norway and Israel(Montanya et al. 2016; Holzworth et al. 2019). The characteristics of winter lightning differ from those of summer lightning in several aspects. Storms with only a few
flashes are often observed and are called “single lightning flash storms” (Michimoto 1993; Hayashi and Marui2016). The positive cloud-to-ground (CG) flash frequently occurs compared with other lightning. The electric current in a single lightning strike is sometimes large,
and the lightning damage is much more serious than that
due to single lighting in summer (Goto and Narita 1995;
Ishii and Saito 2009). The length of the flash is longer
than that of summer lightning, so a wider area of charges
is removed (Yoshida et al. 2019). However, although thesefeatures would affect the performance of explicit BLMs,
only a few previous studies have used them for this typeof winter lightning (Altaratz et al. 2005; Kawano et al.
2019).
For the reality in Japan, the only previous studies to use
BLMs were those by Sato et al. (2022), Hayashi (2006),
and Kawano et al. (2019), but they did not evaluate theperformance of the BLMs. To apply a BLM to operational forecasting, its performance should be evaluatedin various cases including the specific winter lightning
described above. In addition, one-moment cloud microphysics models are often used in operational forecast
models because of their low computational cost, but Satoet al. (2022) used a two-moment cloud microphysicalscheme. The one-moment scheme predicts the mixing
ratio of hydrometeors, while the two-moment scheme
predicts the number concentration in addition to themixing ratio and so offers a more accurate representationof the particle size distribution of hydrometeors but at
higher computational cost. Barth et al. (2005) noted thedependence of BLMs on the choice of cloud microphysical scheme, so it is necessary to investigate the sensitivityof a BLM to that choice.
In this study, we conducted numerical simulations
using an explicit BLM (Sato et al. 2019) implementing thelightning discharge scheme of Fierro et al. (2013) and targeting six typical lightning events in Japan. By evaluating
the BLM for the six events, our aim was to examine its
advantages. We also conducted sensitivity experimentswith another cloud microphysical scheme. To evaluate
the BLM’s performance and indicate its advantages, its

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

results were compared with those of conventional lightning diagnostic methods and observations by groundbase measurements.

2 Data and methods
2.1 Model

The model used in this study was version 5.3.6 of SCALE(Scalable Computing for Advanced Library and Environment) (Nishizawa et al. 2015; Sato et al. 2015), with an
explicit BLM implemented in SCALE (Sato et al. 2019,2022) used to calculate the electrification of hydrometeors and lightning. Often called riming electrification, thenoninductive charge separation was considered basedon a look-up table by Takahashi (1978). The discharge
process was calculated by the parameterization of Fierro
et al. (2013): in this scheme, the initiation point of a discharge is determined as the point where the absolutevalue of the electric field exceeds the threshold value( Eint ), and the charge of the hydrometeors in a cylinder
with the specific radius (rcyl ) centered at the discharge
initiation point is neutralized; for the present work, Eint
and rcyl were set at 110kVm−1 and 5km, respectively,
based on Fierro et al. (2013) and Sato et al. (2022). In the
BLM, the Bi-Conjugate Gradient STABilized method (BiCGSTAB; van der Vorst 1992) is used to solve the inverse
matrix for calculating the electric field from the charge
density, and we used symmetric Gauss–Seidel preconditioning for Bi-CGSTAB to accelerate the calculationof the electric field; see Appendix 1 for the contributionof the preconditioning to reducing the computationalcost. The radiation and turbulence processes were calculated using the MSTRN-X radiation scheme (Sekiguchi
and Nakajima 2008) and the Mellor–Yamada turbulencescheme (Nakanishi and Niino 2006), respectively. The
bucket land model and single-layer urban canopy model(Kusaka et al. 2001) were used for calculating surface
temperature and moisture over the land in nonurban andurban areas, respectively, and the surface flux was calculated using a Beljaars-type bulk model (Beljaars and Holtslag 1991).
In sensitivity experiments, we used two cloud microphysical schemes to examine the sensitivity of the simulated lightning activity to that choice: (i) a two-momentbulk microphysical scheme of Seiki and Nakajima (2014)(hereinafter SN14), and (ii) a one-moment bulk schemeof Tomita (2008) (hereinafter T08).
2.2 Target cases

Our aim was to evaluate our BLM for typical cases oflightning in Japan. To do so, we conducted numerical simulations targeting six cases including both summer (called S1 and S2) and winter (called W1, W2, W3,
and W4) in Japan. S1 and S2 are heavy rainfall events inPage 3 of 22summer, and the others are lightning events in winter.
The characteristics and the motivation behind the choiceof each case are given below.
S1 and S2, which are the same cases discussed by Satoet al. (2022), are heavy rain events in summer, but thelightning frequency differs greatly between these twocases. S1 is a heavy rainfall event in 2017 over the Kyushu
area of Japan in summer; large amounts of rainfall were
triggered by tall convective clouds generated under
unstable conditions with high CAPE, and this event
occurred during a short period of time and with a highflash rate (Tsuji et al. 2020; Kawano et al. 2018; Kawano
and Kawamura 2020). S2 is a heavy rainfall event in 2018
over a wide area of western Japan in summer; low convective clouds generated under moist and stable conditions produced long-lasting heavy rainfall over the wholearea (Tsuji et al. 2020). The flash rate in S2 was lower thanthat in S1 (Kawano et al. 2018; Nakakita et al. 2019), andcomparing the performance of the BLM between the twocases enables us to examine whether it can distinguishthe lightning frequency in lightning events with high andlow flash rates.
The other cases are lightning event in winter. For W1,
we focused on a typical case of lightning in winter, i.e.,that triggered by strong convection such as in a lowpressure system. A low-pressure system passed through
Hokkaido in Japan on November 11, 2019, and rainfall
and lightning from cumulonimbus clouds accompanying
a cold front were observed near Sapporo. Evaluating theBLM for this case enables us to examine whether it can
reproduce the lightning that occurs in a winter thunderstorm with strong convection.
W2 is a snowfall event caused by a northwest monsoon
over Hokkaido on February 9, 2016. In this case, lightning was detected only a few times near Sapporo, andso evaluating the BLM for this case enables us to examine whether it can reproduce lightning with only a few
flashes in winter.
W3 and W4 are winter lightning cases that occurredin the Hokuriku region along the Sea of Japan. In bothcases, lightning and precipitation were caused by low
cumulonimbus clouds generated by a northwest monsoon. Evaluating the BLM for this case enables us to
examine its ability to reproduce winter lightning in the
Hokuriku region.
2.3 Experimental setup

The calculation domain and experimental setup for eachcase are summarized in Fig. 1 and Table 1, respectively.
Numerical simulations were conducted over the shaded
areas in Fig. 1, and the regions surrounded by dashedblack lines indicate the analysis areas. The data fromupwind areas and the first several hours were not usedTomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 4 of 22

Fig. 1 Calculation domain (shaded gray) and analysis areas (surrounded by thick dashed lines) for S1 (a), S2 (b), W1 and W2 (c), W3 (d), and W4 (e)Table 1 Experimental setup for each case. NX, NY, and NZ are the numbers of grid points of the model in the zonal, meridional, andvertical directions, respectively
Calculation period (UTC)

Analysis period (UTC)NX × NY

NZ

S1

2017/7/5 00:00–7/6 00:00

2017/7/5 06:00–7/6 00:00

300 × 300

60

S2

2018/7/6 21:00–7/8 00:00

2017/7/7 03:00–7/8 00:00

720 × 600

57

W1

2019/11/10 15:00–11/11 15:00

2019/11/11 6:00–11/11 15:00

600 × 600

60

W2

2016/2/9 00:00–2/10 00:00

2016/2/9 15:00–2/10 00:00

600 × 600

60

W3

2020/12/22 00:00–12/22 18:00

2020/12/22 09:00–12/22 18:00

900 × 660

60

W4

2019/12/11 00:00–12/11 21:00

2019/12/11 12:00–12/11 21:00

900 × 90060

LIDEN data was gridded to the nearest model gridfor the analyses to avoid the effects of lateral boundary
conditions and spin-up time (see Table 1). The horizontal grid spacing was 1 km. The vertical layer was divided
into 60 layers from 40 m above the ground to 21,545 m
at the model top, with the layer thickness increased from40 to 683 m toward the upper layers except for the simulation of S2; in S2, the vertical layer was divided into 57
layers from 40 m above the ground to 19,528 m at themodel top, with the layer thickness increased from 40
to 651 m toward the upper layers following Sato et al.
(2022), which reported that the cloud top height in S2
was lower than that of S1 and 57 layers is sufficient for
simulating the clouds in S2. The horizontal wind, air density, temperature, specific humidity, soil moisture, soil
temperature, and sea surface temperature (SST) of Mesoscale ANaLysis (MANL) provided by JMA (2022) wereused for the initial and boundary conditions; MANL has

a horizontal resolution of 5 km, 50 vertical layers, andtemporal resolution of 3 h. SST was fixed as the initialvalue, and the time interval of the model output was set
to 5 min.
2.4 Observational data

The 1-km-mesh synthetic-radar grid point values provided by JMA (2009) were used to evaluate the precipitation calculated by the BLM, and the lightning location
data observed by the Lightning Detection Network System (LIDEN; Ishii et al. 2014) operated by the JMA wereused to evaluate the performance of the BLM in that
respect. The LIDEN data was grided to the nearest model
grid, but the radar data was used with the original grid
setup of the observations. The BLM outputs the flash
origin density (FOD) defined by Eq. (9) in Fierro et al.
(2013), which corresponds to the lightning frequency

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

of CG and intra-cloud (IC) lightning combined. Basedon Fierro et al. (2013), the FOD is a good surrogate for
the Earth Networks Total Lightning Network (ENTLN;
https://​w ww.​weath​erbug.​com/​weath​er-​forec​a st/​now/),which is a system similar to LIDEN. Thus, the locationsof the initiation discharge points of CG and IC observedby LIDEN can be compared with the FOD.
It has been suggested that LIDEN can give false detections caused by electromagnetic waves other than lightning, and our pre-analyses showed that the location oflightning estimated by LIDEN can have errors. To reduce
the effects of such errors, JMA conducts quality control of LIDEN data using radar observations, and in thisstudy, we eliminated any LIDEN data with no precipitation observed by radar within a 10 km radius and within
10 min before or after, the same time interval as the radar
data, in the same way as done by JMA (Kasahara 2011).
2.5 Diagnostic scheme for lightning frequency

McCaul et al. (2009, hereafter MC) developed a lightning
threat diagnostic scheme, which is a widely used lightning
diagnostic method (e.g., Dafis et al. 2018; Fierro et al. 2013;
McCaul et al. 2020). In this scheme, lightning frequency(FMC) is calculated by a linear combination of the upward
flux of graupel in the mixed-phase region at −1 5 °C (F1) andthe vertically integrated mass of ice, snow, and graupel (F2):
F1 = k1 wqg ,(1)

F2 = k2

ρ qg + qs + qi dz,

(2)
(3)

FMC = r1 F1 + r2 F2 ,where qg, qs, and qi are the mixing ratios of graupel, snow,
and ice, respectively, w is the vertical velocity, ρ is the air
density, z is the height, k1 (= 0.042), k2 (= 0.2), r1 (= 0.95),
and r2 (= 0.05) are empirical parameters, and F1 and F2
were set to zero to reproduce the regional distribution oflightning if they are less than the threshold values of 0.01
and 0.4, respectively.
chance =Page 5 of 22

Hong and Lim 2006) and is not adjusted for SN14 or T08
of SCALE; adjusting the empirical parameterization forthem may improve the results. To examine reasonable
parameter values of the MC scheme for SCALE, we conducted 16 sensitivity experiments with k2 = 0.2, 0.15, 0.1,
and 0.05 and threshold values of 0.4, 0.6, 0.8, and 1.0, andthe results of those experiments are discussed in Sect. 4.
2.6 Scores for performance evaluation

The lightning model was evaluated in terms of absolutevalue, temporal variation, and geographical distribution.
To evaluate these three elements, the observed lightning
was converted to gridded data, and the following threestatistical scores were used.
To evaluate the errors in the absolute value of the lightning frequency, we calculated the root mean squared
error (RMSE) defined as

2
 
pre
 n
 i=1 Fi − Fiobs(4),
RMSE =nprewhere Fi is the accumulated lightning frequency calculated by the model for each time step for the model output, and Fiobs is the observed one.
To verify that the model can reproduce relative changes
in lightning frequency, we calculated the Pearson correlation coefficient (r) for the temporal relationship between
the model and the observations for hourly area-accumulated lightning frequencies.
The equitable threat score (ETS) was calculated to validate the distribution of precipitation and lightning quantitatively. Traditional ETS calculation methods use a 2 × 2
contingency table at individual grid points, in which theelements of the table are hits (correctly forecast events),
misses (observed but not forecast events), false alarms(forecast but not observed events), and correct negatives(correctly forecast nonevents). The ETS is given as

ETS =hits − chance,hits + misses + false alarms − chance

where(hits + misses)(hits + false alarms).
hits + misses + correct negatives + false alarms

To evaluate the performance of the BLM, we compared
the lightning distribution and its frequency predicted by
the BLM and diagnosed by the MC scheme with LIDEN
observations. Note that the MC scheme was developed
for the Weather Research and Forecasting (WRF) onemoment scheme (WRF Single-Moment Six-Species;

(5)

(6)

To reduce the effect of displacement between simulated and observed values, the ETS was calculated withina neighborhood framework (Clark et al. 2010; Lynn et al.
2015). In this framework, the grid point is defined as a hit
if the event is predicted at any grid point within a certain
radius R from the observation. Similarly, the grid point isdefined as a hit if the event is observed at any grid point

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

within a radius R of the predicted grid. The value of R was
set to 5 km, and the threshold for accumulated precipitation during the entire analysis period was set in increments of 50 mm for S1 and S2, in increments of 10 mm
for W1 and W3, and in increments of 5 mm for W2 andW4. Unlike the ETS for precipitation, to examine the spatiotemporal performance of the BLM, hits were considered if there was a LIDEN observation within 5 km and15 min before or after the FOD simulated by the model at
each time step according to Hayashi (2006). The ETS for
lightning verification is shown hourly as well as for theentire analysis period to capture time evolution. Because
the BLM output can take decimal values, grid points with
FOD values above zero were counted as having lightning.

3 Results
3.1 Evaluation of simulated precipitation

First, the simulated precipitation was compared with theobserved precipitation to assess the validity of each cloud
microphysical model. Figure 2 shows the geographicaldistribution of accumulated surface precipitation during
the analysis period in the analysis area. In all cases, SN14
reproduced the observed accumulated precipitation
approximately well, but there were discrepancies in somecases, such as shifts in the location of the precipitation
area. In S1, a large precipitation area exceeding 450 mm
was observed around 33.4°N, 130.8°E (Fig. 2a). An areawith large precipitation amount was reproduced in SN14,
but its location shifted to the southeast, near 33°N, 131°E(Fig. 2b). In S2 (Fig. 2d and e), SN14 reproduced well thelarge precipitation amount around 35.7°N, 135°E, but it
overestimated the precipitation amount over the Sea ofJapan. In W1 (Fig. 2g and h), SN14 reproduced the location of a large amount of precipitation near Sapporo, but
it overestimated the precipitation amount and did not
reproduce the precipitation areas near 42°N, 142°E. In
W3 and W4 (Fig. 2m, n, p, and q), the model reproduced
the observed precipitation area, but the simulated areawith a large amount of precipitation exceeding 30 mm
showed discrepancies with the observed one.
T08 also reproduced the observed precipitation well in
all cases (Fig. 2c, f, i, l, o, and r). However, T08 overestimated the maximum precipitation amount in S2 (Fig. 2f )and W2 (Fig. 2l), and the large-precipitation area shifted
southward as with SN14 in W4 (Fig. 2r).
As discussed by Sato et al. (2022), the above discrepancies originated from the mismatched meteorological fields and different cloud microphysical schemes in
the child model (i.e., SCALE in this study) and parent
model (i.e., MANL in this study), which often occur insuch dynamical downscaling simulations. Thus, suchdiscrepancies would be reduced by improving the initialand boundary conditions for the meteorological fields. In

Page 6 of 22

spite of such discrepancies, the ETS for the low precipitation threshold (Fig. 3) of both SN14 and T08 is around
0.25–0.5, which is similar to the ETS of the JMA mesoscale model used for operational weather forecasts (e.g.,
Honda and Sawada 2009), and therefore both SN14 andT08 generally reproduced the geographical distributionof observed precipitation. From these results outlined
above, it is concluded that the model reasonably simulated the distribution of the convective clouds that triggered the large amount of precipitation, if we consider
the uncertainty of the initial and boundary conditions
and model errors.
3.2 Evaluation of simulated lightning for summer cases

Figure 4 shows the geographical distribution of the flash
density observed by LIDEN, the FOD simulated by the
BLM, and flashes diagnosed by the MC scheme accumulated during whole analyses period for S1 and S2.
First, differences between the lightning simulated by
the BLM with SN14 and the observed lightning are discussed as a basic performance of the BLM, because theBLM with SN14 has been used previously (Sato et al.
2019, 2021, 2022). In S1, the BLM with SN14 reproduced the lightning distribution approximately, but thesimulated lightning distribution was displaced from themeasured lightning distribution as with the precipitation
area (Fig. 4a and b). The displaced lightning distribution
originated for the same reason as the displaced precipitation area as discussed in Sect. 3.1. The displacement
would be lessened by preparing reasonable meteorological fields and correcting the problem of the mismatchedcloud microphysical models, so we do not discuss these
displacements in this study.
The cumulative number of observed lightning events
exceeded 10 per grid point at the most frequent points,
but the BLM with SN14 underestimated the maximum
number as being about seven per grid point (Fig. 4a andb). Nevertheless, the BLM with SN14 overestimated theamount of accumulated lightning over the whole analysis
area by about 1.4 times the observed amount.
In S2 (Fig. 4f and g), the BLM with SN14 reproduced
well the geographical distribution and the maximum
accumulated lightning frequency per grid point, but it
overestimated the lightning frequency accumulated over
the analysis area by an order of magnitude. The overestimation of area-accumulated lightning frequency in thesetwo cases was originated from the grids with lightning
but small frequency. The lightning frequency with thesegrids are not seen clearly in Fig. 4f and g, but their contributions to the cumulative frequency over the wholeanalysis area were not so small.
These results indicate that the BLM with SN14 reproduced well the lighting with high flash rate (e.g., theTomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 7 of 22

Fig. 2 Geographical distribution of observed precipitation (a, d, g, j, m, p) and precipitation simulated with cloud microphysics schemes of SN14(b, e, h, k, n, q) and T08 (c, f, i, l, o, r) accumulated within the analysis period over the analysis area for each caseTomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 8 of 22

Fig. 3 Equitable threat score (ETS) at each threshold value of accumulated precipitation simulated with cloud microphysics schemes of SN14 (red)and T08 (blue) for each case for the whole analysis period and for the entire analysis area: S1 (a), S2 (b), W1 (c), W2 (d), W3 €, and W4 (f)Fig. 4 Geographical distribution of accumulated flashes detected by LIDEN (a, f), accumulated FOD simulated by BLM with SN14 (b, g) and BLM
with T08 (c, h), and flashes diagnosed by MC with SN14 (d, i) and MC with T08 (e, j) within the analysis period over the analysis area for S1 (a, b, c, d,e) and S2 (f, g, h, i, j)lightning in S1) within a factor of 2 (i.e., a factor of about
1.4 in this study). The results also suggest that the BLM
with SN14 overestimated the lighting with low flash rate
within a factor of 10 (i.e., one order of magnitude). However, the BLM with SN14 reproduced the difference ofabout one order of magnitude in the observed lightning
frequency between the two cases. Thus, we conclude
that the BLM with SN14 can distinguish between summer lightning with high and low flash rates, even though
it overestimates the frequency of lightning with low flashrate.

In both cases, the BLM with T08 simulated lightning
over wider areas than the observations and overestimatedthe lightning frequency per grid point and the area-integrated lightning counts (Fig. 4c and h). Especially in S2(Fig. 4h), a large amount of lightning was calculated evenon the west side of the analysis area, where no lightning
was observed, and the area-integrated lightning amounts
were nearly 100 times higher than the observations. The
amounts of area-integrated lightning in S1 and S2 were of
the same order of magnitude. These results indicate that

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 9 of 22

Fig. 5 a, b Time series of flashes detected by LIDEN (black), FOD simulated by BLM with SN14 (solid red) and T08 (solid blue), and flashes diagnosed
by MC with SN14 (dashed red) and T08 (dashed blue). c, d Time series of hourly root mean squared error (RMSE) of flashes observed by LIDEN
and FOD simulated by BLM with SN14 (solid red) and T08 (solid blue) or flashes diagnosed by MC with SN14 (dashed red) and T08 (dashed blue).
The mean RMSE for the entire analysis time is shown in the legends of (c) and (d). a, c and b, d show the results for S1 and S2, respectively

the BLM with T08 could not reproduce the difference inlightning frequency between summer lightning with highand low flash rates. The reason for the overestimatedlightning frequency by T08 and difference in the electrification in T08 and SN14 are discussed in Sect. 3.5.
Figure 4d, e, i, and j show the geographical distribution of lightning frequencies calculated by the MC
scheme using SN14 and T08. Compared to the observation results (Fig. 4a, and f ), the MC scheme with SN14
overestimated the lightning frequency in both casesand simulated a large amount of lightning even in areas
where lightning was not observed. These overestimations
originated from the overestimation of high clouds constructed by snow, as also discussed by Sato et al. (2022).
Based on Eq. (1), the overestimation of qs results in thelarge lightning frequency. When T08 was used, the MC
scheme showed results similar to those from the BLM
with T08. However, the MC scheme did not reproduce
the difference in lightning frequency between S1 and S2,
with a difference factor of only about 1.5 in the amount ofarea-integrated lightning between the two cases in bothcloud microphysics schemes.
From the results outlined above, we conclude that theBLM with SN14 can distinguish between lightning eventswith high and low flash rates, whereas the other methods(i.e., BLM with T08, and MC) cannot.
Quantitative evaluations of each scheme based on thethree scores support the above analyses. Figure 5 showstime series of the integrated lightning frequency over
the analysis area and its hourly RMSE for S1 and S2. The
RMSE for the entire analysis time is also shown in thelegend of Fig. 5c and d. In both cases, the BLM with SN14
had lower RMSEs than the other methods, indicatingthat it performs better in reproducing lightning frequencies. The large RMSE with the other methods originatedfrom the overestimation of the lightning frequency.
Scatter plots and correlation coefficients of the accumulated lightning frequency during every hour for
LIDEN and BLM or MC are shown in Fig. 6a and b. The
correlation coefficients of all the methods do not differ

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 10 of 22

Fig. 6 Scatter plots of hourly accumulations of flashes observed by LIDEN and FOD simulated by BLM with SN14 (a) and T08 (b) or flashes
diagnosed by MC with SN14 (a) and T08 (b), and time series of hourly equitable threat scores (ETS) for FOD simulated by BLM with SN14 (red)and T08 (blue) and flashes diagnosed by MC with SN14 (dashed red) and T08 (dashed blue) for S1 (c) and S2 (d). In the legends, the correlation
coefficient between observed flashes and FOD or diagnosed flashes for each case is shown as r (a, b), and ETS for the entire analysis time is shown(c, d). The lines in a and b show simulated value equals 10 times and 0.1 times observed value

appreciably in both cases. However, there are more plotsof the BLM with SN14 than those of the other methods
between the two straight lines in Fig. 6a and b showingthat the difference from observation is within a factor of
10. This result also suggests that the BLM with SN14 performs better at simulating lightning frequency.

Figure 6c and d show the ETS for each hour and for theentire analysis time (in the legends) in S1 and S2. Because
there was no lightning during 3:00–8:00 and 21:00–24:00
in S2, the value of ETS for these hours was zero, andthat for the entire analysis time for S2 was also small(Fig. 6d). In S1 using the BLM with SN14, the simulated

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 11 of 22

Fig. 7 Same as Fig. 4 but for W1 (a–e) and W2 (f–j)

precipitation area and lightning area were shifted to thesoutheast in the first half of the analysis period (figure
not shown), so the ETS values for this period are lower
than those with T08. Except for the period in S1, there
is no appreciable difference in ETS values between the
BLM with SN14 and the other methods.
Based on these results, we conclude that the BLM
with SN14 performs better at simulating lightning
frequency in summer events compared with othermethods, and the performances for geographical distribution and relative change with time are equivalent
with those of the other methods. The advantage of theBLM with SN14 is its ability to reproduce the differencein lightning frequency between summer lightning withhigh and low flash rates.
3.3 Evaluation of simulated lightning for winter cases
over Hokkaido

Figure 7 shows the geographical distribution of the
flashes observed by LIDEN, the FOD simulated by BLM,
and the diagnosed lightning frequency by the MC scheme
accumulated during whole analyses period for W1 andW2. In W1, the BLM with SN14 reproduced well theobserved lightning frequency and the geographical distribution (Fig. 7a and b). In W2, lightning was observed only
once during the analysis period (Fig. 7f ), which is similar
to the results of the BLM with SN14 in which no lightning was simulated (Fig. 7g). In both cases, the BLM with
T08 calculated a large amount of lightning and overestimated the observations (Fig. 7c and h). The lightning was
simulated by the BLM with T08 in almost all areas with
precipitation (Fig. 2i and r). The lightning diagnosed by
the MC scheme with both T08 and SN14 overestimatedthe measured lightning in W1 (Fig. 7d and e). In addition,lightning was diagnosed by the MC scheme with both
SN14 and T08 in W2 (Fig. 7i and j), when lightning was

rarely observed. These results indicate that the BLM with
SN14 outperforms the other schemes.
Figure 8 shows time series of the integrated lightning
frequency over the analysis area and its RMSE for W1
and W2. In both cases, the BLM with SN14 had lower
RMSEs than those with the other methods, which indicates that it performs better in reproducing lightning frequencies. The large RMSEs for the MC scheme and theBLM with T08 originated from the overestimation of thelightning. Figure 9a and b show scatter plots and correlation coefficients of the accumulated lightning frequency
during every hour for LIDEN and BLM or MC for W1
and W2, respectively. For W1, all plots due to the BLM
with SN14 are within a factor of 10, whereas those due tothe other methods are mostly within a factor greater than
10. The correlation coefficient by SN14 was higher thanthat by T08 with both BLM and MC (Fig. 9a and b), indicating that SN14 reproduced well the relative change oflightning with time. Note that Fig. 9b contains only one
data point for W2 because lightning was not observed
for most of the analysis period. The ETS is shown in
Fig. 9c, but W2 is not shown because lightning was rarely
observed and the ETS value was almost zero. There was
no appreciable difference in ETS values between the
BLM and MC schemes regardless of the choice of cloud
microphysical scheme for W1.
All the above results indicate that the BLM with SN14
performed better at simulating lighting in both cases andthat it has an advantage for reproducing winter lightning
triggered by strong convection by low-pressure systems,
which has a small flash rate compared with summer
lightning. Another advantage of the BLM with SN14 is
its ability to reproduce events with either no or very little lightning, which are not reproduced by the othermethods.

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 12 of 22

Fig. 8 Same as Fig. 5 but for W1 (a, c) and W2 (b, d)

3.4 Evaluation of simulated lightning for winter casesin Hokuriku

Figure 10 shows the geographical distribution of the
flash density observed by LIDEN, the FOD simulatedby the BLM, and the diagnosed lightning frequency by
the MC scheme accumulated during the whole analysis period for W3 and W4. In both cases, the BLM with
SN14 underestimated the lightning frequency per grid
point by more than one order of magnitude (Fig. 10b
and g), and it underestimated the area-accumulated
lightning frequency by between a fifth and a sixth compared with the observations. In contrast, the BLM with
T08 overestimated in both cases, which is similar to theother four cases (Fig. 10c and h). Overestimation was also
seen in the results of MC, except for MC with T08 in W4(Fig. 10j). From these results, it is difficult to determine
which method is better for W3 and W4, and so the scores
are useful.
The RMSE, scatter plots, correlation coefficients, andETS are shown in Figs. 11 and 12. First, we evaluate theresults for W3. The BLM with SN14 showed smaller
RMSE than the other methods (Fig. 11c). The correlation
coefficient with SN14 was larger than that with T08 forboth the BLM and the MC scheme (Fig. 12a and b), andthe correlation coefficients with the BLM and the MC
scheme were of the same magnitude. These results meanthat the performance of each model in terms of the relative change of the lightning frequency was similar. However, the BLM with T08 and the MC scheme with both
SN14 and T08 had more data points in the lower rightof the scatterplot, which corresponds to overestimationof lightning frequency. In addition, the ETS of all methods was similar, which means that the performance ofeach model in terms of the geographical distribution oflightning did not differ greatly. From these results, weconclude that there was no appreciable difference in theperformance of either method for reproducing the distribution and increase or decrease of lightning. However,the BLM with SN14 was advantageous for simulating thelightning frequency more accurately for W3, albeit with
underestimation (Fig. 10b).
In W4, the MC scheme with T08 and that with
SN14 had almost equivalent RMSEs, but SN14 had asmaller hourly RMSE during the time when lightning
was observed (between 13:00 and 17:00) (Fig. 11d).
The BLM with SN14 underestimated the lightning

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 13 of 22

Fig. 9 Same as Fig. 6 but for W1 (a, b, c) and W2 (a, b)

frequency compared to the MC scheme with T08, resulting in a larger RMSE. The correlation coefficient by the
MC scheme with T08 was larger than that by the othermethods, and the ETS did not differ greatly among all
the methods. Thus, for W4, the MC scheme with T08
outperformed the BLM, which was contrary to the performance evaluation for the other five cases. However,the MC scheme involves empirical parameters and is
not applicable for lightning events that have never been

experienced previously, and the tuning of those parameters based on a specific case does not always improve
the performance for other cases, as discussed in Sect. 4.
Considering these disadvantages of the MC scheme, theadvantages of the BLM with SN14 shown via the analysesof S1–W3 exceed the disadvantages, and we suggest thatthe BLM with SN14 is more suitable for lightning simulation than other methods.

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 14 of 22

Fig. 10 Same as Fig. 4 but for W3 (a–e) and W4 (f–j)

Fig. 11 Same as Fig. 5 but for W3 (a, c) and W4 (b, d)

3.5 Sensitivity of simulated lightning to cloud
microphysical models

Here, we discuss the sensitivity of the BLM and the MC
scheme to the choice of cloud microphysical model basedon the vertical profiles of the mixing ratio of hydrometeors around convective core, where the graupel generates

frequently (Fig. 13) and the amount of graupel charge
production over the cloudy grids (Fig. 14). Compared to
SN14, T08 had a lower mixing ratio of total hydrometeors, and its vertical distribution tended to shift to lower
layers where graupel was more dominant. According toprevious studies (Roh and Satoh 2014; Kondo et al. 2021),

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 15 of 22

Fig. 12 Same as Fig. 6 but for W3 (a, b, c) and W4 (a, b, d)

T08 shows the large conversion rate from snow to graupel and the fast fall velocity of graupel compared with
SN14. This would be because the graupel can be easily created by the accretion process once the graupel isgenerated by the riming according to the prescribed constant parameters such as interception parameter, slope
parameter and so on in the single moment bulk scheme
including T08 as briefly summarized in Appendix 3. Bythese characteristics of T08, hydrometeors in the atmosphere quickly become graupel and fall to the ground as

precipitation. Thus, in T08, graupel in the clouds tends
to be more dominant, especially around the height of− 10 °C where riming effectively occurs (Pflaum andPruppacher 1979), and the hydrometeors are distributed
in lower layers than in SN14 for all cases.
T08 has a higher collision rate between graupel andsnow or ice because the graupel mixing ratio tends
to be larger and the graupel falls faster (Kondo et al.
2021). This higher collision rate of graupel with snow orice causes active charge separation through the riming

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 16 of 22

Fig. 13 Vertical profiles of mixing ratios of hydrometeors (black), graupel (red), and ice (green) simulated with cloud microphysics schemesof SN14 (solid lines) and T08 (dashed lines) averaged over the convective core area in the analysis area during the entire analysis period. Bluedotted and dashed lines indicate the − 10 °C and − 15 °C altitudes, respectively, averaged over the entire clouded grid points in the analysis area
during the entire analysis period. The results of T08 and SN14 were both averaged for calculating the − 10 °C and − 15 °C altitudes. The hydrometeor
mixing ratio (Qhyd in the legends) is the sum of the mixing ratios of cloud, rain, ice, snow, and graupel, and the ice mixing ratio (Qice in the legends)
is the sum of the mixing ratios of ice, snow, and graupel (= qi + qs + qg). A convective core area was defined as a grid where the sum of the liquidwater and ice water paths exceeded 1 kg ­m−2electrification mechanism, and therefore T08 had alarger amount of charge separation and overestimatedthe lightning frequency in the BLM (Fig. 14). Note thatthe peak of the charge production by T08 was similar
to that by SN14 in S2 and W4 as shown in Fig. 14b andf. In these cases, the charge separation in T08 occurred
wider layer than SN14 due to the higher collision rate
graupel and snow or ice. This result supports the conclusion that the charge separation through the riming
electrification mechanism in T08 was more active compare to SN14. In summary, the BLM is highly sensitive
to the choice of cloud microphysical model, which may
cause overestimation of lightning frequency, especially in models in which a large amount of graupel issimulated.

4 Conclusions and discussion
By comparisons with observations for six summer andwinter cases over Japan, this study evaluated a BLM that
explicitly predicts lightning, as well as investigating its
sensitivity to the choice of cloud microphysical model.
The BLM was also compared with the lightning diagnostic method of McCaul et al. (2009). Through these sensitivity experiments and comparisons, the advantages of

the BLM were examined. The evaluation results showedthat the BLM with SN14 had the following advantages:(1) it reproduced well the lightning frequency of summer lightning with high flash rates; (2) it reproduced
the difference between summer lightning with high andlow flash rates; (3) it reproduced winter lightning witha lower flash rate than summer lightning; (4) it reproduced the absence of or only a few cases of lightning in
winter; and (5) it reproduced the lightning frequency in
Hokuriku without an empirical diagnostic scheme.
The results of the sensitivity experiments for the choiceof cloud microphysics model showed that SN14 reproduced the distribution and frequency of lightning qualitatively well. By contrast, T08 overestimated the lightning
frequency and simulated its distribution over wider areas
than the observations; this is because graupel is produced
easily by the large conversion rate from ice and snow
to graupel (Kondo et al. 2021), thus charge separation
through the riming electrification mechanism is caused
frequently by collisions with snow and ice in T08. These
results suggest that the BLM with SN14 could be a more
powerful tool for lightning forecasting than conventional
methods, especially for forecasting the lightning frequency. However, for the BLM to be used in operational

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 17 of 22

Fig. 14 As Fig. 13 but for vertical profiles of the amounts of positive (red) and negative (blue) graupel charge production simulated with the cloud
microphysics schemes of SN14 (solid lines) and T08 (dashed lines) averaged over the entire clouded grid points in the analysis area during the entire
analysis period. Black dotted and dashed lines indicate the − 10 °C and − 15 °C altitudes, respectively, averaged over the entire clouded grid points
in the analysis area during the entire analysis period. The results of T08 and SN14 were both averaged for calculating the − 10 °C and − 15 °C
altitudes. A clouded grid point is defined as a grid where the sum of the liquid water and ice water paths exceeded 1 g ­m−2

forecasting, future work is required involving forecast
experiments conducted over long periods of time.
As well as the advantage of the BLM, we must discuss its uncertainties. In this paper, we assumed thatthe electric field threshold (Eint) and the radius of charge
neutralization (rcyl) were constant when calculating discharges with the BLM, but the seasonal dependence of
these parameters must be discussed. Previous numerical simulations reported that the electric field required
to trigger the breakdown is larger at lower altitudes (e.g.,
Gurevich and Zybin 2005), which indicates that it may be
more appropriate to assume a larger Eint in the BLM in
winter cases. It is also observed that the average lightning
length is longer in winter (Yoshida et al. 2019), and therefore a longer rcyl would be more appropriate in the BLM
in winter cases. However, it is difficult to determine theseparameters correctly based on theoretical consideration and observations. Thus, the experiments conductedin this study were focused on how well the BLM could
perform with the same parameters as in previous studies (Fierro et al. 2013; Sato et al. 2022) and did not varythese parameters. In future work, sensitivity experiments
should be conducted to determine reasonable parametervalues.

In addition to the uncertainties of the BLM, it is necessary to discuss the empirical parameters of the MC
scheme. In this study, the coefficients and threshold values proposed by McCaul (2009) were used without modification in the calculation of the MC scheme, whereas
tuning them for SCALE may improve the performance ofthe MC scheme. Thus, we conducted several sensitivity
experiments to examine suitable parameter values for the
MC scheme in SCALE and discussed whether the advantages of the BLM persist in comparison with the MC
scheme with tuned parameters.
In the MC scheme, the lightning frequency is calculated
by a linear combination of F1 and F2 as shown in Sect. 2.5,
but F2 was dominant for most of the diagnosed lightning
frequency in this study. The coefficients and thresholds ofthe MC scheme are determined empirically, so tuning theseparameters against F2 may improve the overestimation tendency of the MC scheme. Accordingly, for MC with T08,
which scored better in the present study, we conducted
sensitivity experiments on k2 and the threshold values
under which F2 set zero. The results of the sensitivity experiments showed that the RMSE was reduced most effectively when k2 and the threshold were set as 0.1 and 0.6,
respectively (Table 2), although the correlation coefficientTomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 18 of 22

Table 2 Accumulated flash number over the whole analysis area and RMSE for each case, diagnosed with MC with T08 tuned tok2 = 0.1 and a threshold of 0.6. RMSEs of BLM with SN14 and MC with T08 not tuned are also given for comparison
S1

S2

W1

W2

W3

W4

Accumulated flash number

18,856

6954

405

0

713

241

RMSE (MC with T08, tuned)

81.56

31.11

7.13

0.10

11.42

39.03

RMSE (BLM with SN14)

116.86

22.76

3.93

0.10

12.89

35.95

RMSE (MC with T08, not tuned)

195.40137.18

30.39

0.70

45.60

24.70

and ETS did not change appreciably in any experiment. The
RMSE of the tuned MC with T08 was comparable to thatof the BLM with SN14, except for W4, in which the RMSE
became worse and larger than that of the BLM with SN14,
and the amount of accumulated lightning over the wholeanalysis area was underestimated compared to the BLM
with SN14. In addition, the difference in the diagnosed
number of lightning events between S1 and S2 was only
about a factor of three, so it is difficult to reproduce the difference between the summer cases with different lightning
frequencies, even if the parameters are tuned. From theresults shown above, it is concluded that the MC scheme
can be tuned to SCALE, but it is difficult to find suitableparameters for all cases because of the strong dependenceof the parameters on each case. This strong dependence ofthe parameters was also reported in previous observational
studies. Basarab et al. (2015) and Hayashi et al. (2021) proposed relational expressions for the observed microphysical properties in clouds and lightning frequency, similar to
the relational expression reported by McCaul et al. (2009).
They derived the relational expressions from the results ofobservations of several cases, but each relational expression
was highly dependent on the individual cases. Thus, it is
not easy to determine a relational expression that is applicable for all cases because the suitable parameters of an MC
scheme are not always suitable for other cases. The BLM’s
strength is that it can calculate lightning without such tuning for each case and can be used for cases that have neverbeen experienced before and for future experiments, which
are difficult to predict with empirical parameters.
In addition, we must discuss the uncertainty of the lightning observation system for comparing the BLM. According to Fierro et al. (2013), the FOD (i.e., the output of theBLM) is supposed to be a good surrogate for ENTLN, asystem similar to LIDEN, if it detects one or two points per
flash. However, the LIDEN observations seem not to meetthis requirement, especially for the winter lightning in
Hokuriku based on the analyses using the Lightning Location System, which is constructed similarly to LIDEN (e.g.,
Tajiri et al. 2021). In general, lightning observation systems

including LIDEN, which have long base lines between eachobservation station, are less accurate than three-dimensional lightning location data such as the observation dataof the Lightning Mapping Array (Thomas et al. 2001) and
BOLT (Yoshida et al. 2014; Wu et al. 2014). Comparing
the BLM with such three-dimensional lightning location
data would enable us to evaluate it in more detail, and suchcomparisons should be conducted in the near future.

Appendix 1
Application of preprocessing

In our BLM (Sato et al. 2019), Eqs. (7) and (8) are solved
at each time step to calculate the electric field (E) from thecharge density (ρe) with Neuman boundary condition atlateral boundary, and Dirichlet boundary with electricalpotential with zero at the ground and the model top:

∇ 2φ = −ρe,εE = −∇φ,(7)(8)

where φ and ε are the electrical potential and atmospheric permittivity, respectively, and ρe is the charge
density at each grid point, which is defined as the sumof the charge density of each hydrometeor at each grid
point. Equation (7) is Poisson’s equation, and the computational cost of a Poisson solver is large, especially ona state-of-the-art supercomputer using massively parallel computing, this being because global communicationis necessary. A Poisson solver based on the Bi-CGSTAB
method (van der Vorst 1992) is used in SCALE for solving
Eq. (7), and because of the large computational cost, theelapsed time for the lightning component is much larger
than those for the other physical components in SCALE(Table 3). To reduce the computational cost of the Poisson solver, three preconditioning methods (Saad 2003)
were applied for Bi-CGSTAB in SCALE, i.e., the Gauss–
Seidel method (GS), the symmetric Gauss–Seidel method(SGS), and incomplete lower–upper factorization (ILU).

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Table 3 Elapsed times for calculating physical components in
SCALE, evaluated by simulations of five time steps from 03 UTC
on 5 July, 2017 for S1
Component

WO

GS

SGS

ILU

Microphysics

3.58

3.62

3.56

3.60

Turbulence

0.85

0.81

0.84

0.85

Lightning (Mean
iteration numberof Bi-CGSTAB*)

155.50 (2137) 29.81 (193) 21.38 (105) 24.20 (119)Radiation

4.14

4.12

4.10

4.08

Land/Surface

0.052

0. 052

0.053

0.053

Total

171.06

45.42

37.11

40.19

*The mean iteration number means the number of iterations per every callof subroutine for Bi-CGSTAB. For this test, the iteration was finished when therelative residual norm becomes ­10–10 or smaller
WO Without preconditioning; GS With Gauss–Seidel preconditioning; SGS Withsymmetric Gauss–Seidel preconditioning; ILU With incomplete lower–upper
factorization preconditioning. For the calculation, the radiation processes wascalled once, and other physical processes were called five times

Fig. 15 Same as Fig. 5 but with the axis (b) and (d) were changed

Page 19 of 22

Table 3 shows that all three preconditioning methods
helped greatly to reduce the number of iterations done by
Bi-CGSTAB, and fewer iterations resulted in less elapsed
time (i.e., smaller computational cost) of Bi-CGSTAB. Of
the three preconditioning methods, SGS was the mosteffective, and the computational cost of Bi-CGSTAB was
reduced by about 84% compared with that without preconditioning. Based on the results shown above, SGS
preconditioning was applied.

Appendix 2
Enlarged view of Fig. 5 for small value of the lightning
frequency

Figure 15 is the same as Fig. 5, but with the axis (b) and(d) was changed.

Tomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Appendix 3
Charge separation mechanism in T08

In this part, we briefly describe the equation of the charge
separation and relating assumptions in T08 model. The size
distribution function for rain, snow, and graupel assumed by
T08 is Marshall-Palmer distribution (Tomita 2008):
n[r,s,g] (D) = N0,[r,s,g] exp −[r,s,g] D(9)where D, [r,s,g], and N0,[r,s,g] is diameter, slope parameter, and interception parameter for rain, snow, andgraupel, respectively. The monodisperse distribution isassumed for ice and cloud. The interception parameterof snow and graupel N0,s, N0,g , which are directly relating to the charge separation (riming electrification), aregiven as 3 × ­106 ­m−4, 4 × ­106 ­m−4, respectively. [r,s,g] areshown in Eq. (27) of Tomita (2008). These parameters are
fixed during whole calculation time and whole calculation area. The calculation of the collision and coagulation
is based on the stochastic collision equation (SCE; c.f.,Rogers and Yau 1989). Substituting Eq. (9) and the equation of monodisperse size distribution function into the
SCE, we obtain the tendency of the graupel and snow/ice(PGACS and PGACI) by the collision and coagulation as follows (Tomita 2008):
  1/2
π EGI N0,g cg qi Ŵ 3 + dg
ρ0
PGACI =
(10)
3+dg
ρ
4g

PGACS =π as |VG − VS |ESG N0,g N0,s4ρ Ŵ(bs + 1)Ŵ(3)
bg +1 3gs+2

Ŵ(bs + 2)Ŵ(2)(11)

bg +2 2gs+Ŵ(bs + 3)Ŵ(1)
bg +3

s

g

�

where EGI, ESG, is the collection efficiency of graupel andice, the collection efficiency of the graupel and snow. as,bs, bg, cg, and dg are prescribed empirical parameters, andρ0, ρ, and qi density at ground level, density, and mixing
ratio of ice. According to the similarity of the SCE andEqs. (10) and (11), the charge separation rate of graupel by the collision with ice (PGACI,e) and that with snow
(PGACS,e) are calculated as follows:  1/2π(1 − EGI )N0,g cg qi Ŵ 3 + dg
ρ0
PGACI,e = δρe
3+dg
ρ
4g
(12)

Page 20 of 22

πas |VG − VS |(1 − EGS )N0,g N0,s
PGACS,e =δρe4ρ Ŵ(bs + 1)Ŵ(3) + 2 Ŵ(bs + 2)Ŵ(2)
bg +1
bg +2
s 3gs 2g�Ŵ(bs + 3)Ŵ(1)+bg +3
s  g(13)

where δρe is charge separation late by the single collision
between graupel and snow or ice based on look up tableof Takahashi (1978).
According to the results from the previous studies suchas Roh and Satoh (2014) and Kondo et al. (2021), T08
tends to show the large conversion rate from snow/ice tograupel and the fast fall velocity of graupel, and it generates the large amount of the graupel compared with SN14
(Kondo et al. 2021). This large amount of graupel promotes the collision of graupel with snow/ice, and therefore, results in the large electrification rate in T08.
Abbreviations
Bi-CGSTAB	Bi-Conjugate Gradient STABilized method
BLM	Explicit bulk lightning model
CAPE	Convective available potential energy
CG	Cloud-to-ground lightning
ENTLN	Earth networks total lightning network
ETS	Equitable threat score
GS	Gauss–SeidelIC	Intra-cloud lightning
ILU	Incomplete lower–upper factorization
JMA	Japan meteorological agency
LIDEN	Lightning detection network system
MANL	Mesoscale ANaLysis
NWP	Operational numerical weather prediction
RMSE	Root mean squared error
SCALE	Scalable computing for advanced library and environment
SGS	Symmetric gauss–seidel
SN14	Seiki and Nakajima (2014)
SST	Sea surface temperature
T08	Tomita (2008)
WRF	Weather research and forecasting
Author contributions
TT conducted all simulations, analyzed the results, and drafted the manuscript.
YS helped in conducting simulations, interpreting the results, and writing themanuscript. SH provided LIDEN data and discussed the analyses and results.
SY gave useful comments on the analysis and examined the data qualityof LIDEN data. TI applied the preconditioning of the Bi-CGSTAB method to
reduce the computational cost. All authors discussed the results and read andapproved the final manuscript.
Funding
This work was supported by JSPS KAKENHI (Grant Number 20H04196), the JST[Moonshot R&D][Grant Number JPMJMS2286], the Support Program for Next
Generation Supercomputer of the Information Initiative Center, Hokkaido University, the collaborative research program 2022 Information Initiative Center,
Hokkaido University Support Program for Young Researchers, the collaborative
research of Hokkaido University and NTT Space Environment and Energy Laboratories, and the SECOM Science and Technology Foundation. Yousuke Sato
supported by funds from the Research Field of Hokkaido Weather Forecast
and Technology Development (endowed by the Hokkaido Weather Technology Center Co., Ltd.). The SCALE library was developed by Team-SCALE of theTomioka et al. Progress in Earth and Planetary Science(2023) 10:60

Page 21 of 22

RIKEN Center for Computational Sciences. The MANL data and the 1-km-meshsynthetic-radar grid point values were supplied via the Weather Business Consortium and the Research Institute for Sustainable Humanosphere of Kyoto
University (http://database.rish.kyoto-u.ac.jp/index-e.html). We also thankThinkSCIENCE, Inc. (Tokyo, Japan) for English language editing.

cyclone, and a winter storm. Mon Wea Rev 141:2390–2415. https://​doi.​
org/​10.​1175/​MWR-D-​12-​00278.1
Fujibe F (2017) Long-term change in lightning mortality and its relation toannual thunder days in Japan. J Nat Dis Sci 38:17–29. https://​doi.​org/​10.​2328/​jnds.​38.​17
Goto Y, Narita KI (1995) Electrical characteristics of winter lightning. J Atmos
Terr Phys 57:449–458. https://​doi.​org/​10.​1016/​0021-​9169(94)​00072-V
Gurevich AV, Zybin KP (2005) Runaway breakdown and the mysteries of lightning. Phys Today 58:37–43
Hayashi S (2006) Numerical simulation of electrical space charge density
and lightning by using a 3-dimensional cloud-resolving model. SOLA
2:124–127. https://​doi.​org/​10.​2151/​sola.​2006-​032
Hayashi S, Marui C (2016) Frequency and seasonal variation of single lightning
flash storm, “Ippatsurai”, in Japan (“Ippatsurai” no hasseihindo to kisetsu
henka) (in Japanese). J Atmos Electr 36:13–22. https://​doi.​org/​10.​1541/​
jae.​36.​13
Hayashi S, Umehara A, Nagumo N, Ushio T (2021) The relationship between
lightning flash rate and ice-related volume derived from dual-polarization radar. Atmosph Res 248:105166. https://​doi.​org/​10.​1016/j.​atmos​res.​2020.​105166
Helsdon JH, Wojcik WA, Farley RD (2001) An examination of thunderstormcharging mechanisms using a two-dimensional storm electrification
model. J Geophys Res Atmos 106:1165–1192. https://​doi.​org/​10.​1029/​
2000J​D9005​32
Holzworth RH, McCarthy MP, Brundell JB, Jacobson AR, Rodger CJ (2019)
Global distribution of superbolts. J Geophys Res Atmos 124:9996–10005.
https://​doi.​org/​10.​1029/​2019J​D0309​75
Honda Y, Sawada K (2009) Upgrade of the operational mesoscale 4D-Var
system at the Japan Meteorological Agency. CAS/JSC WGNE Res Activ
Atmos Oceanic Modell 39:01.11–01.12 https://​www.​wcrp-​clima​te.​org/​
WGNE/​BlueB​ook/​2009/​indiv​idual-​artic​les/​01_​Honda_​Yuki_​jnova.​pdf
accessed 10 Jan 2023
Hong SY, Lim JOJ (2006) The WRF single-moment 6-class microphysics scheme(WSM6). Asia Pac J Atmos Sci 42:129–151
Ishii M, Saito M (2009) Lightning electric field characteristics associatedwith transmission-line faults in winter. IEEE Trans Electromagn Compat
51:459–465. https://​doi.​org/​10.​1109/​TEMC.​2009.​20254​96
Ishii K, Hayashi S, Fujibe F (2014) Statistical analysis of temporal and spatial
distributions of cloud-to-ground lightning in Japan from 2002 to 2008. J
Atmos Electr 34:79–86. https://​doi.​org/​10.​1541/​jae.​34.​79
JMA (2009) Technical information on distribution materials (Meteorological
Edition) No. 296 (Haishin shiryou ni kansuru gizyutsu jouhou (kisyou hen)dai 296 gou) (in Japanese). https://​www.​data.​jma.​go.​jp/​suish​in/​jyouh​ou/​
pdf/​296.​pdf Accessed 7 Feb 2023
JMA (2022) Outline of the operational numerical weather prediction at thejapan meteorological agency. https://​www.​jma.​go.​jp/​jma/​jma-​eng/​jma-​
center/​nwp/​outli​ne2022-​nwp/​index.​htm Accessed on 2023/02/09
Kasahara S (2010a) Thunder nowcast (Kaminari naukyasuto) (in Japanese).
Textb Forecast Technol Train 2009:116–135
Kasahara S (2010b) Launch of thunder nowcast (Kaminari naukyasuto no
teikyou kaishi) (in Japanese). Tenki 57:847–852
Kasahara S (2011) Lightning analysis and prediction techniques for Thunder
Nowcast and how to use it (Kaminari naukyasuto ni okeru kaminari no
kaiseki ･ yosoku gizyutu to riyou houhou) (in Japanese). Weather Service
Bulletin 78:3 https://​www.​jma.​go.​jp/​jma/​kishou/​books/​sokkou/​78/​vol78​
p095.​pdf accessed 22 Dec 2022
Kawano T, Suzuki K, Kawamura R (2018) Lightning properties in line-shaped
precipitation systems inducing the 2017 Northern Kyushu rainfall event(Heisei 29 nen 7 gatsu gouu wo motarashita senjyoukousuitai no hatsuraitokusei) (in Japanese). Abstract A356 presented at the MSJ Autumn
Meeting 2018, Sendai, Japan, 20 Oct-1 Nov 2019. https://​www.​metsoc.​
jp/​defau​lt/​wp-​conte​nt/​uploa​ds/​2019/​05/​2018_​MSJ_​Autumn_​Meeti​ng_​
180915.​pdf Accessed 22 Dec 2022
Kawano T, Nakao T, Kawamura R, Suzuki K, Sugimoto S, Takahashi T (2019)
Numerical simulation of Hokuriku winter lightning (Hokuriku toukiraino suuchi simyure-syon) (in Japanese). Abstract A305 presented at theMSJ Spring Meeting 2019, Tokyo, Japan, 15–18 May 2019. https://​www.​metsoc.​jp/​defau​lt/​wp-​conte​nt/​uploa​ds/​2019/​11/​2019_​MSJ_​Spring_​
Meeti​ng_​LowRes_​v2.​pdf Accessed 22 Dec 2022
Kawano T, Kawamura R (2020) Genesis and maintenance processes of a quasistationary convective band that produced record-breaking precipitation

Availability of data and materials
The source code of the SCALE library is available from the SCALE website(https://​scale.​riken.​jp/). The data are available upon request from the corresponding author (YS).

Declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interest.
Received: 27 February 2023 Accepted: 25 September 2023

References
Altaratz O, Reisin T, Levin Z (2005) Simulation of the electrification of winterthunderclouds using the three-dimensional Regional Atmospheric
Modeling System (RAMS) model: Single cloud simulations. J Geophys Res
110:D20205. https://​doi.​org/​10.​1029/​2004J​D0056​16
Barthe C, Molinié G, Pinty J-P (2005) Description and first results of an explicit
electrical scheme in a 3D cloud resolving model. Atmos Res 76:95–113.
https://​doi.​org/​10.​1016/j.​atmos​res.​2004.​11.​021
Barthe C, Chong M, Pinty J-P, Bovalo C, Escobar J (2012) CELLS v1.0: updated
and parallelized version of an electrical scheme to simulate multiple
electrified clouds and flashes over large domains. Geosci Model Dev
5:167–184. https://​doi.​org/​10.​5194/​gmd-5-​167-​2012
Barthe C, Hoarau T, Bovalo C (2016) Cloud electrification and lightning activity
in a tropical cyclone-like vortex. Atmos Res 180:297–309. https://​doi.​org/​
10.​1016/j.​atmos​res.​2016.​05.​023
Basarab BM, Rutledge SA, Fuchs BR (2015) An improved lightning flash rate
parameterization developed from Colorado DC3 thunderstorm data foruse in cloud-resolving chemical transport models. J Geophys Res Atmos
120:9481–9499. https://​doi.​org/​10.​1002/​2015J​D0234​70
Beljaars ACM, Holtslag AAM (1991) Flux parameterization over land surfaces
for atmospheric models. J Appl Meteor 30:327–341. https://​doi.​org/​10.​
1175/​1520-​0450(1991)​030%​3c0327:​FPOLSF%​3e2.0.​CO;2
Bovalo C, Barthe C, Pinty J (2019) Examining relationships between cloudresolving model parameters and total flash rates to generate lightning
density maps. Q J R Meteorol Soc 145:1250–1266. https://​doi.​org/​10.​
1002/​qj.​3494
Chen Z, Qie X, Yair Y, Liu D, Xiao X, Wang D, Yuan S (2020) Electrical evolutionof a rapidly developing MCS during its vigorous vertical growth phase.
Atmos Res 246:105201. https://​doi.​org/​10.​1016/j.​atmos​res.​2020.​105201
Clark AJ, Gallus WA Jr, Weisman ML (2010) Neighborhood-based verificationof precipitation forecasts from convection-allowing NCAR WRF Model
simulations and the operational NAM. Wea Forecasting 25:1495–1509.
https://​doi.​org/​10.​1175/​2010W​AF222​2404.1
Curran EB, Holle RL, López RE (2000) Lightning casualties and damages in the
United States from 1959 to 1994. J Clim 13:3448–3464. https://​doi.​org/​10.​
1175/​1520-​0442(2000)​013%​3c3448:​LCADIT%​3e2.0.​CO;2
Dafis S, Fierro A, Giannaros TM, Kotroni V, Lagouvardos K, Mansell E (2018)
Performance evaluation of an explicit lightning forecasting system. J Geophys Res Atmos 123:5130–5148. https://​doi.​org/​10.​1029/​2017J​D0279​30
Fierro AO, Mansell ER, Macgorman DR, Ziegler CL (2013) The implementationof an explicit charging and discharge lightning scheme within the WRFARW model: benchmark simulations of a continental squall line, a tropicalTomioka et al. Progress in Earth and Planetary Science(2023) 10:60

in Northern Kyushu, Japan on 5 July 2017. J Meteorol Soc Japan Series II
98:673–690. https://​doi.​org/​10.​2151/​jmsj.​2020-​033
Kondo M, Sato Y, Inatsu M, Katsuyama Y (2021) Evaluation of cloud microphysical schemes for winter snowfall events in Hokkaido: a case study ofsnowfall by winter monsoon. SOLA 17:74–80. https://​doi.​org/​10.​2151/​
sola.​2021-​012
Kusaka H, Kondo H, Kikegawa Y, Kimura F (2001) A simple single-layer urban
canopy model for atmospheric models: comparison with multi-layerand slab models. Bound-Layer Meteorol 101:329–358. https://​doi.​org/​10.​1023/A:​10192​07923​078
Lopez P (2018) Promising results for lightning predictions. ECMWF Newslett
155:14–19. https://​doi.​org/​10.​21957/​plz73​1tyg2
Lu J, Qie X, Xiao X, Jiang R, Mansell ER, Fierro AO (2022) Effects of convective mergers on the evolution of microphysical and electrical activity in a severe squallline simulated by WRF coupled with explicit electrification scheme. J Geophys
Res Atmos 127:e2021JD036398. https://​doi.​org/​10.​1029/​2021J​D0363​98
Lynn BH, Kelman G, Ellrod G (2015) An evaluation of the efficacy of using observed
lightning to improve convective lightning forecasts. Wea Forecast 30:405–423.
https://​doi.​org/​10.​1175/​WAF-D-​13-​00028.1
MacGorman DR, Straka JM, Ziegler CL (2001) A lightning parameterization for
numerical cloud models. J Appl Meteorol 40:459–478. https://​doi.​org/​10.​
1175/​1520-​0450(2001)​040%​3c0459:​ALPFNC%​3e2.0.​CO;2
Mansell ER, Ziegler CL (2013) Aerosol effects on simulated storm electrificationand precipitation in a two-moment bulk microphysics model. J Atmos Sci
70:2032–2050. https://​doi.​org/​10.​1175/​JAS-D-​12-​0264.1
Mansell ER, MacGorman DR, Ziegler CL, Straka JM (2002) Simulated threedimensional branched lightning in a numerical thunderstorm model. J
Geophys Res Atmos 107:ACL2-1-ACL2-12. https://​doi.​org/​10.​1029/​2000J​
D0002​44
Mansell ER, MacGorman DR, Ziegler CL, Straka JM (2005) Charge structure andlightning sensitivity in a simulated multicell thunderstorm. J Geophys Res
110:D12101. https://​doi.​org/​10.​1029/​2004J​D0052​87
McCaul EW, Goodman SJ, LaCasse KM, Cecil DJ (2009) Forecasting lightning threat using cloud-resolving model simulations. Wea Forecast
24:709–729. https://​doi.​org/​10.​1175/​2008W​AF222​2152.1
McCaul EW, Priftis G, Case JL, Chronis T, Gatlin PN, Goodman SJ, Kong F (2020)
Sensitivities of the WRF lightning forecasting algorithm to parameterized
microphysics and boundary layer schemes. Wea Forecast 35:1545–1560.
https://​doi.​org/​10.​1175/​WAF-D-​19-​0101.1
Michimoto K (1993) A study of radar echoes and their relation to lightning
discharges of thunderclouds in the Hokuriku District, II, Observation
and analysis of “single-flash” thunderclouds in midwinter. J Meteorol Soc
Japan 71:195–204. https://​doi.​org/​10.​2151/​jmsj1​965.​71.2_​195
Montanyà J, Fabró F, van der Velde O, March V, Williams ER, Pineda N, Romero
D, Solà G, Freijo M (2016) Global distribution of winter lightning: a threat
to wind turbines and aircraft. Nat Hazards Earth Syst Sci 16:1465–1472.
https://​doi.​org/​10.​5194/​nhess-​16-​1465-​2016
Nakakita E, Osakada Y, Yamaguchi K (2019) The characteristics of the heavy
rainfall event in July 2018 and global warming (in Japanese). DPRI Annuals 62:1–5
Nakanishi M, Niino H (2006) An improved Mellor-Yamada level-3 model: its
numerical stability and application to a regional prediction of advection fog. Bound-Layer Meteorol 119:397–407. https://​doi.​org/​10.​1007/​
s10546-​005-​9030-8
Nishizawa S, Yashiro H, Sato Y, Miyamoto Y, Tomita H (2015) Influence of grid
aspect ratio on planetary boundary layer turbulence in large-eddy
simulations. Geosci Model Dev 8:3393–3419. https://​doi.​org/​10.​5194/​
gmd-8-​3393-​2015
Pflaum JC, Pruppacher HR (1979) A wind tunnel investigation of the growth ofgraupel initiated from frozen drops. J Atmos Sci 36:680–689
Rogers RR, Yau MK (1989) A short course in cloud physics. Int Ser Nat Philos
113:290
Roh W, Satoh M (2014) Evaluation of precipitating hydrometeor parameterizations in a single-moment bulk microphysics scheme for deep convective systems over the tropical central pacific. J Atmos Sci 71:2654–2673.
https://​doi.​org/​10.​1175/​JAS-D-​13-​0252.1
Saad Y (2003) Iterative methods for sparse linear systems, 2nd edn. SIAM,
Philadelphia
Sato Y, Nishizawa S, Yashiro H, Miyamoto Y, Kajikawa Y, Tomita H (2015) Impactsof cloud microphysics on trade wind cumulus: which cloud microphysics

Page 22 of 22

processes contribute to the diversity in a large eddy simulation? Prog
Earth Planet Sci. https://​doi.​org/​10.​1186/​s40645-​015-​0053-6
Sato Y, Miyamoto Y, Tomita H (2019) Large dependency of charge distribution
in a tropical cyclone inner core upon aerosol number concentration. Prog
Earth Planet Sci 6:62. https://​doi.​org/​10.​1186/​s40645-​019-​0309-7
Sato Y, Miyamoto Y, Tomita H (2021) Lightning frequency in an idealized hurricane-like vortex from initial to steady-state using a coupled meteorological and explicit bulk lightning model. Mon Weather Rev 149:753–771.
https://​doi.​org/​10.​1175/​MWR-D-​20-​0110.1
Sato Y, Hayashi S, Hashimoto A (2022) Difference in the lightning frequency
between the July 2018 heavy rainfall event over central Japan and the2017 northern Kyushu heavy rainfall event in Japan. Atmos Sci Lett
23:e1067. https://​doi.​org/​10.​1002/​asl.​1067
Seiki T, Nakajima T (2014) Aerosol effects of the condensation process on aconvective cloud simulation. J Atmos Sci 71:833–853. https://​doi.​org/​10.​
1175/​JAS-D-​12-​0195.1
Sekiguchi M, Nakajima T (2008) A k-distribution-based radiation code andits computational optimization for an atmospheric general circulation
model. J Quant Spectrosc Radiat Transf 109:2779–2793. https://​doi.​org/​
10.​1016/j.​jqsrt.​2008.​07.​013
Tajiri T, Morimoto T, Nakamura Y, Sakai H, Takayanagi Y, Shimizu M (2021) A
study on progress process of lightning discharge with long lightning
channels (Tyoukyori sinten suru kaminari houden no houden katei nokaiseki) (in Japanese). J Atmos Electr 40:32–36. https://​doi.​org/​10.​1541/​
jae.​40.​32
Takahashi T (1978) Riming electrification as a charge generation mechanism inthunderstorms. J Atmos Sci 35:1536–1548. https://​doi.​org/​10.​1175/​1520-​
0469(1978)​035%​3c1536:​REAACG%​3e2.0.​CO;2
Takahashi T (1984) Thunderstorm electrification—a numerical study. J Atmos
Sci 41:2541–2558. https://​doi.​org/​10.​1175/​1520-​0469(1984)​041%​3c2541:​
TENS%​3e2.0.​CO;2
Thomas RJ, Krehbiel PR, Rison W, Hamlin T, Harlin J, Shown D (2001) Observations of VHF source powers radiated by lightning. Geophys Res Lett
28:143–146. https://​doi.​org/​10.​1029/​2000G​L0114​64
Tomita H (2008) New microphysical schemes with five and six categories bydiagnostic generation of cloud ice. J Meteorol Soc Japan 86A:121–142.
https://​doi.​org/​10.​2151/​jmsj.​86A.​121
Tsuchida (2018) Lightning guidance (in Japanese). JMA Numer Predict Div Rep
64:165–176
Tsuji H, Yokoyama C, Takayabu YN (2020) Contrasting features of the July
2018 heavy rainfall event and the 2017 Northern Kyushu rainfall event inJapan. J Meteorol Soc Japan Series II 98:859–876. https://​doi.​org/​10.​2151/​
jmsj.​2020-​045
van der Vorst HA (1992) Bi-CGSTAB: a fast and smoothly converging variantof BiCG for the solution of nonsymmetric linear systems. SIAM J Sci Stat
Comput 13:631–644. https://​doi.​org/​10.​1137/​09130​35
Wu T, Ushio T, Yoshida S (2014) LF lightning location system (BOLT) in Osaka
area and its applications. In: international conference on lightning protection (ICLP) pp 82–84. https://​doi.​org/​10.​1109/​ICLP.​2014.​69730​99.
Yokoyama S (2002) Facts and research issues on lightning damage problems
in the advanced information society (in Japanese). IEEJ Technical Report
902:13–18
Yoshida S, Wu T, Ushio T, Kusunoki K, Nakamura Y (2014) Initial results of LF
sensor network for lightning observation and characteristics of lightningemission in LF band. J Geophys Res Atmos 119:12,034-12,051. https://​doi.​
org/​10.​1002/​2014J​D0220​65
Yoshida S, Yoshikawa E, Adachi T, Kusunoki K, Hayashi S, Inoue H (2019) Threedimensional radio images of winter lightning in Japan and characteristics
of associated charge structure. IEEJ Trans Electr Electron Eng 14:175–184.
https://​doi.​org/​10.​1002/​tee.​22795

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.