Numerical Simulation of Ionospheric Disturbances Due to Rocket Plume and Its Influence on HF Radio Waves Propagation

Abstract

In this paper, the ionospheric disturbances of CO₂, which is released by rocket exhaust plumes, was simulated. The effect of this disturbance on the propagation of high-frequency (HF) radio waves at different incident frequencies was also simulated by using three-dimensional digital ray tracing technique. The results show that CO₂ can effectively dissipate the background electrons and form ionospheric holes after being released in the ionosphere. At the peak height of ionospheric electron density (about 300 km), the electrons are dissipated fastest and the radius of ionospheric hole is also largest. This is due to the fact that the diffusion coefficient of CO₂ usually increases with height while the electron density just increases before reaching its peak height and then decreases with height, and the chemical reaction rate between ions and CO₂ also becomes largest at the peak height of electron density (about 300 km). Around 100 s after the release of CO₂, when the radio waves at a frequency of 8 MHz pass through the ionosphere with an elevation range of 85~95°, the “secondary focusing effect” can occur, and we believe that this is due to the reflection of HF shortwaves on the walls of the ionospheric holes. With time going on, this phenomenon disappears at 300 s and only one focus is left at this time. For the HF shortwaves at same incident frequency, the focusing effect of waves displays a weakening trend with time increasing, and the height of focus center also ascends gradually. At the same time after CO₂ releasing, with the increasing of radio waves frequency, the focusing effect also becomes weaker and the focus center displays an ascending trend.

1. Introduction

The ionosphere is an important link in the solar-terrestrial space environment, in which a large number of plasmas can affect the propagation trajectory of electromagnetic wave, make it reflect, refract or scatter, and cause special effects such as polarization plane rotation. As a natural medium of radio wave propagation, the ionosphere has an important influence on shortwave communication, navigation and positioning, monitoring and early warning and other technologies. With the expansion of human activities into space, in addition to natural disturbance factors such as solar activity and magnetic storms, products including rocket plume inevitably disturb the ionosphere. Rocket plume is a kind of plasma with weak ionization, high temperature and non-magnetization, which is mainly composed of H₂O, CO₂, H₂ and Al₂O₃ and so on. It can effectively dissipate electron density and produce artificial ionospheric holes. When the Vanguard II spacecraft was launched in 1959, Booker first observed the ionospheric holes by using a ground-based ionosonde [1]; this phenomenon was observed repeatedly in subsequent rocket launches. In 1973, the launch of NASA’s Skylab by a Saturn V rocket caused such a large-scale F region depletion (50% depletions spanning several million square kilometers) that a revised theory for ionospheric holes had to be offered. This disturbance lasted about four hours, and the incident interrupted short-wave communication over the vast area of the Atlantic Ocean, and the lowest reflection frequency was reduced to about 2 MHz [2,3]. Mendillo et al. hypothesized and analyzed the reasons for the formation of the ionospheric hole and gave the possible factors for its gradual disappearance. Their studies show that there are a large number of oxygen ions (O+) in the F region of the background ionosphere, and the oxygen ions react violently with the neutral gas molecules in the rocket exhaust plume to form molecular ions, and then immediately reacts with the electrons in the background ionosphere. This leads to the decrease in electron-ion concentration and the formation of huge ionospheric holes. Zinn et al. [4] reported that when the satellite HEAO-C was launched in 1979, the interaction between the rocket plume and the ionospheric plasma resulted in the emergence of a large electron density cavity with a north-south span of 600–800 km and a thickness of 40 km in the ionospheric. It was observed that the total electron content decreased to the same level as at night. Savastano et al. [5] used the real-time TEC observation data obtained from dual-frequency GEO (Geostationary Orbit) satellites to describe the physical characteristics of ionospheric plasma density depletion caused by the 2017 Vandenberg Air Force Base Falcon 9 rocket launch in California.

The frequency range of high-frequency (HF) radio waves is 3–30 MHz and high frequency communication mainly propagates information in the form of sky wave and propagates forward by reflecting back and forth between the ionosphere and the ground. The change of ionospheric state will directly affect the propagation path of HF radio wave in the ionosphere, and then affect the quality of high-frequency communication. Over the horizon radar (OTHR) is an application of high frequency communication. It utilizes the reflection of electromagnetic waves between the ionosphere and the ground to transmit high-frequency energy in order to detect long-range moving targets below the horizon that cannot be detected by conventional radar, and its detection distance is not limited by the curvature of the earth. The ionosphere has the characteristic of irregular electron density, which will affect the propagation path of HF radio wave to a great extent. The ionospheric hole produced by the rocket plume will have a great impact on the propagation path of the over-the-horizon radar [6]. Three-dimensional ray tracing technology is one of the most effective methods to describe and analyze the echo propagation path of over-the-horizon radar [7,8]. In this paper, we simulate the disturbance effect of CO₂ released by rocket exhaust plume in the ionosphere, and then use three-dimensional ray tracing technique to investigate the impact of the plume disturbance on the propagation path of HF radio wave. Moreover, we hope that it can be used as a reference for the study of the influence of rocket plume on the radio wave propagation of over the horizon radar.

2. Kinetic Model

2.1. Neutral Gas Diffusion

When the rocket plume is released to the early stage of the background ionosphere, the neutral gas molecules in the plume will push the surrounding plasma away, similar to a snow blower under the high pressure of the thrusters [9,10]. This process is carried out at supersonic speeds and takes a very short time, usually only a few seconds. When these chemicals are released into the background ionosphere, the pressure difference decreases sharply. When the density of the released chemical substance is close to that of the background gas, the release is fully mixed with the surrounding plasma and diffuses into the surrounding space. This process takes a long time, and the ion chemical reaction mainly takes place at this stage.

Assuming that the background of the ionosphere and thermosphere is horizontally stratified, ignoring the motion of the background atmosphere, the thermal diffusion effect of neutral molecules and the expansion process of chemical substances, and considering the effects of chemical reaction loss and gravity, the selected Cartesian coordinate system takes the neutral gas release point as the center, the x axis points eastward, the y axis points north, and the z axis goes up vertically, then the approximate solution of the diffusion equation of neutral gas under the condition of point source release is [11]:

$$\begin{array}{ll}{n}_i\left(x,y,z,t\right)=& \frac{N_0}{{\left(4\pi D_{0}t\right)}^{3/2}}exp[-\left(z-z_0\right)\left(\frac{3}{4H_a}+\frac{1}{2H_i}\right)-\alpha t\\ & -\frac{H_a^2{\left\{1-exp\left[-\left(z-z_0\right)/\left(2H_a\right)\right]\right\}}^2}{4D_{0}t}-\frac{\left(x^2+y^2\right)exp\left[-\frac{z-z_0}{\left(2H_a\right)}\right]}{4D_{0}t}\\ & -{\left(\frac{1}{H_a}-\frac{1}{H_i}\right)}^2\frac{D_{0}texp\left[\left(z-z_0\right)/\left(2H_a\right)\right]}{4}]\end{array}$$

(1)

where $z_0$ is the height of the release point, $N_0$ is the total amount of neutral gas molecules released, $H_a=kT/\left(m_{a}g\right)$ is the scale height of atmosphere, $H_i=kT/\left(m_{s}g\right)$ is the scale height of the released gas, $k$ is the Boltzmann constant, $T$ is the background gas temperature, $m_a$ is the average molecular weight of the atmosphere, $m_s$ is the molecular weight of the released neutral gas, $g$ is the gravitational acceleration, and $\alpha t$ is the loss term caused by the chemical reaction. $D_0$ is the diffusion coefficient of the released gas. If only O, O₂ and N₂ in the background atmosphere are considered, the expression of the diffusion coefficient of the chemical CO₂ at the release point is [12]:

$$D_{\mathrm{C}{\mathrm{O}}_2}={\left(\frac{N(\mathrm{O})}{5.87\times {10}^{17}{T}^{0.5}}+\frac{N({\mathrm{N}}_2)}{6.58\times {10}^{16}{T}^{0.752}}+\frac{N({\mathrm{O}}_2)}{5.77\times {10}^{16}{T}^{0.749}}\right)}^{-1}{\mathrm{cm}}^2/\mathrm{s}$$

(2)

Considering that when the rocket is flying, the neutral gas released by the rocket exhaust plume cannot be released at a single point only at a fixed altitude; thus, when we analyze the disturbance effect of the rocket plume to the ionosphere, it is necessary to consider the dynamic process of multi-source release of neutral gas, that is, the variation of neutral gas with time and space. In the case of a point source, the density distribution function of the release is $n\left(x,y,z,t\right)$. Assuming that the same amount of neutral gas is released at any release point on the rocket trajectory, the approximate solution of the diffusion equation of neutral gas in the case of multi-source release can be obtained as [13]:

$$n_D\left(x,y,z,t\right)={{\displaystyle \int }}_0^t{{\displaystyle \iiint }}^{ }n\left(x-\mathsf{\alpha },y-\mathsf{\eta },z-\mathsf{\lambda },t-\mathsf{\gamma }\right) \mathrm{d}\mathsf{\alpha }\mathrm{d}\mathsf{\eta }\mathrm{d}\mathsf{\lambda }\mathrm{d}\mathsf{\gamma }$$

(3)

where $n_D$ is the distribution function of release density under multi-source release, and $\left(\mathsf{\alpha },\mathsf{\eta },\mathsf{\lambda },\mathsf{\gamma }\right)$ is the point on the release path, the upper and lower limits of the integral for $\mathsf{\alpha }$, $\mathsf{\eta }$ and $\mathsf{\lambda }$ are determined by the start and end coordinates of the release path. In numerical simulation, multi-source release can be regarded as composed of multiple discrete release point sources, and different ionospheric hole shapes can be obtained by modulating the temporal and spatial distribution of each release point source. In this paper, a spindle-shaped ionospheric hole is constructed by releasing CO₂ vertically at equal height intervals in 250–350 km. In order to ensure that the point sources in the source term are dense enough, the height release interval is selected as 5 km and the time interval is selected as 1 s.

2.2. Ion Chemical Reaction Process

The release of neutral gas molecules (such as CO₂, H₂O, etc.) in the background ionosphere can reduce the electron concentration in the plume region of the rocket. This is mainly because there are a large number of oxygen atomic ions ${\mathrm{O}}^{+}$ in the F layer of the ionosphere, whose recombination coefficient with electrons is about ${10}^{-12}{\mathrm{cm}}^3·{\mathrm{s}}^{-1}$, whereas neutral gases can easily transform atomic ${\mathrm{O}}^{+}$ into molecules, and the recombination coefficient can reach ${10}^{-7}{\mathrm{cm}}^3·{\mathrm{s}}^{-1}$ or even larger, which is about 3 orders of magnitude larger than the recombination coefficient of ${\mathrm{O}}^{+}$ and electrons [2,14]. Therefore, it can consume a lot of electrons in the ionosphere and form artificial ionospheric holes. The main chemical reaction equation after the release of neutral gas molecule CO₂ is shown in

Table 1. Chemical reaction equation of CO₂ in the ionosphere.

	Reaction Equation	Reaction Coefficient/cm3s−1
1	${\mathrm{CO}}_2+{\mathrm{O}}^{+}\stackrel{{\mathrm{k}}_1}{\to }{\mathrm{O}}_2{}^{+}+\mathrm{CO}+1.20\mathrm{eV}$	${\mathrm{k}}_1=9.4\times {10}^{-10}$
2	${\mathrm{O}}_2{}^{+}+{\mathrm{e}}^{-}\stackrel{{\mathrm{k}}_2}{\to }{\mathrm{O}}^{\ast }+{\mathrm{O}}^{\ast \ast }+6.90\mathrm{eV}$	${\mathrm{k}}_2=1.9\times {10}^{-7}{\left(300/T_e\right)}^{0.5}$

¹ Note: An asterisk indicates that the atom or molecule is in an excited state.

[15].

Among them, $T_e$ is the electron temperature, ${\mathrm{k}}_1$ and ${\mathrm{k}}_2$ are the reaction coefficients. Then, the reduction in electron density is $\Delta {\mathrm{n}}_{\mathrm{e}}={\mathrm{k}}_2\left({\mathrm{k}}_1·{\mathrm{n}}_{{\mathrm{CO}}_2}·{\mathrm{n}}_{{\mathrm{O}}^{+}}·\Delta \mathrm{t}\right)·{\mathrm{n}}_{\mathrm{e}}·\Delta \mathrm{t}$.

2.3. Plasma Diffusion

Due to the disturbance of the rocket plume, the change of the electron density in the release area of the plume chemical destroys the density distribution structure of the original charged particles and breaks the original dynamic balance. Affected by the effects of electron and ion density gradient, magnetic field, electric field, neutral wind, ionospheric disturbance zone collision, etc., plasma will drift, and its detailed force condition is more complicated. In the F layer of the ionosphere, the ion cyclotron frequency is much larger than the collision frequency of ion-ion and ion-neutral molecules; thus, we suppose that plasma is constrained to move along geomagnetic field lines, for the fact that its force is assumed stronger than others in the ionosphere, and ignoring the influence of magnetic deflection angle, the plasma diffusion equation can be obtained [13]:

$$\begin{array}{ll}\frac{\partial n_p}{\partial t}=D{\cos }^{2}I\frac{{\partial }^2{n}_p}{\partial y^2}& +D\frac{\cos I\sin I}{H_p}\frac{\partial n_p}{\partial y}-v_D\cos I\frac{\partial n_p}{\partial y}+D{\sin }^{2}I\frac{{\partial }^2{n}_p}{\partial z^2}\\ & +\left({\sin }^{2}I\frac{\partial D}{\partial z}+D\frac{1}{T_p}{\sin }^{2}I\frac{\partial T_p}{\partial z}+D\frac{1}{H_p}{\sin }^{2}I-v_D\sin I\right)\frac{\partial n_p}{\partial z}\\ & +(\frac{1}{T_p}{\sin }^{2}I\frac{\partial D}{\partial z}\frac{\partial T_p}{\partial z}+\frac{1}{H_p}{\sin }^{2}I\frac{\partial D}{\partial z}+D{\sin }^{2}I\frac{{\partial }^2\ln T_p}{\partial z^2}\\ & +D{\sin }^{2}I\frac{\partial \left(\frac{1}{H_p}\right)}{\partial z})n_p+P-L\end{array}$$

(4)

where $n_p$ is the numerical density of plasma (mainly refers to ${\mathrm{O}}^{+}$), D is the effective bipolar diffusion coefficient, $D=\left(1+\frac{T_e}{T_i}\right)D_i$, where $D_i$ is the ion diffusion coefficient [16], ${\upsilon }_p$ is the external drift velocity and its influence is ignored in this article, $T_p$ is the temperature of the plasma, $T_p=\left(T_e+T_i\right)/2$, $H_p$ is the scale height of plasma, $H_p=2k{T}_p/\left(m_{p}g\right)$, and $I$ is the magnetic inclination, $P$ is the production rate of plasma and $L$ is the loss rate of plasma.

2.4. 3D Ray Tracing

The radio wave propagation is described by Maxwell equations, and the high frequency wave propagation generally satisfies the geometrical optics approximation, so the propagation problem can be studied by ray tracing. The ray tracing method primarily includes analytical ray tracing [17,18,19,20] and digital ray tracing [21,22,23]. Since digital ray tracing technology can take into account the influence of geomagnetic field and other factors, this paper adopts digital ray tracing technology for simulating calculation. The three-dimensional ray tracing technology is one of the effective methods to describe and calculate ionospheric radio wave propagation [23]. The 3D ray tracing model in this paper is combined with Appleton Hartree equation. Considering the effect of the geomagnetic field and the more complex ionospheric structure on the ray path in the simulation, the Haselgrove ray equation must be solved digitally. Therefore, we need to use the Runge–Kutta method to solve the Haselgrove ray equation. Figure 1 is a schematic diagram of its three-dimensional coordinates.

Here, we assume that the earth is approximately a sphere. In Figure 1, OXYZ forms a Cartesian rectangular coordinate system. The coordinate origin O is located at Earth’s centre, OX points to 0° longitude, OY points to 90°E, and OZ points to the North Pole along the Earth’s axis, OXY is the equatorial plane. The coordinate of Q point is (r,$\phi ,\theta$), where r is the distance between Q and the centre of Earth, $\phi$ is the deviation angle from OX to the east to the longitude of Q point, and $\theta$ is the complementary angle of the latitude of Q. The wave vector at the Q point is $\overrightarrow{k}$ = ($k_r,k_{\phi },k_{\theta }$), where $k_r$ is perpendicular to the ground upward, $k_{\phi }$ points to the geographical east direction, and $k_{\theta }$ points to the geographical south direction. The ray tracing equations are solved based on the above coordinate system, then the Haselgrove equations become in the following form:

$$\{\begin{array}{c}\frac{dr}{d{P}^{\prime }}=\frac{c}{\omega }{k}_r\\ \frac{d\theta }{d{P}^{\prime }}=\frac{c}{\omega r}{k}_{\theta }\\ \frac{d\phi }{d{P}^{\prime }}=\frac{c}{r\omega \sin \theta }{k}_{\phi }\\ \frac{d{k}_r}{d{P}^{\prime }}=-\frac{\omega }{2c}\frac{\partial X}{\partial r}+k_{\theta }{}^2\frac{c}{r\omega }+k_{\phi }{}^2\frac{c}{r\omega }\\ \frac{d{k}_{\theta }}{d{P}^{\prime }}=\frac{1}{r}(-\frac{\omega }{2c}\frac{\partial X}{\partial \theta }-k_r{k}_{\theta }\frac{c}{\omega }+k_{\phi }{}^2\frac{c}{\omega }\cot \theta )\\ \frac{d{k}_{\phi }}{d{P}^{\prime }}=\frac{1}{r\sin \theta }(-\frac{\omega }{2c}\frac{\partial X}{\partial \phi }-k_{\theta }{k}_r\frac{c}{\omega }\sin \theta +\frac{c}{\omega }{k}_{\theta }{k}_{\phi }\cos \theta )\end{array}$$

(5)

In the above formula, P’ is the group path, which represents the wave packet distance moving at the speed of light. c is the speed of light, ω represents the angular frequency of the radio wave; ${X=\omega }_{plasma}{}^2{/\omega }^2$, where ${\omega }_{plasma}$ is the plasma angular frequency. After the initial emission conditions of the radio wave are determined: the emission position ($r_0,{\phi }_0,{\theta }_0$), the azimuth $\alpha$, the elevation angle $\beta$ and so on, the initial value of the ray differential equation can be obtained, which is shown in the following formula:

$$\{\begin{array}{c}r=r_0\\ \theta =\frac{\pi }{2}-{\lambda }_0\\ \phi ={\phi }_0\\ k_r=\left|\stackrel{\rightharpoonup }{k}\right|\cos (\frac{\pi }{2}-\beta )=\frac{\omega }{c}\sin \beta \\ k_{\theta }=-\left|\stackrel{\rightharpoonup }{k}\right|\cos \beta \cos \alpha =-\frac{\omega }{c}\cos \beta \cos \alpha \\ k_{\phi }=\left|\stackrel{\rightharpoonup }{k}\right|\cos \beta \cos (\frac{\pi }{2}-\alpha )=\frac{\omega }{c}\cos \beta \sin \alpha \end{array}$$

(6)

This is the ray component equation in the ionosphere that ignores magnetic fields and collisions. It is not only related to the position of the point and the component of the wave vector at this point, but also to the rate of change of plasma parameters at this point. Then the equations are solved according to the step size of the group path P’, and the obtained points are connected. Finally, the ray trajectory in the spherical coordinate system is obtained.

2.5. Experimental Parameter Setting

In this paper, the neutral atmospheric parameters are obtained by the NRLMSISE-00 model. The background ionospheric parameters are obtained by the empirical model IRI-2016, and the geomagnetic field parameters are obtained by the IGRF-13 model. Both the release point and the ray emission points are set in Nanjing (longitude 119° E, latitude 32° N), and the simulated release time is 1 July 2010. The numerical simulation platform is Matlab R2019b. Assuming that the rocket travels vertically at a speed of approximately 10 km/s and 3000 mol CO₂ is released every 0.5 s from 250 km to 350 km in the vertical direction, that is, 3000 mol CO₂ is released every 5 km on the vertical ascent path. After obtaining the spatial-temporal distribution of neutral gas and electron density, the three-dimensional digital ray tracing method can be used to simulate the effect of ionospheric disturbances due to rocket plume on the propagation path of HF radio wave. In this paper, a fast algorithm is used for shortwave ray tracing, and the group path step is adaptively adjusted according to the gradient of plasma frequency. The elevation range of the simulated shortwave is 75~105°, the elevation progression is 2°, and the azimuth is 90°.

3. Simulation Results and Discussion

3.1. CO₂ Releases Simulation Results

The process of ionospheric hole formation is simulated under the condition that the plume releases 3000 mol CO₂ every 5 km on the vertical ascending path of rocket. Figure 2 shows the evolution of the isosurfaces with electron densities of $0.3\times {10}^6{\mathrm{cm}}^{-3}$ and $0.6\times {10}^6{\mathrm{cm}}^{-3},$ respectively, after the CO₂ produced by the rocket plume is released. As shown in Figure 2, at the beginning of the disturbance, the volume of the ionospheric hole is relatively small and disturbance effect is relatively weak at the time of 30 s. With time goes by, the volume of the ionospheric hole gradually increases and the disturbance effect gradually deepens. It can be seen from the figure that the CO₂ produced by the rocket plume can form a large range of electron density holes in the ionosphere in a short time.

Figure 3 shows the vertical profile distribution of electron density at different times after CO₂ release. It can be seen from the figure that the electron density changes rapidly with time before t = 100 s. When t = 100 s, the maximum relative change rate of electron density is 85.0%, but at this time the ionospheric hole is still in the stage of development and has not yet reached saturation. Comparing the vertical profiles of electron density from 0 s to 100 s and from 100 s to 300 s, it is found that with the passage of time, the decreasing rate of electron density decreases greatly after 100 s, and the shape of ionospheric hole is gradually formed. The red dotted line in the figure shows the decrease in electron density at 500 s, and the decrease in electron density reaches the maximum at 296 km, when the maximum relative change rate of electron density has reached 96.4%. Combining the evolution process of t = 30 s − 500 s electron density vertical profile, it is found that the height of the ionospheric hole is decreasing, which is consistent with the result of the continuous decrease in the height of the ionospheric hole center in the CO₂ single-source release simulation experiment carried out by Lyu et al. [24].

Figure 4 shows the spatial distribution of CO₂ concentration after 100 s release of CO₂. Figure 4a is a three−dimensional view of the release CO₂, which is spatially distributed as an ellipsoid with a lower narrow and upper width, because the diffusion coefficient of the release is anisotropic in the horizontal plane [24]. Figure 4b shows the vertical profile of CO₂ concentration with height in the east-west direction of y = 0 km, which is shaped like the flame of a candle, and the release concentration is of the order of 10⁶ cm⁻³ or above. It can be seen from the figure that the diffusion range of CO₂ increases gradually with the increase in height, this is because the background atmospheric temperature increases with height, and the diffusion coefficient of the release is positively correlated with the background atmospheric temperature [25]. Figure 4c is a horizontal section at the height of z = 255 km. The difference in concentration in the horizontal direction is not obvious after 100 s release of CO₂. The horizontal section is an approximate circle, and the diffusion radius in the east-west and north-south directions is about 24 km.

Figure 5 shows the spatial distribution of electron density Ne after 100 s release of CO₂. Similar to the distribution of the release concentration, the spatial distribution of the electron density is also an ellipsoid with a narrow bottom and a wide top, but the difference is that the vertical profile of the electron density in the east-west direction at Figure 5b y = 0 km is approximately spindle−shaped. This is similar to the “pear” −shaped structure with larger top and smaller bottom in the vertical profile of electron density distribution in Hu et al. [6] and the approximate spindle-shaped ionospheric cavity cross section after 100 s release of H₂O [26]. Figure 5b shows that the radius of the ionospheric cavity on 335 km is largest, about 42 km, whereas that of the ionospheric cavity on 237 km is smallest, about 24 km. This is also similar to the simulation experiment carried out by Zhu et al. [26] on the formation of ionospheric holes by combustion and release of H₂O during ballistic missile flight (releasing 2240 mol H₂O per 20 km from 200 km to 360 km in the vertical direction). The maximum radius of the cavity is approximately 40 km and the minimum radius is 20 km. Figure 5c shows the horizontal section of the ionospheric hole at the height z = 300 km, which is similar to the horizontal distribution of CO₂ concentration and is also approximate to a circle. It can be seen from the figure that the radius of the cavity at the height of 300 km is about 36 km.

3.2. Simulation Results of Short-Wave Ray Tracing

When the rocket passes through the ionosphere, the neutral gas released by the rocket plume will cause the ionospheric electron density loss and form the ionospheric hole. In order to further study the influence of the electron density disturbance structure caused by CO₂ release on the information transmission link, shortwave ray tracing is carried out to simulate the change in the propagation path and direction of shortwave signals with different frequencies after passing through the disturbance zone [27]. Figure 6, Figure 7, Figure 8 and Figure 9 shows the disturbance effect of ionospheric holes to shortwave propagation at various frequencies after CO₂ release at 100 s, 300 s, 500 s and 700 s. When the radar transmitted wave frequency f is 5 MHz, the elevation angle is 75~105°, the elevation progressive length is 2°, and the azimuth is 90°. The simulation results of shortwave ray tracing are shown in the following figure.

As shown in Figure 6, when the release time is 100 s, the volume of the ionospheric hole is still very small. However, when the disturbance time reaches 700 s, the volume of the ionospheric hole is already very large and the shape of the disturbance region has changed. With the passage of time, the shape of the ionospheric cavity gradually evolved from the ellipsoid at the beginning of 100 s, in which the upper is wide and the lower is narrow, to the ellipsoid with middle width, narrow top and bottom at 700 s. In addition, it can be found that at the peak height of ionospheric electron density (i.e., ~300 km), the electron density dissipates fastest and the radius of the ionospheric hole is largest. When t = 700 s, the radius of the ionospheric hole near 300 km has reached about 85 km. The reason for this phenomenon may be that because the background atmospheric temperature increases with the increase in height, and the diffusion coefficient of the release is positively correlated with the background atmospheric temperature, the diffusion coefficient of the released CO₂ increases gradually with the increase in height [24]. The electron density increases at first and then decreases with the increase in height, the electron density is highest around 300 km, and the rate of ion chemical reaction with CO₂ is fastest, so the radius of the ionospheric hole is largest at the peak height of ionospheric electron density. Due to plasma diffusion, the horizontal distribution of ionospheric disturbances is uneven. At the peak height of the ionospheric electron density (i.e., ~300 km), the electron density dissipates most quickly and the disturbance radius of the ionospheric hole is largest.

By analyzing the ray tracing results of Figure 6, it is found that the radio waves of 5 MHz are completely reflected by the ionosphere and cannot reach the ionospheric disturbance region. Although with the passage of time, the volume of the disturbance region continues to increase and the lower boundary of the cavity continues to extend downward, the shortwave of 5 MHz is still less than the critical frequency of the lower boundary of the cavity and cannot reach the disturbance zone. This demonstrates that the ionospheric hole produced by the rocket plume has no effect on the propagation path of the 5 MHz shortwave.

Figure 7 shows the simulation results of the propagation path of HF radio wave when the frequency is 6 MHz. As can be seen, the radio waves of 6 MHz are completely reflected by the ionosphere and cannot reach the ionospheric disturbance zone at t = 100 s. This is because the range of ionospheric disturbance is relatively small, the lower boundary of the cavity is still relatively high, and the shortwave rays cannot reach the lower boundary of the cavity at t = 100 s. With the passage of time, when the release time t is 300 s, most of the shortwave rays of 6 MHz are directly reflected by the ionosphere and cannot reach the ionospheric disturbance region. However, as the volume of the cavity increases and extends downward, the shortwave rays with the elevation nearly 90° penetrate through the bottom of the ionospheric hole and are finally returns to the ground after complex refraction and reflection in the ionospheric disturbance zone. The reason for this phenomenon is that the electron density of the upper boundary of the ionospheric cavity is larger than that of the lower boundary, so that the shortwave rays can penetrate the lower boundary of the cavity but not the upper boundary [24]. When the release time reaches 500 s and 700 s, with the continuous development of the ionospheric cavity, more and more rays can penetrate the bottom of the ionospheric cavity and are finally returned to the ground after complex refraction and reflection in the ionospheric disturbance zone.

In the short beam of 8 MHz, the shortwave with the elevation closest to 90° can penetrate the disturbed ionospheric region, whereas the shortwave with a smaller elevation angle is reflected back to the ground at t = 100 s. Among them, rays of approximately 85~95° pass through the ionospheric disturbed region and produce electromagnetic wave focusing effect at the height of about 280 km. This is because the decrease in electron density in the ionospheric disturbance region increases the refractive index of the signal [28]. In addition, we find an interesting phenomenon that the “secondary focusing effect” occurs when the rays with an elevation range of 85~95° pass through the top of the cavity, and there is a “second focal point” at the height of about 400 km, which has not appeared in the previous simulation experiments. However, this phenomenon is accidental. With the passage of time, this phenomenon disappears at 300 s, then there is only one focus. We also find that with the passage of time, the disturbance range of the ionospheric hole gradually increases, and the disturbance effect gradually deepens. When t = 300 s, most of the radio waves penetrate the ionosphere to a higher height, but there are still a few rays that do not reach the ionospheric disturbance region and directly reflect back to the ground. When t = 500 s, with the further increase in the disturbance range of the ionospheric hole and the further deepening of the disturbance degree, all the radio waves pass through the ionospheric disturbance region, but the shortwave with a small elevation still cannot penetrate the upper boundary of the ionospheric hole and passes through the side boundary of the cavity after complex refraction and reflection.

Figure 9 illustrates the paths of 10 MHz short-wave rays. As can be seen, the radio waves of 10 MHz all penetrate the disturbance region and the ionosphere, but the decrease in the electron density in the disturbance region increases the refractive index of the signal, so the signal passing through the disturbance region has an obvious focusing effect. For the same incident frequency, with the increase in the release time of CO₂, the disturbance volume of the ionosphere increases, the effect of disturbance deepens gradually, the focusing effect of radio waves weakens gradually, and the focus becomes higher and higher.

Figure 10 shows the three-dimensional ray path simulation results of HF radio waves with different frequencies after 100 s of CO₂ release. The incident frequencies of shortwaves are 6 MHz, 8 MHz, 10 MHz and 15 MHz, respectively. The elevation angle is 75–105°, one beam is emitted every 2°, and the azimuth is 90°.

According to Figure 10, all the radio waves of 6 MHz are completely reflected by the ionosphere and cannot reach the ionospheric disturbance region. When the incident frequency of HF shortwave increases to 8 MHz, most of the shortwaves do not reach the ionospheric disturbance region and are directly reflected to the ground. While the “secondary focusing effect” occurs when the rays with an elevation range of 85~95° pass through the top of the cavity, and there is a “second focal point” at the height of about 400 km. However, this phenomenon is accidental, it disappears with the increase in CO₂ release time, then there is only one focus. We believe that the “second focal point” is due to the reflection of HF shortwave from the walls of the ionospheric holes, and according to this idea, the formation of the “second focal point” requires the following two conditions to be met, the first condition is the ionospheric hole is as narrow as possible vertically to let the ray refracted past the focal point to reach the reflection point on the wall of the ionospheric cavity, the second condition is the frequency of HF shortwave cannot be too high to make the altitude of the first focal point as low as possible. In this situation, why the “secondary focusing effect” occurs only at the early 100 s with 8 MHz HF shortwave can be explained. As the incident frequency increases to 10 MHz, all the shortwaves penetrate the disturbed region of the ionosphere and are no longer reflected back to the ground. The radio waves of 10 MHz all penetrate the disturbance region and the ionosphere, but the decrease in the electron density in the disturbance region increases the refractive index of the signal, so the signal passing through the disturbance region has an obvious focusing effect. All the radio waves of 12 MHz also penetrate the disturbance region and the ionosphere, but the focusing effect is weaker than that of 10 MHz signal because the refractive index of the disturbance region decreases with the increase in radio signal frequency.

4. Conclusions

In this paper, the process of ionospheric disturbance produced by the plume of a rocket in a vertical ascending path is simulated, with an assumption that 3000 mol CO₂ are released per 5 km in the height range from 250 km to 350 km. The influence of this ionospheric hole on the propagation of HF radio waves at different frequencies is also simulated. In this work, the initial neutral atmospheric parameters, background ionospheric parameters and geomagnetic field parameters are obtained from NRLMSISE-00 model, IRI-2016 model and IGRF-13 model, respectively. Taking the neutral gas diffusion process, ion chemical reaction and plasma diffusion process into account, a three-dimensional kinetic model is constructed, by which we can obtain the ionospheric disturbance effect of CO₂ released from the rocket plume at different times. In addition, the 3-D digital ray tracing technique is used to simulate the effect of the ionospheric hole on the propagation path of HF radio waves with different incident frequencies. A fast algorithm is used in the simulation, and the group path step is adaptively adjusted according to the gradient of plasma frequency. In lower ionospheres, the density of oxygen atomic ions and the recombination rate are much lower than those in higher ionosphere, so the ionospheric hole caused by a rocket in lower ionosphere may not as obvious as that in the higher ionosphere. Thus, we did not consider the lower ionospheric disturbances in this paper. Finally, we draw the following conclusions:

(1)

The spatio-temporal distribution of ionospheric disturbances at 100 s, 300 s, 500 s and 700 s after the release of main chemical CO₂ is simulated. Around 100 s after release, the spatial distribution of CO₂ is a narrow ellipsoid at the bottom and wider at the top, which is due to the anisotropy of diffusion coefficient in the horizontal plane [24]. The vertical section of the east-west direction at y = 0 km is shaped like the flame of a candle, the concentration of the release is on the order of 10⁶ cm⁻³ or above, and the horizontal section at z = 255 km is an approximate circle. After the release of CO₂ for 100 s, the spatial distribution of electron density Ne is also similar to an ellipsoid, in which the upper is wide and the lower is narrow, but the vertical profile of electron density in the east-west direction at y = 0 km is approximately spindle-shaped. The horizontal section at z = 300 km is also close to a circle.

(2)

At the peak height of ionospheric electron density (~300 km), the electron density dissipates fastest and the radius of ionospheric hole is largest. We think that the reason for this phenomenon may be that the background atmospheric temperature increases with height increasing, and the diffusion coefficient of the release is positively correlated with the background atmospheric temperature. Thus, the diffusion coefficient of the released CO₂ increases gradually as height increases [24]. The electron density increases at first and then decreases with the increase in height. The electron density is highest around 300 km, and the rate of ion chemical reaction with CO₂ is fastest, so the radius of the ionospheric hole is largest at the peak height of ionospheric electron density.

(3)

By analyzing the ray tracing results, it can be found that the radio waves of 5 MHz are completely reflected by the ionosphere and cannot reach the ionospheric disturbance region. The radio waves of 6 MHz are completely reflected by the ionosphere and cannot reach the ionospheric disturbance zone at t = 100 s. With time going on, the volume of the cavity gradually increases and extends downward. More and more shortwave rays of 6 MHz with the elevation nearly 90° penetrate through the bottom of the ionospheric hole and finally return to the ground after complex refraction and reflection in the ionospheric disturbance zone. For the shortwave of 8 MHz, we find an interesting phenomenon that the “secondary focusing effect” occurs when the rays with an elevation range of 85~95° pass through the top of the cavity, and there is a “second focal point” at the height of about 400 km. However, this phenomenon is accidental and we believe that the “second focal point” is due to the reflection of HF shortwave on the walls of ionospheric holes. As time progresses, this phenomenon disappears at 300 s, and then there is only one focus. In addition, with the continuous increasing of incident frequency, the shortwaves of 10 MHz and 15 MHz will all penetrate the disturbance region and the ionosphere.

(4)

For the same incident frequency, with the increase in CO₂ release time, the focusing effect of the radio wave is weakened gradually and the focus becomes higher; when the release time of CO₂ remains the same, the focusing effect decreases with the increase in the radio signal frequency and the focus also rises.