© P. Mungufeni et al., Published by EDP Sciences 2022

1 Introduction

The peak electron density in the F2-region (NmF2) and total electron content (TEC) are widely used parameters to characterize the ionosphere (Rishbeth & Garriott, 1969; Gerzen et al., 2013). The NmF2 affects high‐frequency (3–30 MHz) radio wave communication applications. For instance, Geeta and Yadav (2014) stated that for frequencies lower than 30 MHz, the ionosphere acts as a helpful aid for radio wave propagation, but for the frequencies slightly greater than 30 MHz, it may cause attenuation. The TEC affects Global Navigation Satellite System (GNSS) based positioning by introducing ionospheric refraction (Hofmann-Wellenhof et al., 2007) in which the code delay or carrier phase advance occurs, resulting in a pseudo-range measurement instead of a true range measurement.

Several research efforts have tried to model NmF2 and TEC. For instance, the Committee on Space Research and the International Union of Radio Science formed a working group in the late sixties to produce an empirical standard model of the ionosphere, the International Reference Ionosphere (IRI), based on all available data sources (Bilitza et al., 1993). The output of the IRI includes NmF2 and TEC, among other parameters. The NeQuick is another global ionospheric model that can be used to obtain NmF2 and TEC. The NeQuick model and its subsequent modifications (NeQuick G and NeQuick 2) are a three-dimensional, time-dependent ionospheric electron density model developed by the Abdus Salam International Center for Theoretical Physics (ICTP) in Trieste, Italy, and the Institute for Geophysics, Astrophysics, and Meteorology of the University of Graz, Austria (Nava et al., 2008, and references therein).

In order to produce a good empirical model of the ionosphere over a region, there is a need to have extensive observations to be used for constructing the model. It may be difficult to obtain adequate observations that can be used in modeling since maintaining stable and reliable ground-based instruments all over regions may be expensive or not practicable, particularly over the seas and deserts. It is now known that estimations of ionospheric parameters by empirical and theoretical models may be improved by assimilating observations to the model.

Bust et al. (2004) presented the Ionospheric Data Assimilation Three-Dimensional (IDA3D) algorithm, which uses a three-dimensional variational data assimilation technique (3DVAR). The IDA3D is capable of incorporating most electron density-related measurements, including GNSS-TEC measurements, low-Earth-orbiting beacon TEC, and electron density measurements from radars and satellites. Bust and Datta-Barua, (2013) stated that Ionospheric Data Assimilation Four-Dimensional (IDA4D) is an ionospheric data assimilation algorithm that provides global 3-D time-evolving maps of ionospheric electron density. A computationally practical data assimilation technique known as Best Linear Unbiased Estimator (BLUE) has been implemented by Angling and Cannon (2004) for combining RO data with background ionospheric models. Recently, Mengist et al. (2019), investigated the IDA4D technique over Korea and the neighboring areas, considering the IRI model as the background. They showed that assimilation of ground-based global positioning system (GPS) Slant TEC (STEC), NmF2 obtained from Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) radio occultation (RO), and ionosondes yielded the best results. Ssessanga et al. (2019) presented a preliminary study that assessed the capability of their developed four-dimensional (in space and time) data assimilation scheme to more accurately estimate the 3-D picture of the ionosphere over the South African region. Yue et al. (2007) assimilated electron densities observed by the Millstone Hill incoherent scatter radar (ISR) into a one-dimensional midlatitude ionospheric theoretical model by using an ensemble Kalman filter (EnKF) technique.

The current study performed assimilation of TEC derived from ground-based GPS receiver and ionosonde stations as well as COSMIC RO to the NeQuick model. Moreover, the study assimilated the TEC derived from GPS receivers onboard Low Earth Orbiting (LEO) Swarm and COSMIC satellites. The classical Kalman filtering technique applied to assimilate TEC observations to NeQuick model and the successive corrections method used to correct the NeQuick model generated ionospheric parameters at locations and epochs that lacked TEC observations make this study unique. After describing in Section 2 the data used in this work, we present in Sections 3 and 4 the classical Kalman filtering technique and the successive corrections method, respectively. The results and discussions are presented in Section 5, while the conclusions are presented in Section 6.

2 Data sources

During geomagnetic storms, the variations in zonal electric fields and composition of the neutral atmosphere contribute significantly to the occurrence of negative and positive ionospheric storm effects (Rishbeth & Garriott, 1969; Buonsanto, 1999). In order to ascertain the ability of the assimilation results to emulate the complex ionospheric changes during geomagnetic storms, the assimilations were done during March 15–20, June 4, 21–26, October 6–11, 27, and December 18–24, 2015. In these four different months, representing seasons of March equinox, June solstice, September equinox, and December solstice, the recorded minimum Dst were −222, −204, −124, and −155 nT, respectively. For purposes of comparing ionospheric parameters on quiet and disturbed geomagnetic conditions, there was one quiet day (Kp ≤ 3) in each of the selected months. The selected geomagnetically quiet (Kp ≤ 3) days were March 15, June 4, October 27, and December 18. Figure 1 presents the variation of Dst (panel (a)’s) and Kp indices (panel (b)’s) during the period under study. The hourly Dst and 3 hourly Kp indices can be obtained from the World Data Center of Kyoto, Japan (http://swdcwww.kugi.kyoto-u.ac.jp/).

Figure 1 Variation of (a) Dst and (b) Kp during (1) March 15–20, 2015, (2) June 4, 21–26, 2015, (3) October 6–11, 27, 2015, and (4) December 18–24, 2015.

The pattern of the variation of Dst shown in Figure 1 indicates several main and recovery phases of geomagnetic storms during the period mentioned. The high Kp index values (Kp > 5) on March 17–18, June 22–23, October 7–8, and December 20–21 confirm the occurrence of geomagnetic disturbances during the period. Therefore, the ionosphere might vary greatly in the period under study.

Most global climatological ionospheric models such as NeQuick and IRI might not emulate possible rapid variations of the ionosphere due to geomagnetic storms. In order to improve the estimation of electron density up to the altitude of GPS satellites (~20,200 km) during the previously stated disturbed geomagnetic periods, we performed data assimilation to NeQuick model driven by daily F10.7 values. It is important to mention that, unlike NeQuick model, which can estimate electron density up to the altitude of GPS satellites, IRI only estimates it up to an altitude of 2000 km. Moreover, since NeQuick is a quick-run model, it is suitable for the data assimilation process as computation time might be reduced.

Figure 2 shows with blue crosses and red diamonds the locations of the ground-based GPS receivers and ionosonde stations, respectively, used to obtain the TEC data. The stations indicated in the figure had available data during most of the days considered in the study period. Table 1 provides the geophysical parameters (e.g., geographic and geomagnetic coordinates) associated with the stations in Figure 2. The Receiver Independent Exchange (RINEX) data format files of GPS receivers can be obtained from the University NAVSTAR Consortium network (ftp://data‐out.unavco.org). The RINEX files were processed using the software described in Ciraolo et al. (2007), yielding TEC together with other parameters such as time, elevation and azimuth angles, and geographic longitude and latitude of the ionospheric pierce points. The software can produce TEC data at resolutions of 30 s, 1, 5, 10, and 15 min. This study considered 5 minutes’ resolution data to reduce the computational burden of working with high-resolution TEC data.

Figure 2 Locations of GPS (blue crosses) and ionosonde (red diamonds) stations used in the study.

The TEC derived by integrating electron density profiles obtained from ionosonde stations listed in Table 1 can be accessed from the National Oceanic and Atmospheric Administration (NOAA) website via the link, ftp://ftp.ngdc.noaa.gov. The data obtained from the NOAA website was in the form of auto-scaled ionospheric parameters such as peak height in F2-region, foF2, and TEC, which are stored in Standard Archiving Output (SAO) format files. The TEC data provided in SAO files have a resolution of 15 min and are obtained by integrating electron density profiles up to an altitude of ~700 km. Reinisch and Huang (2001) stated that the auto-scaling program (real-time ionogram scaler with true height (ARTIST)) approximates the electron density profile above the F2 layer peak by an α-Chapman function with a constant scale height that is derived from the bottom-side profile shape near the F2 peak. The ionospheric parameters in SAO files are associated with flags showing confidence level, which ranges from low (11) to high (55). The foF2 and TEC data obtained from the SAO files used in this study were those with high confidence levels (≤22).

As mentioned in Section 1, TEC measurements obtained from the GPS receivers onboard Swarm and COSMIC satellites were also used. The Swarm constellation is composed of three identical satellites, namely, Alpha (A), Bravo (B), and Charlie (C). Detailed information about orbital characteristics of the Swarm satellites can be found in Zakharenkova and Astafyeva (2015). Each Swarm spacecraft carries a Precision Orbit Determination (POD) antenna. The GPS signal phase measurements obtained from this antenna can estimate the line of sight TEC between Swarm and GPS satellites. This line of sight TEC data can be freely downloaded from the European Space Agency (ESA) website (http://www.earth.esa.int/swarm). Its resolution was 1 and 10 s after and before July 15, 2014, respectively. The COSMIC TEC data, which is based on POD antenna measurements, is processed at 1-second resolution and archived at the COSMIC Data Analysis and Archive Centre (CDAAC) (http://cosmic-io.cosmic.ucar.edu/cdaac/index.html). The line of sight TEC (STEC) obtained from Swarm, and COSMIC satellites were converted to vertical TEC (VTEC) at the position of the LEO satellites as in Zhong et al. (2016).

The TEC data resulting from the integration of electron density profiles associated with COSMIC RO used in this study were also obtained from CDAAC. The integrated electron density (integration being done up to the altitudes of the COSMIC satellites) can be obtained from ionPrf files. The TEC associated with a particular electron density profile was assigned to the geographic coordinate of NmF2 in the same file. The electron density profiles are obtained by the Abel inversion of RO data, assuming local spherical symmetry of the electron density in a large region (a few thousand kilometers in radius) around the ray path tangent points (Krankowski et al., 2011). This assumption may not always be valid, and horizontal ionospheric gradients may significantly affect the retrieved electron density profiles, in particular below the F-layer. In addition, the geographical location of the ray path tangent points at the top and at the bottom of a profile may differ (horizontal smear) by several hundred kilometers. Several studies (e.g., Krankowski et al., 2011; Mengist et al., 2019) that have used COSMIC data commonly consider measurements with horizontal smear >1500 km prone to errors, and they reject such measurements. Typically, over an area bounded by 5° and 4° longitude and latitude ranges, respectively, the total number of COSMIC TEC data obtained in a day may be ~3. Therefore, rejection of data may make this number reduce further. The current study did not reject COSMIC TEC with horizontal smear >1500 km since Mungufeni et al. (2020) analyzed COSMIC TEC data which were coincident with TEC estimated by ionosonde stations over South Africa, finding that, compared to measurements with horizontal smear >1500 km, some measurements with horizontal smear <1500 km were far from the linear least squares fitting line.

3 The classical Kalman filtering technique

In this study, we considered TEC measurement y at time t_k to improve the NeQuick model estimate of the ionospheric electron density profile x along the path that contains electrons that constitute y. We further considered that y is linearly related to x via the equation (Grewal & Andrews, 2001; Angling & Cannon, 2004),

(1)

where H is the measurement sensitivity matrix and w is the measurement noise. If x consists of p electron density values, H will be a row matrix with dimension p. As justified later in this section, elements of H were considered to be vertical grid spans corresponding to electron density values that constitute x. The improved estimate of electron density profile x_a based on the assimilation of measurement y into the NeQuick model profile x_b is given by (Angling & Cannon, 2004),

(2)

where K is the Kalman gain which can be determined as

(3)

In equation (3), B and R are the background and observation covariance matrices, respectively, while superscripts T and −1 denote transpose and inverse of the matrix, respectively.

In general, for n values of y within a period of 15 min and horizontal grid cell having geographic longitude and latitude spans of 5° and 4°, respectively, equations (2) and (3) can be written as

(4)

and

(5)

respectively, where i = 1, 2, …, n.

As in Yue et al. (2007), measurement noise values w_i can be considered as white noise with zero expectation so that R_i in equation (5) is

(6)

where c_r = 0.01 and y_i is the ith observation which is assimilated.

The background error was considered so that p × p matrix B_i consists of elements determined as (Yue et al., 2007),

(7)

where c_b = 0.001, u = 1, 2, 3, …, p, and v = 1, 2, 3, …, p. Following Bust et al. (2004), the distances L_i(u) for scaling down H_i(u) were considered as 20 km in E- and F-regions and 500 km in the plasmasphere. This implies that

(8)

Equation (4) corresponds to a specific grid cell and assimilation time window. This approach reduces the huge demand for computational resources, especially when all the grid cells are considered at once, as mentioned in Rodgers (2000). This study termed the filtering as indicated in equation (4), as classical Kalman filtering due to its similarity to the application of the same in classical mechanics (e.g., prediction of river floods and tracking/navigation of ships and spacecrafts).

It is important to emphasize that y_i can be obtained from any of the five data sources discussed in Section 2. The y values obtained from ground-based GPS receivers were treated similarly to the TEC measurements of ionosonde and COSMIC RO, which are considered vertical. This was done by limiting y values obtained from GPS receivers to slant TEC (STEC) observed at very high elevation angles (>60°). Moreover, the ionospheric pierce points associated with the STEC were restricted to that of the specific spatial grid cell considered. The first advantage of this procedure is the minimization of multipath effects on TEC observations, and the second is rendering simplicity to the assimilation process (y occupies one horizontal grid cell). In the future, we might adopt the method described in Angling and Cannon (2004) to treat low elevation angle TEC observations that cross several horizontal grid cells.

Since the NeQuick model can yield electron density profile for heights starting from about 60 to 20,200 km (approximate GPS satellite altitude), the altitudinal intervals for computing electron densities ( in Eq. (4)) were varied in this range. These altitudinal intervals, which constitute elements of H_i were set based on the known typical vertical electron density profile as follows; in E- and F-regions (<600 km), electron densities were computed at intervals of 10 km, for altitude region of 600–2000 km, the interval was increased to 50 km and above 2000 km, the interval was further increased to 2000 km. Therefore, for the case of ground-based GPS receiver TEC, all the elements of H_i were non-zero.

In order to account for the lack of consideration of the entire electron density profile within the altitude range of 60–20,200 km during the determination of y_i associated with TEC data from Swarm and COSMIC satellites, COSMIC RO, and ionosondes, some elements of H_i were set to zero. For instance, while considering measurements associated with ionosonde and COSMIC RO, elements of H_i corresponding to heights >700 and >800 km, respectively, were set to zero. For the case of Swarm A and C, which fly at altitude ~460 km, elements of H_i corresponding to heights below this altitude were set to zero. While for Swarm B and COSMIC satellites, elements of H_i corresponding to heights <510 and <800 km, respectively, were set to zero.

It can be deduced from the sentence preceding equation (4) that TEC data assimilations were done at horizontal grid cells, which contained TEC observations. The influence of TEC data assimilation at horizontal grid cells, which lacked TEC data was achieved through successive corrections methods (Bergthorsson & Döös, 1955; Bratseth, 1986; Rodgers, 2000).

4 The Successive Corrections Method

In order to correct NeQuick model generated electron density profile at horizontal grid cell d based on assimilation results (f = 1, 2, 3, …, F) in the nearest neighborhood of d, the successive corrections method described in Bratseth (1986) and Rodgers (2000) was applied after modifying it as,

(9)

where is the corrected electron density profile, is the background electron density profile associated with , and is the background electron density profile at the horizontal grid cell d. In equation (9), the expression in the second term was introduced in order to obtain the average effect of assimilations at F grid cells in the nearest neighborhood of grid d. Overall, simplification of the expression that constitutes the second term in equation (9) yields small quantities. These small quantities either increase or reduce elements (depending on the sign of ) of the NeQuick model generated electron density profile at grid d.

Although the expression in equation (9) contains elements that are much lower than the typical NmF2 value (~10¹² electrons/m²), the elements might still have a factor of 10 raised to a number <12 as a power. Based on this idea, the term 10^−r in equation (9) was introduced to allow the quantity needed to correct background electron density () to reduce as r increases. In fact, the term 10^−r in equation (9) dictates that at large values of r, assimilation results do not (minimally) influence background electron density profile at grid d. On the other hand, when r is small, the assimilation results influence the background electron density profile maximally.

Actually, r quantifies the effect of temporal (d = f) or spatial separation (at a fixed epoch) between grid cells d and f. The consideration of r as spatial or temporal separation was based on the idea that the local time difference between two locations separated by 15° longitude is about 1 h. After performing several trials, this study established that to achieve smooth variations of electron densities spatially and temporarily, r values should vary in steps of 0.25. Moreover, since elements in the expression are expected to be small (≤hundreds) compared to typical NmF2, the first r value was assigned to 1. Equation (10) describes the variation of r as a function of the longitudinal difference between cells d and f (d_lon(f)), the latitudinal difference between cells d and f (d_lat(f)), and the time difference between epoch with observation data and that which lacked observation data (dt(f)).

(10)

where the two vertical bars represent the magnitude, and i is one of the positive integers, which was specified after knowing the spatial or temporal difference between grid cells d and f. For example, when all conditions set in the top row of the right-hand side of equation (10) are satisfied and |d_lat(f)| = 6°, i will take the value of 2.

Then, r = r_f(4 < |d_lat(f)| ≤ 8) = 1 + (2 − 1) × 0.25 = 1.25. Generally, at a particular assimilation time window, all grid cells that lacked observation data were corrected in two ways. In the first case, d_lon(f) was restricted to ≤5°, while r varied with d_lat(f) as in the top expression on the right-hand side of equation (10). In the second case, d_lat(f) was restricted to ≤4°, while r varied with d_lon(f) as in the middle expression on the right-hand side of equation (10). Concerning corrections at a particular grid cell (d = f) for epochs that lacked observation data, the r values varied with dt(f) as in the bottom expression on the right-hand side of equation (10).

5 Results and discussions

5.1 An example of data assimilation process

Figure 3 presents an example of assimilation result () based on equation (4) and TEC observed by ground-based GPS receivers on October 7, 2015. The data assimilated were those at grid cell centered at longitude and latitude 127.5° E and 38° N, respectively as well as time interval 0:00–0:15 UT. In this assimilation window, we considered only 21 TEC observations (n = 21) to be assimilated. For clarity, we only present in Figure 3 and which are associated with the first three and the last two observed TEC values. The red line in Figure 3 represents background electron density () obtained from the NeQuick model. While implementing equation (4), we set . This implies that the 21 observed TEC values were treated as scalars and assimilated in a manner similar to the recursive approach as in Grewal & Andrews (2001).

Figure 3 Electron density profile obtained from NeQuick model (red) and improved profiles (., ×, +, o, *) associated with assimilation of some TEC values within spatial grid cell centered at longitude and latitude 127.5° E and 38° N, respectively. The date and time interval associated with the data plotted are indicated on top of the panel.

Based on the similarity with the recursive approach, we considered associated with the last assimilated observed TEC value as the assimilation result for the specific assimilation window. We need to mention that the peak electron densities computed from and were 6.0 × 10¹¹ and 8.5 × 10¹¹ electrons/m³, respectively. These peak electron densities can be converted to foF2 using the expression (Davies, 1990),

(11)

yielding 6.96 and 7.41 MHz, respectively. These two values of foF2 can be compared with 5.94 MHz obtained from the ionosonde at IC437 (Lon 127°, Lat 36°) during the date and assimilation window indicated in Figure 3. The comparison reveals that the difference between foF2 obtained by assimilation procedure and the observed is smaller than that between foF2 obtained from the NeQuick model and the observed. For generalization purposes, several test scenarios were performed for the entire period under study.

5.2 Assimilation Test Scenarios

The test scenarios presented in this subsection were validated at horizontal grid cells centered at geographic latitude (longitude) (i) 26° (122.5°), (ii) 34° (137.5°), (iii) 38° (127.5°), (iv) 33.4°(126.3°), (v) 26.3°(127.8°), and (vi) 39.6°(115.9°). The validation at grid cells (i) and (ii) were done using TEC data obtained from GPS receiver at TCMS and foF2 data obtained from ionosonde at TO536, respectively, while that at (iii) was done using TEC data obtained from GPS receiver at YONS as well as foF2 data obtained from ionosonde at IC437. Moreover, the validation at grids (iv) and (v) were done using foF2 data obtained from OK426 and JJ433, while the validation at grid (vi) was done using TEC data obtained from BJFS. It should be noted that validation station data were not used during assimilation.

5.2.1 Assimilation of Ground-based GPS receiver TEC

We assimilated ground-based GPS receiver-derived TEC to the NeQuick model using equation (4). The foF2 derived from the assimilation results over IC437 during March 15–20, October 6–11, 27, and December 18–24, 2015, are presented in Figures 4a, 4b, and 4c, respectively. The corresponding foF2 observed by the ionosonde at the station, as well as foF2 obtained from the NeQuick model, are superimposed. The blue, red, and black colors in Figure 4 and later in Figures 5 and 7 represent parameters obtained from the NeQuick model, assimilation, and observation, respectively.

Figure 4 Panels (a), (b), and (c) present variation of foF2 during March 15–20, October 6–11, 27, and December 18–24, 2015, respectively over IC437. The blue, red, and black colors represent NeQuick model generated parameter, assimilation result, and observed parameter.

Figure 5 Panels (a), (b), (c), and (d) present variation of TEC during March 15–20, June 4, 21–26, October 6–11, 27, and December 18–24, 2015, respectively over YONS. The blue, red, and black colors represent NeQuick model generated parameter, assimilation result, and observed parameter.

Figures 5a–5d present TEC over YONS obtained from assimilation, NeQuick model, and GPS receiver station during March 15–20, June 4, 21–26, October 6–11, 27, and December 18–24, 2015, respectively. The increase and reduction in comparison to a background value of observed ionospheric parameters during geomagnetic storms are usually termed as positive and negative ionospheric storm effects, respectively (Buonsanto, 1999). As mentioned in Section 2, this study associated the ionospheric parameters during March 15, June 4, October 27, and December 18, 2015, when Kp ≤ 3 with background values. It can be deduced from Figures 4c, 5b, and 5d that values of ionospheric parameters increased during the main phases of geomagnetic disturbances in June and December solstices when compared to those on quiet days in these seasons. Moreover, Figures 4a, 4b, 5a, and 5c clearly show that values of ionospheric parameters reduced significantly during the main phases of geomagnetic disturbances in March and September equinoxes when compared to those on quiet days in these seasons. Since detailed discussions about the physical mechanisms responsible for the generation of positive and negative ionospheric storm effects are out of the scope of the current study, interested readers may refer to the previous studies (e.g., Buonsanto, 1999; Mendillo & Klobuchar, 2006) for such discussions.

Generally, Figures 4 and 5 show that the assimilation yielded ionospheric parameters, which are in most cases closer to the observed parameter than to those generated from the NeQuick model. Therefore, the assimilation process yields ionospheric parameters, which depict the response of the ionosphere to geomagnetic disturbance. Another general feature of Figures 4 and 5 is that assimilation results and the observed data exhibited daily/diurnal variability, which is almost not present in the climatological NeQuick model values. However, visual inspection of Figure 4 shows some cases (e.g., panel (c) on December 18 and 24) where NeQuick model results are closer to the observed data compared to assimilation results. These isolated cases might result because the data assimilated are obtained using a GPS receiver while the observation data presented in Figure 4 are obtained using ionosonde. The inherent discrepancy in the instruments might manifest in the observed cases of assimilation results not performing well. The overall statistical analysis presented in Figure 6b confirms that these are isolated cases.

Figure 6 Panels (a) and (b) present scatter plots of foF2 obtained from NeQuick model (red color) and assimilation result (black color) as a function of coincident observed foF2 over TO536 and IC437, respectively. Panel (c) presents TEC obtained from NeQuick model and assimilation result as a function coincident observed TEC over YONS. The observations are those that fall within the study period.

Figures 6a–6c present the scatter plots of NeQuick model generated ionospheric parameter (red color) and assimilation results (black color) as a function of a coincident observed parameter over TO536, IC437, and YONS, respectively. It should be noted that the data plotted in Figure 6 is the same as that in Figures 4 and 5. The correlation coefficients r and root mean squared error (RMSE) associated with the data plotted are indicated on the respective panels. It can be seen in Figure 6 that the r-value associated with the assimilation result parameter over a particular station is always higher than that of the NeQuick model. Moreover, the RMSE value associated with the assimilation result parameter over a particular station appears significantly improved (reduced) compared to that associated with the NeQuick model. In fact, the improvement percentage (IP) of the RMSE values were determined as

(12)

where NQ_E and ASM_E represent RMSE values associated with the NeQuick model and assimilation result, respectively. The IP over TO536, IC437, and YONS were 69, 56, and 81%, respectively.

The RMSE improvement percentage values (≥56%) we obtained in this study are higher than the 44% of Mengist et al. (2019) when they investigated the performance of the Ionospheric Data Assimilation Four‐Dimension technique over Korea and adjacent areas, considering International Reference Ionosphere model as the background model. The study by Mengist et al. (2019) was done during both geomagnetic quiet and disturbed days (March 15–18, 2015). The limitation of the vertical grids to an altitude of ~1336 km in Mengist et al. (2019) might contribute to the low RMSE improvement percentage obtained in their study compared to that of the current study. The study on imaging South African regional ionosphere using a 4D-var technique by Ssessanga et al. (2019) might not reasonably estimate TEC up to an altitude of GPS satellites since it limited the vertical grid to an altitude of 1336 km.

The results presented in Figures 4 and 6 signify that assimilation of TEC data obtained from ground-based GPS receiver to NeQuick model can be helpful in determining fairly well foF2 over a location that does not have an ionosonde station.

5.2.2 Assimilation of TEC obtained from ionosonde stations

After assimilating ionospheric TEC data obtained from ionosonde stations to the NeQuick model, the validation of the assimilation exercise was done using TEC data obtained from YONS. It should be noted that most records of ionosonde stations did not have available TEC data with the exception of IC437, which belongs to the same spatial grid cell (ii) as YONS. Figure 7 presents the variations of TEC obtained from assimilation (red) over YONS during the period under study. Superimposed over the figure is the corresponding TEC obtained from the NeQuick model as well as observed TEC data obtained from GPS receiver over YONS. Figure 7 clearly shows similar observations that were deduced from Figures 4 and 5. Overall, the TEC obtained by assimilation closely follows the observed TEC.

Figure 7 Panels (a)–(c) present TEC obtained from NeQuick model (blue), assimilation (red), and observed by ground-based GPS receiver (black) during March 15–20, October 6–11, 27, and December 18–24, 2015, respectively over YONS.

In order to quantify how well the assimilation has improved estimation of TEC over YONS, we present in Figure 8 a scatter plot of TEC obtained from assimilation (black) and NeQuick model (red) as a function of coincident observed TEC over YONS. In Figure 8, the r and RMSE associated with TEC obtained from assimilation and the NeQuick model are indicated. The r values portray that the TEC obtained from assimilation correlates with observed TEC better than that obtained from the background model. By performing assimilation, the RMSE improved by ~50%. This low improvement compared to that associated with assimilation of ground-based GPS receiver data could be due to the lack of inclusion of TEC above 700 km during assimilation. In addition, although electron density can be computed by the NeQuick model up to an altitude of 700 km, the ionosonde observations may not in some cases reach this altitude. This might partly explain why the assimilation TEC data in Figure 7 appears to be noisy.

Figure 8 Scatter plot of TEC obtained from NeQuick model (red color) and assimilation result (black color) as a function of coincident observed TEC over GPS receiver station at YONS. The data plotted are the same as that in Figure 7.

It needs to be noted that the 50% increment signifies those over locations that do not have ground-based GPS receivers, assimilation of TEC obtained from ionosondes into the NeQuick model can be helpful in improving the estimation of TEC data that would be measured by ground-based GPS receivers.

5.2.3 Assimilation of TEC obtained from COSMIC RO

Due to the scarcity of RO TEC data, we validated the results of assimilation of the data into the NeQuick model with foF2 obtained from all the ionosonde stations considered in this study. Moreover, validations were done using TEC observed over Japan (KGNI), South Korea (YONS), Taiwan (TCMS), and China (BJFS). Figure 9a presents a scatter plot of foF2 obtained from assimilation (black) and NeQuick model (red) as a function of coincident observed foF2 over the ionosonde stations, while Figure 9b presents TEC obtained from assimilation and NeQuick model as a function of coincident observed TEC over the IGS stations. The r and RMSE associated with data plotted in panels of Figure 9 are shown on the respective panels. Also indicated on the panels of Figure 9 are the numbers of observed data that are coincident with assimilation results. Despite considering 27 days in the year 2015 and all 4 ionosonde stations, there were only 10 observations of foF2 values, which were coincident with assimilation results. For the case of validation with TEC data, there were only 17 observations that were coincident with assimilation results.

Figure 9 Panel (a) presents scatter plot of foF2 obtained from NeQuick model (red) and assimilation (black) as a function of coincident observed foF2 over the ionosonde stations considered in this study. Panel (b) presents scatter plot of TEC obtained from NeQuick model and assimilation as a function of coincident observed TEC over KGNI, TCMS, YONS, and BJFS GPS receiver stations. The data plotted are those on the days of the study period.

Figure 9 exhibits that the assimilation of COSMIC RO TEC greatly improves the estimation of foF2 and TEC since the r values associated with assimilation results are higher than those associated with the NeQuick model. In fact, the RMSE improvement percentages for foF2 and TEC are 57 and 10%, respectively. The percentage associated with TEC data is much less compared to that of foF2, maybe due to the limitation of COSMIC RO electron density profile altitude of ~800 km.

Based on the observations and discussions associated with Figure 9, the assimilation of COSMIC RO TEC seems to improve the estimation of foF2 and TEC ionospheric parameters.

5.2.4 Assimilation of TEC obtained from GPS receivers onboard Swarm and COSMIC satellites

The results of assimilation of TEC obtained from GPS receivers onboard Swarm, and COSMIC satellites were validated with ground-based GPS receivers over Japan (KGNI), South Korea (YONS), Taiwan (TCMS), and China (BJFS). Similar to Figure 9, Figure 10 presents the scatter plot of NeQuick model TEC (red color) and assimilation TEC (black color) as a function of coincident observed TEC over the ground-based GPS receiver stations. Indicated in Figure 10 are the r and RMSE associated with the data plotted in the figure, as well as the number of coincident assimilation TEC and observed TEC over the GPS receiver stations. Even though 27 days in the year 2015 and TEC data from GPS receiver stations located in 4 different countries were considered, there were still few (35) coincident data, as indicated in Figure 10. This observation may be due to the fact that LEO satellites pass over a particular location reoccurs after several days.

Figure 10 Scatter plot of TEC obtained from NeQuick model and assimilation as a function of coincident observed TEC over GPS receiver stations at KGNI, TCMS, YONS, and BJFS. The data plotted are those on the days of the study period.

Figure 10 reveals that the assimilation of TEC obtained from GPS receivers onboard LEO satellites yields a lower correlation coefficient compared to that associated with the NeQuick model. Actually, the RMSE deteriorated by ~107%.

The poor estimation of ionospheric parameters by the assimilation of TEC obtained from GPS receivers onboard LEO satellites need to be investigated further. However, we tentatively attribute this poor performance to the (i) dynamics of the receiver, which might not have been considered in computing the POD TEC (ii) limitation of the assimilation technique.

In a practical application situation, assimilation of TEC data would depend on the precedence of the data source, which can be set based on the results depicted in the various test scenarios presented above. For instance, TEC data obtained from ground-based GPS receiver would be given the highest precedence, followed by TEC data from ionosonde, and lastly, COSMIC RO TEC.

Based on the results from the test scenarios, for now, we would not recommend assimilation of TEC derived from GPS receivers onboard the LEO satellites to the NeQuick model. If a horizontal grid cell does not have TEC data from any of the three recommended sources, the NeQuick model generated electron density profile at the cell would be corrected as described in Section 4. Examples of assimilation results obtained from data assimilations followed by application of successive corrections method are presented and compared with TEC data processed at the Center for Orbit Determination in Europe (CODE) in Section 5.3. As one of the international GPS services (IGS) for geodynamics analysis centers, CODE provides daily Global Ionospheric TEC data Maps (GIMs) at http://www.aiub.unibe.ch/download/CODE/.

5.3 Validation of successive corrections method

This section presents in Figure 11 results where TEC data from all the ground-based GPS receiver and ionosonde stations indicated in Figure 2 as well as COSMIC RO were assimilated followed by successive corrections method to yield the final assimilation result. Later in this section, we also present in Figure 12 similar results, but where TEC data from ionosonde stations at JJ433 and OK426 were not assimilated. These two stations appear to be suitable for validating further successive corrections method since Figure 2 shows there were no ground-based GPS receivers within the grid cells that contain the stations. Figures 11a–11c (panels in column) present CODE GIMs, TEC obtained from assimilation result, and NeQuick model, respectively. Figures 11 (i)–11 (v) (panels in row) present the TEC data at 5:00 UT (14:00 and 15:00 LT over Korea and Japan, respectively) when considerable ionization is expected during October 6–10, 2015, respectively. These dates were chosen to reveal the ionospheric changes before and during the main phase (October 6–7, 2015) of the storm as well as during and after the recovery phase (October 8–10).

Figure 11 Panels in columns (a)–(c) present TEC obtained from CODE GIMs, assimilation, and NeQuick model, respectively. Panels in rows (i)–(v) present the TEC data at 05:00 UT during October 6–10, 2015, respectively.

Figure 12 Panels 1 (a) and 2 (a) show foF2 obtained from NeQuick model (blue line), assimilation results (red line), and ionosonde stations (black dots) at JJ433 and OK426, respectively. Panels 1 (b) and 2 (b) present the magnitudes of the differences between foF2 obtained from NeQuick model (blue bars) and assismilation results (red bars) and that observed at ionosonde stations at JJ433 and OK426, respectively.

It can be seen in Figure 11 that on October 7, 2015, TEC data from CODE and assimilation were the highest and reduced significantly on October 8, 2015. This observation is consistent with the results that were presented and discussed in Section 5.2 (see Figs. 4–6), where the TEC values during the main phase of the storm on October 7, 2015 (see Fig. 1) were found to be higher than that during the recovery phase of the storm on October 8, 2015. As expected, these variations in TEC data during the main and recovery phases of the storm are not strongly reflected in the NeQuick model presented in Figure 11c. Visual inspection of panels in a row (iii) of Figure 11 shows that at longitude >132°, TEC data from CODE GIMs are higher than that of assimilation and NeQuick. Since CODE GIMs are constructed using a series of spherical harmonics functions whose coefficients are determined using available TEC data from IGS stations (Schaer, 1999), CODE GIMs TEC data over locations that lacked IGS stations may contain the high error value. Therefore, although TEC data from both CODE and assimilation seem to respond to TEC changes due to the occurrence of a geomagnetic storm, TEC data from the two sources may not perfectly correlate. Furthermore, although the resolution of the CODE GIMs was changed from 2 h to 1 h on October 19, 2014, this hourly averaging might still prevent the capturing of fine structures in the maps.

As mentioned before in this subsection, the successive corrections method was validated further using foF2 obtained from ionosonde stations at JJ433 and OK426. Panels 1 (a) and 2 (a) in Figure 12 present the foF2 obtained from the NeQuick model (blue line), assimilation results (red line), and ionosonde stations (black dots). The magnitudes of the differences between foF2 obtained from (i) NeQuick model and (ii) assimilation results and the observed, denoted as ΔfoF2, are presented in panels 1 (b) and 2 (b) of the Figure 12. It is important to mention that the data plotted in Figure 12 are for the months and days in the month indicated on the horizontal axis of panel 2b. For each day, the data were sampled at 4 h intervals. Particularly, the data corresponding to 02:00, 06:00, 10:00, 14:00, 18:00, and 22:00 LT are plotted in Figure 12.

It can be seen in panels 1 (b) and 2 (b) of Figure 12 that the ΔfoF2 associated with assimilation results are mostly smaller than those corresponding to the NeQuick model. The average ΔfoF2 associated with assimilation results over JJ433 and OK426 was established as 0.79 and 1.30 MHz, respectively, while the average ΔfoF2 associated with NeQuick results over JJ433 and OK426 were 1.04 and 1.40 MHz, respectively. These results imply that over a particular station, the average ΔfoF2 associated with assimilation reduces significantly compared to that associated with the NeQuick model. The high average ΔfoF2 and foF2 values over OK426, as depicted in Figure 12 might be associated with the closeness of the station to the equatorial region where high ionization and electro-dynamic processes occur. The sparse availability of ground-based GPS receivers within the vicinity of OK426, as shown in Figure 2 might be another reason for the high average ΔfoF2 observed over the station. This is expected since the effectiveness of the successive corrections method decreases as the distance from locations of data increases.

6 Conclusions

We made an effort to improve the estimation of ionospheric parameters over Korea and adjacent areas by employing the classical Kalman filtering technique to assimilate TEC data from various sources into the NeQuick model. Successive corrections method was applied to spread the effect of TEC data assimilation at a given location to others that lacked TEC observations. The results from different assimilation scenarios showed that data assimilation of ground-based GPS-derived TEC data could improve RMSE associated with the model estimation by ≥56%. Assimilation of TEC measured by ionosonde stations can improve RMSE associated with the model estimation of TEC data by ~50%. The assimilation of TEC obtained from COSMIC RO revealed RMSE improvement of ~10%. Assimilation of TEC measured by GPS receivers’ onboard LEO satellites degraded the RMSE associated with the model estimation by ~107%, probably due to either the dynamics of the receivers or limitation of the assimilation technique. Validation of our assimilation results with global ionosphere TEC data maps processed at CODE revealed that both reproduced similar TEC changes, showing response to a geomagnetic storm. However, TEC data from the two sources may not perfectly correlate.

For practical applications, we propose the assimilation of TEC data into the NeQuick model depending on the precedence of the data source, which can be set based on the results presented in this study. That is, TEC data obtained from ground-based GPS receiver would be given the highest precedence, followed by TEC data from ionosonde, and lastly, COSMIC RO TEC.

Acknowledgments

Part of the work presented in this article was done when the first Author, Patrick Mungufeni, was a post-doc researcher in Chungnam National University, South Korea, under the supervision of Prof. Young-Ha Kim. The first author, Patrick Mungufeni, greatly appreciates Dr. Sripathi Samireddipalle of the Indian Institute of Geomagnetism for the discussions about the manuscript. We thank the developers of the NeQuick model for making their model available. Dst data is provided by the World Data Center for Geomagnetism at Kyoto (http://swdcwww.kugi.kyoto-u.ac.jp/), and the Kp data is provided by GFZ Potsdam (ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/). The ionosonde observations were obtained from the website of the National Oceanic and Atmospheric Administration (http://www.swpc.noaa.gov/), while ionPrf files used to derive COSMIC RO TEC were obtained from COSMIC Data Analysis and Archive Center (http://cosmic-io.cosmic.ucar.edu/cdaac/index.html). The RINEX files were used in this study were obtained from the University NAVSTAR Consortium website (ftp://data-out.unavco.org). The global ionosphere TEC data maps processed at CODE were obtained from http://www.aiub.unibe.ch/download/CODE/. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

	Figure 1 Variation of (a) Dst and (b) Kp during (1) March 15–20, 2015, (2) June 4, 21–26, 2015, (3) October 6–11, 27, 2015, and (4) December 18–24, 2015.
In the text

	Figure 2 Locations of GPS (blue crosses) and ionosonde (red diamonds) stations used in the study.
In the text

	Figure 3 Electron density profile obtained from NeQuick model (red) and improved profiles (., ×, +, o, *) associated with assimilation of some TEC values within spatial grid cell centered at longitude and latitude 127.5° E and 38° N, respectively. The date and time interval associated with the data plotted are indicated on top of the panel.
In the text

	Figure 4 Panels (a), (b), and (c) present variation of foF2 during March 15–20, October 6–11, 27, and December 18–24, 2015, respectively over IC437. The blue, red, and black colors represent NeQuick model generated parameter, assimilation result, and observed parameter.
In the text

	Figure 5 Panels (a), (b), (c), and (d) present variation of TEC during March 15–20, June 4, 21–26, October 6–11, 27, and December 18–24, 2015, respectively over YONS. The blue, red, and black colors represent NeQuick model generated parameter, assimilation result, and observed parameter.
In the text

	Figure 6 Panels (a) and (b) present scatter plots of foF2 obtained from NeQuick model (red color) and assimilation result (black color) as a function of coincident observed foF2 over TO536 and IC437, respectively. Panel (c) presents TEC obtained from NeQuick model and assimilation result as a function coincident observed TEC over YONS. The observations are those that fall within the study period.
In the text

	Figure 7 Panels (a)–(c) present TEC obtained from NeQuick model (blue), assimilation (red), and observed by ground-based GPS receiver (black) during March 15–20, October 6–11, 27, and December 18–24, 2015, respectively over YONS.
In the text

	Figure 8 Scatter plot of TEC obtained from NeQuick model (red color) and assimilation result (black color) as a function of coincident observed TEC over GPS receiver station at YONS. The data plotted are the same as that in Figure 7.
In the text

Figure 9 Panel (a) presents scatter plot of foF2 obtained from NeQuick model (red) and assimilation (black) as a function of coincident observed foF2 over the ionosonde stations considered in this study. Panel (b) presents scatter plot of TEC obtained from NeQuick model and assimilation as a function of coincident observed TEC over KGNI, TCMS, YONS, and BJFS GPS receiver stations. The data plotted are those on the days of the study period.

In the text

	Figure 10 Scatter plot of TEC obtained from NeQuick model and assimilation as a function of coincident observed TEC over GPS receiver stations at KGNI, TCMS, YONS, and BJFS. The data plotted are those on the days of the study period.
In the text

	Figure 11 Panels in columns (a)–(c) present TEC obtained from CODE GIMs, assimilation, and NeQuick model, respectively. Panels in rows (i)–(v) present the TEC data at 05:00 UT during October 6–10, 2015, respectively.
In the text

Figure 12 Panels 1 (a) and 2 (a) show foF2 obtained from NeQuick model (blue line), assimilation results (red line), and ionosonde stations (black dots) at JJ433 and OK426, respectively. Panels 1 (b) and 2 (b) present the magnitudes of the differences between foF2 obtained from NeQuick model (blue bars) and assismilation results (red bars) and that observed at ionosonde stations at JJ433 and OK426, respectively.

In the text