The Propagation Effects of Lightning Electromagnetic Fields Over Mountainous Terrain in the Earth-Ionosphere Waveguide.

In this paper, a full-wave two-dimensional Finite-Difference-Time-Domain model is developed to evaluate the propagation effects of lightning electromagnetic fields over mountainous terrain in the Earth-ionosphere waveguide. In the model, we investigate the effect of the Earth-ionosphere waveguide structure and medium parameters, including the effect of the ionospheric cold plasma characteristics, the effect of the Earth curvature, and the propagation effects over mountainous terrain. For the first time, the obtained results are validated against simultaneous experimental data consisting of lightning currents measured at the Säntis Tower and electric fields measured in Neudorf, Austria, located at 380-km distance from the tower. It is shown that both the time delays and amplitudes of the lightning electromagnetic fields at 380-km distance can be strongly affected by the ionospheric electron density profile, the mountainous terrain, and the Earth curvature. After taking into account the effect of the irregular terrain between the Säntis Tower and the field measurement station, the vertical electric fields calculated by using our model are found to be in good agreement with the corresponding measured cases occurred in both daytime and nighttime. The ideal approximation used in either the classical solutions or the simplified models might lead to inaccuracies in the estimated reflection height. Furthermore, we discuss the sensitivity of our results by considering different return stroke models, as well as different typical values of the return stroke speed and of the ground conductivity. ©2019. The Authors.

Introduction

Very low frequency (VLF) waves generated by lightning discharges propagate in the Earth‐ionosphere waveguide (EIWG) through a process of multiple reflections occurring between the lower region of the ionosphere (60–90 km) and the Earth ground surface. Since this altitude range is too low for satellites and too high for balloons, ground‐based recorded VLF signals that involve the information of the ionosphere and its variability are well known as an efficient tool to probe the localized variation of the ionospheric region parameters (Inan et al., 2010).

In the past 30 years, a number of models and methods have been developed to study the lightning discharge interaction with the lower ionosphere. Different full‐wave theoretical approaches have been used, such as the wave‐hop (ray theory) method (Jacobson et al., 2009, 2018; Qin et al., 2017; Shao and Jacobson, 2009), the waveguide mode theory (Budden, 1961; Cheng and Cummer, 2005; Cummer et al., 1998), or numerical methods such as the Finite‐Difference‐Time‐Domain (FDTD) method (Azadifar et al., 2017; Han et al., 2011; Hu and Cummer, 2006; Marshall, 2012; Tran et al., 2017) and the full‐wave Finite Element Method (Lehtinen and Inan, 2008; Lehtinen and Inan, 2009). Although empirical models such as the International Reference Ionosphere (IRI) model (Bilitza and Reinisch, 2008) and the semiempirical FIRI model based on radio wave propagation data from rocket soundings (Friedrich and Torkar, 2001) could provide the specification of the ionosphere parameters at a given time and location, to simplify the calculation, most VLF remote sensing studies make use of the two‐parameter profile introduced by Wait and Spies (1964) consisting of the reference height and the steepness to infer the ionospheric parameters. The ionospheric parameters could change from day to night at different locations even during a magnetically quite periods (Wait and Spies, 1964; Wait and Walters, 1963).

Recent studies show that the amplitude and phase perturbation for VLF signals have a complicated relationship with the ionospheric region parameters. The propagation of the VLF signal between the Earth ground surface and the lower ionosphere can be affected by many factors, such as the propagation distances (Aoki et al., 2015; Azadifar et al., 2017; Tran et al., 2017), the attenuation due to the ground conductivity (Aoki et al., 2015; Tran et al., 2017), the electron and neutral particle densities (Lay and Shao, 2011; Lay et al., 2014; Marshall, 2012), the Earth curvature (Tran et al., 2018) and the presence of the geomagnetic field (Lay et al., 2014; Shao et al., 2013).

However, for simplicity and computational efficiency, most models and methods have been simplified. Some studies ignore the effect of the Earth curvature (Aoki et al., 2015; Tran et al., 2017); others consider the effect of the Earth curvature based on the correction algorithm and full‐wave ray theory over the spherical Earth (Hu and Cummer, 2006; Jacobson et al., 2018, 2009; Qin et al., 2017; Shao and Jacobson, 2009) or use a staircase approximation of the curved surface of the Earth (Tran et al., 2018). Most studies considered the effect of the ground conductivity by assuming a homogeneous finitely conducting smooth ground (Aoki et al., 2015; Hu and Cummer, 2006; Qin et al., 2017; Shao et al., 2013; Tran et al., 2018). More recently, the effect of the propagation of lightning electromagnetic fields over a rough surface has been investigated in many studies (Cooray and Ming, 1994; Ming and Cooray, 1994; Li et al., 2016a, 2016b, 2017, 2014; Liu et al., 2016; Schulz and Diendorfer, 2000; Zhang, Jing, et al., 2012; Zhang, Yang, et al., 2012). It was noted that the waveshape, peak value, and time delay of the lightning‐radiated fields could be strongly affected when the propagation proceeds along a mountainous terrain or a nonflat ground. The assumption of a finitely conducting smooth ground might result in a significant underestimation of the peak of the electric fields in mountainous areas (Li et al., 2017, 2016b). Moreover, previous studies presenting comparison between simulations and experimental data are based on relative values because the incident sources were unknown. In many cases, the lightning current waveform is calculated based on analytical expressions, and the lightning channel is treated unrealistically as a vertical electric dipole (Bérenger, 2005; Hu and Cummer, 2006; Qin et al., 2017; Schmitter, 2014; Shao and Jacobson, 2009), or only one of the transmission‐line‐based models is used to discuss the influences of the lightning source characteristics (Aoki et al., 2015; Azadifar et al., 2017; Tran et al., 2017). To the best of our knowledge, none of the previous studies has considered the propagation effect of lightning electromagnetic fields over mountainous terrain along the Earth's curved surface based on the simultaneously measured channel‐base lightning current waveform and its radiated electromagnetic fields. The reason for this is that, as discussed by Thomson (2010), the measurements at intermediate locations were not generally practicable due to factors such as the presence of the sea, lack of ready access via roads, or low, uncertain, rapidly varying, ground conductivity.

Therefore, to better investigate the propagation effect of lightning‐radiated electromagnetic fields in the EIWG, a full‐wave two‐dimensional (2‐D) spherical Finite‐Difference‐Time‐Domain (FDTD) model is developed in this paper including the effect of the ionospheric cold plasma characteristics, the effect of the Earth curvature, and the propagation effects over a mountainous terrain. The simulation results are validated against simultaneous experimental data obtained from the Säntis Tower in Switzerland and the 380‐km distance electric field station at Neudorf in Northern Austria. Furthermore, we discuss the sensitivity of the obtained results by considering different return stroke models, as well as different typical values of the return stroke speed and of the ground conductivity.

Experimental Data

Simultaneous experimental data were obtained from the Säntis Tower in Switzerland and from the 380‐km distance electric field station at Neudorf in Northern Austria. The 124‐m‐tall Säntis Tower is located on the top of Mount Säntis in Switzerland. It has been instrumented since 2010 to measure the lightning current waveforms and their time derivatives at two different heights (24 and 82 m) along the Tower. More details on the measurement sensors and instrumentation system can be found in (Romero et al., 2012; Romero et al., 2013; Azadifar et al., 2014). Figure 1 shows the locations and pictures for the Säntis Tower (red triangle) and the electric field station in Neudorf (red dot).

Figure 1

The topographic map of the selected region including the locations of the Säntis Tower (red triangle) and the 380‐km Electric field station in Neudorf (red dot).

The data used in this study were recorded on May 7, 2014. They consist of 10 flashes that occurred during the period from 9:00 to 23:00 local time (LT) in Switzerland. Table 1 summarizes the information on all the events captured on May 7, 2014: occurrence time, peak of the largest return strokes, number of return strokes, and number of negative and positive return strokes. Note that, throughout the paper, we only focus on the highest return stroke current in each flash. As shown in Table 1, all the peak values of the largest return stroke currents are negative, except the one corresponding to the flash that occurred at 20:42:26 LT. Unfortunately, the electric field corresponding to this particular positive return stroke current is too noisy to identify it in the background signal at the 380‐km station. Therefore, in this study, we only focus on the negative return strokes. Figures 2a and 2c give, respectively, the channel‐base current of the flash that occurred at 09:02:39 LT during daytime and the corresponding 380‐km electric field. Figures 2b and 2d further present the highest stroke current in this flash and the associated 380‐km electric field. Figure 3 is similar to Figure 2 but for the flash that occurred at 23:26:59 LT during nighttime. It should be noted that our 380‐km field measuring system exhibited a high noise level at frequencies of 200, 300, 400, and 500 kHz. In order to reduce the noise, 20‐kHz bandwidth notch filters centered at each of the abovementioned frequencies were applied to all recorded signals (see Figure 5 in Azadifar et al., 2017 for more information). In addition, a 60‐kHz ringing can be seen in all the measured field waveforms. The origin of this ringing is unknown and currently under investigation. It might be due to a malfunction of the integrator at the measuring station or electromagnetic noise coupling to the digitizer card caused by a bad PC power supply. In order to reduce the effect of the 60‐kHz ringing noise, throughout this study, we applied a third‐order Butterworth filter that stops the frequencies between 40 and 200 kHz based on the MATLAB digital Butterworth Bandstop filter toolbox. We selected this filter since the ringing noise in our case being relatively well separated in frequency from the signal without ringing noise. In Appendix A, we discuss the difference between the used filter and the other available filters (such as Bessel filter, Gauss filter, and a simple brick‐wall filter; see Appendix A for details). In order to assess the effect of the used filter, we applied it to the FDTD simulated field waveform (see Appendix A). The results of the comparison between the original and filtered waveforms support the fact that the ringing noise in our case is sufficiently separated in frequency from the received signal, and it can be effectively removed after applying the filter.

Table 1

Occurrence Time (LT), Peak of the Highest Return Stroke Current (kA), Total Number of Return Strokes (Total_RS), and Number of Negative (Neg_RS) and Positive Return Strokes (Pos_RS) for All the Flashes Captured on May 7, 2014

Time (LT)	Peak of the largest RS current (kA)	Total_ RS	Neg_RS	Pos_RS
2014/5/7 09:02:39a	‐5.5	10	10	0
2014/5/7, 09:08:55	‐5.0	5	5	0
2014/5/7, 18:11:48	‐9.5	12	12	0
2014/5/7, 20:42:26	33.6	19	16	3
2014/5/7, 20:46:00	‐7.2	7	6	1
2014/5/7, 21:26:05	‐6.4	15	15	0
2014/5/7, 21:28:25	‐8.4	18	18	0
2014/5/7, 21:29:56	‐9.0	8	7	1
2014/5/7, 21:32:02	‐14.8	13	13	0
2014/5/7, 23:26:59b	‐13.7	17	17	0

Corresponding to Figure 2.

Corresponding to Figure 3.

Figure 2

Simultaneously recorded experimental data associated with the lightning flash occurred on May 7, 2014 at 09:02:39 LT. The flash current (a) and its associated field waveform (c) and the largest stroke current in this flash (b) and its associated field waveform (d) (measured and filtered field data are shown in red and black lines, respectively).

Figure 3

Similar to Figure 2 but for the simultaneous experimental data associated with the lightning flash occurred on May 7, 2014 at 23:26:59 LT (measured and filtered field data are shown in red and black lines, respectively).

The final filtered data for field waveforms are shown in black lines in Figures 2c, 2d, 3c, and 3d; the detailed method for the filter used to minimize the effect of the 60‐kHz noise is presented in Appendix A.

FDTD Modeling

In this section, a full‐wave FDTD model including the effect of the ionospheric reflections and the propagation effects of a mountainous terrain is presented. In our simulation, the measured return stroke currents obtained from the Säntis Tower were directly used as input, and the electric field waveforms at 380‐km distance from the Säntis Tower were calculated to compare with the corresponding measured electric field data at the Neudorf station. As shown in Figure 4, the FDTD modeling 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). Note that the model itself naturally includes the effect of the Earth curvature by considering the irregular terrain along the Earth's curved surface. The lightning return stroke channel is modeled using transmission‐line‐type models assumed as a vertical phased current‐source array (Baba and Rakov, 2003) in our FDTD model. Mount Säntis with the Tower on the top and the Neudorf field measurement point are also shown in Figure 4.

Figure 4

The geometry of the computational domain of the FDTD model.

Ionospheric Model

As mentioned in section 1, lightning VLF wave propagation between the Earth's ground surface and the lower ionospheric region can be affected by the electron and neutral particle densities and the geomagnetic field (Lay and Shao, 2011; Lay et al., 2014; Marshall, 2012). The FDTD model used in this study solves the fundamental equations as suggested by Marshall (2012) including the Maxwell's equations and Langevin equation. The effect of the ionospheric parameters is included in the Langevin equation as shown in Equation (1). In the simulation, only the effects of electrons are considered since ion mobility can be neglected in the lower region ionosphere (Han and Cummer, 2010). where 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, is the collision frequency between electrons and the neutral air, and , , and are the charge, mass, and number density of electrons, respectively. The electron mobility is considered as a function of the local electric field and the background gas number density based on the profile proposed by Pasko et al. (1997). According to Marshall (2012) and Liu et al. (2017), the propagation of the lightning electromagnetic field in the lower ionosphere is mainly determined by the interaction of electrons and neutrals; the geomagnetic field can be neglected when the electron‐neutral collision frequency is much higher than the electron gyrofrequency. In our case, considering a geomagnetic field of about 50,000 nT based on the U.S./U.K. World Magnetic Model 2019, the calculated electron gyrofrequency  rad/s is much lower than the electron‐neutral collision frequency, namely,  rad/s at 60 km altitude and  rad/s at 70 km altitude. Therefore, it is reasonable to neglect the effect of the geomagnetic field in our case (the gyrofrequency vector in Equation (1) is assumed to be zero).

Figure 5 shows a comparison of the electron density profile between the IRI 2007 model (Bilitza and Reinisch, 2008) and the two‐parameter profile introduced by Wait and Spies (1964) consisting of the reference height and the steepness . It can be seen that the traditional Wait and Spies' exponential profile gives a reasonable approximation for the electron density profile below 90‐km altitude compared to the IRI model. This two‐parameter approach has been widely used in many VLF remote sensing studies to infer the ionospheric parameters, and it has been found to give good agreement between observed and calculated region ionosphere characteristics (Cummer et al., 1998; McRae and Thomson, 2000, 2004; Shao et al., 2013; Thomson, 2010).

Figure 5

Distribution of electron density profile as a function of the altitude taken from the IRI model 2007 (blue line) and the two‐parameter electron density profile of Wait and Spies (1964) (red line) consisting of the reference height  km and the steepness  km .

The neutral density profile of gas in the ionosphere is considered to be composed of 78% N and 22% O according to the MSIS‐E‐90 model (Hedin, 1991), which is valid below 100‐km altitude (Marshall, 2012). Figure 6 gives the distribution of gas number density profile as a function of the altitude by assuming the mentioned 78% N and 22% O atmosphere composition.

Figure 6

Distribution of gas number density profile as a function of the altitude taken from the MSIS‐E‐90 model assumed atmosphere of 78% N (blue line) and 22% O (red line).

In our study, since we only consider altitudes below 90 km, the profile of electron density ( 90 km) is obtained from the two‐parameter electron density equation below (Shao et al., 2013; Volland, 1995) where is a constant and and are the steepness and the reference reflection height of the profile, respectively. The effective reflection height corresponds to the height where most of the VLF energy is reflected.

The Effect of Mountainous Terrain

To take into account the topography between the Säntis Tower and the 380‐km distance sensor, the Global Digital Elevation Model Version 2 data from the Advanced Spaceborne Thermal Emission and Reflection radiometer with a horizontal grid spacing of 1 arc sec (approximately 30 m) were utilized. More details for the ground boundary can be found in (Li et al., 2016b, 2017). Figure 7 gives the top view of the terrain map (panel a) and 2‐D cross section of the topographic profile (panel b) along the direct path between the Säntis Tower and the 380‐km Neudorf distant sensor site (see the white dashed line in Figure 7a). Note that, in our case, the Säntis Tower is located on the highest point (the top of the Säntis Mountain) in the considered region (see Figure 7b), which allows us to take advantage of the mirror and rotational features by using the 2‐D symmetric assumption to reduce the computation time and storage capability required for the simulation. However, the 2‐D symmetric assumption can only be considered as reasonable as long as the terrain profile does not exhibit very intense variations along its path. A detailed analysis of this problem was presented in Li et al. (2017). Note that the effect of the 124‐m‐tall Säntis Tower is neglected in our study since the shortwave propagation round‐trip time along the tower relative to the rise time of the current waveforms does not result in any significant effect on the radiated fields (Baba and Rakov, 2007; Visacro and Silveira, 2005). Since the 380‐km electric field antenna in Neudorf was installed on the roof of a building with an estimated enhancement factor of about 2.6 (Pichler et al., 2010), all the simulated field values were multiplied by the factor 2.6.

Figure 7

The terrain map (a) and 2‐D cross section (b) of the topographic profiles along the direct path between the Säntis Tower and the 380‐km distant Neudorf sensor site.

Modeling Parameters

The working space of the FDTD model is 400 km   110 km, which is divided into 100‐m quadrilateral cells. Considering 10 steps per wavelength in the FDTD simulation to avoid numerical dispersion, the maximum frequency which can be simulated by using our full‐wave FDTD model is about 300 kHz. The current distribution along the lightning channel was specified according to the transmission‐line‐type model with exponential current decay (MTLE) (Nucci, 1988; Rachidi and Nucci, 1990). In the MTLE model, the temporal and spatial variation of the return stroke current is given by , , where is the current at ground level, is the return stroke speed, and is a constant describing the current decay with height, inferred from sets of simultaneously measured electric and magnetic fields at two different distances (Nucci and Rachidi, 1989). Note that was recently found to vary in the range of hundreds meters to about 2 km in two triggered lightning events (Shao et al., 2012). The channel‐base return stroke currents measured on the Säntis Tower were directly imported as into the FDTD model. The lightning channel height was assumed to be  km, and the return stroke speed was set to  m/s. The Convolutional Perfectly Matched Layer (CPML) absorbing boundary (Roden and Gedney, 2000) is used to eliminate reflections in the upper, lower, and right directions of the computational domain. The ground conductivity for the entire path was considered to be homogeneous with a conductivity  S/m and a relative permittivity . The depth of the ground layer was 10 km. The developed FDTD simulation code has been validated against results obtained using the Finite Element Method (Li et al., 2015; Paknahad et al., 2015) and using as reference for both theoretical equations and experimental data (Azadifar et al., 2017; Li et al., 2016b). Note that the influence of different types of return stroke models, different return stroke speeds and ground conductivities will be analyzed in section 6.

Comparison Between the FDTD Simulation and the Experimental Data

In this section, we consider the two sets of simultaneously recorded return stroke current and associated field mentioned above and shown in Figures 2b, 2d, 3b, and 3d, occurred during daytime at 09:02:39 LT and during nighttime at 23:26:59 LT, respectively. Based on the FDTD model described in section 3, we calculated the electric fields at 380 km over the mountainous terrain along the Earth's curved surface by varying the ionospheric parameters, namely, the steepness (0.25–3.0  with step ) and the reference reflection height (60–90 km with step  km) in Equation (1). In order to find the best pair of and values to match the experimentally observed field waveforms, the peak value ratios ( ) and time delays ( ) between the ground wave and the first reflected skywave (see Figures 2 and 3) were calculated for both measured and simulated electric fields. Figure 8 shows the squared difference maps calculated by using the peak value ratios ( ) and time delays ( ) corresponding to different pairs of and based on the method proposed by Shao et al. (2013). The white star markers in Figure 8 represent the local minimum of the squared difference map, which corresponds to the optimal pair of ‐ values. It was found that the best‐matched ‐ values (marked by the white star) are , and , for the lightning return strokes during the daytime (see Figure 2b) and nighttime (see Figure 3b), respectively. These values agree well with previous studies in which the reference height of the ionosphere varied in the range of about 60 to 90 km as a function of local time (Smith et al., 2004). The reference height of the ionosphere in the nighttime is higher than that in the daytime since the solar radiation increases the total ionization in the region (Han et al., 2011; Smith et al., 2004). By using the best‐matched ‐ values, Figure 9 further shows the comparison between the measurement and the FDTD modeling results. It can be seen in the figure, after taking into account the effect of the irregular terrain along the Earth's curved surface, the vertical electric fields calculated by using our model are in generally good agreement with the measurements obtained from the 380‐km sensor in Neudorf for both the daytime and nighttime cases. Discrepancies are observed in the overshooting hump of the skywave in the daytime case. The disagreement is probably due to the inadequacy of the used simplified two‐parameter exponential model for the description of the daytime ionosphere (Han et al., 2011). It is noted that the effect of the consecutive skywaves in our case might be affected by the presence of the small‐angle mountainous profile at close distance.

Figure 8

Squared difference map of between the FDTD simulation results and the obtained simultaneous experimental data associated with the pairs of ‐ values.

Figure 9

Comparison between the measurement and the FDTD modeling results considering the best‐matched ‐ values. (a) Daytime case:  km,  km and (b) nighttime case:  km,  km .

In addition, we compared our results with classical solutions presented by Schonland et al. (1940) and Smith et al. (2004). Both methods are based on the simplified ray theory of reflection within the EIWG considering the effect of the Earth's curvature. However, in order to simplify the calculation for VLF/LF waves, they assumed the lightning source to be a vertical dipole on the ground or at a certain altitude and considered the entire ionosphere to be ideal. The raypath is assumed to be straight, following the rules of geometric optics and ideally getting reflected from the ground at an effective (or virtual) ionospheric reflection height. In our case, the calculated reference reflection height using Schonland's method is 67 km for the daytime case and 81 km for the nighttime case. By assuming that the altitude of the lightning source varied from 0 to 8 km, the calculated reference reflection height using Smith's method is 66–70 km for the daytime case and 81–85 km for the nighttime case. It is found that both simplified ray theory methods overestimate the reference reflection height by a few kilometers in the range from 2 to 8 km, which confirms the conclusion that the effect of medium parameters in the EIWG could have a significant effect on the evaluated ionosphere height.

In order to investigate the effect of the medium parameters in the EIWG, we considered the Earth's curved surface as opposed to a simple flat surface. As mentioned in section 1, the effect of the Earth curvature has been considered either using correction algorithms, using full‐wave ray theory over a spherical Earth (Hu and Cummer, 2006; Jacobson et al., 2018, 2009; Qin et al., 2017; Shao and Jacobson, 2009) or using a staircase approximation of the curved surface of the Earth (Tran et al., 2018). In our study, to provide a more accurate boundary approximation, the 2‐D spherical coordinate FDTD model presented in section 3 was developed, while for the purpose of comparison, a 2‐D cylindrical coordinate FDTD model (Li et al., 2016b) was also used to represent the case of the flat ground surface.

In the analysis, we considered four different cases for both the cylindrical and the spherical coordinate FDTD models: (i) a perfectly conducting smooth ground without electron density profile in the ionosphere, (ii) a perfectly conducting smooth ground with electron density profile in the ionosphere, (iii) a finitely conducting smooth ground (electrical parameters  S/m and ) with electron density profile in the ionosphere, and (iv) a finitely conducting ground considering the mountainous terrain profile (electrical parameters  S/m and ) with electron density profile in the ionosphere. Figure 10 shows the simulation results for the four different cases calculated by using the cylindrical and spherical FDTD models for both the daytime and the nighttime cases.

Figure 10

The simulation results for four different cases calculated by using 2‐D cylindrical (dashed line) and spherical coordinate (solid line) FDTD models for both the daytime (a, May 7, 2014 at 09:02:39 LT) and nighttime (b, May 7, 2014 at 23:26:59 LT) cases.

It is shown that both the groundwaves and ionospheric reflected skywaves of the lightning electromagnetic fields at 380‐km distance can be affected by the presence of the medium parameters in the EIWG. Tables 2 and 3 further summarize the peak value ratios ( ) and time delays ( ) between the ground wave and first reflected skywave of the simulation results considering the four different cases. Note that the and for the measurement data are also shown in the tables. As expected, the parameters and considering the effect of the irregular terrain along the Earth's curved surface and the electron density in the ionosphere show the best consistency with the measurements for both the daytime and nighttime cases. As can be seen, both peak value ratios ( ) and time delays ( ) can be significantly affected by the presence of the electron density in the ionosphere. The results associated with the perfectly conducting smooth ground surface without the electron density profile in the ionosphere always underestimate and compared with the measurements since the evaluated and are not fitted for the ideal cases. It is interesting to note that, after considering the electron density in the ionosphere, for both the cylindrical (without Earth curvature) and spherical (with Earth curvature) coordinates FDTD models, the peak value ratios seem not to be affected so much, but, in contrast, the time delays ( ) changed by a few microseconds due to the effect of the Earth curvature and the mountainous terrain. The Earth curvature can be neglected for short distance propagation with small incident angles (Hu and Cummer, 2006), but in our case, it seems that, even for the case of a perfectly conducting smooth ground, the effect of the Earth curvature is still significant when the observation distance is 380 km from Säntis Tower. The peak values of both the ground wave and the first skywave can be enhanced by the Earth's curvature relative to a flat‐ground case. The enhancement of the groundwave peak is due to the definition of the observation distance in the cylindrical and spherical coordinate systems. In order to make the comparison, we considered the observation distance in the cylindrical and spherical coordinate systems to be the same. However, for the spherical coordinate system, is the propagation distance between the source and the receiver considering a spherical Earth geometry; the direct line distance in the spherical coordinate system is smaller than the propagation distance in the cylindrical coordinate system. Thus, the amplitude of the signal corresponding to the Earth curvature case is enhanced. Additionally, as discussed in Li et al. (2016b), the presence of the tall mountain on the top of which the tower is sitting can also result in a significant enhancement of the ground wave field peak. Note that, for more complex situations, for example, when high mountains are present (higher than the location of the lightning source), the effect of the mountainous terrain can be much stronger at both short and long distances (Li et al., 2017, 2016a).

Table 2

The Peak Value Ratios ( ) and Time Delays ( ) Between the Ground and the First Reflected Skywave for the Four Considered Cases Associated With the Lightning Return Stroke Occurred at 09:02:39 LT on May 7, 2014

	Without Earth curvature		With Earth curvature
Different cases	R	τ(μs)	R	τ(μs)
(i)	0.7	72.1	0.6	76.3
(ii)	1.6	82.8	1.4	90.8
(iii)	1.5	82.2	1.4	86.5
(iv)	1.6	85.5	1.5	89.7
Measurement	R	τ(μs)
	2.3	89.2

Note. (i) A perfectly conducting smooth ground without electron density profile in the ionosphere, (ii) a perfectly conducting smooth ground with electron density profile in the ionosphere, (iii) a finitely conducting smooth ground with electron density profile in the ionosphere, and (iv) a finitely conducting ground considering the mountainous terrain profile with electron density profile in the ionosphere.

Table 3

The Peak Value Ratios ( ) and Time Delays ( ) Between the Ground and the First Reflected Skywave for the Four Considered Cases Associated With the Lightning Return Stroke Occurred at 23:26:59 LT on May 7, 2014

	Without Earth curvature		With Earth curvature
Different cases	R	τ(μs)	R	τ(μs)
(i)	0.7	91.2	0.6	106.6
(ii)	2.1	107.5	1.9	118.4
(iii)	2.0	107.0	1.8	116.3
(iv)	2.0	109.4	1.9	119.2
Measurement	R	τ(μs)
	1.8	120.2

Note. (i) A perfectly conducting smooth ground without electron density profile in the ionosphere, (ii) a perfectly conducting smooth ground with electron density profile in the ionosphere, (iii) a finitely conducting smooth ground with electron density profile in the ionosphere, and (iv) a finitely conducting ground considering the mountainous terrain profile with electron density profile in the ionosphere.

Additionally, as seen in Figure 10, the waveform corresponding to Case (i) perfectly conducting smooth ground without electron density profile in the ionosphere looks very different from the other cases with the peak value of the first reflected skywave bigger than the groundwave and of same polarity as the ground wave. This is due to the fact that the results for Case (i) are obtained by assuming a Perfect Electrical Conductor boundary in the upper ionosphere region considering 65 and 77.5 km for daytime and nighttime, respectively. We included this case that assumes the ionosphere as a perfect conductor because many simplified theoretical approaches, such as the mentioned simplified ray theory method (Smith et al., 2004; Schonland et al., 1940) and the simplified waveguide model (Wait, 2013) assumed a sharply bounded perfectly conducting ionosphere.

The amplitude difference between the groundwave and the first reflected skywave can be significantly affected by the geometry for the reflection of the lightning electromagnetic waves propagating in the EIWG with different incident angles (see further discussion in Krider, 1992 and Shao et al., 2005).

As mentioned in section 5, the classical solutions based on the simplified ray theory of reflection that ignores all the medium parameters in the EIWG have been widely used in many VLF remote sensing studies to infer the ionospheric parameters based on the time delays ( ) between the ground wave and the reflected skywave of the signals. It turns out that both time delays and amplitudes of the lightning‐radiated electromagnetic fields might be affected by the medium parameters in the EIWG, such as the electron density profile, the Earth curvature, and the presence of the mountainous terrain. The ideal approximation in the calculation might therefore lead to inaccuracies in the evaluation of the ionospheric parameters.

Discussion

In this section, we discuss the sensitivity of the obtained results by considering different return stroke models, different typical values of the return stroke speed, and the ground conductivity.

Influence of Return Stroke Models

Engineering models of the lightning return stroke have been widely used for lightning electromagnetic pulse calculations in numerous studies including lightning electromagnetic coupling with power and communication lines (Rachidi, 2012; Rakov and Rachidi, 2009; Rakov and Uman, 1998). For most modeling studies involving the lightning discharge interaction with the region ionosphere, the lightning return stroke channel is always assumed to be a vertical electric dipole, which means the return stroke current waveform has been assumed to be the uniform along the entire lightning channel. However, the engineering models of the lightning channel account for the temporal and spatial variation of the return stroke current. In this section, three different engineering (transmission‐line type) models, (i) the transmission line model (TL) (Uman and McLain, 1969), (ii) the modified transmission‐line model with linear current decay with height (MTLL) (Rakov and Dulzon, 1991), and (iii) the modified transmission line model with exponential current decay with height (MTLE) (Nucci, 1988; Rachidi and Nucci, 1990), are used. In the Transmission Line (TL) model, the temporal and spatial variation of the return stroke current is given by , , where is the current at ground level and is return stroke speed. In the MTLE model, the spatial and temporal variation of the return stroke current is given by , , where  km is a constant describing the current decay with height. According to the MTLL model, the return stroke current at any level along the channel is given by , , where is the height of the return stroke channel.

In order to compare with different return stroke models, we assumed the channel length for all the cases. Figure 11 shows the lightning electric fields at 380‐km distance over mountainous terrain with an electron density profile for different return stroke models. Although the lightning channel itself is much smaller than the propagation distance (380 km), both the amplitude and waveform of the groundwave and skywave can be affected by using different models. The TL model exhibits disagreement with respect to the other transmission‐line‐type models due to the absence of current attenuation along the lightning channel. It can also be seen that, by assuming a dipole model for the lightning source in our case, the electric field peaks are much higher than those using transmission‐line‐type models.

Figure 11

Lightning electric fields at 380‐km distance over mountainous terrain by considering three different transmission‐line‐type models: (i) the transmission line model (TL), (ii) the modified transmission line model with linear current decay with height (MTLL); (iii) the modified transmission line model with exponential current decay with height (MTLE); and (iv) the Dipole model.

According to Table 4, the time delays seem not to be affected so much by the different return stroke models, although they do change slightly. However, the peak value ratios for the dipole model are found to be much bigger than the other transmission‐line‐type models. For the considered cases, the MTLE model is found to have the best agreement with the measurements for both the daytime and the nighttime waveforms.

Table 4

The Peak Value Ratios ( ) and Time Delays ( ) Between the Ground Wave and the First Reflected Skywave of the Simulation Results Obtained With Three Different Transmission‐Line‐Type and Dipole Models for Two Lreturn Strokes Occurred at 09:02:39 LT and 23:26:59 LT on 7 May 2014

	09:02:39 LT		23:26:59 LT
Different return stroke models	R	τ(μs)	R	τ(μs)
TL	1.4	81.8	1.7	117.8
MTLL	1.3	86.6	1.6	116.8
MTLE	1.5	89.7	1.9	119.2
Dipole	2.4	87.0	2.5	123.5
Measured data	2.3	89.2	1.8	120.2

Influence of the Return Stroke Speed

In this section, we consider typical values of the lightning return stroke speed (Rakov and Uman, 2003), namely, , , and  m/s. The considered values for the return stroke speed are in agreement with the range of the observed experimental data (CIGRE, 2014; Shao et al., 2012). The adopted return stroke model is MTLE. As shown in Figure 12, the peak amplitude of the lightning electric field increases by 30% when the return stroke speed increases from to  m/s. The return stroke speed has an impact on the effective lightning channel length (Blaes et al., 2014), leading to different effective current moments in the lightning electromagnetic radiation. Table 5 further gives peak value ratios ( ) and time delays ( ) associated with different return stroke speeds. The time delays change slightly as the return stroke speeds increase, but the peak value ratios remain the same since the relative amplitudes between the ground wave and the first reflected skywave are unchanged for different return stroke speeds.

Figure 12

Lightning electric fields at 380 km over mountainous terrain considering three different return stroke speeds: (i) , (ii) , and (iii)  m/s.

Table 5

The Peak Value Ratios ( ) and Time Delays ( ) Between the Ground Wave and the First Reflected Skywave of the Simulation Results for Different Return Stroke Speeds for Two Lightning Return Strokes Occurred at 09:02:39 LT and 23:26:59 LT on 7 May 2014

	09:02:39 LT		23:26:59 LT
Different return stroke speeds (m/s)	R	τ(μs)	R	τ(μs)
1.0×10^8	1.4	80.8	1.8	119.4
1.5×10^8	1.5	89.7	1.9	119.2
2.0×10^8	1.4	82.0	1.8	120.7
Measured data	2.3	89.2	1.8	120.2

Influence of the Ground Conductivity

The finite ground conductivity affects essentially the early time response of lightning‐radiated electromagnetic fields by attenuating its peak and slowing down its rise time (Cooray, 2009; Cooray et al., 2000; Delfino, Procopio, & Rossi, 2008; Delfino, Procopio, Rossi, Rachidi, & Nucci, 2008). In Figure 13, the lightning electric fields at 380 km are evaluated by considering three different values for the ground conductivity, namely, 0.1, 0.01, and 0.001 S/m. For comparison, the results obtained assuming a perfectly conducting ground ( ) are also shown. The obtained results are very similar for both the ground waveform and sky waveform parts. As the ground conductivity decreases, the rise times of the ground and the first reflection pulses increase, and their peak amplitudes decrease. Table 6 further gives peak value ratios ( ) and time delays ( ) associated with different ground conductivities. It is noted that the ground conductivities have comparatively more effect on the time delays than on the peak value ratios . The time delays change a few microseconds as the ground conductivity decreases, but the peak value ratios are almost unchanged for both the daytime and the nighttime cases.

Figure 13

Lightning electric fields at 380 km over mountainous terrain by considering different ground conductivities: (i) , (ii)  S/m, (iii)  S/m, and (iv)  S/m. The relative permittivity was set to in all cases.

Table 6

The Peak Value Ratios ( ) and Time Delays ( ) Between the Ground Wave and the First Reflected Skywave of the Simulation Results Associated With Different Ground Conductivities for Two Lightning Return Strokes Occurred at 09:02:39 LT and 23:26:59 LT on 7 May 2014

	09:02:39 LT		23:26:59 LT
Different ground conductivities ( σ_g, S/m)	R	τ(μs)	R	τ(μs)
∞	1.3	83.9	1.9	123.0
0.1	1.3	83.9	1.9	123.0
0.01	1.3	83.7	1.9	122.8
0.001	1.5	89.7	1.9	119.2
Measured data	2.3	89.2	1.8	120.2

Summary and Conclusion

In this paper, we developed a full‐wave 2‐D spherical FDTD model to analyze the propagation effects of lightning electromagnetic fields over mountainous terrain in the EIWG. The obtained results were validated against simultaneous experimental data consisting of lightning currents measured at the Säntis Tower and distant electric fields measured at 380 km in Northern Austria. After taking into account the effect of the irregular terrain between the Säntis Tower and the field measurement station, as well as the electron density profile in the ionosphere, the vertical electric fields calculated by using our model are in good agreement with the corresponding measured waveforms, for both daytime and nighttime. It was shown that the time delays and amplitudes of the lightning‐radiated electromagnetic fields can be affected by the medium parameters in the EIWG including the effect of the ionospheric cold plasma characteristics, the effect of the Earth curvature, and the propagation effects over mountainous terrain. Classical solutions based on the simplified ray theory of reflection by using an ideal approximation in the EIWG have been widely used in many VLF remote sensing studies to infer the ionospheric parameters based on the time delays ( ) between the ground wave and the reflected skywave of the signals. It turns out that both the time delays and the amplitudes of the lightning‐radiated electromagnetic fields might be affected by the medium parameters in the EIWG, such as the electron density profile, the Earth curvature, and the presence of the mountainous terrain. The ideal approximation in the calculation might therefore cause some uncertainties in the evaluation of the ionospheric parameters. Furthermore, although the lightning channel is much smaller than the propagation distance (380 km), both the amplitude and the waveform of the groundwave and skywave can be affected by using different models. The MTLE model is found to have the best agreement with the measurements for both the daytime and the nighttime cases. Both ground conductivities and return stroke speeds can affect the waveshape of the lightning‐radiated field. The peak amplitude of the lightning electric field increases by about 30% when the assumed return stroke speed increases from to  m/s. As the ground conductivity decreases, the lightning electric fields are found to be influenced by an increase in the rise time and a decrease of the amplitude.