Optimal Estimation Inversion of Ionospheric Electron Density from GNSS-POD Limb Measurements: Part I-Algorithm and Morphology

Abstract

GNSS-LEO radio links from Precise Orbital Determination (POD) and Radio Occultation (RO) antennas have been used increasingly in characterizing the global 3D distribution and variability of ionospheric electron density (N_e). In this study, we developed an optimal estimation (OE) method to retrieve N_e profiles from the slant total electron content (hTEC) measurements acquired by the GNSS-POD links at negative elevation angles (ε < 0°). Although both OE and onion-peeling (OP) methods use the Abel weighting function in the N_e inversion, they are significantly different in terms of performance in the lower ionosphere. The new OE results can overcome the large N_e oscillations, sometimes negative values, seen in the OP retrievals in the E-region ionosphere. In the companion paper in this Special Issue, the HmF2 and NmF2 from the OE retrieval are validated against ground-based ionosondes and radar observations, showing generally good agreements in NmF2 from all sites. Nighttime hmF2 measurements tend to agree better than the daytime when the ionosonde heights tend to be slightly lower. The OE algorithm has been applied to all GNSS-POD data acquired from the COSMIC-1 (2006–2019), COSMIC-2 (2019–present), and Spire (2019–present) constellations, showing a consistent ionospheric N_e morphology. The unprecedented spatiotemporal sampling of the ionosphere from these constellations now allows a detailed analysis of the frequency–wavenumber spectra for the N_e variability at different heights. In the lower ionosphere (~150 km), we found significant spectral power in DE1, DW6, DW4, SW5, and SE4 wave components, in addition to well-known DW1, SW2, and DE3 waves. In the upper ionosphere (~450 km), additional wave components are still present, including DE4, DW4, DW6, SE4, and SW4. The co-existence of eastward- and westward-propagating wave4 components implies the presence of a stationary wave4 (SPW4), as suggested by other earlier studies. Further improvements to the OE method are proposed, including a tomographic inversion technique that leverages the asymmetric sampling about the tangent point associated with GNSS-LEO links.

1. Introduction

Ionospheric electron density (N_e) and its variability are critical for positioning, navigation, and timing (PNT) in the global navigation satellite system (GNSS) technology. The ionosphere also plays an important role in radio communication, navigation, and radar applications where space plasma can strongly influence electromagnetic wave propagation. Because of large ionospheric variability, with scales ranging from minutes to years in time and from sub-kilometer to thousands of kilometers in three-dimensional (3D) space, it is very challenging to characterize its structural changes with sufficient spatiotemporal sampling from observations on a global basis.

Empirical models, such as International Reference Ionosphere (IRI-2016) [1] and NeQuick [2], have a basic description of global 3D ionospheric N_e profiles and Total Electron Content (TEC). Near-real-time data analysis such as global ionosphere maps (GIMs) from the International GNSS Service (IGS) have been frequently used to validate the vertical TEC (vTEC) from empirical models [3,4,5] and other forecast models [6,7]. However, lack of detailed vertical information on N_e profiles leaves the modeled plasma structural variabilities poorly constrained within the ionosphere. As a result, the primary characterization of the ionosphere is currently limited to the peak electron density (NmF2), peak electron density height (HmF2), and TEC.

Until the recent proliferation of constellations of SmallSat/CubeSat low-Earth orbit (LEO) GNSS radio occultation (RO) satellites, spaceborne N_e observations have not been able to achieve sufficient spatiotemporal coverage for ionosphere observations. In particular, the constellations from the Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) [8], COSMIC-2 [9], and Spire [10] have achieved over 20,000 profiles/day for the lower ionosphere and ~12,000 profiles/day for the upper in 2022. The GNSS-RO system is recognized as a key technique for sounding atmospheric temperature and water vapor through atmosphere-induced bending [11]. In addition, the ionosphere-induced delay in RO L-band signals can be also used to measure the Ne profile with a good vertical resolution [12,13]. The GNSS-RO satellites are often equipped with two receiver antennas: one located at the top of the satellite, called the Precise Orbital Determination (POD) antenna, and the other on the side, called the RO antenna. As shown in Figure 1, for example, COSMIC-2 satellites have two sets of these antenna/receiver systems, situated in the fore and aft directions with respect to the satellite flight direction. The POD antenna, usually a low-gain antenna with a wide field of view (FOV), is used to provide satellite position and navigation with 1-Hz sampling through GNSS-POD links, whereas the GNSS-RO links are obtained from a high-gain and narrow-FOV antenna with a 50/100-Hz sampling rate. Other antenna configurations have been developed for modern GNSS receivers on LEO satellites to allow a wider range of applications, including GNSS-reflection (GNSS-R) for ocean surface winds [14] and flood inundation [15], and grazing angle reflection (GNSS-GR) for sea ice detection [16].

Both GNSS-POD and GNSS-RO links can be used to retrieve ionospheric N_e profiles from the phase delay associated with radio wave propagation between the GNSS transmitter and the LEO receiver. For the N_e profile retrieval, these links need to have measurements from negative elevation angles (ε < 0°) such as the radio wave propagation that transverses the bulk of the ionosphere below the satellite altitude [Figure 1]. The GNSS-RO links, by design, acquire most of their RO measurements at ε < 0°; but it is not always the case for the GNSS-POD links because their main purpose is to provide PNT for the LEO satellite by tracking GNSS signals above LEO (i.e., ε > 0°). However, the wide FOV of POD antennas often allows the extension of GNSS tracking to negative ε, which enables the N_e retrieval. As the POD-tracked GNSS signals occur at an angle away from its antenna boresight (i.e., ε < 0°), the signal strength begins to weaken, and is eventually lost at a low ε < 0°. As a result, not all GNSS-POD data are useful for the N_e retrieval. In fact, there are not many GNSS-POD links that have data acquired below the tangent height (h_t) of 100 km. On the other hand, the high-gain RO antenna can track GNSS signals even where the GNSS satellites are completely occulted behind the Earth’s shadow. However, most of the GNSS-RO operations have been focused on the neutral atmosphere, not acquiring data at h_t > 90 km. Although the GNSS-RO profiles yield little sampling for the F-region, they were found to be scientifically useful for inferring the D/E-region N_e [16]. In summary, both GNSS-POD and GNSS-RO measurements at ε < 0° can be used for the limb sounding of N_e, but the GNSS-POD links have a more complete coverage of the F-region ionosphere, whereas the GNSS-RO links are limited to the D/E-region N_e.

To retrieve the N_e profile from the GNSS-POD measurements, the onion-peeling (OP) inversion technique is commonly used in level-2 data processing [13]. The algorithm assumes spherical symmetry of ionospheric N_e around the tangent point (TP) of each limb sounding profile and uses the Abel vertical weighting functions to invert N_e from the horizontal TEC (hTEC) profile observed at ε < 0°. Sequentially carried out from the top to the bottom of spherical shells, the OP method solves the top layer first and then uses it as the a priori to solve the next layer below. The OP-derived N_e appears to work well for inferring HmF2 and NmF2 [17,18], but in the lower ionosphere, the OP-derived N_e profile often exhibits a large oscillatory and sometimes negative artifact [19,20].

To advance the N_e retrieval technique, in this study, we developed an optimal estimation (OE) method and applied it to all GNSS-POD data that have the ε < 0° hTEC measurements. The study has two main objectives. First, it seeks to use the OE algorithm to reduce the N_e oscillations in the lower ionosphere such that this portion of N_e profile is scientifically useful. The OE algorithm has the advantage of optimizing the error reduction for the entire F-region in the inversion process. Second, this study establishes an OE framework that will enable tomographic inversion of the N_e profiles in the presence of large ionospheric inhomogeneity. The paper is organized to detail the method in Section 2, and its performance from selected results in Section 3. Section 4 discusses impacts of ionospheric inhomogeneity on the current OE algorithm, which assumes the spherical symmetry. Finally, further improvements for the OE algorithm and future developments are proposed based on the new findings from this study.

2. Data and Methods

2.1. GNSS-POD and GNSS-RO hTEC Profiles

The hTEC measured at ε < 0° is a function of h_t, which is defined as the TP height along the straight line between GNSS and LEO satellites. Here, hTEC is used for the TEC from limb viewing and is differentiated from slant TEC (sTEC), which is often referred to as a slant path which is associated with an uplooking or downlooking view angle near the zenith/nadir. The hTEC(h_t) from GNSS-POD and GNSS-RO links can be derived from the excess phase profiles (${\varphi }_{exL1}$ and ${\varphi }_{exL2}$) at L1 ($f_1$ = 1.57542 GHz) and L2 ($f_2$ = 1.22760 GHz) frequencies by

$$\mathrm{h}\mathrm{T}\mathrm{E}\mathrm{C}(h_t)=({\varphi }_{exL1}(h_t)-{\varphi }_{exL2}(h_t))\frac{f_1^2{f}_2^2}{40.3\times {10}^{16}(f_2^2-f_1^2)}$$

(1)

Ionospheric bending is neglected in this calculation because it is a secondary effect compared to the phase advance induced the ionospheric plasma at h_t > 60 km. As the LEO satellite moves, the limb scan of the ionosphere produces a profile of hTEC(h_t) with respect to h_t [Figure 2], which is an integral along the GNSS-LEO link. hTEC(h_t) is a key input variable for the N_e inversion and can be expressed by

$$hTEC\left(h_t\right)={\int }_{LEO}^{GNSS}{N}_e·ds$$

$$={\int }_{h_t}^{Z_{GNSS}}{N}_e\left(z\right)·z/\sqrt{z^2-h_t^2}dz+{\int }_{Z_{LEO}}^{h_t}{N}_e\left(z\right)·z/\sqrt{z^2-h_t^2}dz$$

(2)

where Z_LEO is the LEO satellite altitude, and Z_GNSS is the GNSS satellite altitude. Z_GNSS is so high above the ionosphere that $Z_{GNSS}\approx \infty$ is a good approximation. From the perspective of the receiver in LEO, the first integral is from the far side of the TP, whereas the second integral is from the near side. The term, $z·dz/\sqrt{z^2-h_t^2}$, is known as the Abel weighting function. Within each GNSS-LEO limb sounding, the TP location does not move much, although the LEO satellite position may change significantly. This creates an asymmetric limb sounding between the far side and the near side, which can be used as the basis of a tomographic inversion as discussed later in Section 4.

Spherical symmetry in the N_e distribution and in the Abel weighting functions are assumed in all 1D retrieval algorithms. However, the LEO satellite motion can affect the pathlength (i.e., Abel weighting function) of the second integral to a significant extent. The OP method assumes the Abel weighting functions from the first and second interval are the same, which can potentially induce a large error and unstable N_e inversion in the lower ionosphere. In this study, although the near and far sides of the Abel weighting functions are treated differently, the same spherical symmetry is assumed for the N_e profile (i.e., a 1D N_e distribution). The OE algorithm can be modified to treat them differently in future improvements. As the first step, it is imperative to estimate the OE framework for the 1D N_e inversion before improving the algorithm for the 2D inversion.

2.2. hTEC Calibration

Absolute calibration of the hTEC measurement is required for the N_e retrieval. For example, the Spire hTEC products archived under NASA’s Commercial Smallsat Data Acquisition (CSDA) are not calibrated. Forsythe et al. [21] proposed a mitigation solution, using the differential hTEC from the far- and near-side hTEC, namely ΔhTEC = hTEC(-ε) − hTEC(ε), and assuming the same hTEC(ε > 0) values at both the far and near sides. This assumption may induce another error for hTEC calibration in the situations where the hTEC(ε > 0) values are different due to horizontal inhomogeneity.

The podTEC files produced by CDAAC contain vTEC from the measurements at ε > 0, which is the TEC above the satellite altitude converted from hTEC using sin(ε). An improved calibration scheme for hTEC measurements derived from the pseudorange of GNSS-LEO links is required, which has a typical navigation error of 1–3 m or 6–18 TECu [22].

In this study, we developed an empirical method for hTEC calibration, using the vertical derivative of the hTEC profile. The hTEC vertical gradient does not require absolute calibration, but it is linearly correlated with the absolute hTEC nearly universally at ε between −2° and −10°. As shown in Figure 3, this linear relationship between the limb-viewing hTEC and its vertical gradient at the top side of the ionosphere provides a valuable means for fast empirical calibration of hTEC on a profile-by-profile basis. For a satellite flying inside the ionosphere with the GNSS-LEO link for limb sounding, the observed hTEC vertical derivatives at elevation angles ε between −2° and −10° can be used to determine the absolute value of hTEC. Since the low flying LEO satellites often observe part of the ionosphere, it is a challenge to obtain very accurate calibration for both hTEC (or pseudorange) and PNT from a single link. Solving both variables accurately would require a least-squared solution to the measurements from multiple links in real time, which involves larger computing and data acquisition resources. Hence, the empirical model based on the linear relationships [Figure 3] has a practical advantage as a fast, ad-hoc method for hTEC calibration.

The empirical hTEC calibration model contains a set of tabulated coefficients that determine the slope of the relationship between hTEC and its derivative with respect to ε. These coefficients are also a function of the LEO satellite altitude. To generate a full lookup table for the hTEC calibration, we computed these derivative-to-hTEC coefficients using the IRI-2016 data. The residual error of this empirical model is <2 TECu. In reality, additional errors such as spatial inhomogeneity may contribute the final calibration uncertainty. The empirical model has 21 LEO satellite altitude bins between 400 and 800 km and 9 elevation angle bins between −2° and −10°. Because the GNSS-POD observations sometimes do not have all ε angles, this empirical calibration method would still work as long as there are hTEC measurements at ε between −2° and −10°. More importantly, it does not require the spherical symmetry assumption for hTEC(-ε) to be same as hTEC(ε).

The empirical calibration method was applied to all GNSS-POD data that contain hTEC measurements with ε < 0, including the data in COSMIC-1 and COSMIC-2 podTEC files from CDAAC. Although the hTEC data distributed by CDAAC are calibrated, we found some significant differences between two calibration methods, as shown in Figure 4. There appears to be an ~10 TECu offset between the two results, with the CDAAC values being higher. A similar high CDAAC bias (~4 TECu) was found by [23] in a comparison of the vTEC above 800 km between NeQuick-2 and COSMIC-1. We further grouped the bias between our empirical calibration and CDAAC hTEC according to LEO satellite altitudes but did not find any significant differences, which is the case for both COSMIC-1 and COSMIC-2 during the constellation formation period. We also sampled the data from an early and a later (mature) mission operation day, but no significant difference was found for either. This hTEC bias would affect the top of the retrieved N_e(z) profile, and the error impact reduces sharply at heights 10 km below the satellite altitude. For consistent F-region N_e retrievals, this study applied the same calibration method to all the COSMIC-1 and COSMIC-2 hTEC data, to compare with the Spire retrievals in the following sections. In addition to the recalibration for CDAAC hTEC data, we also need to correct an error of the reported LEO and GNSS satellite positions in the podTEC file, to obtain the right TP location and local solar time (LST).

2.3. OE Inversion Algorithm

The N_e inversion problem for the F-region is formulated in the same way as in the OE method used in the D/E-region retrieval [24], except that in this study, the input measurement vector y is the $\mathrm{h}\mathrm{T}\mathrm{E}\mathrm{C}$ derived from Equation (2), i.e., y $\equiv \left\{hTEC\left(h_t\right);h_t=h_1\dots h_{max}\right\}$, where $h_{max}$ is the maximum h_t below Z_sat. The state vector x is the N_e profile, i.e., x $\equiv \left\{N_e\left(z\right);z=100\mathrm{k}\mathrm{m},\dots ,z_{max}\right\}$. For the F-region N_e retrieval, $z_{max}$ is set at 800 km and the inverted N_e(z) profile above Z_sat is essentially an extrapolation with the information from the upper tail of the Abel weighting functions.

We used a linearized forward model for the measurement vector ($\mathbf{y}$) and the state vector ($\mathbf{x}$) that are related by the weighting function K, or the Jacobians as in $\mathbf{K}=\mathbf{\partial }\mathbf{y}/\mathbf{\partial }\mathbf{x}$, as follows:

$$\mathbf{y}={\mathbf{y}}_0+\mathbf{K}·\left(\mathbf{x}-{\mathbf{x}}_0\right)+{\mathit{\epsilon }}_{\mathit{y}}$$

(3)

where the linearization state vector is ${\mathbf{x}}_0$ (N_e in el/m³) and the associated measurement vector is ${\mathbf{y}}_0$ ($\mathrm{h}\mathrm{T}\mathrm{E}\mathrm{C}$ in TECu). The model is linearized on a single state vector ${\mathbf{x}}_0$, which is the annual mean of IRI-2016 in 2008 for the F-region [1] and the mean FIRI (Faraday-International Reference Ionosphere) E-region profile for the solar minimum condition [25].

The current version (V6p) of the OE algorithm assumes spherical symmetry and homogeneity in the N_e distribution, in which the inverted profile $\mathbf{x}$ depends only on altitude and is registered at the TP. The algorithm parameters used for the F-region N_e retrieval are similar to those in the similar version (V4) of the E-region retrieval [17], with a measurement uncertainty (2 TECu) added to measurement vector y, or ${\mathit{\epsilon }}_{\mathit{y}}$. The $\mathbf{x}$ inversion in the optimal estimation method can be expressed by

$$\widehat{\mathbf{x}}={\left[{\mathbf{S}}_{\mathbf{a}}^{-1}+{\mathbf{K}}^{\mathit{T}}{\mathbf{S}}_{\mathit{y}}^{-1}\mathbf{K}\right]}^{-1}\left[{\mathbf{S}}_{\mathbf{a}}^{-1}\mathbf{a}+{\mathbf{K}}^{\mathit{T}}{\mathbf{S}}_{\mathit{y}}^{-1}\mathbf{y}\right]$$

(4)

$${\mathbf{S}}_{\widehat{\mathbf{x}}}={\left[{\mathbf{S}}_{\mathbf{a}}^{-1}+{\mathbf{K}}^{\mathit{T}}{\mathbf{S}}_{\mathit{y}}^{-1}\mathbf{K}\right]}^{-1}$$

(5)

where $\widehat{\mathbf{x}}$ is an optimal solution to x in Equation (3) with its uncertainty determined by ${\mathbf{S}}_{\widehat{\mathbf{x}}}$; a is the a priori of x, which is set to be same as ${\mathbf{x}}_{\mathbf{0}}$ in this study; and ${\mathbf{S}}_{\mathbf{a}}={\mathit{\epsilon }}_{\mathit{x}\mathbf{0}}^{\mathbf{2}}·\mathit{I}$ and ${\mathbf{S}}_{\mathbf{y}}={\mathit{\epsilon }}_{\mathit{y}}^{\mathbf{2}}·\mathit{I}$ are covariance matrices for a and y, respectively. As in the D/E-region retrieval, ${\mathbf{S}}_{\mathbf{a}}$ is empirically determined to allow a stable inversion for the large N_e dynamic range and variability with respect to altitude.

3. Results

3.1. N_e Retrievals from OE/V6p and OP/CDAAC

The valid range of the OE/V6p N_e retrievals is between the LEO satellite altitude and the lowest h_t, or min(h_t), where the hTEC(h_t) measurements are available. Outside this range, the N_e values come largely from the extrapolation and their uncertainty is assigned to the negative value. In addition, the N_e data should be discarded if they are too close (within 10 km) to the satellite altitude. Thus, the negative N_e uncertainty and the height clearance may be used to screen the OE/V6p N_e data as a quality control. To ensure a successful N_e inversion, the OE/V6p algorithm imposed additional screening criteria for the hTEC measurements. No N_e retrieval is produced, if the min(h_t) of hTEC is greater than 110 km or no data exist at elevation angles between −2° and −10° for hTEC calibration.

The OE/V6p N_e retrieval algorithm has been applied to all COSMIC-1, COSMIC-2, FY3-C/D/E, and Spire GNSS-POD data that have hTEC measurements at ε < 0. Most of this dataset was obtained from the CDAAC distribution, except for the Spire data from NASA’s CSDA program and the FY3 data from the Fengyun Satellite Data Center. The hTEC data from CDAAC have their pseudorange calibration, which was processed by the CDAAC algorithm. In addition, CDAAC also provided the F-region N_e retrieval in the COSMIC-1 and COSMIC-2 level-2 data (ionPrf). Comparisons of the V6p and CDAAC COMSIC-1 N_e profiles are shown in Figure 5, whereas more can be found in a companion paper for the F2 peak height (HmF2) and density (NmF2) [26].

The OE/V6p N_e and OP/CDAAC N_e (aka. ionPrf) profiles generally agree with each other at the top side of the F-layer. However, significantly large differences are evident in the bottom side and in the E-region. As seen in Figure 5, the CDAAC N_e profiles in the lower ionosphere are often oscillatory, which can sometimes result in large negative N_e values in the E-region. This is an artifact of the OP retrieval because the method is particularly sensitive to the residual errors induced by the N_e retrievals from higher altitudes. On the other hand, the OE/V6p method is designed to balance the hTEC measurement uncertainties of all levels such that the bottom side of the F-layer is less prone to the top side error. To prevent the N_e retrievals in the E-region from oscillating, the OE algorithm allows the covariance values to be comparable to the variability of the a priori state vector in the E-region. As a result, the retrieved N_e profile is less spiky in the E-region compared with the OP results, although in some cases, the negative N_e values still exist.

Both OP and OE methods use the Abel weighting function in the N_e inversion process and assume a spherically homogeneous N_e profile. As discussed later in Section 4, this assumption can be problematic in an inhomogeneous ionosphere, as the radio waves trespass the ionosphere over a large vertical and horizontal spatial extent in limb sounding. The horizontal gradients are particularly large in the equatorial region, causing significant errors in the E-region N_e profiles [20,27]. Climatologically-constrained [21,28] and model-aided inversion schemes [29] were proposed to mitigate the artifact in the lower ionosphere, but these retrievals are subject to the a priori N_e climatology that might not be suitable for the highly-variable ionosphere on a daily basis.

3.2. Sensitivity to Ionospheric Variability

Diurnal N_e variations provide a valuable test of the OE algorithm since it is based on the forward model linearized on a single state vector, to cover a full dynamic range of N_e. It would be a problem if the single linearization model in Equation (3) were not sufficient to cope with the large diurnal and E-to-F variations. In that case, one would need to consider multiple linearization models from different state vectors, which could be computationally expensive. To evaluate this sensitivity concern, we compared the OE/V6p-retrieved N_e profiles over the full diurnal period, during which the ionosphere can change from zero to a large daytime maximum. If the OE algorithm were constrained too much to the a priori or its linearization state vector, the retrieval results would return the a priori with little variability. As seen in Figure 6, the OE algorithm is able to produce large N_e diurnal variations in both the E- and F-regions and the results agree generally well with the IRI-2016 model. The daytime N_e retrievals show a value varying by an order of amplitude greater than the linearization point, suggesting that the linearization model based on a single state vector is sufficient for capturing the ionospheric diurnal variations at all altitudes.

The equatorial diurnal variations also compare reasonably well between the OE and OP N_e retrievals in terms of hourly means [Figure 7]. Given that the retrievals are obtained from very different inversion and calibration schemes, this consistency confirms the robust characteristics seen in the N_e diurnal variations at these altitudes. For example, in the early morning hours, the enhancement at 200–300 km is weakly present in the OE and OP retrievals of COSMIC-2 N_e, both of which show a less pronounced morning overshoot than IRI-2016. The differences between the COSMIC-2 and IRI-2016 diurnal variations stress the importance of global GNSS-POD observations. Unlike the empirical IRI-2016 model, which is driven mostly by land-based ionosonde observations, new GNSS-POD data have a global coverage of lands and oceans, which help to reduce the sampling-induced bias significantly.

Quality control criteria used in the OE/V6p and OP/CDAAC algorithms may cause the differences seen in the daily number of N_e retrievals (e.g., 2465 at 250 km in Figure 6 vs. 1464 at 250 km in Figure 7). The OE/V6p algorithm can handle the hTEC data with acquisition slips or large ionospheric variability, which might cause a problem for the OP inversion. In addition, the OE method does not require uniform sampling in the hTEC profile. Also, as revealed in the scatters around each hourly bin, the OE/V6p data show fewer outliers than the OP/CDAAC data, especially in the lower ionosphere. Below 250 km, the OP/CDAAC N_e values at nighttime are lower than those of OE/V6p and spread more into the negative. The daytime CDAAC values at 150 km, on the other hand, are higher with a wider spread.

Compared with COSMIC-2 (covering latitudes between 44°S–44°N), the number of Spire N_e retrievals is much less in the tropics, because Spire is a constellation of polar-orbiting CubeSats with a full latitude coverage (90°S–90°N). As a result, gaps in LST sampling are evident for this day [Figure 8], showing missing measurements around the noontime and sparse sampling in the late afternoon hours. Nevertheless, the hourly mean N_e values from Spire agree well with the COSMIC-2 retrievals, suggesting that the hTEC calibration scheme developed in this study (Section 2.2) is consistent with the calibration provided at CDAAC. The new hTEC calibration method is readily applied to other GNSS-POD data that have the L1 and L2 excess phase profiles at negative elevation angles. making all GNSS-POD data available from these LEO receivers will help to greatly improve the spatiotemporal coverage of F-region N_e observations. Thus, it is recommended for all operational centers to archive and release the excess phase data from the GNSS-POD links, so that they can be processed for ionospheric N_e research and further improve the global spatiotemporal sampling.

3.3. Spire and COSMIC-2 N_e Sampling

This section highlights the new Spire and COSMIC-2 N_e results retrieved with the OE/V6p algorithm. Since the COSMIC-1 and COSMIC-2 N_e data have been analyzed intensively in previous studies [4,30], here we focus more on the Spire results. Spire Global (or Spire) operates a constellation of more than 100 CubeSats that acquire over 20,000 GNSS-RO and GNSS-POD profiles per day for Earth’s atmosphere [31] and ionosphere observations [10]. Under NASA’s CSDA agreement, these Spire data are archived for government-sponsored research to augment/complement the Earth observations acquired by NASA and other U.S. and international government agencies. Unlike the COSMIC-2 constellation that covers only latitudes between 44°S–44°N, Spire has a pole-to-pole coverage, which is critical for studying polar processes such as auroral electron precipitation and ionosphere-magnetosphere couplings.

The combined COSMIC-2 and Spire constellations have revolutionized the ionosphere observations with unprecedented global coverage and spatiotemporal sampling. Because the RO technique has been targeted mainly towards the atmospheric applications, there are more measurements for the atmosphere (GNSS-RO) than for the ionosphere (GNSS-POD). However, Wu [32] showed that the GNSS-RO data can also be used to derive the E-region N_e profile. As illustrated in

Table 1. Typical daily numbers of GNSS-RO and GNSS-POD profiles from COSMIC-1, COSMIC-2, and Spire constellations.

Number of Daily Profiles	GNSS-RO (Atmos&D/E-Region)		GNSS-POD (F-Region)
	Total L1B	V4 Ne	Total L1B	V6p Ne	CDAAC Ne
COSMIC-1 (1 September 2006)	2221	1697	3469	2570	2868
COSMIC-1 (31 December 2013)	1555	1307	1416	608	764
COSMIC-2 (31 December 2021)	6264	6160	9578	7225	4268
Spire (1 January 2022)	15,900	15,756	18,433	5977	N/A

, the number of daily measurements increases by ~7× and ~10× respectively in the GNSS-RO and GNSS-POD Level-1B (L1B) podTEC files, from the COSMIC-1 to the combined COSMIC-2 and Spire constellations. In the Level-2 products, the number of daily N_e profiles increases by ~8× and ~9× for the E-region (V4) and F-region (V6p) retrievals, respectively.

The improved GNSS-RO sampling from the COSMIC-2 and Spire constellations helps to overcome the sparse geographic distribution of ionospheric measurements from COSMIC-1. The new constellations can now produce 2-hourly statistics of the D/E-region N_e on a monthly basis for a 2° × 2° longitude-latitude grid [17], which is important for disentangling the seasonal variation from the diurnal variation. For the F-region, the number of GNSS-POD observations can yield the 2-hourly statistics on a 4° × 4° grid. High spatiotemporal sampling is also critically needed to better understand the ionospheric space-weather processes that often occur at a very short time scale and are associated with a strong coupling to the magnetospheric and atmospheric forcings. As shown in Figure 9, the number of hTEC measurements has been steadily increasing. More useful F-region N_e retrievals were produced after the development of COSMIC-2 and Spire constellations began to mature. Since 2021, the COSMIC-2 yield has been maintained at 70–80% while that of Spire is around 25–35%.

The maturation of the COSMIC-2 and Spire constellations is perhaps reflected in the time series of the daily number of GNSS-POD measurements, showing that the percentage of useful N_e retrievals increases gradually with time [Figure 9]. Compared to the Spire data in NASA’s CSDA archive, COSMIC-2 has a relatively higher percentage of L1B data that can be used for the N_e retrieval. To develop a constellation that spans over different LSTs, six COSMIC-2 satellites needed to lower their orbital altitude one by one from ~715 km to ~545 km. This was achieved over a period of one year. The lower orbital altitude also helps to increase the number of OE/V6p N_e retrievals, because the GNSS-POD links can obtain more observations from a negative elevation angle. The orbital maneuver was very different in developing the COSMIC-1 constellation. Instead of lowering the satellite orbit, it was decided to increase orbital altitude of the COSMIC-1 satellites from the initial ~525 km to ~810 km when the constellation was completed. Thus, the COSMIC-1 GNSS-POD links had fewer hTEC profiles in later years for the F-region observations.

3.3.1. Zonal Mean Morphology

The zonal mean N_e distribution for January 2022 is an example of the 2-hourly climatology that the Spire constellation can produce for characterizing the background ionosphere and its variabilities [Figure 10]. The Spire constellation now has enough sampling within a month to capture the global E-to-F-region transition of N_e distribution across the meridional plane from different LSTs. The OE/V6p algorithm is able to produce a reasonable F1-layer (140–200 km) that appears during daytime soon after the Sun rises. This layer is the most interesting for long-distance communications and has a similar behavior to the E-layer in terms of ionization processes. The F1-layer does not disappear rapidly in the evening hours after sunset, as it remains ionized till midnight before merging with the F2-layer. The latitude-dependent F1-layer N_e distribution is evident in the morning and evening hours.

The so-called equatorial “fountain effect” is clearly seen in the daytime zonal means, where the vertical ExB drift lifts the plasma to a higher altitude in the magnetic tropics [Figure 10]. The double-peak structure (e.g., LT = 21 h) is indicative of the “fountain effect”, whereby the plasma is redistributed along the magnetic field lines from the equator to higher latitudes in the F-region. Atmospheric tides and planetary waves can modulate the ExB drift considerably, creating significant longitudinal and LST variations to this redistribution [33,34].

3.3.2. Diurnal Variations

As expected for the sunlight-driven photoionization, ionospheric N_e variations are dominated by the diurnal availability of solar irradiance, which is often characterized by describing variables as a function of solar zenith angle (SZA) or LST. However, as shown in Figure 11, the diurnal variation is not symmetric about the local noon, suggesting that additional processes other than photolysis need to be considered, such as ion chemistry, cosmic rays, and lightning storms. In the upper ionosphere (z > 300 km), the N_e distribution is significantly skewed towards the evening hours, with a higher concentration in the magnetic tropics, compared with the middle and high latitudes. In the lower ionosphere, more symmetry with respect to LT is found in the N_e distribution, and the bulk of the N_e enhancement is confined within SZA < 90°. The so-called winter anomaly of daytime N_e is evident near 30 °N latitude in the lower ionosphere (110–210 km). A larger molecular-to-atomic ratio of the neutral atmosphere is thought to cause more ion loss in the summer hemisphere. In the lower ionosphere, auroral energetic electron precipitation (EEP) at high latitudes becomes increasingly important and shows significant diurnal variability. Significant differences can be found between two polar regions. The summer polar N_e from the auroral EEP is generally stronger than the winter, with a peak in the morning hours from January 2022. Similar variations are found for July 2022. There are two maxima from January 2022, peaking at LT = 11 h and 15 h, in the daytime N_e at the summer midlatitudes in the lower ionosphere.

Analyzing the COSMIC NmF2 and HmF2 data, Burns et al. [35] found that NmF2 is symmetric about the magnetic equator, which is consistent with the OE/V6p results for the upper-ionospheric N_e distribution. No winter anomaly was found in the NmF2 variation, which suggests that the primary control of NmF2 and upper-ionospheric N_e distribution is magnetic rather than geographic. An annual variation in HmF2 and a semiannual variation in NmF2 were found in their analysis [36], which is consistent with the OE/V6p observations, indicating that other forcings may play a significant role in regulating the N_e distribution.

3.3.3. Longitudinal Variations

To highlight the geographic N_e variability, we generated the N_e maps separately for daytime and nighttime at three selected heights (150, 270, and 510 km) for the lower, middle, and upper ionosphere [Figure 12, Figure 13 and Figure 14]. In these day and night maps, the N_e distributions observed by COSMIC-2 and Spire are generally consistent with each other, although different horizontal smearing in the two constellations may cause under/over-estimation of N_e amplitudes in various places. Because most of the Spire satellites are on a polar orbit, their GNSS-POD links tend to have a line-of-sight (LOS) toward the south-north direction. For the low-inclination COSMIC-2 satellites, their LOS directions are biased toward the east-west direction. As a result, limb smearing effects on a heterogeneous ionosphere would manifest themselves differently in these two constellations. More detailed discussions and quantitative error analysis about the spherical symmetry assumption can be found in Section 4.

Despite sampling and LOS differences in the Spire and COSMIC-2 constellations, they reveal a consistent longitudinal variation, suggesting that other processes can produce strong perturbations to the background N_e. Processes such as the geomagnetic field, atmospheric waves, and magnetospheric forcings can alter the N_e distribution with an amplitude similar to that of the photoionization process. Although the N_e distributions follow the magnetic field, as depicted by the dip angle contours, the large longitudinal variations along the same dip angle are comparable to, and, in some cases, greater than the amplitude modulated by the magnetic field. For example, at 150 and 270 km, the daytime N_e appears to be stronger in the America and Atlantic sectors at the northern midlatitudes for January 2022, but in the Pacific and Africa sectors at the southern midlatitudes. A wavenumber 4 (WN4) pattern is clearly seen in the daytime N_e at 510 km along the magnetic tropics.

Longitudinal variations of the nighttime N_e differ substantially from those observed during the daytime. At 270 km, there are large anomalies over the tropical Atlantic and Asia for January 2022, which resembles the distribution of equatorial plasma bubbles (EPBs) derived from COSMIC-1 S4 measurements [37,38]. Also evident in the middle and upper ionosphere is the so-called Weddell Sea Anomaly (WSA), where the summer nighttime N_e enhancement over the southern Pacific exceeds the daytime magnitude substantially. The primary driving mechanism of the WSA is thought to be the combined influences from neutral wind advection and magnetic declination [39].

In the polar regions, the EEP in the auroral oval has become a prominent feature in the lower ionosphere, especially in the E-region (<150 km), where the N_e from the solar photoionization diminishes. Note that the polar N_e enhancements more closely follow the magnetic dip angles of those derived from one Earth radius than those from 100 km, as expected, since the EEP process originates from the magnetosphere. Much stronger N_e enhancements at 150 km are observed in the summer pole than in the winter pole, which is consistent with the seasonal variation of auroral electron precipitation flux inferred from nitric oxide measurements [40], the field-aligned current (FAC) observed at the 780 km orbital altitude [41], and the D/E-region N_e measurements from GNSS-RO [17].

3.4. Frequency–Wavenumber Spectra of EIA

The unprecedented spatiotemporal sampling from the Spire and COSMIC-2 constellations allows a detailed analysis of the N_e variability over a wide range of frequency–wavenumber spectra. For example, the equatorial ionization anomaly (EIA) is a salient ionospheric structure along the geomagnetic equatorial band (15°S–15°N) [42]. It shows a persistent wavenumber-4 (WN4) feature in zonal winds [43], UV emissions [33], and ionospheric peak density, NmF2 [44]. The longitudinal variations at different spatial scales are thought to be due to the dynamo interaction between atmospheric tides with the daytime lower ionosphere (E-layer). The WN3 eastward-propagating non-migrating diurnal tide (DE3), induced by tropospheric forcing, has been suggested as a primary tidal component that is responsible for the WN4 variation [45]. However, other wave components can also contribute the apparent WN4 structure, including the semidiurnal eastward-propagating WN2 semidiurnal tide (SE2) and the stationary WN4 planetary wave (sPW4) [44,46,47,48].

In the following analysis, we followed the same naming convention as in [0], to express ionospheric oscillations as

$$A_{n,s}·\cos (n\mathsf{\Omega }t+s\lambda -{\varphi }_{n,s})$$

(6)

where t is universal time, Ω is the Earth’s rotation rate, λ is longitude, n is a subharmonic of a day, and s is zonal wavenumber with negative values indicating the eastward propagation. $A_{n,s}$ and ${\varphi }_{n,s}$ are the amplitude and phase of wave component (n,s), respectively. Wave components (n,s) are also named by letters with ‘D’, ‘S’, and ‘T’ for diurnal, semidiurnal and terdiurnal oscillations, and ‘E’ and ‘W’ as eastward- and westward-propagating waves. For example, the migrating diurnal tide is (1,1) or DW1, and the non-migrating eastward-propagating WN3 diurnal tide is (1,−3) or DE3.

To better characterize the spatiotemporal variability of EIA, we performed a frequency–wavenumber spectral analysis of the N_e data from 2021–2022. Since the EIA also varies with the magnetic field, the spectral analysis was carried out along the geomagnetic tropics (10°S–10°N), instead of the geographic tropics, to focus on longitudinal variations at the constant magnetic dip angles. A least-squares spectral analysis method, developed by Wu et al. [49], was used to process the data series from a non-uniform spatiotemporal sampling, which is the case for most GNSS-LEO constellations. The GNSS-POD sampling from Spire and COSMIC-2 constellations is nearly random in longitude and time. As discussed in [49], the random sampling in a time series helps to resolve a broad range of wave spectra without much aliasing to a distinct set of wave components associated with a regular spatiotemporal sampling.

Each frequency–wavenumber spectrum in Figure 15 represents the N_e variations in January 2022 at an altitude along the magnetic tropics. The spectrum was computed by scanning the N_e measurements over a period of one month using the least-squares method. The scan was carried out for one wave component each time at frequency intervals of 0.05 day⁻¹ in each wavenumber. As shown in Figure 15, the wave power spectra from Spire and COSMIC-2 data reveal a similar distribution for major wave components. Despite very different spatiotemporal samplings with the two constellations, the consistent spectra reveal fundamental N_e variability in the magnetic tropics. It also implies that GNSS-POD constellations and the receiver sensitivity are adequate for capturing these salient ionospheric variations. It is worth noting that the Spire spectra appear to be slightly noisier than those of COSMIC-2 at the same altitude, likely because of a maller number of samples in the tropics from Spire. Compared to Spire (21,899), COSMIC-2 has nearly three times more measurements (77,298) from Jan 2022, which helps to reduce spectral fluctuations in the frequency–wavenumber spectra.

The frequency–wavenumber spectrum of N_e varies considerably with altitude because of impacts of complex processes from F-region dynamics and forcings from the atmosphere below and the magnetosphere above. To illustrate the impacts of the neutral atmosphere, we computed the frequency–wavenumber spectrum with a similar approach to that for the N_e, but for a modeled atmospheric density from the Whole Atmosphere Community Climate Model with Ionosphere/Thermosphere eXtension (WACCM-X). This first-principles physics-based model simulates the whole atmosphere from the surface to the ionosphere/thermosphere. We analyzed the density data at 150 and 450 km from Specified Dynamics WACCM-x (SD-WACCM-X), a version that is nudged by the Modern-Era Retrospective Analysis for Research and Applications (MERRA) Reanalysis dataset from the surface to around 50 km [50,51], to incorporate a realistic forcing from the lower atmosphere. Salinas et al. [52] applied the same least-squares method in a SW2 study. The frequency–wavenumber spectrum of the modeled neutral density from January 2019 shows similar enhancements at (1,1), (2,2) and (3,3) components, as seen in the observed N_e, suggesting the importance of atmospheric forcings coupled to the ionosphere [Figure 16]. While the (1,−3) or DE3 component in modeled density is not pronounced in the observed N_e, the enhanced power at (2,−4) agrees with the observation. The modeled N_e is not used to compare with the observations, because it lacks these wave components. Further improvements are needed to realistically model ionospheric variabilities.

The N_e wave spectra exhibit significant seasonal and interannual variations, as expected due to the impacts from atmospheric/magnetospheric forcings and solar irradiance changes. Figure 17 illustrates such variations in the diurnal and semidiurnal components during 2021–2022. At 150 km, although the (1,1) component is relatively steady over seasons, the semidiurnal wave components exhibit significant seasonal variations with the (2,2) peaking near the equinoxes and the (2,−4) near the solstices. As the solar activity increases from 2021 to 2022, the other diurnal components appear to have more power in 2022 compared with 2021. The increase in the 2022 wave spectral power is evident at 300 and 450 km. At these higher altitudes, the (1,−3) or DE3 component begins to play an increasingly important role, showing peaks near the equinoxes, which is in phase with the (1,1) and (2,2) waves. Finally, it is worth noting that the WN4 components are present in both eastward and westward semidiurnal waves, although they are associated with different seasonal variations. In the case when they are equally important, it may indicate the presence of the stationary WN4 planetary wave (sPW4), as discussed in [48], since the linear combination of these two waves in the opposite propagating directions can yield a stationary wave.

4. Discussions

4.1. Impacts of Horizontal Inhomogeneity

Limb sounding is known for its coarse horizontal resolution along the LOS, which, with a strongly inhomogeneous ionosphere, would impact the OE or OP N_e retrievals. These algorithms assume spherical homogeneity (or 1D) of the N_e profile, i.e., N_e = N_e(z), and uses the Abel weighting functions to invert N_e(z) from the GNSS-POD hTEC(h_t). The retrieved N_e(z) profile is assigned to the TP location. Hence, the retrieved N_e(z) and the associated HmF2 and NmF2 variables are subject to errors induced by any 2D ionospheric inhomogeneity, since the limb sounding can span over 4000 km in the F-region. As shown in Figure 18, N_e may change substantially and the 1D retrieval sometimes is not a good representation of this variable ionosphere.

To quantify the impacts of ionospheric inhomogeneity, we simulated hTEC(h_t) profiles by integrating the 3D IRI-2016 N_e distribution along the GNSS-POD ray path at each h_t and used the simulated hTEC(h_t) as the input to the OE/V6p retrieval. Figure 18 illustrates the integrating calculation of hTEC(h_t) along the limb ray paths. Each hTEC(h_t) integration is carried out from the LEO retriever to the top of ionosphere. As in the observed GNSS-POD link, the integration excludes the near-side N_e component above the LEO height but includes the full vertical extent of the far-side N_e. Due to the fast (slow) motion of LEO (GNSS) satellites, the occultation ray paths are not symmetric about the TP, with a wider spread on the near side of the TP than on the far side.

The IRI-2016 N_e data from 12Z of December 2008 are used in the simulations, and two ray-path integrations were carried out in the simulation to mimic horizontal smearing in different directions. One set of integration is in the latitude direction whereas the other is in the longitude direction. The latitudinal smearing is an analogy to Spire GNSS-POD observations from high-inclination orbits, whereas the longitudinal smearing is an analogy to COSMIC-2. Different 2D effects are expected for the longitudinal and latitudinal integrations because the ionospheric horizontal N_e gradients can vary significantly with location and direction. Errors induced by the inhomogeneity are estimated by comparing the true N_e profile at the TP with the OE/V6p retrieval from the simulated hTEC(h_t).

The HmF2 and NmF2 errors induced by the ionospheric inhomogeneity are shown in Figure 19 and Figure 20. The cross-latitude gradients are largely determined by the F-region N_e distribution from the so-called fountain effect, whereas the cross-longitude gradients are driven primarily by the diurnal variation. Although the HmF2 and NmF2 retrieved from OE/V6p are mostly along the 1:1 line, non-negligible systematic errors are evident, which appear to be associated with the F-region N_e structure. In Figure 19, the HmF2 and NmF2 errors, which are induced by horizontal gradients across latitude, are generally larger during daytime, when the N_e in the F-region ionosphere is higher. As expected, the stronger the horizontal N_e gradient, the greater the errors in the OE/V6p retrieved HmF2 and NmF2. The 1D retrieval from OE/V6p tends to produce a higher (lower) daytime F2 peak in the daytime NH (SH) where HmF2 is lower (higher), and a lower (higher) NmF2 where it is higher (lower). These systematic errors occur from horizonal smearing by the limb sounding and by neglecting horizonal N_e gradient with the 1D algorithm.

Compared with the retrieval errors due to cross-latitude gradients, systematic HmF2 and NmF2 errors from cross-longitude gradients are relatively smaller [Figure 20]. This is expected for the amplitude of horizontal inhomogeneity gradients viewed from different directions. The cross-longitude errors occur in the region where the morning transition or horizontal N_e gradient in the ionosphere is the largest. For low-inclination orbiting satellites such as COSMIC-2, which have an occultation direction similar to the simulated cross-longitude smearing, the impact of ionospheric inhomogeneity is less serious than that for the polar-orbiting sensors (e.g., COSMIC-1 and Spire) that have a smearing in the cross-latitude direction. Thus, one would expect a larger systematic error in the HmF2 and f0F2 from COSMIC-1 and Spire.

In the lower ionosphere, the OE/V6p N_e is constrained to the linearization vector, as described by Equations (3) and (4), to prevent the retrieved profile from oscillating unstably in the presence of large horizontal inhomogeneity. As seen in Figure 5, Figure 6, Figure 7 and Figure 8, such a constraint appears to work effectively to reduce the unrealistic oscillations as well as the number of negative values in the lower ionosphere. However, negative values are still sometimes present in the OE/V6p algorithm, but occurring less frequently compared with the OP/CDAAC algorithm. Relaxing the constraint of the OE algorithm to the linearization vector would increase its response to lower-ionospheric variability, including occasional sporadic-E (Es) layers. These improvements can be achieved by further optimizing the inversion parameters in the OE algorithm.

In summary, impacts of ionospheric inhomogeneity on the 1D GNSS-POD N_e retrieval depend on the horizontal gradient of N_e distribution along the limb path. Thus, the COSMIC-1 and Spire constellations, of which the satellites have a higher inclination angle, would suffer from greater errors, compared with the low-inclination COSMIC-2 constellation. These simulations were based on the IRI-2016 climatology of a quiet background ionosphere, whereas ionospheric inhomogeneity is likely great in reality.

4.2. LST Sampling Differences Between COSMIC-2 and Spire

In addition to spatial inhomogeneity effects, COSMIC-2 and Spire sampling differences in LST may also contribute to some of the discrepancies seen in Figure 12, Figure 13 and Figure 14, especially near the turning latitudes of COSMIC-2. As shown in Figure 21, the LST sampling of monthly COSMIC-2 data is not uniform at the 35°–44° latitudes, and the Spire monthly statistics still have gaps in LST-latitude sampling, despite the constellation having a large number of CubeSats. The uneven sampling of COSMIC-2 at the tuning latitudes was likely due to attitude-dependent variations of GNSS-POD links. As previously mentioned for the quality control in OE/V6p data processing, the algorithm requires that the min(h_t) of hTEC profiles needs to be greater than 110 km. The requirement would lead to systematic discrimination on the GNSS-POD links at latitudes that are not favorable for COSMIC-2 POD antennas to reach the limb h_t of 110 km.

In the Spire case, the average daily number of OE/V6p N_e retrievals is still less than that from COSMIC-2, because of its low yield rate. In other words, the Spire sampling still has gaps in the 2h-2° LST-latitude coverage on a monthly basis. The gaps in the current Spire constellation are mostly at mid-latitudes near the terminator, but they appear to be consistent from month to month. To overcome the LST sampling limitation, a merged N_e dataset is being generated using the COSMIC2, Spire, and FY3-C/D/E constellations, which will eliminate most of the data gaps in the monthly LST-latitude coverage.

4.3. Tomographic Inversion for Inhomogeneous Ionospheres

To minimize the 2D inhomogeneity effects from the limb sounding, auxiliary data sets, such as a 3D climatology data set [28] and an ionosphere model data set [29], were introduced as the a priori of the N_e distribution. This a priori distribution would determine the partition between the near and far sides of N_e profiles, so that the retrieved 1D profile is a better representation of the N_e at the TP. These correction approaches, however, are static and biased to the climatology/model data set used, and therefore lack the skill of measuring the real N_e variability at the TP.

Here, we propose a non-conventional tomographic inversion technique for mitigating inhomogeneity effects in the GNSS-POD limb sounding. As shown in Figure 18, the GNSS-POD ray paths are asymmetric about the TP, which provides an opportunity to disentangle the horizonal inhomogeneous ionosphere to some extent. The 2D information about the ionosphere is mostly in the hTEC(h_t) profile below the F2 peak (HmF2). As shown in Figure 22, the simulated hTEC(h_t) profiles differ significantly if we flip the horizonal N_e gradient as used in Figure 18. In this flip, the N_e profile at the TP is identical, but the near and far sides of the N_e profiles are different, and so are the hTEC(h_t) profiles in these flipped cases.

These differences in the hTEC(h_t) profiles provide the basis for a potential tomographic retrieval of N_e along the observing plane. In particular, the hTEC(h_t) differences below HmF2 are large enough to create the sensitivity to the N_e distribution in the near and far sides of the TP, which can be used to form a 2-element tomographic inversion. This technique appears to produce a very promising result for retrieving two TECs at the near and far sides of the TP (not shown). To carry out this tomographic inversion for N_e, one would need to evaluate the weighting function of hTEC(h_t) carefully to make sure that there is enough degrees of freedom for the extended dimensions.

It is one of the advantages of the OE algorithm to extend the hTEC-to-N_e inversion from a 1D-N_e case, which assumes the spherical homogeneity, to a 2D-N_e inversion using the tomographic technique. Unlike the conventional tomographic inversion problem, the proposed double-N_e-profile scheme (at the near and far sides of the TP) would make use of the hTEC(h_t) measurements from a single limb scan in the 2D plane. The conventional tomographic techniques have multiple consecutive scans from different locations and angles, which are combined to constrain the 2D variable field in the retrieval. In the GNSS-POD link case, the vertical sampling densities by hTEC(h_t) are different between N_e profiles from the near and far sides of the TP. A fast-moving LEO satellite makes the near-side N_e sampling coarser than that of the far side, which is closer to a GNSS satellite. As a result, small differences in the hTEC(h_t) profile can be used to infer the horizonal gradient of the N_e profile.

Because the OE algorithm is readily scalable under the same mathematical framework [i.e., Equations (3) and (4)], it can include both GNSS-POD and GNSS-RO links in the measurement vector if these are observed closely in time. The GNSS-RO measurements, as described by Wu et al. (2022), provide a good constraint to the E-region N_e. Combining GNSS-POD and GNSS-RO measurements for retrieving a more complete ionospheric N_e profile would be very valuable in future algorithm development. Because the E-region N_e is prone to the residual errors from the F-region N_e result, a sequential OE method would be desired to first retrieve the E-region N_e and then the F-region N_e with a tight constraint on the E-region N_e.

Finally, the OE/V6P algorithm does not use the hTEC data from ε > 0, which can provide the vTEC information above the LEO satellite altitude. Combining the vTEC with the near-side partial TEC from the tomographic retrieval will produce a complete TEC at the near side. Instead of correcting the TEC or N_e at the TP, the future OE algorithm is aimed at producing these measurements at the near and far sides of the TP, which may help to improve the horizontal resolution from ~3000 km to ~1500 km.

5. Conclusions and Future Work

In this study, we successfully developed and applied an optimal estimation (OE) inversion algorithm to GNSS-POD measurements to retrieve N_e profiles at 100–500 km from three SmallSat/CubeSat constellations (COSMIC-1, COSMIC-2, and Spire). An empirical hTEC calibration scheme was also developed, using its own vertical gradients at elevation angles between −2° and −10°, to yield a quick, self-service absolute calibration of hTEC on a profile-by-profile basis. The empirical method was applied to both Spire and COSMIC-1/2 data. A ~10 TECu bias was found between the empirically calibrated and CDAAC-calibrated COSMIC-2 hTEC data. The hTEC differences resulting from these calibration methods warrant further evaluation using independent measurements.

The OE/V6p algorithm demonstrated that the N_e retrieval has a quality similar to those from the onion-peeling (OP) method in the upper ionosphere. In the lower ionosphere, the OE/V6p retrieval showed fewer negative N_e values and oscillations, which make the data more reliable and scientifically useful. The Spire and COSMIC-2 N_e data from 2021–2022, retrieved with the OE/V6p algorithm, showed an unprecedented spatiotemporal coverage of the ionosphere, which allows characterization of its morphology in great detail for the diurnal, latitudinal, and longitudinal variations. For the first time, a full frequency–wavenumber spectrum of the EIA was obtained to unambiguously resolve both eastward- and westward-propagating wave components with eight wavenumbers and periods as short as the terdiurnal.

Impacts of the spherical symmetry assumption were evaluated globally by comparing the OE/V6p N_e profiles retrieved from the simulated hTEC data using a realistic heterogeneous ionosphere. As expected for the limb smearing effect, the largest impacts occur in the regions where the limb LOS passes across a strong horizontal gradient. As a result, different satellite constellations would produce different error amplitudes even though the ionospheric inhomogeneity is the same.

The new algorithm under the OE architecture allows further improvements for the N_e retrieval, for example, incorporating the tomographic inversion technique. Because the ray paths of the GNSS-POD link are asymmetric about the TP, this limb sounding has sensitivity to the N_e profiles on both the near and far sides of the TP. Thus, the OE algorithm can be extended to simultaneously retrieve two N_e profiles from each GNSS-POD hTEC profile, so as to improve the horizontal resolution N_e measurements from ~3000 km to ~1500 km.

Lastly, most operational GNSS-POD observations contain measurements from negative elevation angles, but these data have not been made available except for COSMIC-1 and COSMIC-2 L1b data from CDAAC, and the Spire L1b data delivered to NASA’s CSDA program. In light of the great potential of the POD limb data for monitoring the global ionosphere and its variability, as demonstrated in this study, it is recommended that GNSS-POD data providers ought to make all hTEC data available for scientific research and further algorithm/product developments.