On the Accuracy of Ray-Theory Methods to Determine the Altitudes of Intracloud Electric Discharges and Ionospheric Reflections: Application to Narrow Bipolar Events.

Narrow bipolar events (NBEs) (also called narrow bipolar pulses [NBPs] or compact intracloud discharges [CIDs]) are energetic intracloud discharges characterized by narrow bipolar electromagnetic waveforms identified from ground-based very low frequency (VLF)/low-frequency (LF) observations. The simplified ray-theory method proposed by Smith et al. (1999, https://doi.org/10.1029/1998JD200045; 2004, https://doi.org/10.1029/2002RS002790) is widely used to infer the altitude of intracloud lightning and the effective (or virtual) reflection height of the ionosphere from VLF/LF signals. However, due to the large amount of high-frequency components in NBEs, the propagation effect of the electromagnetic fields for NBEs at large distance depends nontrivially on the geometry and the effective conductivity of the Earth-ionosphere waveguide (EIWG). In this study, we investigate the propagation of NBEs by using a full-wave Finite-Difference Time-Domain (FDTD) approach. The simulated results are compared with ground-based measurements at different distances in Southern China, and we assess the accuracy of the simplified ray-theory method in estimating the altitude of the NBE source and the effective reflection height of the ionosphere. It is noted that the evaluated NBE altitudes have a slight difference of about ±1 km when compared with the full-wave FDTD results, while the evaluated ionospheric reflection heights are found to be bigger than those obtained from FDTD model by about 5 km. ©2020. The Authors.

Introduction

Narrow bipolar events (NBEs), also called narrow bipolar pulses (NBPs) or compact intracloud discharges (CIDs), are intracloud electrical discharges characterized by narrow bipolar electromagnetic waveforms in ground‐based very low frequency (VLF)/low‐frequency (LF) observations. Several properties of NBEs stand out when compared to other lightning processes: They exhibit short time durations (10–20 m/s), fast propagation speeds (107–108 m/s), and strong very high frequency (VHF) radiation, and they typically occur isolated from other electrical activities within a time frame of a few milliseconds, although in many cases NBEs are precursors of a leader process inside the thunderstorm (e.g., Gurevich et al., 2003; Le Vine, 1980; Smith et al., 1999, 2002; Willett et al., 1989). These unique properties have turned NBEs into an active research topic that promises to illuminate the initiation mechanism of lightning in a thunderstorm (Rison et al., 2016). Furthermore, the physical mechanism of NBEs is as yet poorly understood (Arabshahi et al., 2014; Gurevich et al., 2002).

Due to the great difficulty in directly measuring the NBE discharge channel and the intracloud current inside the thundercloud, most studies use the identified groundwave and their reflected skywaves in a ground‐based VLF/LF signal to determine the altitude of the intracloud NBE source and the effective (or virtual) reflection height in the D region of the ionosphere. Smith et al. (1999) expanded the previous techniques to propose a simplified ray‐theory method that considers an elevated intracloud source at a finite altitude. Later, Smith et al. (2004) improved the method by accounting for the Earth's curvature. They found that the evaluated ionosphere heights agree well with the expected diurnal variation from 60 to 95 km and with the evaluated NBEs source heights ranging from 7 to 20 km above ground level for positive narrow bipolar events (PNBEs) and for negative narrow bipolar events (NNBEs) occurring above 15 km. This simplified ray‐theory method has been applied to many ground‐based VLF/LF studies about the altitudes of the PNBEs/NNBEs and the variation of the ionospheric heights in different regions around the world (e.g., Karunarathne et al., 2015; Liu, Zhu et al., 2018; Liu, Zilong et al., 2018; Wu et al., 2012, 2014; Zhang et al., 2016; Zhu et al., 2010). More recently, the altitude of Energetic In‐cloud Pulses (EIPs), which are linked to Terrestrial Gamma‐ray Flashes (TGFs), has also been evaluated based on this simplified ray‐theory method (Lyu et al., 2015, 2016).

Other approaches based on full‐wave methods such as wave‐hop (ray‐theory) method, waveguide mode theory, and the finite‐difference time‐domain (FDTD) method have also been used to determine the ionospheric height as revealed by ground or intracloud lightning discharges in both daytime and nighttime (e.g., Azadifar et al., 2017; Cummer et al., 1998; Han et al., 2011; Han & Cummer, 2010a, 2010b; Jacobson et al., 2009; Qin et al., 2017; Shao & Jacobson, 2009; Shao et al., 2013; Tran et al., 2017). The propagation effect of NBE electromagnetic pulses or the radiation fields generated by CIDs over finitely conducting ground is mentioned by Rison et al. (2016) and Watson and Marshall (2007), and further investigated by Cooray et al. (2018), but their studies were limited to the groundwave part of the NBE signal. The results show that, due to a large amount of high‐frequency components for the narrow features of NBEs, the peak amplitude, the zero‐crossing time, the full width at half maximum (FWHM), and the time derivatives of the groundwave of the NBEs were significantly distorted by propagation effects over the finitely conducting ground. However, the propagation effect of the NBEs at far distance in the Earth‐Ionosphere WaveGuide (EIWG) will not only depend on the ground conductivity, but it is also a complex function of the geometry and the effective conductivity in the Earth‐ionosphere cavity (Smith et al., 1999).

In this paper, we investigate the propagation effect of NBEs by using a full‐wave FDTD approach. The simulated results are compared with the ground‐based VLF/LF measurements of NBEs at different distances. Moreover, we assess the accuracy of the simplified ray‐theory model in evaluating the ionospheric height and the altitude of the NBE source due to the propagation effects of NBEs in the EIWG.

Theory and Method

In this section we present a full‐wave FDTD model including the propagation effect of the ionospheric reflections and the finitely conducting ground. Note that our FDTD model is very memory and CPU‐time intensive due to the large amount of high‐frequency components stemming from the narrow features of NBEs. We faced three important challenges: (1) the spatial discretization must be fine enough in comparison with the wavelength (usually 10–20 grid cells per wavelength) in order to make the numerical dispersion error negligible; (2) the time step must be small enough to satisfy the FDTD stability condition; and (3) we are interested in distant fields measured up to several hundred kilometers away from the source. Even though we used the open parallel message passing interface (MPI) technique (Guiffaut & Mahdjoubi, 2001) to improve the FDTD computation efficiency, we were unable to run three‐dimensional (3‐D) FDTD simulations to solve this large‐scale problem. Luckily, the whole FDTD domain is symmetric around the z axis, which allows us to reduce our problem to a two‐dimensional (2‐D) one by means of symmetry, which vastly reduces the computation and storage required for the simulation. As shown in Figure 1, the geometry of the full‐wave FDTD model is based on 2‐D spherical coordinates with the sector region defined by the center of the Earth (Bérenger, 2002; Thèvenot et al., 1999).

Figure 1

Geometry of the full‐wave FDTD model.

The FDTD model used in this study solves three fundamental equations as suggested by Marshall (2012) including two Maxwell's curl equations and a modified Ohm equation (sometimes called Langevin equation in this context). In the simulation, only electrons are mobile, since ion motion can be neglected in the lower D region ionosphere (Han & Cummer, 2010a) within our time scale of interest. Our system of equations is where J is the current density of the source, is the self‐consistent conduction current driven by the electric field, is the gyrofrequency vector associated with the geomagnetic field vector , is the plasma frequency, ν =e/(μ m ) is the collision frequency between electrons and the neutral air, and e, m , and N are, respectively, the charge, mass, and number density of electrons. The electron mobility μ is considered as a function of the local electric field and the background gas number density N based on the profile proposed by Pasko et al. (1997). The background gas number density N is taken from the MSIS‐E‐90 model (Hedin, 1991) by assuming a mixed air with 78% N2 and 22% O2.

Because we are limited to 2‐D symmetrical conditions, the effect of the geomagnetic field is neglected. This is justified according to (Marshall et al., 2010; Marshall, 2012), which claim that the effect of the geomagnetic field can be neglected when its dip angle lies between 0° and 45° in the lower ionosphere, but it can be important for high geomagnetic latitudes. In our case the geomagnetic fields in the campaign region have a dip angle lower than 45°. Note that solar disturbances may influence the geomagnetic field (Han & Cummer, 2010a) but we did not consider that possibility. With this approximation, equation (3) reduces to

Following Wait and Spies (1964), our electron density has an exponential height profile given by where β is the steepness of the profile and = 1e9 m−3 is the electron density where most VLF energy can be reflected at the effective reflected height H , which can be evaluated based on the plasma frequency equation and the detected bandwidth of the lightning signal.

FDTD Modeling Parameters

In the FDTD simulation, the computational domain is 500 km×110 km, with spatial steps of 50 m in both r and z directions. Considering 10 steps per wavelength to avoid the numerical dispersion, the maximum frequency that can be simulated by using our full‐wave FDTD model is about 600 kHz, which is enough to cover the VLF/LF frequency band. As shown in Figure 1, we assume a vertical, straight NBE channel above ground and consider the current distribution along the channel according to the modified transmission line model with a linear decay of current with height (MTLL) (Rakov & Dulzon, 1991). The parameters of the NBE current used in this study are based on the VHF observation results from Rison et al. (2016) as shown in Table 1. The NBE current pulse starts at an altitude H 2 and propagates downward at a speed υ a distance L to reach the final, lower altitude H 1. At its uppermost point, the NBE current is assumed to follow a double‐exponential waveform proposed by Rison et al. (2016): where α=1/τ 1 and β=1/τ 2 are the rise and fall time constants of the current and the peak current is normalized to one by setting . Figure 2 shows the distribution of the normalized current at different altitude along the NBE channel.

Table 1

The NBE Parameters Used for the FDTD Model

	τ 1	τ 2	L	h	υ	σ g	ϵ r
Parameters	(μs)	(μs)	(m)	(km)	(m/s)	(S/m)	—
Values	1	6	300	10	3.5e7	0.01	10

Figure 2

The distribution of the normalized current at different altitude along the NBE channel.

The observation point p is located at ground level at a distance d from the NBE source. The ground is assumed to be homogeneous with conductivity σ and relative permittivity ϵ . The ionosphere is considered to be vertically stratified, that is, horizontally homogeneous in the r direction. To avoid artificial reflections at the upper, lower, and outer boundaries of the computational domain, we used absorbing convolutional perfectly matched layers (CPML) (Roden & Gedney, 2000).

Observation Data

The measured data used in this study were recorded by the Jianghuai Area Sferic Array (JASA) in China. The JASA was established in 2009 to investigate lightning‐induced effects in the near‐Earth space. It consists of more than 10 LF electric field sensors that continuously record VLF/LF broadband (between 300 Hz and 300 kHz) lightning signals at a sampling rate of 5 MHz with a detection range of up to several thousand kilometers. By using a GPS time receiver with timing accuracy better than 0.1 m/s for each station, the location of an NBE source is obtained by means of the time of arrival (ToA) technique for at least three stations. The location error in the considered area is estimated to better than 2 km (Qin et al., 2015). More details about the measured sensors and instrumentation system can be found in Zhu et al. (2015).

Results

During the 2012 campaign, 1,539 NBE events were recorded by JASA including 1,395 positive NBEs and 144 negative NBEs from 4 July to 7 July. Further details about the statistical analysis of all NBE events can be found in Liu et al. (2016). In this study, we focus on one typical positive NBE event captured by JASA at local time 08:57:29 on 7 July 2012. Figure 3 shows a comparison between the simulated waveforms from our full‐wave FDTD model and the VLF/LF waveforms for the considered positive NBE event measured at different distances. The typical VLF/LF waveform of NBEs is clearly appreciated in Figure 3, with a groundwave and two ionospheric reflected sky waves marked as G, 1S, and 2S. Notice that, by using the appropriate physical parameters for the Earth and ionosphere, all the simulated results obtained from the full‐wave FDTD model agree perfectly with the VLF/LF measurements at distances ranging from about 100 km to about 400 km, thus validating our full‐wave FDTD results. The bipolar shape of the electric field waveform is due to the lightning radiation field that is proportional to the time derivative of the current (Uman et al., 1975). Since our measured electric field signal is uncalibrated and the currents of NBEs are not available, the electric fields are normalized to the maximum amplitude measured at different sensors from JASA.

Figure 3

Comparison of normalized electric fields between the simulated FDTD results and the VLF/LF waveforms measured at different sensors from JASA: (a) d = 169 km, (b) d = 245 km, (c) d = 340 km, and (d) d = 413 km. The effective reflected height H and the steepness β of the profile are 66 km and 0.5 km−1, respectively. The ground conductivity σ = 0.01 S/m and relative permittivity ϵ =10.

After comparing our FDTD model to the observational results at different distances, let us now compare the FDTD results with the method proposed by Smith et al. (1999, 2004) (henceforth referred to as Smith method) to assess the accuracy of this simplified ray‐theory method. As mentioned above, this method has been widely used to evaluate the altitude of the intracloud source and the effective (or virtual) reflection height of the D region ionosphere. In the Smith method, the relative times of arrival of the simplified signal reflections between the ground and ionosphere are used to determine both the altitude of the NBE source (h) and the effective (or virtual) reflection height of the ionosphere (H) (see the appendix for the details of the Smith method). As shown in Figure 3a there are several possibilities to calculate the time delays between the groundwave and the first and second sky waves. By and , we denote the time delay between the positive peaks of both groundwave and reflected skywaves; these are the time delays used in many studies based on the method proposed by Smith et al. (2004). In order to minimize the effect of the noises in the measured data, we applied a moving average filter with a window length of 30 data points to find the negative and positive peaks to determine the time delays. Table 2 lists the altitude of the NBE source (h) and the effective (or virtual) reflection height of the ionosphere (H) evaluated by Smith method based on the time delays and at different distances as in Figure 3. The table also includes the effective reflected height H and the altitude of the NBE source h in the FDTD model. The effective (or virtual) reflection height of the ionosphere (H) calculated by using Smith method overestimates the effective reflected height in the FDTD simulation (H ) by about 5 km. The Smith method also overestimates the altitude of NBE source h by about 1 km. This agrees with the comparison results from the VHF radio emission of NBEs observed by the FORTE satellite (Smith et al., 2004). Note that the ground conductivity is assumed to be homogeneous in the model. However, the ground conductivity might be inhomogeneous along the propagation paths to different stations. Figure 4 further shows the influence of the ground conductivity on the electric field at different distances. We see that the finite ground conductivity affects essentially the early time response of the groundwave part of the electric fields by attenuating its peak and slowing down its risetime. Meanwhile, the effect of both groundwave and skywave increases with longer distances. Note that the electric fields for each observation distance are normalized to the maximum amplitude of the case corresponding to perfectly conducting ground (σ =∞).

Table 2

Altitudes of the NBE Source (h) and the Effective (or Virtual) Reflection Heights of the Ionosphere (H) Evaluated by Smith Method Based on the Different Time Delays of the Measured NBE Waveforms at Different Distances

	( Tpp1, Tpp2)		( Tpn1, Tpn2)		( Tnp1, Tnp2)
Distance (km)	H (km)	h (km)	H (km)	h (km)	H (km)	h (km)
169	71.75	10.35	67.05	10.93	66.55	9.93
245	71.87	8.94	67.61	11.26	67.55	10.53
340	71.98	9.62	67.64	10.49	67.50	10.84
413	70.70	9.51	66.51	10.78	67.28	11.02
FDTD	H ref (km)	h (km)	H ref (km)	h (km)	H ref (km)	h (km)
—	66	10	66	10	66	10

Figure 4

The simulated electric fields by considering different ground conductivities: (i) σ =∞, (ii) σ = 0.01 S/m, and (iii) σ = 0.001 S/m at different observation distances: (a) d = 169 km, (b) d = 245 km, (c) d = 340 km, and (d) d = 413 km. The relative permittivity was set to ϵ =10 in all cases. The effective reflected height H and the steepness β of the profile are 66 km and 0.5 km−1, respectively.

Figure 5 further compares the full‐wave FDTD results with the geometric ray path of Smith method (see equations (A1) and (A2) in Appendix A) by using cross‐sectional views of the normalized electric fields at different time steps. In the Smith method, the signal is assumed to follow the rule of geometric optics by propagating in straight lines and reflecting ideally on the ground and on an effective (or virtual) ionospheric reflection height (H) (Smith et al., 1999). We see in Figure 5 that these assumptions yield the maximum height reached without significant attenuation by the highest‐frequency component. However, the FDTD results show that the full‐wave electromagnetic pulses from an NBEs contain different frequency components, which are reflected at different heights in the ionosphere. Actually, from Figure 3, we also see that the reflected skywaves are wider and weaker than the groundwaves. Due to dispersion and attenuation in the ionosphere, the higher‐frequency components penetrate deeper than the lower ones and are then rapidly attenuated by the higher conductivity.

Figure 5

Cross‐sectional views of the normalized electric fields for FDTD model associated with the calculated geometric raypath of Smith method at different time steps (a) 0.35002 ms, (b) 0.42002 ms, (c) 0.52503 ms, and (d) 0.6417 ms. The effective (or virtual) reflection height of the ionosphere (H) is calculated by using Smith method, and the effective reflected height (H ) is obtained from the full‐wave FDTD model.

Furthermore, Figure 3 shows that waveform of the skywave undergoes a polarity reversal comparing to the groundwave. This is because different frequency components in the VLF/LF broadband (300 Hz to 300 kHz) are reflected at different altitudes due to the dependency of the plasma frequency with the electron density. It is therefore more reasonable to measure time delays between peaks of opposite polarity. We propose two other time delay definitions as shown in Figure 3a: and are the time delays between the negative peak of groundwave and the positive peaks of the reflected skywaves; and are the time delays between the positive peak of groundwave and the negative peaks of the reflected skywaves. Table 2 further lists the altitude of the NBE source (h) and the effective (or virtual) reflection height of the ionosphere (H) evaluated by Smith method based on the time delays / and / of the simulated NBE waveforms at different distances. Note that in our case by using the new time delays, both altitudes match the full‐wave FDTD results well within ±2 km. Comparing the full‐wave FDTD model with the Smith method that considers the new time delays in reverse polarity improves the accuracy of estimated effective (or virtual) reflection height of the ionosphere.

Summary and Conclusion

We have investigated the propagation effect of NBEs in the EIWG by using a full‐wave FDTD approach. The simulated FDTD results were compared with the ground‐based VLF/LF measurements of NBEs from different sensors in Southern China. With the appropriate parameters for the EIWG, the full‐wave FDTD results agree perfectly with measurements at distances ranging from about 100 km to about 400 km. Moreover, we assessed the accuracy of the simplified ray‐theory model when evaluating the altitude of NBE sources and the ionospheric reflection height. We also proposed different definitions of the time delay between the groundwave and the reflected skywaves in Smith method. We found that considering the new time delays in reverse polarity improves the accuracy of the effective (or virtual) reflection heights of the ionosphere (H); however, by applying the old time delay in the same polarity for Smith method, the evaluated ionospheric reflection heights are found to be larger than the effective reflected heights (H ) obtained from FDTD model about 5 km. For the altitude of NBE source h, the values from the Smith method have a slight difference of about ± 1 km when compared with the full‐wave FDTD results, which is in agreement with the previous results calculated from the VHF radio emission of NBEs observed by FORTE satellite (Smith et al., 2004).