Investigating the prediction ability of the ionospheric continuity equation during the geomagnetic storm on May 8, 2016

Hany Mahbuby; Yazdan Amerian

doi:10.1515/jogs-2022-0175

Article Open Access

Investigating the prediction ability of the ionospheric continuity equation during the geomagnetic storm on May 8, 2016

Hany Mahbuby and Yazdan Amerian

Published/Copyright: January 24, 2025

Published by

Become an author with De Gruyter Brill

Author Information

From the journal Journal of Geodetic Science Volume 15 Issue 1

Abstract

Ionospheric electron density (IED) and vertical total electron content (VTEC) are two of the most important characteristics that interpret the ionosphere. Today, satellite geodesy plays an important role in providing these data. Apart from monitoring the ionosphere by means of Global Navigation Satellite System observations, predictions of these ionospheric characteristics have also been taken into consideration. The majority of prediction methods in the ionosphere are focused on VTEC prediction. The continuity equation in the ionosphere is a spatiotemporal partial differential equation that can be used to predict both IED and VTEC. In this study, we address this issue during May 8, 2016, which was a day of high geomagnetic activity. For this purpose, first, high-accuracy IED grids with a spatial resolution of 0.5° × 0.5° in longitude and latitude, 50 km in altitude, and a temporal resolution of 10 min are prepared. These IED grids are called analysis grids and are derived by assimilation of GPS-derived VTECs into the background IED grids provided by international reference ionosphere. Second, the analysis grids can be utilized as the initial and boundary values in the continuity equation to predict IED for the next 3 h with a step of 10 min. The analysis grids are constructed over the Iran region, and the performance of the continuity equation in the ionosphere prediction is investigated during the strongest geomagnetic storm of 2016. Finally, the accuracy of the predictions has been evaluated in two ways. First, analysis grids with high accuracy are available for each epoch in which the prediction is presented. On average, the root mean square (RMS) of the differences between the predicted IEDs and the analysis IEDs is 1.66 × 10¹⁰ el/m³ for 1-h predictions and 5.33 × 10¹⁰ el/m³ for 3-h predictions. Second, by numerical integration of the predicted grids in the vertical direction, VTECs are estimated and compared with VTEC values of test data that were already known. On average, the RMS of their difference was 1.29 total electron content unit (TECU) for 1-h predictions and 2.57 TECU for 3-h predictions. An acceptable accuracy for both IED and VTEC prediction by a continuity equation has been obtained up to 3 h on the most active ionospheric day.

Graphical abstract

Keywords: ionospheric continuity equation; ionospheric data assimilation; ionospheric electron density; ionospheric prediction; vertical total electron content

MSC 2010: 86-08; 86-10; 86A10; 86A22; 86A32; 65F05; 65F20; 65N06; 65N22; 65N50; 65N55

1 Introduction

Ionosphere is the ionized layer of Earth’s atmosphere, which starts from a height of about 50 km above the topography and has a remarkable concentration of ions and electrons up to about 1,000 km. The concentration of electrons causes this environment to affect the propagation of electromagnetic signals. Therefore, the distribution and concentration of electrons in the ionosphere and their total number are the parameters that describe the physical characteristics of the ionosphere. Today, Global Navigation Satellite System (GNSS) receivers are one of the most important tools for monitoring the ionosphere. The geometry-free code measurements that are smoothed by geometry-free phase measurements are one of the most common linear combinations of GNSS observations for ionosphere reconstruction. This observation leads to the finding of the slant total electron content (STEC), which actually shows the number of electrons in a cylinder with a unit cross-sectional area between the satellite and the receiver along the line of sight. The point of intersection between the line of sight and the single layer of the ionosphere, which is usually considered at the height of the ionospheric electron density (IED) peak, is called the ionospheric pierce point (IPP). STEC is divided by a mapping function in terms of the satellite’s elevation angle and imaged in the zenith direction of IPP. In this way, the so-called vertical total electron content (VTEC) will be obtained. VTEC is one of the most important parameters to show the characteristics of the ionosphere in a geographical region.

To estimate instantaneous ionospheric corrections for the electromagnetic signal’s propagation, especially in real-time GNSS positioning and also for radio-communication purposes, ionospheric VTEC predictions are of importance. In this respect, the majority of the studies have utilized machine learning (ML) techniques and achieved accuracies from 1 total electron content unit (TECU) to 5 TECU measured by root mean square (RMS) of the differences between the predicted and known VTEC values (Chen et al. 2022, Lei et al. 2022, Liu et al. 2022, Wen et al. 2021). The accuracy of the predicted values in these studies is better than the empirical international reference ionosphere (IRI) model. Besides, other regional ionospheric studies have proven that ML techniques are also powerful tools for regional VTEC prediction (Huang et al. 2022, Razin and Voosoghi 2016, Shi et al. 2022). Utilizing the multi-factor neural prophet method for regional ionospheric prediction, Huang et al. (2022) achieved an accuracy of 2.33 TECU and 3.12 TECU in stable and geomagnetic storm periods, respectively. Shi et al. (2022) achieved an accuracy of 1.12 TECU and 1.68 TECU for 1 and 2 h ahead predictions on the test data using the bidirectional long short-term memory method. Also, Razin and Voosoghi (2016) obtained errors in the range of 1.8 TECU to 3 TECU on the test data in a regional study using wavelet neural networks.

An advantage of IED prediction is that in addition to VTEC, the critical frequency of the ionosphere can also be predicted, which is important for telecommunications. It is clear that in ML techniques, the prediction interval in which the accuracy is reasonable strongly depends on the training period and the number of input data. Hence, if the prediction of IED in different layers of the ionosphere is desired, ML techniques may not be very effective in terms of complexity and cost. This is due to the significant increase in unknown parameters, which requires more input data and a longer training period. Dutta and Cohen (2022) predicted IED in the topside ionosphere using ML on defense metrological satellite program data with mean absolute error between 0.5 × 10¹⁰ el/m³ and 1.5 × 10¹⁰ el/m³ for 2 years of test data. Chen et al. (2023) predicted global IED grids using a deep neural network on FORMOSAT-3/COSMIC electron density profiles and reported RMS of the differences with known values between 0.5 × 10¹¹ el/m³ and 1 × 10¹¹ el/m³. They concluded that the accuracy of predictions is better at lower latitudes. Yang and Fang (2023) utilized IRI as a prior knowledge to develop an artificial neural network for IED prediction. They obtained an RMS of 0.9 × 10¹¹ el/m³ for the difference of predictions with 2 years of test data. They showed that a prior knowledge can improve the accuracy of ML techniques for ionospheric prediction.

In addition to ML techniques, data-based empirical modeling and data assimilation methods are also utilized for IED prediction. However, a very efficient method for IED prediction that has a long history is physics-based modeling. Zhang et al. (2018) evaluated the assimilation of some physical parameters on IED predictions. They showed using proper physical parameters in data assimilation can improve the accuracy by 40%. Jee (2023) described the fundamental aspects of physics-based ionospheric modeling compared to data-based empirical and assimilative models. In this study, it is concluded that the physics-based model has a fundamental importance since it can describe the physical processes occurring in the ionosphere and also provides suitable forecasts of the state of the ionosphere in space weather applications. To study the spatial and temporal changes of IED, the ionospheric continuity equation (ICE) has been used for many years (Hsu et al. 2014, Li et al. 2020, Shimazaki 1965). ICE considers the photochemical activities in the ionosphere that lead to electron production and disappearance and also describes the electron transport. Therefore, physics-based models that use ICE can be really useful in studying physical phenomena related to the ionosphere. In this regard, the variations of IED profiles of lower ionosphere obtained from ICE with respect to solar flares and solar irradiation are investigated by Chakraborty and Basak (2020). Also, Žigman et al. (2023) presented a new model based on ICE and space and ground observations to predict IED and recombination coefficient of the lower ionosphere under solar flare conditions.

Clearly, investigation of IED variations can also lead to understanding VTEC changes, while the opposite is not possible. Hence, being able to predict electron densities and then VTEC is an advantage over studies that only use ML techniques to predict VTEC.

ICE is a nonlinear partial differential equation, which can only be solved using numerical schemes if not simplified. Usually, this equation is simplified into an ordinary differential equation, and an analytical solution is found for it (Chapman 1931, Memarzadeh 2009, Prölss 2012). In this study, the average rate of electron production by ionospheric gases is considered a multiple of the solar radiation flux received on the horizon plane. The rate of electron disappearance is assumed to be a quadratic parabola in terms of IED in different layers, and the electron transport rate is assumed to be a constant value in a small interval, e.g., 10 min (Li et al. 2020). Under the stated conditions, ICE is utilized to predict IED values for 10-min intervals.

The vital requirement to provide an accurate prediction by ICE is to have accurate initial conditions with a proper geometric distribution, i.e., in the form of a 3D regular grid, which is called an analysis grid. In fact, before the prediction process, the IED modeling process must be implemented. In order to make an analysis grid, it is necessary to have enough observations that are properly distributed in the study region. By data assimilation process, it is possible to combine two categories of data, which are called observation and background, and produce a proper analysis grid. In this study, the background is IED values obtained from the global empirical international reference ionosphere (IRI2016) model (Bilitza et al. 2017), and the observations are obtained VTECs from permanent dual-frequency GPS receivers in Iran. The period and region of the study have been chosen quite challenging. Because in the study region, the density of permanent GPS stations and consequently IPPs is low, and the study date, i.e., May 8, 2016, is also the most active ionosphere day with the highest values of Kp and Ap indices.

The method described in Mahbuby and Amerian (2021) is used to estimate STECs along with the receivers’ differential code biases (DCBs). Thereafter, STECs are converted to VTECs, and the first category of data, i.e., observations, are prepared for data assimilation. The second category, i.e., background, is IRI electron density grids. The accuracy of background data in the study region is low, and the obtained IED values from IRI on May 8, 2016, are biased in the region (Mahbuby and Amerian 2023). The analysis grid is obtained from the assimilation of GPS-derived VTECs into the electron densities obtained from IRI via the presented data assimilation method by Mahbuby and Amerian (2023). This analysis grid provides an accurate initial condition for ICE.

Considering electron production, disappearance and transport processes that ICE implies, this equation has parameters that must be correctly determined for accurate prediction. In fact, these parameters are production, attachment, and recombination coefficients, and the electron transport rate, which have been quantified in numerous studies. In this respect, various experimental relationships have been introduced for these parameters with regard to the physical characteristics of the ionosphere, e.g., solar radiation flux, the diameter of ionospheric particles, the temperature of the electrical charges, and the pressure (Biondi 1969, Dogariu et al. 2010, Seaton 1947).

In this study, since the initial conditions (analysis grids) are prepared with appropriate accuracy, they are also used to estimate the required parameters of ICE. That means, given the initial analysis grids at the first two consecutive epochs, the coefficients are estimated with regression, and then the estimated coefficients and the initial conditions are used for prediction. The advantage of this method is that it can be easily used for the estimation of ICE parameters and prediction, even if no enough data are available for the computation of parameters.

2 Methodology

In this section, we first describe ICE and the related parameters. To solve the equation, we need initial conditions. In the following, the method of preparing the initial conditions with data assimilation will be explained. Also, the allocation of parameters will be investigated, and finally, the total processing scheme used to predict IED and VTEC will be reviewed.

2.1 ICEs and the related parameters

There are a significant number of free electrons in the ionosphere. Free electrons and ions are created by the ionization of neutral particles by the energy of solar radiation. The concentration of electrons in the ionosphere depends on the altitude, and in this respect, the highest concentration of electrons is between 90 and 600 km above the earth. This range is divided from bottom to top into three areas, E, F1, and F2 (Figure 1). The temporal changes of the ionized substance in a unit volume are expressed by the continuity equation. This equation considers the photochemical activities in the ionosphere that lead to production, disappearance, and also transport. This equation is expressed as follows (Ivanov-Kholodny and Mikhailov 2012):

(1) ∂ n i ∂ t + div ( n i V i ) = q i − L i ,

where n i stands for the ith ion-type density in the plasma and q i , L i represents the production and disappearance rates per unit volume, respectively. V i shows the transfer speed of the ith ion species and the term div ( n i V i ) expresses the density changes of the ith ion species due to the transport action. The quantity n i V i stands for the electric charge flux, also known as the electric current density, i.e., the amount of electric charge that flows in a unit cross section per second.

Figure 1

Various regions of the ionosphere and electron density distribution.

The number of produced or disappeared electrons in the plasma is equal to the number of ions that are produced or disappeared; therefore, if n i displays an ionic species density and the ionosphere has a total K ion ionic species, the electron density N e follows the rule of the same amount of positive and negative charges in plasma. Hence, N e reads

(2) N e = ∑ i = 1 K ion n i .

The rate of ion production by ionospheric gases is totally dependent on the amount of solar energy absorbed in the upper atmosphere, i.e., the average electron production rate q ¯ is directly related to the amount of solar radiation flux which is projected on the horizon plane (Li et al. 2020):

(3) q ¯ = k 1 I ,

where I stands for the received solar radiation flux on the horizon plane. If the solar radiation flux in the plane perpendicular to the radiation path is shown by I ′ , the following relationship exists between I and I ′ (Memarzadeh 2009):

(4) I = I ′ cos χ ,

where χ shows the zenith angle of solar radiation, which is shown in Figure 2. This angle is a function of local time, declination of the sun δ , and the latitude ϕ :

(5) cos χ = sin φ sin δ + cos ϕ cos δ cos h s .

Figure 2

Solar constant, solar radiation flux on the horizon plane, the zenith angle of solar radiation, and declination of the sun.

In this equation, h _s shows the hour angle of the sun, which can be obtained from the local time t as follows:

(6) h s = 2 π ( t − 14 ) 24 .

The declination of the sun varies from −23.5° to +23.5° during a year. It can be calculated for each day of year (DOY) based on the following relationship:

(7) sin δ = 0.398 cos 0.986 ( DOY - 173 ) π 180 .

The received solar radiation flux in the plane perpendicular to the radiation path at the top of the ionosphere is almost a constant value I ′ = 1 , 366 W m 2 , which is called the solar constant (Figure 2). In addition, it should be noted that in the hemisphere where it is nighttime, i.e., the sun’s hour angle is h s ≤ − π 2 or h s ≥ + π 2 , the radiation flux is considered zero. Figure 3 shows the received solar radiation flux at the top of the ionosphere projected on the horizon plane over the Iran region on May 8, 2016, at 10 UT. The projected solar radiation flux on the horizon plane is obtained from equation (4) and around the local noon, it is decreasing from south to north.

Figure 3

The solar radiation flux projected on the horizon plane over the Iran region on May 8, 2016, at 10 UT.

In general, the ionosphere can be divided into five regions D, E, F1, F2, and the topside of the ionosphere above F2, as can be seen in Figure 1. The rate of electron disappearance varies in different regions of the ionosphere. The lowest region of the ionosphere is D. The altitude of this region is less than 90 km from the earth’s surface and has a low electron density. This region is formed immediately after sunrise and disappears immediately after sunset. Due to the low electron density in this region, it can be inferred that the role of this region in forming VTEC is insignificant. The main cause of the electron disappearance in the E region is the process of dissociative recombination (DR) of electrons with ionized molecules of NO + and O 2 + :

(8) O 2 + + e → O + O NO + + e → N + O .

If L E represents the electron disappearance rate in the E region, it can be defined as follows (Prölss 2012):

(9) L E ≃ K NO + DR n NO + N e + K O 2 + DR n O 2 + N e ,

where K NO + DR and K O 2 + DR are the coefficients that represent the rate of electron disappearance as a result of DR with the ionized molecule written in the subscript. These coefficients are almost equal, and if we consider the relationship between the density of ionized molecules and IED in the E region:

(10) n NO + + n O 2 + = N e .

It can be simply deduced that the disappearance rate in the E region is a multiple of the square of IED:

(11) L E = κ 2 N e 2 .

In the F region, the most important factor in the loss of electrons is charge exchange (CE):

(12) O + + N 2 → NO + + N O + + O 2 → O 2 + + O .

In this region, the disappearance rate L F can be written as follows:

(13) L F ≃ K O + , N 2 CE n N 2 n O + + K O + , O 2 CE n O 2 n O + ,

where K O + ,N 2 CE and K O + ,O 2 CE are the coefficients that represent the rate of electron disappearance as a result of CE between ions and molecules written in the subscript. The dominant ion in the F region is O⁺, whose density is almost equal to the electron density. Hence, it can be inferred that the rate of electron disappearance in this region is a multiple of IED:

(14) L F = κ 3 N e .

According to equations (11) and (14), it can be stated that the electron disappearance rate in the ionosphere is generally a quadratic function of the IED value (Li et al. 2020):

(15) L ( N e ) = κ 2 N e 2 + κ 3 N e .

If ICE in equation (1) is written for free electrons in the ionosphere, we obtain

(16) ∂ N e ∂ t = q ¯ − L ( N e ) − div ( N e V ) ,

where q ¯ and L ( N e ) represent the average electron production rate and electron disappearance rate in the ionosphere, which can be replaced by equations (3) and (15), respectively. The term div ( N e V ) expresses the changes created in IEDs due to electron transport. In a short period of time, this term can be assumed to be a constant value κ 4 in the continuity equation (Li et al. 2020). Therefore, ICE for prediction in short time steps can be simplified to

(17) ∂ N e ∂ t = κ 1 I − κ 2 N e 2 − κ 3 N e − κ 4 .

2.2 Providing initial conditions

To provide an accurate IED analysis grid as initial conditions for equation (17), two sets of data are utilized. First, STECs are obtained from code and phase observations of dual-frequency GPS receivers at each epoch ( n ) . For this purpose, phase-smoothed geometry-free code measurements are utilized (Arikan et al. 2008, Ciraolo et al. 2007):

(18) STEC r s ( n ) = 1 40.3 f 1 2 f 2 2 f 1 2 − f 2 2 ( B + L 4,r s ( n ) + c ( DCB r + DCB s ) ) ,

where f 1 and f 2 indicate carrier phase frequencies. Subscript r and superscript s stand for the receiver and satellite indices, respectively. DCBs for receivers and satellites are shown by DCB r and DCB s . B in equation (18) stands for the mean of the differences between geometry-free linear combination of code measurements P 4,r s and geometry-free linear combination of phase measurements L 4,r s and is computed in a continuous arc of phase measurements by

(19) B = 1 τ ∑ i =1 τ ( P 4,r s ( i ) − L 4,r s ( i ) ) .

Here, τ shows the number of observations in a time interval without cycle-slip. As is obvious in equation (18), STEC values at each epoch are obtained using phase-smoothed geometry-free code measurements, i.e., B + L 4,r s ( n ) . Each satellite’s DCB value is provided by the Center for Orbit Determination in Europe (CODE) with an accuracy of about hundredths of a nanosecond. The receivers’ DCB values are estimated hourly using Global Ionospheric Maps (GIMs) and a constrained least-squares method represented by Mahbuby and Amerian (2021). Thereafter, STECs are divided into the following mapping function M ( ϵ s ( n ) ) and converted to VTECs:

(20) M ( ϵ s ( n ) ) = 1 − R E cos ( ϵ s ( n ) ) R E + h 2 − 1 / 2 ,

where ϵ s ( n ) represents the local elevation angle of the satellite (s) at each epoch. R E and h indicate the radius of Earth and the applied spherical shell height for computation of IPP locations. These retrieved VTEC values are used as observations set in the data assimilation process to provide the analysis grid. The second set of data used to provide the initial conditions is a 3D regular grid of electron densities obtained from the IRI2016 online service. This dataset is the background. Analysis grids of electron densities, which are considered as initial conditions for ICE, are provided throughout a data assimilation process that assimilates GPS-derived VTECs as observations into the IRI electron density girds as background. The following reference epochs are considered for providing initial conditions:

(21) t 0 = { 00 : 10 , 03 : 10 , 06 : 10 , 09 : 10 , 12 : 10 , 15 : 10 , 18 : 10 , 21 : 10 } UT .

Data assimilation is done in 10-min subintervals. In the middle epoch of each subinterval, a 3D regular grid of electron densities is received from the IRI2016 online service. The grid has a resolution of 0.5° × 0.5° in longitude and latitude and consists of 39 height levels from an altitude of 100–2,000 km with a fixed step of 50 km.

For data assimilation, the optimal interpolation method has been used. The background grid is represented by N B and the analysis grid resulted from data assimilation is denoted by N A and obtained as follows:

(22) N A ( λ i , j , ϕ i , j , h k ) = N B ( λ i , j , ϕ i , j , h k ) + G ( V R ( λ 0 , ϕ 0 ) − H N B ( λ i , j , ϕ i , j , h k ) G = B H T ( R + H B H T ) − 1 ,

where i , j stand for the indices of the regular spherical grid points in the directions of longitude λ and latitude ϕ and h k denotes the altitude of the regular spherical grids. V R is the vector of observations, i.e., retrieved VTECs at IPP coordinates ( λ 0 , ϕ 0 ) . Here, H represents the observation operator matrix, which maps the electron densities of the background to VTECs. The optimal data assimilation weight G is a function of the covariance matrices of observations and background, i.e., R and B .

Integration in the zenith direction over the weighted average of IED values obtained from the nearest four grid points to the IPP location in each layer results in VTEC (Figure 4). Hence, weighted averaging and integration map the electron densities of the background to VTEC at IPP (Mahbuby and Amerian 2023)

(23) V B ( λ 0 , ϕ 0 ) = ∑ k = 1 K ∑ t < 4 > w t , k N B ( λ t , ϕ t , h k ) Δ h k , w t , k = 1 g t ∑ t = 1 4 1 g t ,

where subscript t refers to the indices of the closest four grid points to the IPP coordinates ( λ 0 , ϕ 0 ) in each height level k . The integration step between consecutive layers is denoted by Δ h k , and the weights are shown by w t , k , which is a function of the distances between four closest grid points and IPP, i.e., g t . Figure 4 illustrates the IPP coordinates, the nearest grid points, and the distances between them in each height level.

$Figure 4 The closest four grid points to the IPP coordinates ( λ 0 , ϕ 0 ) ({\lambda }_{0},{\phi }_{0}) in each height level (Mahbuby and Amerian 2023).$

Figure 4

The closest four grid points to the IPP coordinates ( λ 0 , ϕ 0 ) in each height level (Mahbuby and Amerian 2023).

In each subinterval, the covariance matrix of the observations R is obtained by inversion of the weight matrix P , which is constructed according to the elevation angle of the satellites ε r s (Holschneider et al. 2003):

(24) P = diag ( sin 2 ( ε r s ) ) , R = P − 1 .

The weights of the background data obtained from IRI are considered to be the same in each subinterval. However, a proper variance factor σ B is needed to establish an appropriate connection between background and observation variances in data assimilation, as follows:

(25) B = σ B I .

In order to find the optimal variance factor σ B ⁎ and proper covariance matrix of background in each subinterval, the following optimization problem is solved:

(26) σ B ⁎ = arg min ( ∥ V R ( λ o , ϕ o ) − H N A ( λ i , j , ϕ i , j , h k ) ∥ 2 + ∥ N A ( λ i , j , ϕ i , j , h k ) − N B ( λ i , j , ϕ i , j , h k ) ∥ 1 ) s .t . σ B > 0 .

The constructed analysis grids via the data assimilation process are divided into two categories. The first category is obtained from data assimilation at reference epochs and is applied as the initial conditions for prediction by ICE. The second category includes the analysis grids at other epochs. These grids are utilized to evaluate the quality and accuracy of the predictions.

2.3 Parameter estimation and solution

To solve ICE, some simplifications and assumptions were considered as can be seen in equation (17). For this purpose, the following assumptions, which are consistent with the reality of the ionosphere, have been considered:

At each reference epoch and a step before that, accurate IED values on a 3D regular spherical grid (analysis grid) are known. These values are important for two reasons. First, they are used as the initial condition for integration, and second, due to their proper accuracy, they can also be used to estimate the parameters of ICE.
In a short period of time, the electric and magnetic fields in plasma are stable, and the divergence term, which is related to the electron transport in ICE, can be assumed to be a constant (Li et al. 2020).
In different regions of the ionosphere, electron transport occurs in the horizontal direction, and there is no electron transport between regions in the vertical direction. Also, the constant value κ 4 , which indicates the electron transport rate, is different in various regions because the friction is different.
Due to the insignificance of IED values in the D region, we avoid its prediction and consider that the constants κ 1 , κ 2 , κ 3 , and κ 4 are different in four regions E, F1, F2, and topside. However, their values along each region remain constant in each 10-min interval.

Now, considering equation (17) and assumption Nos 3 and 4, ICE for each of the four regions E, F1, F2, and topside above F2, at the moment t can be written as follows:

(27) ∂ N e ( λ i , j , ϕ i , j , h k , t ) ∂ t = κ 1 I ( λ i , j , ϕ i , j , t ) − κ 2 N e 2 ( λ i , j , ϕ i , j , h k , t ) − κ 3 N e ( λ i , j , ϕ i , j , h k , t ) − κ 4 .

For each region, the loose rate parameters ( κ 2 , κ 3 ) are estimated once by linear regression for a period of 3 h. However, the production parameter ( κ 1 ) and the transport parameter ( κ 4 ) are constant for each region of the ionosphere in 10-min intervals. In other words, in order to predict in each time step with equation (27), these parameters must be estimated for each of the four regions. Considering assumption No. 1 and knowing the variations of IED in two consecutive epochs from the analysis grids, i.e., the reference epoch t 0 and a step before t 0 − 1 , the following equation can be written to estimate the parameters in each of the stated regions in assumption No. 4:

(28) N e ( λ i , j , ϕ i , j , h k , t 0 ) − N e ( λ i , j , ϕ i , j , h k , t 0 − 1 ) Δ t = κ 1 I ( λ i , j , ϕ i , j , t 0 − 1 ) − κ 2 N e 2 ( λ i , j , ϕ i , j , h k , t 0 − 1 ) − κ 3 N e ( λ i , j , ϕ i , j , h k , t 0 − 1 ) − κ 4 .

Considering the known values of electron densities from analysis grids and solar radiation flux, equation (28) becomes a system of linear equations in each of the four regions, which has four unknowns in the day and three unknowns in the night. The unknown parameters κ 1 , κ 2 , κ 3 , and κ 4 can be estimated via a least-squares regression. It should be noted that in each 10-min subinterval, κ 1 and κ 4 values should be upgraded because the solar radiation flux on the horizon plane and electron transport both vary in short intervals and must be estimated again by predicted values. It is necessary to remember that κ 1 is zero in the night and only κ 4 should be upgraded.

Now, the necessary parameters for prediction are obtained, and utilizing the initial conditions for equation (27), it is possible to predict electron densities for the next steps. In order to estimate the unknown parameters, a constrained least-squares process is required. Because κ 1 is a percentage of the solar radiation flux and there is a condition for it to be positive:

(29) κ = ( κ 1 , κ 2 , κ 3 , κ 4 ) T = arg min ( ∥ C κ κ − d κ ∥ 2 2 ) s .t . κ 1 ≥ 0 .

In this equation, d κ stands for the vector of constant values, which contains the differences between IED values in the two initial epochs divided by Δ t (the left-hand size of equation (28)). Here, C κ indicates the design matrix constructed by electron densities and solar radiation flux on the horizon plane at t 0 − 1 (the right-hand size of equation (28)):

(30) C k = [ I − N e 2 − N e − 1 ] Day [ − N e 2 − N e − 1 ] Night .

The constrained optimization problem in equation (29), which has an inequality constraint, is solved using the interior point algorithm, i.e., first, the inequality constraint is added to the main objective function in the form of a logarithmic barrier function, and then the problem becomes unconstrained, and we solve the equivalent unconstrained problem using the gradient and hessian of its objective function with Newton’s method (Boyd and Vandenberghe 2004).

After the estimation of the parameters by solving the optimization problem stated in equation (29), the obtained parameters in each region are supposed to be used for the prediction of IED values in that region by equation (27). This ordinary differential equation is solved by the Runge–Kutta method of the fourth order. For this purpose, equation (27) is considered in the following form:

(31) ∂ N e ∂ t = f ( t , N e ) .

Also, initial conditions are known from the analysis grid at t 0 , i.e., N e ( t 0 ) = N e 0 . Hence, equation (31) can be solved for each 10-min time step Δ t to predict N e ( t n + 1 ) = N e n + 1 knowing N e ( t n ) = N e n . Therefore, a Runge–Kutta method of the fourth order is applied as follows (Press 2007, Süli and Mayers 2003):

(32) N e n + 1 = N e n + Δ t 6 ( μ 1 + 2 μ 2 + 2 μ 3 + μ 4 ) , μ 1 = f ( t n , N e n ) , μ 2 = f t n + Δ t 2 , N e n + Δ t μ 1 2 , μ 3 = f t n + Δ t 2 , N e n + Δ t μ 2 2 , μ 4 = f ( t n + Δ t , N e n + Δ t μ 3 ) .

2.4 Data processing scheme

Figure 5 represents the whole data processing stages for prediction of IED using ICE.

Figure 5

Data processing scheme.

3 Validation

If electron densities in different regions of the ionosphere are predicted, VTEC can also be estimated

(33) VTEC = ∫ h = 0 ∞ N e d h .

In practice, the integration limits can be considered the same as the analysis grids, i.e., from 100 to 2,000 km above the Earth’s surface:

(34) VTEC ( λ , ϕ , t ) = ∑ k = 1 K N e ( λ , ϕ , h k , t ) Δ h k .

In this study, the predicted IED grids have been evaluated in two ways. First, the analysis grid is pre-made and available with high accuracy, and it can be used to evaluate the accuracy of the predictions in 1, 2, and 3-h intervals.

Second, according to equation (34), it is possible to obtain the predicted VTECs from the predicted electron densities. Therefore, a group of GPS-derived VTECs from the main dataset is separated in each 1-h interval, and we call that test data. These data are not used in data assimilation and construction of the analysis grids. The test data are only used to evaluate the predicted VTECs, which are referred to in equation (34).

4 Numerical experiment

The represented method is applied for IED and VTEC prediction for May 8 (DOY 129), 2016, which was the day with the most ionospheric activity and severe geomagnetic storm in that year.

As stated in Section ‎2.2, two types of data have been used in this study: VTECs obtained from permanent dual-frequency GPS receivers in the study region and IED grids obtained from IRI. The study region and the distribution of GPS stations are depicted in Figure 6 (left). Also, in Figure 6 (right), the horizontal distribution of the received IED grids from IRI, which are applied in data assimilation, is shown by red dots. The location of IPPs at the sample interval of [4, 5) UT for which VTEC is retrieved is depicted by blue dots. Data assimilation has been done in 10-min subintervals, and GPS-derived VTECs in each subinterval are considered observations and assimilated into the IED grids of the background. The outputs of data assimilation are 3D IED grids, which are called analysis and are supposed to be used as the initial conditions in prediction by ICE. Also, the initial conditions obtained from analysis grids, as mentioned in Section ‎2.3, are utilized for the estimation of the necessary parameters for prediction by ICE. The accuracy of the constructed analysis grids is investigated by integration over the zenith direction and comparison with test data that contain VTECs. In the study conducted by Mahbuby and Amerian (2023), the accuracy of the resulted VTECs from the analysis grids is claimed to be 0.9 TECU, and also, it is shown that the accuracy of single-frequency positioning by corrected observations for ionospheric delay via the constructed analysis grids is improved about several decimeters at the control station indicated by the blue star on the left panel of Figure 6.

Figure 6

The distribution of permanent dual-frequency GPS receivers and control station (blue star) in the study region (left), and observations and background used for data assimilation to construct the analysis grids (right).

In this study, according to assumption number 2 in Section ‎2.3, ICE is converted to ODE. It is clear that using smaller time steps for prediction results in better accuracy. However, it should be noted that this requires having the analysis grid in the two initial moments with the same selected step. The requirement for providing such an analysis grid is to have more permanent GNSS stations in the region, and also, the data assimilation process should be done in these smaller intervals. This issue can lead to more computational complexity and cost. Therefore, the following argument has been made to investigate the appropriateness of the 10-min step for discretization interval. Here, the analysis grids are applied, and the relative variations of IDE are computed at the level where peak analysis values occur. The function of relative variations is defined as follows:

(35) ρ ( Δ t ) = Δ N e N e = N e ( t 0 + Δ t ) − N e ( t 0 ) N e ( t 0 ) ,

and various time steps are considered as Δ t = 10 , 20 , 30 , 40 , 50 min. The function is drawn for these time steps between 10 UT and 19 UT when IED values are higher than at other times in the study region (Figure 7).

Figure 7

The relative variations of IED for various time steps.

As can be seen from the figure, the amplitude of the relative variations decreases as the time step increases. The amplitude values of variations for time steps of 20 and 10 min are more similar. According to this argument, the time step of 10 min includes the dominant fluctuations, and there is no need to reduce it further. Also, the energy of ρ ( Δ t ) is plotted for the mentioned time steps in Figure 8. The energy is defined as the sum of the squares of the relative variations for each step. We can see that although the energy function decreases as the time step is reduced, there is no significant difference in energy between 20- and 10-min steps, and the slope of the graph drops. With this analysis, it can be said that a 10-min step is appropriate and further reduction may not be cost-effective due to the increase in complexity and time of calculating the analysis grid.

Figure 8

The energy of the relative variations of IED for various time steps.

The estimated parameters of ICE for disappearance rate ( κ 2 , κ 3 ) are considered to be valid for 3 h, but the parameters related to electron production ( κ 1 ) and electron transport rate ( κ 4 ) are upgraded every 10 min. Utilizing these parameters, electron densities are predicted for 1, 2, and 3-h intervals. For sample epochs t 0 = 00 : 10 UT and t 0 = 10 : 10 UT , initial conditions are applied, and the predicted values are compared with the analysis values. The results of this comparison are shown in Figures 9 and 10. These figures contain three rows and three columns. The left panel shows the analysis values for 49,569 grid points, and the middle panel shows the predicted values with ICE for the same grid points. In the right panel, the predicted and analysis values are compared, and the correlation coefficients are indicated in the figures of this panel. The first row of these figures compares 1-h predictions with the analysis, and the second and third rows compare 2 and 3-h predictions with their corresponding analysis values at the same moment. As expected, the accuracy of predictions has decreased with the increase in the prediction period. In Figure 9, the first row, an RMS of 1.20 × 10¹⁰ el/m³, is obtained for the differences between predicted and analysis values at 1 UT, while in the second and third rows, RMS of the differences are 3.86 × 10¹⁰ el/m³ and 7.77 × 10¹⁰ el/m³ at 2 UT and 3 UT, respectively. Similarly, in Figure 10, the first, second, and third rows, RMS values of 2.64 × 10¹⁰ el/m³, 4.75 × 10¹⁰ el/m³, and 7.67 × 10¹⁰ el/m³ are obtained for the differences between predictions and analysis values at 11 UT, 12 UT, and 13 UT. As it is clear, in the right panel of both figures, the correlation coefficients between predicted and analysis values decrease from top to bottom as the prediction period increases.

Figure 9

The electron densities of the analysis grid (left), the predicted electron densities (middle), and the correlation between analysis and predicted values (right) with initial conditions at t ₀ = 00:10 UT for 1-h (first row), 2-h (second row), and 3-h predictions (third row).

Figure 10

The electron densities of the analysis grid (left), the predicted electron densities (middle), and the correlation between analysis and predicted values (right) with initial conditions at t ₀ = 10:10 UT for 1-h (first row), 2-h (second row), and 3-h predictions (third row).

As can be seen from Figure 9, the predicted values for 1 UT, 2 UT, and 3 UT by applying the initial conditions at 0:10 UT have a slightly negative bias and are underestimated, but compared to them, the predicted IED values for 11 UT, 12 UT, and 13 UT by applying the initial conditions at 10:10 UT are more accurate and without any significant bias. Accordingly, the slope of the trend line is close to one, and the correlation coefficients between prediction and analysis are higher (Figure 10).

Applying the initial conditions at reference epochs in equation (21), for various intervals predictions are provided, and RMS values of the differences between predicted and analysis values during the study period are investigated in Figure 11. In this figure, RMS values of the differences for 20, 30, and 40-min intervals are represented in the left panel and RMS values corresponding to 1, 2, and 3-h predictions are shown in the right panel. On average, the RMS values of the differences for 20, 30, and 40-min prediction intervals during the study day are 0.87 × 10¹⁰ el/m³, 1.15 × 10¹⁰ el/m³, and 1.41 × 10¹⁰ el/m³. These daily average values for 1, 2, and 3-h predictions are 1.66 × 10¹⁰ el/m³, 3.08 × 10¹⁰ el/m³, and 5.33 × 10¹⁰ el/m³. If real-time data assimilation is performed and there is an analysis grid at each epoch, the RMS of the differences for 10-min predictions would be smaller, and the daily average reaches 0.66 × 10¹⁰ el/m³.

Figure 11

Average RMS of the differences between predicted IED and analysis in 20, 30, and 40-min intervals (left), and average RMS of the differences for 1, 2, and 3-h intervals (right).

As stated in Section ‎3, interpolation and integration of the predicted electron densities results in predicted VTEC values. Hence, another evaluation method is considered for predictions. For this purpose, about 10% of the GPS-derived VTECs at IPPs are excluded from the observation dataset as test data and did not play a role in data assimilation and construction of the analysis grids. According to equation (23), predicted electron densities can be spatially interpolated and integrated over the zenith direction of IPP to provide predicted VTEC values. RMS of the differences between these predicted VTECs and test data is considered another criterion to evaluate the accuracy of predictions. These values are shown in Figure 12 for the same intervals as mentioned in Figure 11. On average, RMS values of the differences between 20, 30, and 40-min predicted VTECs and test data during the study day are 0.77 TECU, 0.99 TECU, and 1.15 TECU. These daily average values for 1, 2, and 3-h predictions are 1.29 TECU, 1.95 TECU, and 2.57 TECU. As it is clear from these average values and Figures 11 and 12, the shorter the time intervals of the predictions, the better their accuracy. This issue is due to being closer to the reference epoch, in which data assimilation has been done, and the exact initial conditions have been introduced to the continuity equation. If real-time data assimilation is performed and prediction is done using an analysis grid at each epoch for the next 10 min, RMS of the differences with test data would be smaller, and the daily average would reach 0.58 TECU. The maximum and minimum of RMS of the differences between 10-min predictions and test data are 1.2 TECU and 0.2 TECU, respectively.

Figure 12

Average RMS of the differences between predicted VTECs and test data in 20, 30, and 40-min intervals (left), and average RMS of the differences for 1, 2, and 3-h intervals (right).

One important factor that affects the accuracy of predictions is the magnitude of ionospheric activities, which can be described with the 3h-Kp index. The study period is selected on the most active ionosphere day of 2016, and 3h-Kp indices fluctuate between 4.7 and 6.3. In Figure 13, it is possible to see the values of 3h-Kp indices during the study period. Also, in order to investigate the relationship between the accuracy of the predictions and geomagnetic activities of the ionosphere, the RMS of the differences for 1, 2, and 3-h predictions have also been plotted. RMS of the differences between predicted VTECs and obtained VTECs from analysis are shown in the left panel, and RMS values of the differences between predicted IEDs and analysis are shown in the right panel. As can be seen in the figure, between 9 UT to 11 UT, where, according to Kp indices, the ionospheric activities have decreased, RMS values also experience a sudden drop and the accuracy of the predictions has increased. This is true for both IED and VTEC predictions.

Figure 13

3h-Kp indices’ variations during the study period compared to RMS of the differences between predicted VTECs and obtained VTECs from analysis (left) and RMS of the differences between predicted IED and analysis values (right).

5 Discussion and conclusion remarks

Prediction of IED and VTEC variations is important for ionospheric studies. In this regard, physics-based models can play a prominent role. In this cooperation, ICE that interprets the processes of production, disappearance, and transfer of electrons was used to predict IED and subsequently VTEC. Electron production rate during the day was considered as a percentage of the solar radiation flux received on the horizon plane and zero at night. The rate of electron disappearance due to the CE, attachment, and recombination processes was assumed to be a quadratic function of the electron density value. The predictions were made with a step of 10 min. Hence, the electron transfer rate was considered constant in this interval due to the stable electric and magnetic fields of plasma. In this study, the ionosphere is divided into four regions from bottom to top, and according to the assumptions, there are four unknown parameters in each region. In order to make a prediction by ICE, a 3D regular and accurate grid of electron densities is required.

GPS-derived VTECs are 2D at each epoch and allocated to a specific altitude. They do not show the distribution of electrons in the ionosphere. Also, they do not have a regular geometric distribution. It is possible to make a regular 3D grid of IED in any region from empirical ionospheric models such as IRI. Although it is not expected to be accurate everywhere, it has the ability to show the distribution of electrons in the layers of the ionosphere and is a 3D data source at each epoch. These IED values are used as the background in data assimilation, but its relative weight is less than GPS-derived VTECs. The optimal combination of these two sets of data can lead to 3D grids with higher accuracy than the empirical model. This grid is called analysis.

It is clear that the initial conditions play a very important role in the numerical solution of ICE in the form of equation (27). The accuracy of the background and GPS-derived VTECs and the weights used in the data assimilation process have an impressive effect on the accuracy of the initial conditions. More accurate initial conditions will lead to higher accuracy of the numerical predictions.

We used the information from the analysis grid at reference epochs and a time step before that to estimate the unknown parameters in each region. This task requires solving an inequality constrained linear least-squares problem, which was done by the interior point method. Also, the analysis grid at each reference epoch provides the initial conditions of the differential equation. Therefore, having accurate initial conditions from analysis makes it possible to make predictions of IED and VTEC with acceptable accuracy for up to 3 h with the stated method. For predicted VTECs, accuracies between 0.2 TECU and 4 TECU and, for predicted IED values, accuracies between 0.3 × 10¹⁰ el/m³ and 5.8 × 10¹⁰ el/m³ are obtained, and these accuracies are dependent on the duration of the predictions and the intensity of geomagnetic activities. As the prediction interval decreases, the accuracy increases. This is due to the fact that being closer to the epoch of initial conditions in which data assimilation is done will improve the prediction accuracy. In this regard, if data assimilation is done in real-time and the analysis grid as accurate initial conditions is prepared at each epoch, high accuracy 10-min predictions can be provided. Also, the intensity of the geomagnetic activities of the earth almost inversely affects the accuracy of the predictions. The obtained accuracies in this study are within the range of the reported accuracies by other methods, e.g., other physics-based and ML techniques. In comparison, it can be said that in terms of computational complexity, cost, and time, the presented method is efficient because having the analysis grids in only the first two epochs, it can predict both IED and VTEC. However, it should be noted that ML methods are able to achieve longer-term predictions by increasing the number of inputs and training period. Therefore, a suggestion can be made to use ICE as a loss function in ML algorithms and take advantage of the benefits of both methods.

Acknowledgments

The authors acknowledge the National Cartographic Center of Iran (NCC) for providing Iranian Permanent GPS Network (IPGN) data, the Center for Orbit Determination in Europe (CODE) for providing GIMs and satellites’ DCBs.

Funding information: No funds, grants, or other supports were received for this study.
Author contributions: Design the research, H.M., and Y.A.; Analysis, H.M.; Writing, H.M.; Review and editing, Y.A. All authors have read and agreed to the published version of the manuscript.
Conflict of interest: The authors have no competing interests to declare that are relevant to the content of this article.
Data availability statement: The data that support the findings of this study are available from the corresponding author upon reasonable request. Satellites’ DCBs, GIMs in IONEX format are all published by CODE and publicly available from: http://ftp.aiub.unibe.ch/CODE/.

References

Arikan, F., H. Nayir, U. Sezen, and O. Arikan. 2008. “Estimation of single station interfrequency receiver bias using GPS-TEC.” Radio Science 43(4), 1–13.10.1029/2007RS003785Search in Google Scholar

Bilitza, D., D. Altadill, V. Truhlik, V. Shubin, I. Galkin, B. Reinisch, and X. Huang. 2017. “International Reference Ionosphere 2016: From ionospheric climate to real‐time weather predictions.” Space Weather 15(2), 418–29.10.1002/2016SW001593Search in Google Scholar

Biondi, M. A. 1969. “Atmospheric electron–ion and ion–ion recombination processes.” Canadian Journal of Chemistry 47(10), 1711–9.10.1139/v69-282Search in Google Scholar

Boyd, S. and L. Vandenberghe. 2004. Convex optimization. Cambridge, United Kingdom: Cambridge University Press.10.1017/CBO9780511804441Search in Google Scholar

Chakraborty, S. and T. Basak. 2020. Numerical analysis of electron density and response time delay during solar flares in mid-latitudinal lower ionosphere. Astrophysics and Space Science 365(12), 184.10.1007/s10509-020-03903-5Search in Google Scholar

Chapman, S. 1931. “The absorption and dissociative or ionizing effect of monochromatic radiation in an atmosphere on a rotating earth part II. Grazing incidence.” Proceedings of the Physical Society 43(5), 483.10.1088/0959-5309/43/5/302Search in Google Scholar

Chen, Z., B. An, W. Liao, Y. Wang, R. Tang, J. Wang, and X. Deng. 2023. “Ionospheric electron density model by electron density grid deep neural network (EDG-DNN).” Atmosphere 14(5), 810.10.3390/atmos14050810Search in Google Scholar

Chen, Z., W. Liao, H. Li, J. Wang, X. Deng, and S. Hong. 2022. “Prediction of global ionospheric TEC based on deep learning.” Space Weather 20(4), e2021SW002854.10.1029/2021SW002854Search in Google Scholar

Ciraolo, L., F. Azpilicueta, C. Brunini, A. Meza, and S. M. Radicella. 2007. “Calibration errors on experimental slant total electron content (TEC) determined with GPS.” Journal of Geodesy 81(2), 111–20.10.1007/s00190-006-0093-1Search in Google Scholar

Dogariu, A., M. N. Shneider, and R. B. Miles. 2010. “Direct measurement of electron loss rate in air.” In Proceedings CLEO/QELS: 2010 Laser Science to Photonic Applications, IEEE, pp. 1–2.10.1364/QELS.2010.JTuD2Search in Google Scholar

Dutta, S. and M. Cohen. 2022. “Improving electron density predictions in the topside of the ionosphere using machine learning on in situ satellite data.” Space Weather 20(9), e2022SW003134.10.1029/2022SW003134Search in Google Scholar

Holschneider, M., A. Chambodut, and M. Mandea. 2003. “From global to regional analysis of the magnetic field on the sphere using wavelet frames.” Physics of the Earth and Planetary Interiors 135(2), 107–24.10.1016/S0031-9201(02)00210-8Search in Google Scholar

Hsu, C. T., T. Matsuo, W. Wang, and J. Y. Liu. 2014. “Effects of inferring unobserved thermospheric and ionospheric state variables by using an Ensemble Kalman Filter on global ionospheric specification and forecasting.” Journal of Geophysical Research: Space Physics 119(11), 9256–67.10.1002/2014JA020390Search in Google Scholar

Huang, L., H. Wu, Y. Lou, H. Zhang, L. Liu, and L. Huang. 2022. “Spatiotemporal analysis of regional ionospheric TEC prediction using multi-factor neuralprophet model under disturbed conditions.” Remote Sensing 15(1), 195.10.3390/rs15010195Search in Google Scholar

Ivanov-Kholodny, G. S. and A. V. Mikhailov. 2012. The prediction of ionospheric conditions. Dordrecht, Holland: Springer Science & Business Media.Search in Google Scholar

Jee, G. 2023. “Fundamentals of numerical modeling of the mid-latitude ionosphere.” JASS 40(1), 11–18.10.5140/JASS.2023.40.1.11Search in Google Scholar

Lei, D., H. Liu, H. Le, J. Huang, J. Yuan, L. Li, and Y. Wang. 2022. “Ionospheric TEC prediction base on attentional BiGRU.” Atmosphere 13(7), 1039.10.3390/atmos13071039Search in Google Scholar

Li, J., D. Huang, Y. Wang, Y. Zhao, and A. Hassan. 2020. “A new model for total electron content based on ionospheric continuity equation.” ASR 66(4), 911–31.10.1016/j.asr.2020.04.048Search in Google Scholar

Liu, L., Y. J. Morton, and Y. Liu. 2022. “ML prediction of global ionospheric TEC maps.” Space Weather 20(9), e2022SW003135.10.1029/2022SW003135Search in Google Scholar

Mahbuby, H. and Y. Amerian. 2021. “Regional assimilation of GPS-derived TEC into GIMs.” Pure and Applied Geophysics 178, 1317–37.10.1007/s00024-021-02681-7Search in Google Scholar

Mahbuby, H. and Y. Amerian. 2023. “Regional ionospheric electron density modeling by assimilation of GPS-derived TEC into IRI-provided grids on May 8, 2016.” Advances in Space Research 72(6), 2377–90.10.1016/j.asr.2023.05.043Search in Google Scholar

Memarzadeh, Y. 2009. Ionospheric modeling for precise GNSS applications. Delft, Netherlands: Citeseer.10.54419/z48womSearch in Google Scholar

Press, W. H. 2007. Numerical recipes 3rd edition: The art of scientific computing. Cambridge, United Kingdom: Cambridge University Press.Search in Google Scholar

Prölss, G. 2012. Physics of the Earth’s space environment: an introduction. Berlin, Germany: Springer Science & Business Media.Search in Google Scholar

Razin, M. R. G. and B. Voosoghi. 2016. “Wavelet neural networks using particle swarm optimization training in modeling regional ionospheric total electron content.” JASTP 149, 21–30.10.1016/j.jastp.2016.09.005Search in Google Scholar

Seaton, S. 1947. “Temperature and recombination coefficient in the ionosphere.” Journal of the Atmospheric Sciences 4(6), 197–200.10.1175/1520-0469(1947)004<0197:TARCIT>2.0.CO;2Search in Google Scholar

Shi, S., K. Zhang, S. Wu, J. Shi, A. Hu, H. Wu, and Y. Li. 2022. “An investigation of ionospheric TEC prediction maps over China using bidirectional long short‐term memory method.” Space Weather 20(6), e2022SW003103.10.1029/2022SW003103Search in Google Scholar

Shimazaki, T. 1965. “On the continuity equation for electron density in the ionosphere.” JASTP 27(5), 593–604.10.1016/0021-9169(65)90127-3Search in Google Scholar

Süli, E., and D. F. Mayers. 2003. An introduction to numerical analysis. Cambridge, United Kingdom: Cambridge University Press.10.1017/CBO9780511801181Search in Google Scholar

Wen, Z., S. Li, L. Li, B. Wu, and J. Fu. 2021. “Ionospheric TEC prediction using Long Short-Term Memory deep learning network.” Astrophysics and Space Science 366, 1–11.10.1007/s10509-020-03907-1Search in Google Scholar

Yang, D. and H. Fang. 2023. “A low‐latitude three‐dimensional ionospheric electron density model based on radio occultation data using artificial neural networks with prior knowledge.” Space Weather 21(1), e2022SW003299.10.1029/2022SW003299Search in Google Scholar

Zhang, Y., X. Wu, and X. Hu. 2018. “Effects of estimating the ionospheric and thermospheric parameters on electron density forecasts.” Science China Earth Sciences 61, 1875–87.10.1007/s11430-017-9251-4Search in Google Scholar

Žigman, V., M. Dominique, D. Grubor, C. J. Rodger, and M. A. Clilverd. 2023. “Lower-ionosphere electron density and effective recombination coefficients from multi-instrument space observations and ground VLF measurements during solar flares.” Journal of Atmospheric and Solar-Terrestrial Physics 247, 106074.10.1016/j.jastp.2023.106074Search in Google Scholar

Received: 2023-12-05

Revised: 2024-03-02

Accepted: 2024-04-17

Published Online: 2025-01-24

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/jogs-2022-0175

Keywords for this article

ionospheric continuity equation; ionospheric data assimilation; ionospheric electron density; ionospheric prediction; vertical total electron content

Creative Commons

BY 4.0