Assessment of GRAS Ionospheric Measurements for Ionospheric Model Assimilation

Abstract

We conducted a study to assess the GNSS (Global Navigation Satellite System) Receiver for Atmospheric Sounding (GRAS) ionospheric data quality by processing Radio Occultation (RO) observations of ionospheric products. The main objective of the study is to validate ionospheric data generated at EUMETSAT, such as ionospheric bending angle profiles, amplitude and phase scintillations, topside Total Electron Content (TEC) from MetOp-A GRAS instrument as well as generating and validating new ionospheric products derived from GRAS RO observations such as the TEC, rate of TEC and vertical electron density profiles. The assessment is conducted by comparing and evaluating the systematic differences between similar products from other Low Earth Orbit (LEO) satellite missions or from ground-based ionospheric measurements. The study confirms that the GNSS topside and RO observations recorded by the GRAS instrument onboard MetOp satellites are of good quality and are a valuable source of data for ionospheric research.

1. Introduction

The GNSS (Global Navigation Satellite System) Receiver for Atmospheric Sounding (GRAS) instrument on-board the MetOp satellites of the EUMETSAT (European Organisation for the Exploitation of Meteorological Satellites) Polar System (EPS) programme [1] uses Radio Occultation (RO) technique to obtain information about the temperature and humidity in the atmosphere. Therefore, GRAS instruments were designed in such a way that their measurements covered a vertical range from the lower troposphere up to a height of approximately 80 km allowing for neutral atmospheric RO soundings only. Theoretically, if the measurements can cover a height range of up to 600 km (which includes a significant part of the topside ionosphere above the peak ionization), the RO sounding should be capable of providing information on the vertical structure of the ionospheric electron density. Given the increased interest in EUMETSAT’s user community for the provision of space weather and ionospheric data, the Radio Occultation instrument on board the future EPS Second Generation (EPS-SG) EUMETSAT satellites flying in a similar polar orbit at about 800 km height will track GNSS signals up to 600 km. In preparation for this, EUMETSAT conducted a MetOp-A end-of-life testing campaign during the summer of 2020, which offered an opportunity to test an updated configuration of the GRAS instrument, extending its vertical measurement range up to 300 km and 600 km into the lower and middle of the ionosphere. The test campaign lasted for three months, and a large set of GPS radio occultation measurements were recorded by the GRAS instrument during the two different experiments up to 600 km and 300 km heights. It is noted that MetOp-A was flying at an orbit height of about 800 km during the test campaign. To assess the GRAS ionospheric data quality, we have accomplished a study called GIMA (Assessment of GRAS Ionospheric Measurements for Ionospheric Model Assimilation) recently, which includes the following analysis:

• Quality assessment of the dual-frequency bending angle, amplitude, and phase scintillation profiles generated at EUMETSAT by comparing against independent measurements;
• Quality assessment of the topside Total Electron Content (TEC) data generated at EUMETSAT by comparing against independent measurements;
• Processing of RO observations truncated at 600 km and 300 km and generation of ionospheric products, such as slant TEC profiles, vertical electron density values, and extension of electron density profile up to the MetOp-A orbit height;
• Validation of the GRAS RO electron density profiles by comparing them against independent measurements;
• Evaluation of the potential impact of assimilating GRAS data into ionospheric data assimilation models and assessment of the suitability of topside TEC data for data assimilation.

It is noted that MetOp-A successfully provided valuable meteorological data for fifteen years to global users until its deorbiting process started in November 2021. However, MetOp-B and -C are in operation and flying at a high inclination near Polar orbit at an altitude of about 817 km. The MetOp-A end-of-life extension campaign provides an early opportunity of testing truncated data for RO retrievals and confirms that the GRAS instrument is capable of performing valuable ionospheric sounding. Being motivated by the successful test campaign, currently, the GRAS receivers on board MetOp-B and -C are also configured for the remote sensing of the atmosphere up to 300 km altitude for facilitating tropospheric as well as ionospheric sounding. The present paper summarizes the main findings of the GIMA study. The used database and data sources for GRAS data quality assessment are described in Section 2. Section 3 briefly describes the RO inversion method developed for processing GRAS observations to ionospheric products. Section 4 summarizes the GRAS data quality assessment results, including the assessment of bending angle, scintillations profiles, and topside TEC data quality. Section 4 also includes the assessment of the suitability of assimilating topside TEC data into ionospheric data assimilation models. Section 5 concludes the GIMA study results.

2. Database and Data Sources

As already mentioned, the original MetOp-A GRAS measurements only covered a vertical range up to a height of approximately 80 km, suitable for neutral atmospheric RO soundings only. Thanks to the MetOp-A end-of-life extension campaign, which enables the GRAS instrument to extend its vertical measurement range up to 300 km and 600 km height, in two separate experiments, for about a three months period during 2020. The experiment datasets released by EUMETSAT refer to the day of the year (DOY) 176 to 253 of 2020 (i.e., from 24 June until 9 September 2020). The datasets can be divided into two major groups. The first group consists of GRAS RO observations along GPS-MetOp-A links obtained by the limb-sounding antenna. The second group contains the GPS topside TEC data recorded by the zenith antenna.

Table 1. MetOp-A end-of-life experiment period.

600 km Extension Campaign	300 km Extension Campaign
DOY 197–239, 2020 (15 July–26 August 2020)	DOY 175–196 and DOY 240–253, 2020 (23 June–14 July and 27 August–9 September 2020)

gives the specific dates for which the 600 km and 300 km campaigns were conducted.

In the case of the 600 km extension campaign, the GRAS instrument recorded GPS RO signals for the link geometries having a tangential height (equal to or) below 600 km. Similarly, the GRAS instrument recorded RO measurements below 300 km tangential height for the 300 km extension campaign. It is noted that MetOp-A was orbiting at about 800 km altitude during both experiment periods. Since RO observations are not recorded up to MetOp-A orbit height, the datasets are named as truncated data here afterward.

For an independent assessment of the GRAS data products derived during the experiment, similar data products from other satellite missions and sources are considered.

Table 2. Satellite missions are listed in left-most column, and available data products and their sources are listed in the right columns.

Satellite Mission	Inclination [deg]	Height [km]	Product Name	Source
MetOp-A, B, C from UCAR	~98.7	~800	podTec	https://0-doi-org.brum.beds.ac.uk/10.5065/789w-m137 (accessed on 15 January 2023)
COSMIC-1	~72	~800	scnLv1 (s4 index)	https://data.cosmic.ucar.edu/gnss-ro/cosmic1/ (accessed on 15 January 2023)
COSMIC-2	~24	~550	scnLv1, podTc2, ionPrf	https://data.cosmic.ucar.edu/gnss-ro/cosmic2/provisional/spaceWeather/ (accessed on 15 January 2023)
Fengyun-3D	~98.8	~836	Gnosx (bending angle, Ne, TEC)	http://satellite.nsmc.org.cn/portalsite/ (accessed on 15 January 2023)
Swarm (A, B, C)	~87	~460 (A, C) ~510 (B)	Level2/TECATMS	ftp://swarm-diss.eo.esa.int/ (accessed on 15 January 2023)
DMSP	~99	~800	Ne in situ	http://cedar.openmadrigal.org/ (accessed on 15 January 2023)
GIRO (ionoson)		--	NmF2, hmF2	http://giro.uml.edu/didbase/scaled.php (accessed on 15 January 2023)
MetOp-A iono 1d-var	~98.7	~800	Ne profile	via EUMETSAT
SARAL	~98.5	~800	DORIS measurements in RINEXv3	https://cddis.gsfc.nasa.gov/archive/doris/data (accessed on 15 January 2023)
Sentinel-3A	~98.7	814.5	DORIS measurements in RINEXv3	https://cddis.gsfc.nasa.gov/archive/doris/data (accessed on 15 January 2023)
Sentinel-3B	~98.7	814.5	DORIS measurements in RINEXv3	https://cddis.gsfc.nasa.gov/archive/doris/data (accessed on 15 January 2023)

lists the satellite missions and observations used in the GRAS data assessment and also gives the corresponding data sources.

The COSMIC-1 (S4 index) and COSMIC-2 data (topside TEC, S4, RO electron density profile) processed by the University Corporation for Atmospheric Research (UCAR) were obtained through the COSMIC Data Analysis and Archive Center (CDAAC, see Table 2). Additionally, the topside GPS TEC measurements onboard MetOp-A, -B, and -C satellites processed by UCAR were downloaded from the CDAAC site. The low orbit inclination of the COSMIC-2 mission restricts the availability of its products to low and middle-latitude regions. Considering this, the Fengyun-3D (high inclination orbit) data with global coverage was obtained through the China Meteorological Administration National Satellite Meteorological Center. The Fengyun-3D electron density and bending angle profiles are used for the GRAS data assessment. The Global Ionosphere Radio Observatory (GIRO) network provides worldwide vertical-sounding data such as ionosonde measurements. The GIRO data (e.g., F2 layer peak density NmF2 and height hmF2) are obtained from the University of Massachusetts Lowell (UML) portal. Swarm mission has been operating three satellites, namely Swarm A, Swarm B, and Swarm C, since November 2013 at orbit heights of about 460–510 km. The topside GPS data onboard Swarm satellites are obtained from the ESA site (see Table 2) and used for the GRAS data quality assessment. The used in situ electron density data from the Defense Meteorological Satellite Program (DMSP) are originally processed by the Center for Space Sciences at the University of Texas (Dallas) and are obtained through the Madrigal portal (see link in

Table 3. Mean and STD of the differences between the bending angles from MetOp-A and Fengyun-3D using the co-location of 2 h.

Latitude Range (deg)	Mean Difference (murad)	STD of Differences (murad)
−90 to −45	0.73	10.31
−45 to 0	0.64	12.73
0 to 45	−0.10	22.14
45 to 90	3.58	14.41

). MetOp-A electron density profiles retrieved by a new approach recently developed at the ROM SAF (Radio Occultation Meteorology Satellite Application Facility), based on the dual-frequency bending angle (BA) 1d-var data assimilation (see [2,3]), are used for validating the GRAS RO electron density profiles. The BA-derived electron density profiles are provided by EUMETSAT. Dual-frequency DORIS measurements from SARAL, Sentinel-3A, and Sentinel-3B satellites were obtained from Crustal Dynamics Data Information System (CDDIS) site, as mentioned in Table 2. DORIS data was an independent source of ionospheric TEC used in ionosphere model assimilation (see Section 4).

As part of GRAS ionospheric data quality assessment, EUMETSAT made available MetOp-A onboard topside and RO GNSS data (i.e., dual-frequency carrier-phase and code-pseudo ranges) as well as derived products such as the amplitude and phase scintillation indices and bending angle data. Within the scope of the GIMA study, RO observations are processed to generate ionospheric products such as the slant TEC and vertical electron density profiles. The following section briefly describes how the GRAS RO observations were processed for ionospheric products.

3. Process GRAS RO Observations to Ionospheric Products

The GRAS RO observations (e.g., dual-frequency carrier phases) are processed, and the following ionospheric products are generated.

3.1. Profiles of slant Total Electron Content (sTEC)

The slant Total Electron Content sTEC is derived from the dual-frequency carrier-phase measurements by

$$sTEC=\frac{f_1^2{f}_2^2}{40.3\left(f_1^2-f_2^2\right)}\left({\Phi }_1-{\Phi }_2\right)+B_{ambiguity}-{bias}^{satellite,receiver}$$

(1)

where Φ₁ and Φ₂ denote the carrier-phase observations at f₁ and f₂ frequencies, respectively, B_ambiguity is the carrier-phase ambiguity term. The computed sTEC is biased by a constant ambiguity and an inter-frequency phase bias term. Other terms, such as multipath and wind-up terms, are neglected for simplicity. The sTEC data obtained by raw observations (carrier phases) provided at 50 Hz are down-sampled to 1 Hz taking the first measurement value during each 1 s period, and are used as inputs to the RO inversion procedure for retrieving electron density profiles. The RO inversion procedure is discussed in detail in Section 3.3.

3.2. Profiles of sTEC Rate of Change

The sTEC rate of change (ROT) denoted by dsTEC/dt is computed from the dual-frequency carrier-phase data as follows.

$$\frac{dsTEC}{dt}=\frac{s{TEC}^{i+1}-s{TEC}^i}{t^{i+1}-t^i}$$

(2)

$${sTEC}^{i+1}-{sTEC}^i=\frac{f_1^2{f}_2^2}{40.3\left(f_1^2-f_2^2\right)}\left[\left({\Phi }_1^{i+1}-{\Phi }_2^{i+1}\right)-\left({\Phi }_1^i-{\Phi }_2^i\right)\right]$$

(3)

where the superscripts i and i + 1 denote the measurements for consecutive time t epochs. Since the receiver onboard the LEO satellite tracks the same GPS satellite during one occultation event, the ambiguity term and the satellite and receiver phase biases remain the same for consecutive measurements and are cancelled out in Equation (2).

Panels a and b of Figure 1 show an example of the sTEC profile and the corresponding ROT obtained from the same RO event, respectively. Comparing the original 50 Hz (blue curve) and down-sampled 1 Hz (red curve) data in panel a, we see that the estimation of ROT varies depending on the measurement sampling rate. As expected, the ROT estimated from the high rate data (e.g., 50 Hz) experiences more fluctuations compared to the ROT estimated from the low rate (1 Hz) data.

The ROT and Signal to Noise Ratio (SNR) profiles show that, at ionospheric E regions (~90–140 km tangential height), there are a few sudden high fluctuations in both 50 Hz and 1 Hz data. In general, the magnitude of these fluctuations in the 1 Hz data is less than that in the 50 Hz data. One reason may be that the ROT computed from 50 Hz data shows ionospheric gradients on a much smaller scale (50 times smaller) than that from 1 Hz data. Comparing ROT and SNR plots, we see a high correlation between fluctuations in the ROT and SNR data.

Figure 1 shows similar variations for SNR and ROT in the E-layer, whereas they differ above 500 km tangential height. A typical ionosphere profile shows an exponential decay of the electron density at the upper height above the peak density height. Thus, above 500 km height, the electron density value is small. Due to this, in the upper height above 500 km, the change in sTEC due to a small change in RO measurement geometry (especially for 50 Hz sampling rate) is small. We found that the sTEC difference is in the level of measurement accuracy (subject to errors due to phase noise etc.). This causes high fluctuations (positive and negative values) in ROT for 50 Hz data above 500 km height, although corresponding SNR data does not show significant fluctuations.

The decrease in the SNR profile at higher altitudes (>300 km) is due to the Radio Occultation antenna’s pattern. The GRAS antenna was designed to have its best performances (highest and constant gain) for tracking neutral atmospheric occultations. Above 300 km altitude, the gain starts to decrease.

3.3. Profiles of Electron Density up to the MetOp-A Orbit

GNSS RO measurements onboard a LEO satellite is commonly used to compute vertical electron density profiles from LEO orbit down to the Earth’s surface [4,5,6,7,8]. However, MetOp-A experimental setup only provides RO data up to 600 km or 300 km height instead of providing data up to satellite height. On the one hand, the missing topside ionosphere and plasmasphere contribution above 600 or 300 km up to GPS orbit height needs to be correctly modelled. On the other hand, the retrieved electron density needs to be extrapolated up to MetOp-A orbit height. Therefore, the electron density reconstruction from MetOp-A data becomes challenging.

Hernández-Pajares et al. [4] were the first to investigate and successfully reconstruct the electron density from truncated measurements using COSMIC 1 mission data. In this respect, Hernández-Pajares et al. [4] developed the foundations of AVHIRO (Abel-VaryChap Hybrid density profile from topside Incomplete RO data) technique based on Vary-Chapman extrapolation in the topside, which was fully developed at Lyu et al. [5], proposing as well a simplified technique called SEEIRO (Simple Estimation of Electron density profile from topside Incomplete RO data) for electron density retrieval from the truncated data. A comparison with the AVHIRO technique shows that although the SEEIRO technique runs faster than the AVHIRO technique in terms of computational time per occultation, the AVHIRO technique performs better in terms of accuracy [5].

Within the scope of the GIMA study, a new model-assisted RO inversion technique is developed for electron density retrieval from MetOp-A truncated data [6]. The topside ionosphere and plasmasphere above the LEO orbit height are modelled by a Chapman layer function (see Equation (4)) superposed with an exponential decay function (see Equation (5)) representing the plasmasphere [7].

$$N_{eF}=Nm·exp\left(0.5\left(1-z-exp\left(-z\right)\right)\right)$$

(4)

$$N_{ep}=np·exp\left(-z_p\right)$$

(5)

$$z=\frac{h-hm}{Hs}, z_p=\frac{h-hm}{Hsp}$$

where Hs is the atmospheric scale height, Nm is the peak ionization, and hm is the corresponding height of the peak ionization of the Chapman layer. The quantity Hsp is the mean scale height of the plasma density. The quantity np is the plasmaspheric basic density of electrons. The ionosphere and plasmasphere model parameters are derived by an iterative approach by fitting the model functions to the data. For details of the method, we refer to the paper [6]. Figure 2 shows several examples of MetOp-A (data truncated at 600 km height) and COSMIC-2 collocated electron density (Ne) profiles. As a definition of co-located profile-pair, the tangent points of peak parameters within 300 km and 2 h of each other are used [8].

Figure 2 shows similar altitude variations for the collocated MetOp-A and COSMIC-2 Ne profiles. Six COSMIC-2 satellites are orbiting the Earth at an altitude of about 550 km in a low-inclination orbit. Therefore, we found multiple collocated COSMIC-2 Ne profiles for a single MetOp-A Ne profile. Although the collocation criteria are the same (tangential points of peak parameters within 300 km and 2 h) for those profiles, the geometry of the occultation plane (e.g., azimuth) is not the same in each case. Due to this, we see that the collocated COSMIC-2 Ne profiles are not the same but rather differ from one another.

Our investigation shows that using the developed method, a single electron density profile can be reconstructed within less than 1 min. Therefore, the computational cost/demand for the adaptive model approach is low. This makes the implementation of the RO inversion algorithm suitable for operational use in terms of memory use and computational power. The RO inversion method is applied to the MetOp-A truncated data, and numerous electron density (Ne) profiles are generated for the two test campaigns. A comprehensive validation study is accomplished by comparing MetOp-A reconstructed Ne profiles with other satellite data, such as Ne profiles from COSMIC-2, Fengyun-3D missions, and Ne data from many ground ionosonde stations. The quality of the MetOp-A Ne profiles is assessed, and the results are summarized in the next section after the assessment results of bending angle, scintillation indices, and topside TEC data quality.

4. GRAS Data Validation

4.1. Validation of Ionospheric Data Generated at EUMETSAT

The quality of the bending angle profiles, amplitude, and phase scintillation indices profiles generated at EUMETSAT is assessed by studying the systematic differences between similar data products from other LEO satellite missions. For this purpose, we determined pairs of co-located profiles from different LEO satellite missions. It is noted that the temporal, spatial, and azimuthal angle (of the occultation plane) differences between the two missions’ data may limit the use of co-located soundings to compute the accuracy of the relevant RO products. However, their mean differences will still be useful for probing stability, i.e., the re-producibility of the RO products from different instruments [8].

4.1.1. Assessment of Bending Angle Profiles

As Table 2 shows, during the MetOp-A experiment period, we could use bending angle data only from Fengyun-3D missions as validation datasets. Figure 3 shows the mean bending angle profile for the GPS L1 signal determined from the entire dataset of the co-located measurements of MetOp-A and Fengyun-3D during DOYs 176 to 254 of 2020. For determining the co-located pairs, an area of 2.5° × 2.5° in latitude by longitude and a time window of two hours are considered. The profiles in Figure 3 are generated by computing the average bending angle values binned with 1 km in altitude. The averaged bending angle profiles show expected distributions. Positive bending angles refer to the electron density increase associated with the ionospheric bottom side, and negative bending angles represent electron density exponential decay in the topside. The average peak height (bending angle close to zero) is, therefore, around 150 km, which is quite low since most of the colocations occurred in the polar region. The visual comparison shows that the MetOp-A and Fengyun-3D bending angle profiles have a similar vertical structure when the co-location period is limited to 2 h. The largest discrepancy to be highlighted occurs around the E-region peak hmE. It is worth noting that the hmE value mapped by FY3D is slightly lower than that of Metop-A, with a difference of around 10 km. This discrepancy is likely attributed to the broad area of colocation between the missions. Our investigation shows that if we increase the period of the co-location for 4 h, the similarity gets worse, as expected.

Since the 2-h co-location gives reasonable basis for the comparison, we created an analysis for distinct latitude regions considering the 2-h resolution. Figure 4 shows the mean and standard deviation (STD) of the differences between MetOp-A and Fengyun-3D bending angle profiles at GPS L1, considering the entire datasets of both extension campaigns. Comparing top and bottom panel plots, we see that the mean and STD differences are larger in the Northern hemisphere compared to the Southern hemisphere. The reason may be attributed to the seasonal influence being the experiment conducted during the Northern summer, which causes the ionospheric ionization to be higher in the Northern hemisphere. The occurrence of large differences is predominantly observed at low latitudes, where fewer collocated profiles were found. In contrast, at high latitudes, there were more collocations available, allowing for a better understanding of the smaller differences. Table 3 shows the mean and STD values of the differences between the bending angle observations obtained by MetOp-A and Fengyun-3D satellites. As we can see, the mean differences between MetOp-A and Fengyun-3D bending angles are larger in the Northern hemisphere, mainly represented by the standard deviation of the differences.

A comprehensive quality assessment of the bending angle, amplitude, and phase scintillation profiles generated at EUMETSAT is conducted by comparing MetOp-A data against independent measurements. This analysis has been recently summarized and published [9], and a details analysis of the results is out of the scope of the current paper.

4.1.2. Assessment of Scintillation Indices Profiles

As Table 2 shows, the GRAS ionospheric scintillations profiles are evaluated using measurements obtained from FORMOSAT-3/COSMIC-1 and FORMOSAT-7/COSMIC-2 satellites. The panels (a) and (b) of Figure 5 show the MetOp-A S4 profiles observed as a function of altitude and local time, respectively. The averages are computed using bins of 0.5 h in local time and 5 km in altitude. In the MetOp-A profiles, we can see the highest values around 100 km in the E-layer region. Additionally, we can see two main higher concentrations of S4 at two distinct local times. The first peak occurs around 12 LT, which is related to the daily maximum of ionization. The second peak of the highest S4 values occurs after 18 LT. It is noted that the MetOp-A satellite passes the Northern middle latitudes at around 6:00 and 18:00 local time (LT) hours, where and when the sporadic E-layers (Es) are more evident [9]. A literature study shows that the mid-latitude Es is visible only during summertime (being the highest during local noon hours) ([9] and references therein). Such scintillation occurs only at E-layer height (80–140 km).

In the case of COSMIC-2, however, typical daily distributions are not evident (see Figure 5b), most probably due to different orbit inclinations of COSMIC-2 satellites. The LT and latitude (see Figure 6) patterns of the COSMIC-2 S4 index are, therefore, different compared to MetOp-A S4 patterns.

In Figure 6, the S4 averages are computed using bins of 1 degree in latitude and 5 km in altitude. As expected, we see a higher concentration of S4 values around 100 km and close to 40° North again. The COSMIC-2 S4 data are also showing an increasing trend when approaching the Northern hemisphere. However, COSMIC-2 shows lower S4 values in comparison to MetOp-A S4 values.

As already explained, in the case of MetOp-A, the larger S4 around 20–40°N latitude (see Figure 6a) is caused by the sporadic E-layer at the tangential height of about 100 km. In the case of COSMIC-2 (Figure 6a), we do not see the signature of strong S4. The orbit inclination of MetOp-A is about 98.7°, whereas the COSMIC-2 has an inclination of about 24°. Since COSMIC-2 is flying at a very low inclination orbit, it makes the main difference between the ray path geometry received by the GRAS instrument onboard MetOp-A and the GNSS receiver onboard COSMIC-2 satellites. As described by Wu [10], the ray path geometry has a significant role in the occurrence of RO scintillations. The line of sight (LOS) from MetOp-A ROs tends to be more in the meridional direction, whereas the COSMIC-2 ray paths are more in between the meridional and zonal directions. It seems that COSMIC-2 shows weaker scintillation due to the favorable ray path geometry (neither meridional direction nor zonal direction rather falls between the two). For more details of the assessment results regarding the quality of MetOp-A scintillations products, we refer to [9].

4.1.3. GRAS Topside TEC Data Quality Indirect Assessment

The GRAS topside TEC data (denoted as tTEC hereafter) from POD measurements are calibrated in a simultaneous procedure of removing carrier-phase ambiguity terms and satellite and receiver biases (see Equation (1)) and provided by EUMETSAT. Since GPS satellite DCBs are taken from International GNSS Service (IGS) sources, the receiver DCBs and tTEC become highly correlated. Considering this, the quality of tTEC data is assessed by comparing MetOp receiver DCBs estimated by EUMETSAT and the same DCBs estimated by UPC-IonSAT (Universitat Politècnica de Catalunya) using a tomographic approach exclusively based on carrier phase data onboard LEOs implemented in the TOMION software [11].

The first comparison (not shown in Figure 7) of the receiver DCBs between EUMETSAT and TOMION (with either 2- or 4-layer model consisting of one layer or two layers above MetOp orbit height at 790 km) values during DOY 224 to 252 2020 shows (1) a systematic positive bias of 4–6 TECU of EUMETSAT DCBs compared to the TOMION DCBs and (2) the peak-to-peak variability of EUMETSAT DCBs is at the level of 5 TECU whereas the corresponding TOMION DCBs variability is at the level of up to 1 TECU which is expected [12]. Discussions within the study revealed that the main reason for such a large bias of EUMETSAT estimated DCBs was due to the presence of outliers or cycle slips in the GNSS data, which were not successfully removed. After applying a cycle slip detection and removal technique [13], the estimated DCB accuracy significantly improved, as shown in Figure 7a-c (see red plots in comparison to blue and cyan plots).

Daily EUMETSAT DCBs are estimated by applying the zero TEC technique as discussed by Zhong et al. [14], considering all the minimum sTEC computed by zenith data over the previous 7 days. Only data at high latitudes (>60 deg North/South), during night local times (from 18:00 to 6:00), and taken above 40 deg elevation are considered [12]. For more details about the GRAS topside TEC data quality assessment results, we refer to [12].

4.2. Validation of RO Electron Density Profiles

The quality of RO electron density profiles derived by using the Adaptive topside ionosphere/plasmasphere model technique briefly described in Section 3.3 (also in [6]) is assessed by investigating the systematic differences between similar products from other LEO satellite missions. For this purpose, we identified pairs of co-located electron density profiles from different LEO satellite missions. As a definition of co-located profile-pair, the occurrence of peak parameters within 300 km distance and 2 h of each other is used [8]. For example, Figure 8 shows the distribution of F2-layer peak electron density NmF2 and corresponding height hmF2 differences between MetOp-A and COSMIC-2 co-located profiles for the 600 km extension campaign. The mean and STD of differences are given in the plots as well.

Similarly, MetOp-A F2-layer peak density parameters are compared with Fengyun-3D, ionosonde, and also with the electron density data derived using the 1d-var ionospheric retrieval from MetOp-A bending angles [2,3], and statistical estimates are determined. The validation study is also repeated for the 300 km extension campaign.

Table 4. Mean and STD statistics of NmF2 and hmF2 differences when MetOp-A RO data are compared with COSMIC-2, Fengyun-3D, ionosonde, and MetOp-A 1d-var (using bending angle BA) data.

	600 km Truncation				300 km Truncation
	NmF2 (1011 el/m3)		hmF2 (km)		NmF2 (1011 el/m3)		hmF2 (km)
MetOp-A—Ref.	Mean	STD	Mean	STD	Mean	STD	Mean	STD
COSMIC2	0.07 (15%)	1.5 (56%)	0.04 (0.6%)	28.5 (11%)	−0.36 (1.5%)	1.76 (53%)	−3.1 (−0.4%)	26.6 (10%)
Fengyun−3D	0.03 (3%)	0.32 (16%)	12 (6%)	16 (7%)	−0.28 (8.9%)	0.49 (21.6%)	16.1 (7.6%)	19.3 (8.7%)
Ionosonde	0.05 (9%)	1.1 (41%)	−2.4 (−0.03%)	36 (13%)	−0.24 (−1.6%)	1.12 (41.8%)	−4.0 (−0.5%)	33.9 (12.4%)
MetOp-A BA	−0.06 (−1%)	0.7 (20%)	−0.5 (−0.1%)	11 (4%)	−0.4 (−2.5%)	1.34 (39%)	2.5 (1.2%)	16.8 (6.7%)

gives the summary of NmF2 and hmF2 comparisons in terms of residual statistics.

We found that the MetOp-A-RO NmF2 data deviated from the corresponding COSMIC-2 data by an average of about 15%. The mean deviation of MetOp-A NmF2 data from the Fengyun-3D data is 3%. As expected, the mean deviation of MetOp-A-RO NmF2 data from the MetOp-A-BA data is very small, which is 1%. However, the highest STD value is found when compared with COSMIC-2 data (1.5 × 10¹¹ el/m³). Comparison with Fengyun-3D data gives the minimum STD value of about 0.32 × 10¹¹ el/m³. The mean hmF2 deviation of MetOp-A-RO data from COSMIC-2, ionosonde, and MetOp-A-BA data is negligible, which is about −2.4 to 0.1 km. When compared with Fengyun-3D data, the mean deviation is found as 12 km. The largest STD values (e.g., 36 and 28.5 km) are found when compared with ionosonde and COSMIC-2 data. The ionosonde data are vertical-sounding measurements and are not affected by horizontal gradients, but they are not manually checked and scaled. Such factors may cause a bigger STD value. In both NmF2 and hmF2 cases, the STD values are much higher than the mean values. This indicates a large variability between profiles of different data types. The reason for large variability may be (1) the measurement geometry of the occultation plane differs, (2) the ionosphere may change within co-location criteria (i.e., 300 km apart and 2 h period), (3) truncation of MetOp-A data may contribute additional error when estimating the topside ionosphere as well as the peak density parameters.

Table 4 shows that the absolute mean and STD residuals for NmF2 are larger for the 300 km extension experiment compared to the 600 km extension campaign, as expected. However, hmF2 residuals are approximately similar for both experiments. So, it may happen that for most profiles, hmF2 values lie below the 300 km level, and the RO inversion can still retrieve the hmF2 height.

In Figure 9, the mean electron density differences and corresponding STD values along altitude between co-located profiles are plotted. The Ne profile data are pre-processed to compute electron density values from 90 km up to 550 km with a step size of 20 km. This is performed using spline interpolation of the original Ne profile data for MetOp-A (MTA), COSMIC-2 (FM7), and Fengyun-3D (FY3D) profiles.

The top panel of Figure 9 shows that the largest mean deviation and STD values are found as −0.12 × 10¹¹ el/m³ and 1.47 × 10¹¹ el/m³ at about 210 and 250 km height, respectively, when MetOp-A Ne profiles are compared with COSMIC-2 Ne profiles. Negative mean deviation at almost all altitudes indicates that MetOp-A Ne values are slightly larger than the COSMIC-2 Ne values. The maximum mean deviation is found to be similar, whereas the STD value (about 0.48 × 10¹¹ el/m³) is smaller (with respect to the COSMIC-2 case) when comparing co-located MetOp-A Ne profiles with the corresponding Fengyun-3D profiles. The reason for smaller STD values may be that the MetOp-A and Fengyun-3D Ne profiles are distributed over all latitudes instead of concentrated only at low latitude regions, which is the case for the COSMIC-2 data. Panel c of Figure 9 shows that for the 300 km campaign, the mean deviation is much higher compared to the 600 km campaign. Positive mean deviation almost at all altitudes indicates that MetOp-A Ne values are smaller than those of COSMIC-2 Ne profiles. Comparing panels c and d plots, we see that the residual statistics for Fendyun-3D data are smaller than those of COSMIC-2 data, as was the case for the 600 km campaign.

Moreover, the peak density parameters validation is performed for different local times (daytime, nighttime) and rising/setting occultations. The MetOp-A peak density parameters at E-layer height are also compared with corresponding data from COSMIC-2, Fengyun-3D, and MetOp-A 1d-var data. It is noted that the validation study provides consistent results in all these cases.

4.3. Ionospheric Data Assimilation

The quality of GRAS topside TEC data (tTEC) derived at EUMETSAT from POD measurements is assessed in two ways: (1) data assimilation runs into UPC’s multiTOMION [11] and impact evaluation and (2) data assimilation runs using DLR’s tomographic reconstruction technique [15,16] and impact evaluation.

The benefit of multiTOMION is its capability to combine different kinds of measurements (different GNSS systems in particular, different space geodetic systems in general) with different geometries (receivers –and transmitters- at different heights) and with different kinematics in a simple, precise way (based on carrier-phase measurements only). The tool multiTOMION is the recent modernization of TOMION software, also developed at UPC-IonSAT, and which is typically providing one of the best performing Global Ionospheric Maps (UQRG) in the IGS [17]. A potential drawback in scenarios with poor coverage of data and or poor geometries is the lack of a background model in multiTOMION. Knowing this potential drawback, the GRAS tTEC data are further assessed by data assimilation runs using DLR’s tomographic reconstruction technique, which uses a background topside ionosphere/plasmasphere model. In this approach, only the topside tTEC data are assimilated to the background model for topside ionosphere/plasmasphere reconstruction. The bottom part below the satellite orbit is not considered. The main benefit here is that the fitting procedure explicitly works on the topside data, and we have better chances to extract changes in the topside ionosphere and plasmasphere due to space weather impact. Our former investigation showed that such an approach could even detect plasmapause location as well as plasma blobs in the reconstruction. The investigation [15] using MetOp-A data during four ionospheric/geomagnetic storms shows that the approach can be used for space weather monitoring (e.g., detecting storm-related electron density enhancement and depletion afterward) during ionospheric storms. The main drawback is that the bottom side ionosphere is excluded in this approach. However, the multiTOMION approach would complement this drawback.

4.3.1. Data Assimilation Runs into TOMION and Impact Evaluation

The complete electron content above the LEO orbit is considered, described with two layers (see [12] for details), consistently with the limited vertical resolution associated with the predominantly vertical distribution of the POD GNSS line-of-sights. In this way, the topside vertical TEC is estimated in the tomographic forward Kalman filter, exclusively based on the dual-frequency carrier phase measurements provided by the POD GNSS receivers among the GNSS transmitter and receiver orbital information.

We can see representative plots of the electron density estimation provided by TOMION in Figure 10 for the Vertical Upper Electron Content (VUEC in TECU with heights above 790 km) as a function of latitude. The results are shown for the following combinations of input datasets:

–

Ground GPS input data only;

–

MetOp POD GPS data only;

–

MetOp POD and ground GPS data;

–

MetOp POD, ground GPS and DORIS geodetic dual-frequency input data (for details, see [12]).

Comparing plots in Figure 10, we see that the inclusion of DORIS measurements practically does not change the results. We can conclude that the use of MetOp POD measurements can provide by themselves a reasonable VUEC estimate with similar values to the corresponding combined ground GPS and MetOp POD measurements.

4.3.2. Assimilation of Topside TEC Using the DLR Tomographic Reconstruction Technique

The suitability of the topside TEC (tTEC) data for data assimilation is assessed by assimilating the EUMETSAT-provided tTEC for the MetOp-A satellite into the Neustrelitz Electron Density Model (NEDM2020, [18]). DLR routinely reconstructed the topside ionosphere and plasmasphere electron density distribution using CHAMP and GRACE satellites data (see [19]). Recently, DLR has adapted a tomographic reconstruction technique for assimilating COSMIC and MetOp satellite data into the background ionization (see [15,16]). The tomographic reconstruction is performed using a 3D grid specifically built to fit into the orbital geometry of the satellite. The tomographic algorithm is applied with the Algebraic Reconstruction Technique (ART) to invert tTEC measurements into electron density profiles in the topside ionosphere and plasmasphere. In the GIMA study, the GRAS tTEC data are assimilated to the background ionization created by NEDM2020 using the above-mentioned tomographic algorithm.

Figure 11, for instance, shows the VTEC values obtained after applying tomography to update the background representations. The top panel is related to the background, the middle panel to tomography results, and the bottom panels to the input TEC observations from the MetOp-A satellite. As can be seen, tomography VTEC results are closer to the MetOp-A values in comparison to the background, as expected. It is important to notice that, despite the MetOp-A distributions being rather dispersed, we kept a certain level of smoothness on the results to keep the large-scale patterns of the ionosphere rather than trying to represent small-scale variabilities. This gives robustness to our model to outliers since tomography is an ill-conditioned problem, which, in turn, is severely impacted by any residual noise in the input TEC values.

In order to conduct the assessment of the tomography results in comparison to DMSP in situ electron density data, we have run the developed method for the entire dataset using MetOp-A tTEC data as input and NEDM2020 as the background model. We have selected satellite number 17 from the DMSP mission for comparison. Figure 12 shows the comparison between the electron density latitudinal distribution at the height of the DMSP, as obtained directly by the DMSP-17 in situ observations (top panel), the tomography results assimilating Metop-A tTEC (middle panel) and the background model (bottom panel). It is clear that tomography has updated the background to represent the electron density with higher values, especially at low latitudes. At the same time, we see that DMSP electron densities are higher than the background values. This gives a good indication that the input tTEC values from MetOp-A were good enough to improve the background representations.

Table 5. Average error and standard deviation of the error between tomography and DMSP, as well as the background NEDM model and DMSP. Results are separated by latitudinal regions.

Latitude Range (deg)	Average TOMO (el/m3) × 1010	Std. TOMO (el/m3) × 1010	Average NEDM (el/m3) × 1010	Std. NEDM (el/m3) × 1010
−90 to −75	0.32	0.12	0.23	0.14
−75 to −60	0.34	0.14	0.23	0.14
−60 to −45	0.22	0.21	0.04	0.17
−45 to −30	0.00	0.27	−0.32	0.18
−30 to −15	−0.08	0.31	−0.43	0.16
−15 to 0	−0.34	0.37	−0.55	0.28
0 to 15	−0.79	0.39	−1.00	0.34
15 to 30	−0.57	0.39	−0.90	0.33
30 to 45	−0.45	0.32	−0.87	0.23
45 to 60	−0.29	0.21	−0.63	0.17
60 to 75	−0.19	0.17	−0.44	0.19
75 to 90	−0.37	0.23	−0.56	0.23

provides a statistical view of the tomography performance considering all analyzed days averaged over distinct latitudinal regions. As can be seen, the model performance is worse closer to the equatorial region and improves at mid- and high-latitudes. This is an expected pattern due to the higher level of ionization in the equatorial region. As we can see, there is a clear negative bias by the NEDM background model of around 0.45 × 10¹⁰ el/m³. This negative bias has been reduced in tomography by half, up to an average of 0.19 × 10¹⁰ el/m³. The standard deviation has slightly reduced from 0.45 to 0.44 × 10¹⁰ el/m³. The lower standard deviation of the background in a few cases mainly occurs due to the low level of Ne variabilities at the locations with high electron density values, as NEDM represents climatological patterns of the ionosphere.

MetOp-A assessment was extended by comparing the assimilation results of multiple MetOp satellites with Swarm-B topside TEC data. Our investigation shows clear gains when combining MetOp-A, B, and C data in the reconstructions. We found the mean error of about 0.92 and 0.32 TECU for the single satellite case and multi-satellite case, respectively. The important point observed in these results is the bias improvement by tomography, which gives a clear indication that the input TEC data was good enough to be used in 3D assimilation models aiming at the bias improvement of climatological models.

5. Summary

GIMA study results are summarized in the following. We found that different orbits/local time sampling made the validation against bending angles/scintillation index profiles from other RO instruments (on COSMIC-2 and Fengyun 3D) quite challenging. We found good agreement for collocated measurements and observed differences are explainable due to the different line of sight/local time sampling. We found a better capability of GRAS to monitor scintillations induced by the lower E and sporadic E layers thanks to better orbit inclination covering low and mid-latitude regions at sunset/sunrise and higher raw measurements sampling rates.

The investigation shows good agreement between the topside TEC data from MetOp-A and independent products (e.g., topside TEC measurements processed by UCAR, in situ electron densities from Swarm-B and DMSP satellites).

A new algorithm for inverting truncated ionospheric measurements well below the LEO orbit into electron density profiles is developed. In general, the validation against independent products shows a good agreement when the ionospheric peak density is covered by the measurements. So, the validation results show better performance for the 600 km extension than for the 300 km case. We found a very good agreement also with 1-d var results obtained from other inversion algorithms proposed in the meanwhile by ROM-SAF (Radio Occultation Meteorology Satellite Application Facility).

GRAS topside POD data from MetOp-A/B/C have been assimilated into a tomographic algorithm developed by UPC, showing very good performances and improving the results of the assimilation of ground-based data alone. Topside TEC from MetOp-A/B/C has been assimilated alone (or all together) into another tomographic algorithm developed by DLR. Fully consistent results have been found among different assimilation trials, all improving quite well the background model.

The study reveals that the GNSS topside and RO measurements obtained by the GRAS instrument onboard MetOp satellites are of good quality and become a valuable source of data for ionospheric research.