Finite-difference time-domain analysis of ELF radio wave propagation in the spherical Earth–ionosphere waveguide and its validation based on analytical solutions

The finite-difference time-domain (FDTD) model of electromagnetic wave propagation in the Earth–ionosphere cavity was developed under assumption of an axisymmetric system, solving the reduced Maxwell equations in a 2D spherical coordinate system. The model was validated on different conductivity profiles for the electric and magnetic field components for various locations on Earth along the meridian. The characteristic electric and magnetic altitudes, phase velocity, and attenuation rate were calculated. We compared the results of numerical and analytical calculations and found good agreement between them. The undertaken FDTD modeling enables us to analyze the Schumann resonances and the propagation of individual lightning discharges occurring at various distances from the receiver. The developed model is particularly useful when analyzing ELF measurements.

1 Introduction

The finite-difference time-domain (FDTD) model is a numerical analysis technique based on time-dependent differential Maxwell equations. It was originally developed for Cartesian coordinate system, but after elaboration of the code for spherical coordinates, it found applications in studies of extremely low frequency (ELF) and very low frequency (VLF) radio wave propagation in the Earth–ionosphere waveguide (Holland, 1983; Hayakawa and Otsuyama, 2002; Simpson and Taflove, 2002; Otsuyama et al., 2003; Yang and Pasko, 2005; Yu et al., 2012; Samimi et al., 2015). Similar analysis can be performed for Mars and other planets (Soriano et al., 2007; Navarro et al., 2007; Yang et al., 2006; Navarro et al., 2008).

When a small part of the Earth–ionosphere cavity needs to be analyzed, a local volume can be divided into FDTD grid in 2D cylindrical (Cummer, 2000; Hu and Cummer, 2006; Marshall, 2012; Qin et al., 2019), or 3D Cartesian coordinate systems (Araki et al., 2018; Suzuki et al., 2016). It facilitates taking into account some complex inhomogeneities and ionospheric anisotropy in the analysis of ELF–VLF radio waves.

FDTD modeling in 3D Cartesian coordinates system was used for verification of Wait's and Cooray–Rubinstein's analytical formulas, describing lightning-radiated electric and magnetic fields for a mixed propagation path (vertically stratified conductivity), and for the fractal rough ground surface (Zhang et al., 2012a, b; Li et al., 2013, 2014).

The analysis of propagation over a mountainous terrain was performed in a 2D axial symmetric model (Li et al., 2016), which was further developed into an FDTD model in a 2D spherical axisymmetric coordinate system (Li et al., 2019). In such modeling the authors have investigated 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's curvature, and the propagation effects over a mountainous terrain.

A full-wave finite element method was used to calculate electromagnetic fields in a horizontally stratified ionosphere treated as a magnetized plasma (Lehtinen and Inan, 2008). The authors have implemented the source currents with arbitrary vertical and horizontal distributions. The electromagnetic field was calculated both in the Earth–ionosphere waveguide and in the ionosphere as an upward propagating whistler mode wave. This model was used to simulate trans-ionospheric propagation of VLF electromagnetic waves from ground-based transmitters up to satellite altitudes (Lehtinen and Inan, 2009).

In a recent paper by Nickolaenko et al. (2021) the propagation of natural wide-band pulsed radio signals in the spherical Earth–ionosphere cavity is numerically simulated. Applying the realistic vertical profiles of ionospheric conductivity, those authors have obtained ELF–VLF propagation parameters by the full wave solution in the form of Riccati differential equation. Using the day and night ionosphere models, the correspondent parameters were calculated for different source–observer distances.

The influence of ionospheric disturbances on the Schumann resonances was analyzed using the 3D FDTD model in Navarro et al. (2008). The authors estimated the role of day–night asymmetry, polar non-uniformities associated to solar proton events, and X-ray bursts.

Few other numerical approaches were used for the estimation of the parameters of the Earth–ionosphere resonator, like the two-dimensional telegraph equation (TDTE) in a two-dimensional spherical transmission line model of the Earth–ionosphere cavity (Kulak et al., 2003; Prácser et al., 2019; Bozóki et al., 2019) and the transmission line matrix (TLM) numerical method in Cartesian and spherical coordinate systems (Morente et al., 2003; Toledo-Redondo et al., 2016).

Besides the Schumann resonances, another important aspect of ELF studies concerns propagation of ELF waves generated by individual lightning discharges, including Q-burst phenomena (Ogawa et al., 1967). Their waveforms observed at large distance from the source are significantly influenced by the dispersive properties of the Earth–ionosphere waveguide and the round-the-world propagation. Inverse solutions in such cases should take into account full non-uniform solutions of the Maxwell equations. This is particularly important for lighting discharges that have a long continuing current phase or are associated with transient luminous events (TLE) when they occur at large distance from the receiver.

The development of an FDTD model was motivated by our ELF systems: the World ELF Radiolocation Array (WERA) (Kulak et al., 2014) and a new European ELF radiolocation system (EERS) (Mlynarczyk et al., 2018), which are operating in the ELF range, as well as for possible stations on Mars.

In this paper, we present an FDTD uniform model in 2D spherical coordinates and introduce a new approach for validation the FDTD models by introducing computation of complex altitudes of the Earth–ionosphere waveguide, which allows us to compare numerical results with two-scale-height analytical solutions. We also infer the resonance frequencies of Earth–ionosphere waveguides for various conductivity profiles using the decomposition method (Kulak et al., 2006).

2 FDTD model

We have created the FDTD model following the ideas which were originally proposed by Holland (1983) and developed in further studies (Hayakawa and Otsuyama, 2002; Otsuyama et al., 2003; Yang and Pasko, 2005; Navarro et al., 2007; Yang et al., 2006; Navarro et al., 2008). We chose a spherical coordinate system to be able to study the wave's propagation in the Earth–ionosphere cavity and analyze the Schumann resonances.

Since in the present work we did not intend to study the azimuthal dependence of propagation parameters, we reduced a 3D system of Maxwell equations to a 2D axisymmetric system. It allowed us to significantly decrease the required computational power and enabled longer simulation times, which leads to a better frequency resolution (up to df = 0.01 Hz).

2.1 Updated equations

In the case of spherical systems, assuming no dependence on the ϕ coordinate, the system of six Maxwell equations in 3D spherical coordinates (r, θ, ϕ) can be reduced to 2D spherical systems (r, θ) with three equations for E_r, E_θ, and H_ϕ field components (Holland, 1983; Inan and Marshall, 2011):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\epsilon }}_{\mathrm{0}}{\displaystyle \frac{\partial E_r}{\partial t}}+\mathit{\sigma }{E}_r+J_r={\displaystyle \frac{\mathrm{1}}{r\sin \mathit{\theta }}}\left[{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(H_{\mathit{\varphi }}\sin \mathit{\theta }\right)\right],\text{(2)}& {\displaystyle }{\mathit{\epsilon }}_{\mathrm{0}}{\displaystyle \frac{\partial E_{\mathit{\theta }}}{\partial t}}+\mathit{\sigma }{E}_{\mathit{\theta }}=-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{H}_{\mathit{\varphi }}\right),\text{(3)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{0}}{\displaystyle \frac{\partial H_{\mathit{\varphi }}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{r}}\left[{\displaystyle \frac{\partial }{\partial r}}\left(r{E}_{\mathit{\theta }}\right)-{\displaystyle \frac{\partial E_r}{\partial \mathit{\theta }}}\right],\end{array}$$

where σ=σ(r) is the conductivity profile of Earth–ionosphere cavity, and ε₀ and μ₀ are the vacuum permittivity and permeability, respectively.

These equations were discretized using central-difference approximations to the space and time partial derivatives. The resulting finite-difference equations are solved in a leapfrog manner: the electric field vector components in a volume of space are solved at a given instant in time; then the magnetic field vector components in the same spatial volume are solved at the next instant in time and the process is repeated over and over again until the desired transient or steady-state electromagnetic field behavior is fully evolved (Inan and Marshall, 2011). So the update equations for E_r, E_θ, and H_ϕ are the following (Holland, 1983):

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{E}_r{\mathrm{|}}_{i,j}^{n+\mathrm{1}/\mathrm{2}}& ={\displaystyle \frac{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t-{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}}{E}_r{\mathrm{|}}_{i,j}^{n-\mathrm{1}/\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{1}}{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}}\left(-J_r{\mathrm{|}}_{i,j}^n+{\displaystyle \frac{\mathrm{1}}{r_{i+\mathrm{1}/\mathrm{2}}\sin {\mathit{\theta }}_j}}\right.\\ & \cdot \left.\left[{\displaystyle \frac{\sin {\mathit{\theta }}_{j+\mathrm{1}/\mathrm{2}}{H}_{\mathit{\varphi }}{\mathrm{|}}_{i,j}^n-\sin {\mathit{\theta }}_{j-\mathrm{1}/\mathrm{2}}{H}_{\mathit{\varphi }}{\mathrm{|}}_{i,j-\mathrm{1}}^n}{\mathrm{\Delta }\mathit{\theta }}}\right]\right),\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}{E}_{\mathit{\theta }}{\mathrm{|}}_{i,j}^{n+\mathrm{1}/\mathrm{2}}& ={\displaystyle \frac{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t-{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}}{E}_{\mathit{\theta }}{\mathrm{|}}_{i,j}^{n-\mathrm{1}/\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{1}}{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)r_i}}\\ & \cdot \left[{\displaystyle \frac{r_{i+\mathrm{1}/\mathrm{2}}{H}_{\mathit{\varphi }}{\mathrm{|}}_{i,j}^n-r_{i-\mathrm{1}/\mathrm{2}}{H}_{\mathit{\varphi }}{\mathrm{|}}_{i-\mathrm{1},j}^n}{\mathrm{\Delta }r}}\right],\end{array}\text{(6)}& {\displaystyle }\begin{array}{rl}{H}_{\mathit{\varphi }}{\mathrm{|}}_{i,j}^{n+\mathrm{1}}& =H_{\mathit{\varphi }}{\mathrm{|}}_{i,j}^n-{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\mu }}_{\mathrm{0}}{r}_{i+\mathrm{1}/\mathrm{2}}}}\\ & \cdot \left[{\displaystyle \frac{r_{i+\mathrm{1}}{E}_{\mathit{\theta }}{\mathrm{|}}_{i+\mathrm{1},j}^{n+\mathrm{1}/\mathrm{2}}-r_i{E}_{\mathit{\theta }}{\mathrm{|}}_{i,j}^{n+\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }r}}\right.\\ & \left.-{\displaystyle \frac{E_r{\mathrm{|}}_{i,j+\mathrm{1}}^{n+\mathrm{1}/\mathrm{2}}-E_r{\mathrm{|}}_{i,j}^{n+\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }\mathit{\theta }}}\right],\end{array}\end{array}$$

where Δt is the time step, Δr and Δθ are the sizes of grid cell in r and θ coordinates, respectively, $r_i=R_{\mathrm{0}}+i\mathrm{\Delta }r$, $r_{i\pm \mathrm{1}/\mathrm{2}}=R+(i\pm \mathrm{1}/\mathrm{2})\mathrm{\Delta }r$, θ_j=jΔθ, ${\mathit{\theta }}_{j\pm \mathrm{1}/\mathrm{2}}=(j\pm \mathrm{1}/\mathrm{2})\mathrm{\Delta }\mathit{\theta }$, and R₀ is the mean of Earth's radius. Superscript n signifies that the quantities are to be evaluated at time t=nΔt, and i,j represent the point (r=iΔr, θ=jΔθ) in the spherical grid. The half time steps indicate that the electric and magnetic fields are calculated alternately.

The updated equation for E_r cannot be applied for poles (where θ=0 and θ=π) because of sin θ in the denominator. To solve this singularity problem the Holland's approach was used (Holland, 1983). Namely, for the poles the integral form of the curl equation with a small contour around the poles was applied, which leads to updated equations for poles without singularities. After that, the updated equation for θ=0^∘ takes the following form (Holland, 1983):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{E}_r{\mathrm{|}}_{i,j}^{n+\mathrm{1}/\mathrm{2}}& ={\displaystyle \frac{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t-{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}}{E}_r{\mathrm{|}}_{i,j}^{n-\mathrm{1}/\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}}}\\ & \cdot \left[-J_r{\mathrm{|}}_{i,j}^n+{\displaystyle \frac{\sin {\mathit{\theta }}_{\mathrm{1}/\mathrm{2}}{H}_{\mathit{\varphi }}{\mathrm{|}}_{i,\mathrm{0}}^n}{r_{i+\mathrm{1}/\mathrm{2}}(\mathrm{1}-\cos {\mathit{\theta }}_{\mathrm{1}/\mathrm{2}})}}\right].\end{array}\end{array}$$

And for θ=180^∘:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{E}_r{\mathrm{|}}_{i,j}^{n+\mathrm{1}/\mathrm{2}}& ={\displaystyle \frac{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t-{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}{\left({\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}\right)}}{E}_r{\mathrm{|}}_{i,j}^{n-\mathrm{0.5}}\\ & +{\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}_{\mathrm{0}}/\mathrm{\Delta }t+{\mathit{\sigma }}_{i,j}/\mathrm{2}}}\\ & \cdot \left[-J_r{\mathrm{|}}_{i,j}^n-{\displaystyle \frac{\sin {\mathit{\theta }}_{N_{\mathit{\theta }}-\mathrm{1}/\mathrm{2}}{H}_{\mathit{\varphi }}{\mathrm{|}}_{i,N_{\mathit{\theta }}-\mathrm{1}}^n}{r_{i+\mathrm{1}/\mathrm{2}}(\mathrm{1}+\cos {\mathit{\theta }}_{N_{\mathit{\theta }}-\mathrm{1}/\mathrm{2}})}}\right],\end{array}\end{array}$$

where ${\mathit{\theta }}_{\mathrm{1}/\mathrm{2}}=\mathrm{\Delta }\mathit{\theta }/\mathrm{2}$, ${\mathit{\theta }}_{N_{\mathit{\theta }}-\mathrm{1}/\mathrm{2}}=(N_{\mathit{\theta }}-\mathrm{1}/\mathrm{2})\mathrm{\Delta }\mathit{\theta }$, and N_θ is the number of grid cells in the θ direction.

2.2 Mesh

The size of grid cell (in the r and θ directions) is defined by the desired range of frequencies, which are going to be analyzed. For the present analysis we used Δr=1 km and Δθ=0.1^∘, which lets us analyze frequencies up to 1 kHz. The ground was modeled as a perfect electric conductor. Also, the upper boundary R_max for the model is a perfect conductor, but it was placed at high enough altitude to make sure that there would be no reflection of the waves from this boundary.

The total simulation time is defined by the desired frequency resolution of FFT ($\mathrm{\Delta }f=\mathrm{1}/t_{max}$ Hz). The time step is defined by the stability of Courant's criterion $\mathrm{\Delta }t\le \mathrm{1}/c\sqrt{(\mathrm{\Delta }r{)}^{-\mathrm{2}}+(a\mathrm{\Delta }\mathit{\theta }{)}^{-\mathrm{2}}}$, where c is the speed of light and a is Earth's radius (Inan and Marshall, 2011). In order to get sufficient frequency resolution (in our case it is 0.01 Hz, to be able to study the differences between different models), we have conducted the FDTD simulation up to 100 s.

2.3 Source

In order to analyze the propagation of electromagnetic waves originating from lightning discharges, one has to use the source model with similar characteristics to real lightning discharge. We used the time profile of lightning discharge proposed in Kulak et al. (2010) and Rakov (2007), which has a flat spectrum in a wide frequency range. But taking into account the restrictions connected with the FDTD computational grid, we had to modify the source: the highest frequency (or smallest wavelength) is defined by the cell size, and according to the generally accepted rule the smallest wavelength should be at least 10–20 times larger than the cell size. Therefore, the source has to be filtered with a loss-pass filter in order to remove frequency components that do not fit the size of the computational cell.

Also, it should be noted that if the source contains the direct current (DC) in its spectrum, it can introduce artifacts, which are not physical (Li et al., 2013). Therefore, additional filtering has to be used for the lowest frequencies by a high-pass filter. For such filtering we have used a fifth-order Butterworth low-pass filter with a cutoff frequency of 300 Hz and a third-order Butterworth high-pass filter with a cutoff frequency of 2 Hz. The final spectrum of this modified source is almost flat in the range 5–100 Hz but decreases rapidly for f<3 and f>1000 Hz.

We placed the source in the FDTD grid at the pole (θ=0^∘) in radial direction in nodes with i = 0–6, which corresponds to the source length L=6 km. Also, for the implementation of the source into the FDTD grid we took into account the lightning stroke speed, and for our analysis we used v=10⁸ m s⁻¹ (Rakov, 2007). We assumed a cloud-to-ground (CG) discharge and implemented it according to the time delay between adjacent nodes in the source.

3 Validation of the model

The main purpose of the current work is the analysis of resonance phenomena in the Earth–ionosphere cavity. In order to validate the developed FDTD model, we analyzed several configurations of spherical waveguide with known analytical solutions and compared them with our FDTD results. We compared the resonance frequencies for them in relative and absolute units.

3.1 Lossless spherical waveguide

The first model for testing was a lossless spherical waveguide with perfect electric conductors at the ground and the upper boundary. The distance between the conductors was set to h=74 km. The precise theoretical resonance frequencies for such a waveguide for a given height are presented by the equation $f_n=(c/\mathrm{2}\mathit{\pi }a)\sqrt{n(n+\mathrm{1})}\sqrt{\mathrm{1}-h/a}$, where c is the speed of light, a is Earth's radius, and n is the mode of resonance frequency (Bliokh et al., 1977). A comparison of analytical and numerical results for such a model is presented in Table 1. The results of the FDTD analysis for this model are shown at Fig. 1.

Table 1 Comparison of resonance frequencies obtained analytically and numerically for a lossless spherical waveguide.

3.2 Spherical waveguide with a conducting layer

A more complicated model for testing had a perfect conductor at the ground and a constant conductivity of $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{6}}$ S m⁻¹ above 70 km. The analytical solutions for this model were obtained by Kulak and Mlynarczyk (2013). A comparison of analytical and numerical results for such a model is presented in Table 2. The results of the FDTD analysis for this model are shown at Fig. 1.

Table 2 Comparison of resonance frequencies obtained analytically and numerically for a spherical waveguide with a conducting layer above 70 km.

3.3 Two-layered spherical waveguide

The third model for testing was a two-layered model with a perfect conductor at the ground, a constant conductivity of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ S m⁻¹ at altitudes between 70 and 110 km, and a constant conductivity of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ S m⁻¹ above 110 km. The analytical results for such a model can be found in Kulak et al. (2013), which deals with multi-layer waveguides. A comparison of analytical and numerical results for this model is presented in Table 3. The results of the FDTD analysis for this model are shown at Fig. 1.

Table 3 Comparison of resonance frequencies obtained analytically and numerically for a two-layered spherical waveguide.

Figure 1 The FDTD spectrum for spherical waveguide with a lossless perfect cavity (a, d), one conducting layer (b, e), and two layers (c, f) for an electric component at 90^∘ and a magnetic component at 180^∘.

4 Application of a realistic atmospheric conductivity profile

The next validation of our FDTD model was done for a realistic atmospheric conductivity profile (Kudintseva et al., 2016; Nickolaenko et al., 2016). We did not take into account the influence of Earth's magnetic field. The resulting anisotropy of the conductivity has little influence in the analyzed frequency range (i.e., 0–100 Hz). As shown in Yu et al. (2012), the anisotropy has a more significant influence at higher frequencies.

Figure 2 A realistic conductivity profile up to 200 km obtained by combining the profile from Kudintseva et al. (2016) (red part) and an IRI profile (blue part).

For the lower part of the atmosphere (0–100 km) we used a conductivity profile recently proposed by Kudintseva et al. (2016). Since the upper boundary of this profile has the conductivity much below 1 S m⁻¹, which is not enough to attenuate the ELF waves, reflections in the FDTD grid would occur from the perfect electric conductor (PEC) at its upper boundary. To avoid reflections and let the waves attenuate gradually at high altitudes, we extended this profile up to 200 km using the IRI (International Reference Ionosphere) model. We chose it in such a way that the combined profile would be smooth (Fig. 2). The required IRI profile was found on 22 January 2006, at 06:00 UT, for the location typical for mid-latitudes (49^∘ N, 23^∘ E). The results of the FDTD analysis for this model are shown at Fig. 3.

Figure 3 The results of the FDTD analysis are shown, which were obtained using the realistic conductivity profile described in Fig. 2. The solutions for the electric component E_r (a, b) and the magnetic component H_ϕ (c, d) are shown in the time domain (a, c) and frequency domain (b, d). Each plot contains the results for two probe positions (60 and 120^∘).

4.1 Characteristic electric and magnetic altitudes

To be able to compare the numerical results with the analytical solutions, we have extracted the complex characteristic electric and magnetic altitudes from the FDTD model. The corresponding altitudes can be extracted using their definition (Mushtak and Williams, 2002; Kirillov, 1993):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\stackrel{\mathrm{‾}}{h}}_{\mathrm{e}}\left(f\right)={\displaystyle \frac{\mathrm{1}}{{\stackrel{\mathrm{‾}}{E}}_r^{\mathrm{0}}}}\sum _{i=\mathrm{0}}^N{\stackrel{\mathrm{‾}}{E}}_r^i\mathrm{\Delta }r,\text{(10)}& {\displaystyle }{\stackrel{\mathrm{‾}}{h}}_{\mathrm{m}}\left(f\right)={\displaystyle \frac{\mathrm{1}}{{\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }}^{\mathrm{0}}}}\sum _{i=\mathrm{0}}^N{\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }}^i\mathrm{\Delta }r,\end{array}$$

where N is the number of radial nodes, ${\stackrel{\mathrm{‾}}{E}}_r^i$ and ${\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }}^i$ are the complex field values at the radial node i for a given frequency, and ${\stackrel{\mathrm{‾}}{E}}_r^{\mathrm{0}}$, ${\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }}^{\mathrm{0}}$ are the values of these fields at the ground level.

Figure 4 Real and imaginary parts of characteristic electric and magnetic altitudes obtained analytically (denoted by “A”) and from the FDTD model (denoted by “N”).

The complex electric altitude can also be expressed using the conductivity profile by the analytical equation in a normalized form (Mushtak and Williams, 2002):

$$\begin{array}{}\text{(11)}& {\stackrel{\mathrm{‾}}{h}}_{\mathrm{e}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}r}{\mathrm{1}-i{\mathit{\sigma }}_{\mathrm{e}}}},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{e}}=\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{e}}^{\ast }$ and the characteristic conductivity value for electric altitude is ${\mathit{\sigma }}_{\mathrm{e}}^{\ast }=\mathit{\omega }{\mathit{\epsilon }}_{\mathrm{0}}$, and ω=2πf.

Assuming a similar dependence for complex magnetic altitude and taking into account that the characteristic conductivity value for magnetic altitude is (Greifinger and Greifinger, 1978)

$${\mathit{\sigma }}_{\mathrm{m}}^{\ast }={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathit{\mu }}_{\mathrm{0}}\mathit{\omega }{\mathit{\zeta }}^{\mathrm{2}}\left(r\right)}},$$

where the scale height $\mathit{\zeta }\left(r\right)=\mathit{\sigma }\left(r\right)/(\mathrm{d}\mathit{\sigma }/\mathrm{d}r)$, we can write a similar equation for the complex magnetic altitude:

$$\begin{array}{}\text{(12)}& {\stackrel{\mathrm{‾}}{h}}_{\mathrm{m}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}r}{\mathrm{1}+i{\mathit{\sigma }}_{\mathrm{m}}}},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{m}}=\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{m}}^{\ast }$. We compared the characteristic electric and magnetic altitudes calculated analytically from Eqs. (11) and (12) with the FDTD results calculated from Eqs. (9) and (10). The results are presented in Fig. 4.

4.2 Propagation parameters

We calculated the phase velocity V_ph and the attenuation rate α using two different methods: numerical and analytical. Analytically the propagation parameters can be calculated using the conductivity profile and electric and magnetic altitudes given by Eqs. (11) and (12), through the following relationship (Kulak and Mlynarczyk, 2011):

$$\begin{array}{}\text{(13)}& {\displaystyle }{V}_{\mathrm{ph}}={\displaystyle \frac{c}{\mathrm{Re}\stackrel{\mathrm{‾}}{S}}},\text{(14)}& {\displaystyle }\mathit{\alpha }={\displaystyle \frac{\mathit{\omega }}{c}}\mathrm{Im}\stackrel{\mathrm{‾}}{S},\end{array}$$

where ${\stackrel{\mathrm{‾}}{S}}^{\mathrm{2}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{m}}/{\stackrel{\mathrm{‾}}{h}}_{\mathrm{e}}$.

In the case of numerical calculations those parameters were obtained by two approaches: (a) similarly to analytical as described above but using electric and magnetic altitudes from FDTD model by Eqs. (9) and (10); and (b) directly from the spectra of electric and magnetic field components. Assuming that in our coordinate system the wave propagates in the θ direction, the following relationship can be written in the frequency domain:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{{\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }\mathrm{1}}}{{\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }\mathrm{2}}}}\sqrt{{\displaystyle \frac{\sin ({\mathit{\rho }}_{\mathrm{1}}/a)}{\sin ({\mathit{\rho }}_{\mathrm{2}}/a)}}}=\exp \left(\stackrel{\mathrm{‾}}{\mathit{\gamma }}\right({\mathit{\rho }}_{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{1}}\left)\right),\end{array}$$

where the ratio of magnetic field complex amplitudes ${\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }}\left(f\right)$ are calculated for two probes “1” and “2” which are located at distances ρ₁ and ρ₂ from the source, respectively. In this equation the propagation constant $\stackrel{\mathrm{‾}}{\mathit{\gamma }}=\mathit{\alpha }+i\mathit{\beta }$, where α is the attenuation rate (in units Np m⁻¹) and β is the phase constant, so the phase velocity is $V_{\mathrm{ph}}=\mathit{\omega }/\mathit{\beta }$. For further convenient usage we convert the units of attenuation rate to dB Mm⁻¹. A similar relationship to Eq. (15) but for electric field components ${\stackrel{\mathrm{‾}}{E}}_r\left(f\right)$ can be written as well.

Figure 5 The attenuation rate and phase velocity that was calculated analytically and numerically. Using numerical FDTD calculations those parameters were obtained by two approaches: (a) using electric and magnetic altitudes from FDTD by Eqs. (9) and (10); and (b) directly from the spectra of electric and magnetic field components following Eq. (15), where we used probes located at 80 and 90^∘.

We should note that Eq. (15) for ${\stackrel{\mathrm{‾}}{E}}_r$ and ${\stackrel{\mathrm{‾}}{H}}_{\mathit{\varphi }}$ spectra calculated in a closed cavity gives non-monotonic behavior of propagation parameters. This is caused by superposition of wave attenuation along meridian and different amplitudes of the Schumann resonances for different source–observer distances. One possible solution for removing the influence of Schumann resonances is to implement a perfectly matched layer (PML) along the θ direction. However, we solved this problem in a different way, transforming the Maxwell equations at angle θ>90^∘ from spherical coordinates to plane equations, changing the behavior of equations from “closed” to “open”. In this case the waves are unable to propagate around Earth and therefore the Schumann resonance does not occur. The calculated propagation parameters from all the methods listed above are presented in Fig. 5.

4.3 Spectral decomposition method

As an additional validation we applied the spectral decomposition technique to Schumann resonance power spectra. Figure 6 presents the resonance frequencies obtained from the spectra measured at different distance from the source, and the decomposed frequencies that do not depend on the source–observer distance. Following the decomposition algorithm (Kulak et al., 2006; Dyrda et al., 2014), the spectrum is approximated by the function

$$\begin{array}{}\text{(16)}& W\left(f\right)=s+{\displaystyle \frac{z}{f^m}}+\sum _{k=\mathrm{1}}^{\mathrm{3}}{\displaystyle \frac{p_k\left[\mathrm{1}+e_k\left(f-f_k\right)\right]}{{\left(f-f_k\right)}^{\mathrm{2}}+{\left(g_k/\mathrm{2}\right)}^{\mathrm{2}}}},\end{array}$$

where W(f) is the signal power spectrum, s is the white noise component, $z/f^m$ is the color noise term, p_k is the power parameter of the kth resonance peak, e_k is the peak asymmetry parameter, f_k is the resonance frequency, and g_k is the resonant mode half-width parameter. Fitting this function to the FDTD spectra allows us to extract the resonance frequencies f_k of the cavity, which are not equal to Schumann resonance frequencies. The Schumann resonance frequencies obtained from the spectra depend on the distance from the source and represent the superposition of standing and traveling waves (Kulak et al., 2006). To analyze the standing waves separately and reveal the properties of the resonant cavity, we use the spectral decomposition method described in Kulak et al. (2006) and Dyrda et al. (2015). After applying this method, the resonance frequencies become independent of the source–observer distance (see details in Kulak et al., 2006 and Dyrda et al., 2015). The decomposition method shows that the solutions for the electric field are symmetric at θ=180^∘ because the traveling waves cancel each other out and only the standing waves remain. Therefore, the resonant peaks measured from the spectrum at θ=180^∘ for the electric field component ($E_r^{\mathrm{180}}$) represent the resonance frequency of the cavity, and they are in agreement with the decomposed frequencies.

Figure 6 The Schumann resonance frequencies of the first three modes obtained at different source–observer distances and the resonance frequencies of the cavity obtained by the decomposition method (Kulak et al., 2006).

4.4 Resonance frequencies of the Earth–ionosphere cavity

The resonance frequencies of the Earth–ionosphere cavity for a given conductivity profile can be calculated by different approaches, and the consistency between them can be considered as an additional validation of the model. We have considered the followings ways to calculate the resonant frequencies:

1. Analytically using the conductivity profile. The resonance frequencies in this approach are calculated by solving the following equation (Mushtak and Williams, 2002; Galejs, 1972):

   $$\begin{array}{}\text{(17)}& f_c^n={\displaystyle \frac{c}{\mathrm{2}\mathit{\pi }a}}\sqrt{n(n+\mathrm{1})}{\displaystyle \frac{\mathrm{Re}\stackrel{\mathrm{‾}}{S}\left(f_c^n\right)}{\mathrm{|}\stackrel{\mathrm{‾}}{S}\left(f_c^n\right)\mathrm{|}}},\end{array}$$

   where ${\stackrel{\mathrm{‾}}{S}}^{\mathrm{2}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{m}}/{\stackrel{\mathrm{‾}}{h}}_{\mathrm{e}}$ depending on the characteristic complex altitudes, which are discussed in Sect. 4.1;

2. Numerically from the FDTD model using Eq. (17). We denote these resonant frequencies by $f_h^n$;

3. Using FDTD spectra for $E_r^{\mathrm{180}}$ (see Sect. 4.3 for details). We denote these frequencies by $f_{E_r}^n$;

4. Using the spectral decomposition method (see Sect. 4.3 for a description). We denote these resonant frequencies by $f_d^n$.

The obtained resonance frequencies are presented in Table 4. We compared these frequencies with the analytical results, and the absolute and relative differences are presented in the corresponding columns.

Table 4 Resonance frequencies for the conductivity profile shown in Fig. 2 calculated by different approaches described in Sect. 4, and the relative error obtained by comparison with the analytically obtained resonance frequency.

5 Discussion

In this study, we used new methods for validation of numerical simulation:

1. We compared the complex electric and magnetic altitudes of the Earth–ionosphere waveguide, referring to two-dimensional formalism of electromagnetic wave propagation in the Earth–ionosphere cavity. The two complex altitudes were calculated numerically directly from their definitions as given in Eqs. (9) and (10), making use of radial field solutions E(r) and H(r). These altitudes were directly compared with the altitudes obtained with analytical formulas as given in Eqs. (11) and (12) for the same conductivity profile of the atmosphere. This allowed us to fully validate the simulation results.

2. We determined the resonance frequencies of the cavity using the decomposition of power spectra described by Eq. (16). The resonance frequencies cannot be directly determined from the spectra obtained using FDTD, because the field generated by the source in each point of the cavity is a superposition of traveling waves propagating directly from the source and the resonance field resulting from the interference of waves propagating around the world. The resulting Schumann resonance frequencies depend on the source–observer distance (Kulak et al., 2006). Close to the source, where the amplitudes of the traveling waves are significant, the resonance frequencies differ by several percent from the Schumann resonance frequencies obtained by FDTD. With the use of the decomposition method we determined the intrinsic resonance frequencies of the cavity, which are the same at each location, independently from the distance to the source. This enabled us to compare the resonance frequencies obtained from the numerical simulation with the frequencies determined directly from analytical Eq. (9). This method should be recommended as a reference for validation of numerical models.

3. We determined the propagation parameters of the Earth–ionosphere waveguide using the FDTD results in two different points of the great circle. Equation (15) allowed us to determine the phase velocity and the attenuation rate of the Earth–ionosphere cavity. We compared them with the results obtained from analytical Eqs. (13) and (14).

6 Summary and conclusions

In this paper, we analyzed the solutions of Maxwell's equations obtained by the FDTD method for an axisymmetric uniform Earth–ionosphere cavity. We analyzed the propagation of radio waves generated by a short current impulse.

We constructed an FDTD model in axisymmetric spherical coordinate system with the source implemented at the pole. We took into account the finite speed of lightning discharge and implemented the time delay between adjacent nodes in the source.

We validated the model thoroughly, comparing the resonance frequencies, propagation parameters, and electric and magnetic characteristic altitudes. Since the conductivity profile of the atmosphere has a significant influence on radio wave propagation and resonance frequencies, we validated our model for various conductivity profiles.

We paid close attention to the verification of accuracy of the FDTD computations and used a new approach for that purpose. It is based on 2D formalism of wave propagation in the Earth–ionosphere cavity and it allowed us to compare the numerical and the analytical solutions. In this approach, the propagation parameters of the Earth–ionosphere waveguide are defined using the electric and magnetic altitudes. These altitudes can be calculated directly when the vertical conductivity profile of the atmosphere is known.

The obtained analytical solutions were used as reference and compared with the numerical results. First, we compared the solutions for three cases: a perfect cavity, a cavity formed by a perfect ground and a homogenous conducing layer, and a cavity formed by a perfect ground and two-layered upper boundary of the waveguide. As a measure of error between the models, we took the difference between the first three resonance frequencies. The analytical and numerical solutions were in agreement (Tables 1–3). Next, we compared the results for a continuous conductivity profile. We built a realistic conductivity model of which the lower part was based on a recently proposed conductivity profile and the upper part was based on an IRI model. We obtained good agreement between the resonance frequencies of the cavity and the observed Schumann resonance frequencies (Table 4).

Using our FDTD model, we also calculated the spectral dependence of the phase velocity and the attenuation rate. We showed that the analytical and numerical models are in agreement.

The presented model can be used for studying the propagation of ELF electromagnetic waves generated by lighting discharges of various types with the round-the-world propagation taken into account.