Numerical simulation of the propagation of electromagnetic waves in ionospheric irregularities

1.   Introduction

The characteristics of HF radio wave propagation in the ionosphere can be affected by free electrons. Ionospheric irregularities can have a serious impact on the propagation of electromagnetic (EM) waves. Therefore, research on the propagation of radio waves in ionospheric irregularities is of current importance for space weather and short-wave communication. Spread F (ionospheric irregularities) was first reported by Booker and Wells (1938) at nighttime in the equatorial regions. Dungey (1956) then suggested Rayleigh–Taylor instability to explain the formation of ionospheric irregularities (spread F or plasma bubbles). Later, numerous researchers (Basu and Kelley, 1979; Fejer and Kelley, 1980; Ossakow, 1981; Abdu, 2001; Krall et al., 2013) proposed the generalized Rayleigh–Taylor (GRT) instability, which includes the terms of the electric field and winds apart from gravitation, to explain nighttime ionospheric irregularities, or spread F, in the equatorial and low-latitude regions.

A ray-tracing method is a useful tool for investigating the propagation of radio waves in the ionosphere. Many studies (Paul et al., 1968; Mathur and Pandey, 1977; Huang XQ and Reinisch, 2006; Cervera and Harris, 2014; Sokolov et al., 2016; Psiaki, 2019) have used ray-tracing methods to carry out studies on the propagation of radio waves in the ionosphere. Cervera and Harris (2014) summarized the ray-tracing methods that had been applied in modeling studies of ionospheric disturbances. In terms of the propagation of radio waves in ionospheric irregularities, Paul et al. (1968) carried out a ray-tracing synthesis of ionograms for a severely perturbed ionosphere induced by the release of sulfur hexafluoride (SF6). Mathur and Pandey (1977) likewise used a ray-tracing method to study the effects of ionospheric irregularities, as given by an artificial equation, on HF and very high frequency radio waves.

In this study, we conducted a numerical simulation of the propagation of HF radio waves in ionospheric irregularities. In addition, the modeling results were compared with experimental ionograms recorded by an ionosonde. The method is described in Section 2, the model of ionospheric irregularities is introduced in Section 3, and the results and discussion are presented in Section 4.

2.   Methodologies

It is well known that the ordinary traces of ionograms can be synthesized by applying Equation (1) when the magnetic field and collision are not considered and there is no horizontal gradient in the electron density:

where ${h_0}$ is the base height of the ionosphere, ${h_r}$ is the height of reflection of radio waves in the ionosphere, $N$ is the electron density, and $f$ is the frequency of radio waves.

However, Equation (1) is not suitable for synthesizing ionograms when ionospheric irregularities are present in the ionosphere. Kudeki (2010) proposed a ray-tracing method, represented by Equation (2), to simulate the propagation of radio waves and synthesize ionograms in inhomogeneous media as long as the medium is isotropic:

where $k$ is the wave number, $\omega$ is the wave angular frequency, $n$ is the refractive index, $c$ is the speed of light in a vacuum, $r$ is the wave vector, and $\tau$ is the group path of radio waves. In this study, Equation (2) was used to reconstruct vertical ionograms. For simplicity, we ignored the magnetic field and focused only on the propagation of ordinary radio waves in the ionosphere.

To verify the performance of the present ray-tracing method, we reconstructed an unperturbed electron density profile. We then used Equations (1) and (2), respectively, to synthesize the corresponding ionograms. For the ray-tracing method, the direction of radio rays was set to vertical incidence to synthesize the vertical ionograms. Figure 1 shows the traces of synthesized ionograms by the two methods. The red line indicates the electron density profile (which is represented by quasi-parabolic segments (Croft and Hoogansian, 1968; Dyson and Bennett, 1988): f_oE = 2.5 MHz, h_mE = 105 km, y_bE = 15 km; f_oF₁ = 5 MHz, h_mF₁ = 220 km, y_bF₁ = 105 km; f_oF₂ = 12 MHz, h_mF₂ = 300 km, y_mF₂ = 80 km), the green line is the trace calculated by Equation (1), and the blue dots represent the trace calculated by the ray-tracing method. As shown in Figure 1, the ionogram synthesized by the ray-tracing method agreed well with that calculated by Equation (1). Figure 1 shows that the green line calculated by Equation (1) was not smooth because of the error in numerical computation (Scotto et al., 2012). However, the ray-tracing method was able to overcome the error in numerical computation in Equation (1). Therefore, Figure 1 shows that the ray-tracing method presented in this study performed well in simulating the propagation of EM waves and synthesizing vertical ionograms.

Figure  1.  Ionograms synthesized by Equations (1) and (2). The red line represents the electron density profile by the quasi-parabolic segment model (f_oE = 2.5 MHz, h_mE = 105 km, y_bE = 15 km; f_oF₁ = 5 MHz, h_mF₁ = 220 km, y_bF₁ = 105 km; f_oF₂ = 12 MHz, h_mF₂ = 300 km, y_mF₂ = 80 km), and the green line and blue dots indicate the ionograms synthesized by Equations (1) and (2), respectively.

DownLoad: Full-Size Img PowerPoint

3.   The Model of Ionospheric Irregularities

In this study, a two-dimensional model of plasma bubbles, developed by Jiang CH and Zhao ZY (2019), was used to represent the background ionosphere with ionospheric irregularities. The following equations (Kelley, 2009; Yokoyama, 2017; Jiang CH and Zhao ZY, 2019) were utilized to simulate ionospheric irregularities/plasma bubbles in this study:

where $N$ represents the density of electrons or ions, $N = {N_{\rm{e}}} = {N_{\rm{i}}}$ under the quasi-neutrality condition; $t$ is the running time in the simulation. ${V}$ represents the velocity of electrons or ions; ${\nu _R}$ represents the recombination coefficient; ${N_0}$ is the initial electron density; ${M_{\rm{e}}}$ and ${M_{\rm{i}}}$ represent the mass of electrons and ions, respectively; ${{V}_{\rm{e}}}$ and ${{V}_{\rm{i}}}$ represent, respectively, the velocity of electrons and ions; ${E}$ is the electric field, ${E} = {{E}_0} - \nabla \varphi$; ${{E}_0}$ is the background electric field; $\varphi$ is the electrostatic polarization potential; ${B}$ is the geomagnetic field; ${U}$ is the neutral wind; ${\upsilon _{\rm{en}}}$ and ${\upsilon _{\rm{in}}}$ represent, respectively, the collision frequency of electrons and ions with neutrals; ${g}$ is the gravitational acceleration; $q$ is the unit charge; and ${J}$ is the total current density.

The numerical calculations to be presented were performed on a two-dimensional Cartesian coordinate system (X, Z). The number of grid points in the zonal (Z) and vertical (X) directions was 201. The grid spacing was 2 km in both the X and Z directions. Periodic boundary conditions were imposed in the X direction, whereas in the Z direction, transmissive boundary conditions were imposed. The model domain covered a 250- to 650-km altitude range. The running time step was 0.5 s. In the simulation of ionospheric irregularities, the initial seeding source was the traveling ionospheric disturbances (TIDs), with a time period of 40 min and a horizontal wavelength of 300 km. In the numerical simulation of EM wave propagation, the transmitter position, the elevation angle of radio rays, and the working frequency of radio rays needed to be specified in ray tracing.

4.   Results and Discussion

Figure 2 shows typical radio rays (working frequency at 8 MHz) propagated in the ionosphere with ionospheric irregularities (the elevation angle of radio rays was between 80° and 100° with a step of 1°). The solid lines in Figure 2 represent radio rays propagated in the ionosphere. The transmitter site in Figure 2 was set to the original position (the zonal distance was 0 km and the altitude was 0 km). Figure 2 shows that irregularities in the ionosphere had a strong impact on the propagation of radio rays. Some radio rays were able to pass through the topside of the ionosphere or be reflected from the ionosphere to the ground in the distance. However, some radio rays could be reflected from the ionosphere to the ground near the transmitter. Thus, the smaller the elevation angle of the radio rays, the more radio rays the transmitter and adjacent area could receive. These radio rays could then be used to synthesize the vertical ionograms.

Figure  2.  Typical radio rays (solid lines) propagated in the ionosphere with ionospheric irregularities/plasma bubbles. The frequency of radio rays is 8 MHz, and the angle of elevation is between 80° and 100° with a step of 1°.

DownLoad: Full-Size Img PowerPoint

In a numerical simulation of ionospheric irregularities/plasma bubbles, Zalesak et al. (1982) ignored the chemical recombination rate because the recombination was slight (Kelley, 2009). However, Zalesak and Ossakow (1980) introduced a recombination rate to carry out the simulation of plasma bubbles. Thereafter, Jiang CH and Zhao ZY (2019) introduced a factor between the zero value (where the recombination rate was set to 0) and the specified recombination rate to study the effect of the recombination rate on the simulation of plasma bubbles. Figure 3 shows the evolution of ionospheric irregularities/plasma bubbles caused by TIDs in the ionosphere, where the factor was set to 0.01 in this study. The neutral wind was ignored in this simulation. The first plot in Figure 3 is at t = 1800 s, and the time step was 50 s. As shown in Figure 3, ionospheric irregularities/plasma bubbles began to generate at t = 1900 s. As the time progressed, the ionospheric irregularities/plasma bubbles penetrated into the higher altitudes of the ionosphere (see t = 2350 s in Figure 3). Even though Figure 3 illustrates a two-dimensional simulation of ionospheric irregularities/plasma bubbles, the bifurcation and pinching of ionospheric irregularities/plasma bubbles could be reproduced or induced by two-dimensional TIDs (Jiang et al., 2019).

Figure  3.  Evolution of ionospheric irregularities/plasma bubbles in the numerical simulation.

DownLoad: Full-Size Img PowerPoint

Figure 4 shows the ray-tracing synthesis of ionograms, which corresponds to the ionosphere in Figure 3 with ionospheric irregularities/plasma bubbles. As shown in Figure 4, spread F occurred at t = 2150 s at the bottom of the ionograms. At this stage, this kind of spread F was named a range spread F. When ionospheric irregularities/plasma bubbles developed up to the higher altitudes in the ionosphere, the spread F on ionograms developed into a mixed spread F (see t = 2300 s and t = 2350 s in Figure 4). Figure 3 shows that ionospheric irregularities/plasma bubbles occurred at t = 1950 s at the bottomside of the ionosphere. However, Figure 4 shows that spread F occurred at t = 2150 s. We know that in the ray-tracing method, the next step of radio rays is calculated through the refractive index, which depends on the electron density and frequency of the radio rays. When the frequency of radio rays is higher and the gradient of the electron density is smaller, then plasma bubbles have less impact on the propagation of radio rays. Before t = 2100 s, the gradient of the electron density in the ionosphere was not large enough; thus, spread F did not occur on the ionograms. Figure 4 shows that only when the gradient of the electron density in the ionosphere was large enough would spread F appear on the ionograms.

Figure  4.  Ray-tracing synthesis of ionograms corresponding to Figure 3.

DownLoad: Full-Size Img PowerPoint

Figure 5 shows ionograms with spread F that were recorded at the Puer station (PUR, 22.7°N, 101.05°E) on 19 October 2015. As shown in the figure, satellite traces (in which satellite traces (Tsunoda, 2008) on the ionograms appear as one or more quasi replicas of the main trace of the F layer at displaced virtual heights before the onset of spread F) occurred at 19:55 Local Time (LT) on the ionogram. Spread F then appeared at the bottom of the ionogram at 20:00 LT. After 20:05 LT, spread F on the ionograms developed into a mixed spread F. Compared with Figure 4, the ionosonde observations in Figure 5 were consistent with the numerical simulation of the evolution of spread F on the ionograms. Note that the time resolution differed between the observations (Figure 5) and the numerical simulation (Figure 4).

Figure  5.  Ionograms recorded by the ionosonde at the Puer station (22.7°N, 101.05°E) on 19 October 2015.

DownLoad: Full-Size Img PowerPoint

Tsunoda (2008) suggested that satellite traces could be a precursor to equatorial spread F. However, satellite traces could not be reproduced before the occurrence of spread F in Figure 4. In the present simulation of ionospheric irregularities/plasma bubbles, we ignored the zonal movement of the ionosphere with ionospheric irregularities/plasma bubbles, which can be caused by zonal neutral winds or the vertical electric field. The zonal movement of the ionosphere with ionospheric irregularities/plasma bubbles will be considered in a future study. In that study, it may be possible to reproduce satellite traces before the occurrence of spread F on the ionograms, although it is not necessary for satellite traces to occur before the occurrence of spread F according to ionosonde observations.

5.   Conclusions

In this study, we carried out a numerical simulation of the propagation of HF radio waves in ionospheric irregularities and conducted a ray-tracing synthesis of ionograms. The major findings are summarized below:

(1) A ray-tracing method can be utilized to reproduce spread F on ionograms when ionospheric irregularities occur in the ionosphere.

(2) When ionospheric irregularities/plasma bubbles developed from the bottomside to higher altitudes in the ionosphere, results showed that a range spread F occurred first and then developed into a mixed spread F on the ionograms. Both the simulations and the observations confirmed the evolution of ionospheric irregularities in this study.

(3) Compared with ionosonde observations, the results of the numerical simulation were consistent with the evolution of spread F on the ionograms. However, the satellite traces on the ionograms were not able to be reproduced in this simulation. This will be completed in a future study.

Acknowledgments

We acknowledge the Institute of Earthquake Forecasting for providing the ionosonde data. This work was supported by the National Natural Science Foundation of China (no. 41604133).