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

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 (Ne). In this study, we developed an optimal estimation (OE) method to retrieve Ne 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 Ne inversion, they are significantly different in terms of performance in the lower ionosphere. The new OE results can overcome the large Ne 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 Ne morphology. The unprecedented spatiotemporal sampling of the ionosphere from these constellations now allows a detailed analysis of the frequency– wavenumber spectra for the Ne 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 eastwardand 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.


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.
urements. The study has two main objectives. First, it seeks to use the OE algorithm to reduce the Ne oscillations in the lower ionosphere such that this portion of Ne 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 Ne 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.

Figure 1.
A schematic showing two sets of the POD (red) and RO (cyan) L-band receiver antennas onboard a LEO satellite (center black box) that have the fore and aft views in the flight direction (arrow). Their antenna patterns are drawn by circular (red) and triangle (cyan) lines to indicate a wide and narrow FOV respectively. Both POD and RO antennas can be used for limb sounding of the ionosphere from the measurements with negative elevation angles (ε < 0°). The horizontal dashed line, denoted for ε = 0°, separates the uplooking and downlooking (i.e., limb sounding) scenarios. The RO or modified antennas are used to receive GNSS signals from slant angle reflection (GNSS-R) and grazing angle reflection (GNSS-GR). 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.

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 (φ exL1 and φ exL2 ) at L1 ( f 1 = 1.57542 GHz) and L2 ( f 2 = 1.22760 GHz) frequencies by 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 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 ≈ ∞ 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/ z 2 − h 2 t , 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.

hTEC Calibration
Absolute calibration of the hTEC measurement is required for the Ne retrieval. For example, the Spire hTEC products archived under NASA's Commercial Smallsat Data Acquisition (CSDA) are not calibrated. Forsythe et al. [Error! Reference source not found.] proposed a mitigation solution, using the differential hTEC from the far-and near-side hTEC, namely hTEC = hTEC(-ε) − hTEC(ε), and assuming the same hTEC(ε  ) values 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.

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 limbviewing 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(ε). 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 [24] 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 Ne(z) profile, and the error impact reduces sharply at heights 10 km below the satellite altitude. For consistent F-region Ne 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).  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).

Figure 4.
Calibration differences between the empirical algorithm in this study and CDAAC hTEC at the top of hTEC profiles. Four different days from COSMIC-1 and COSMIC-2 are selected to compare potential impacts from satellite altitude and operation. The data with satellite altitudes <600 km are colored in red. The data from 1 January 2007 and 1 January 2020 reflects the early period of COSMIC-1 and COSMIC-2 operation, respectively, when they have satellites in very different altitudes in the constellation.

OE Inversion Algorithm
The Ne inversion problem for the F-region is formulated in the same way as in the OE method used in the D/E-region retrieval [17], except that in this study, the input measurement vector y is the hTEC derived from Equation (2), i.e., y ≡ {ℎ (ℎ ); ℎ = ℎ 1 ••• ℎ }, where ℎ is the maximum ht below Zsat. The state vector x is the Ne profile, i.e., x ≡ { ( ); = 100 km,•••, }. For the F-region Ne retrieval, is set at 800 km and the inverted Ne(z) profile above Zsat 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 ( ) and the state vector ( ) that are related by the weighting function K, or the Jacobians as in = / , as follows: where the linearization state vector is (Ne in el/m 3 ) and the associated measurement vector is (hTEC in TECu). The model is linearized on a single state vector , which is the annual mean of IRI-2016 in 2008 for the F-region [Error! Reference source not found.] . Calibration differences between the empirical algorithm in this study and CDAAC hTEC at the top of hTEC profiles. Four different days from COSMIC-1 and COSMIC-2 are selected to compare potential impacts from satellite altitude and operation. The data with satellite altitudes <600 km are colored in red. The data from 1 January 2007 and 1 January 2020 reflects the early period of COSMIC-1 and COSMIC-2 operation, respectively, when they have satellites in very different altitudes in the constellation.

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 hTEC derived from Equation (2) where h max is the maximum h t below Z sat . The state vector x is the N e profile, i.e., x ≡ {N e (z); z = 100km, . . . , z max }. 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 (y) and the state vector (x) that are related by the weighting function K, or the Jacobians as in K = ∂y/∂x, as follows: where the linearization state vector is x 0 (N e in el/m 3 ) and the associated measurement vector is y 0 (hTEC in TECu). The model is linearized on a single state vector 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 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 ε y . The x inversion in the optimal estimation method can be expressed bŷ wherex is an optimal solution to x in Equation (3) with its uncertainty determined by Sx; a is the a priori of x, which is set to be same as x 0 in this study; and S a = ε 2 x0 ·I and S y = ε 2 y ·I are covariance matrices for a and y, respectively. As in the D/E-region retrieval, S a is empirically determined to allow a stable inversion for the large N e dynamic range and variability with respect to altitude.

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.
pass 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 Ne profiles [Error! Reference source not found.,27]. Climatologically-constrained [21,28] and model-aided inversion schemes [Error! Reference source not found.] were proposed to mitigate the artifact in the lower ionosphere, but these retrievals are subject to the a priori Ne climatology that might not be suitable for the highlyvariable ionosphere on a daily basis. Figure 5. Example of the COMSIC-1 Ne profiles from 2006d244 when both OE/V6p (red) and OP/CDAAC (black) produced a retrieval. The OE/V6p results outside the recommended vertical range are colored in blue (see text). The dashed line and shaded width are, respectively, the linearization state vector and its uncertainty as described in Equation (3). The OE/V6p ̂ retrieval has a 2 km vertical resolution.

Sensitivity to Ionospheric Variability
Diurnal Ne 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 Ne. 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 Ne 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 Ne diurnal variations in both the E-and F-regions and the results agree generally well with the IRI-2016 model. The daytime Ne retrievals show a value varying by an order of amplitude greater than the linearization point, suggesting that the linearization model based on 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 modelaided 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.

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 Eand 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 Ne 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 Ne 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 Ne, 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 Ne 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 Ne values at nighttime are lower than those of OE/V6p and spread 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 Ne retrievals is much less in the tropics, because Spire is a constellation of polarorbiting 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 Ne 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 Ne 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 Ne research and further improve the global spatiotemporal sampling. 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 polarorbiting 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.

Spire and COSMIC-2 Ne Sampling
This section highlights the new Spire and COSMIC-2 Ne results retrieved with the OE/V6p algorithm. Since the COSMIC-1 and COSMIC-2 Ne 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 [Error! Reference source not found.] and ionosphere observations [Error! Reference source not found.]. 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 ionospheremagnetosphere 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 [Error! Reference source not found.] showed that the GNSS-RO data can also be used to derive the E-region Ne profile. As illustrated in Table 1, 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 COS-MIC-2 and Spire constellations. In the Level-2 products, the number of daily Ne profiles increases by ~8× and ~9× for the E-region (V4) and F-region (V6p) retrievals, respectively. Table 1. Typical daily numbers of GNSS-RO and GNSS-POD profiles from COSMIC-1, COSMIC-2, and Spire constellations.

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, 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%.

Zonal Mean Morphology
The zonal mean Ne 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 Ne 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 F2layer. The latitude-dependent F1-layer Ne distribution is evident in the morning and evening hours. 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 tõ 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 tõ 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.

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.
Remote Sens. 2023, 15, x FOR PEER REVIEW 15 of 33 Figure 10. Monthly zonal mean Ne from Spire as a function of magnetic latitude and LST from January 2022. A 2-hourly LT bin is needed for the Spire constellation to produce good statistics. The retrievals are mostly valid up to ~500 km because Spire 3U CubeSats were mostly operated at an orbital altitude of ~525 km.

Diurnal Variations
As expected for the sunlight-driven photoionization, ionospheric Ne 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 Ne 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 Ne distribution, and the bulk of the Ne Figure 10. Monthly zonal mean N e from Spire as a function of magnetic latitude and LST from January 2022. A 2-hourly LT bin is needed for the Spire constellation to produce good statistics. The retrievals are mostly valid up to~500 km because Spire 3U CubeSats were mostly operated at an orbital altitude of~525 km.
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].

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 molecularto-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. NmF2 and upper-ionospheric Ne distribution is magnetic rather than geographic. An annual variation in HmF2 and a semiannual variation in NmF2 were found in their analysis [Error! Reference source not found.], which is consistent with the OE/V6p observations, indicating that other forcings may play a significant role in regulating the Ne distribution.

Longitudinal Variations
To highlight the geographic Ne variability, we generated the Ne maps separately for daytime and nighttime at three selected heights (150, 270, and 510 km) for the lower, middle, and upper ionosphere [ Figures 12-14]. In these day and night maps, the Ne distributions observed by COSMIC-2 and Spire are generally consistent with each other, although 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.

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 [ Figures 12-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 eastwest 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.    Figure 12, but for 270 km. Figure 13. As in Figure 12, but for 270 km.
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]. Remote Sens. 2023, 15, x FOR PEER REVIEW 20 of 33 Figure 14. As in Figure 12, but for 510 km.

Frequency-Wavenumber Spectra of EIA
The unprecedented spatiotemporal sampling from the Spire and COSMIC-2 constellations allows a detailed analysis of the Ne variability over a wide range of frequencywavenumber spectra. For example, the equatorial ionization anomaly (EIA) is a salient ionospheric structure along the geomagnetic equatorial band (15°S-15°N) [Error! Reference source not found.]. It shows a persistent wavenumber-4 (WN4) feature in zonal winds [Error! Reference source not found.], UV emissions [Error! Reference source not found.], and ionospheric peak density, NmF2 [Error! Reference source not found.]. 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 [Error! Reference source not found.]. 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 [Error! Reference source not found.], to express ionospheric oscillations as , • cos ( Ω + − , ) (6) Figure 14. As in Figure 12, but for 510 km.
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].

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 frequencywavenumber 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Ωt + sλ − φ 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 φ 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 eastwardpropagating WN3 diurnal tide is (1,−3) or DE3.
To better characterize the spatiotemporal variability of EIA, we performed a frequencywavenumber 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 −1 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.
sons, 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.   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.

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 Ne retrievals. These algorithms assume spherical homogeneity (or 1D) of the Ne profile, i.e., Ne = Ne(z), and uses the Abel weighting functions to invert Ne(z) from the GNSS-POD hTEC(ht). The retrieved Ne(z) profile is assigned to the TP location. Hence, the retrieved Ne(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, Ne 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(ht) profiles by integrating the 3D IRI-2016 Ne distribution along the GNSS-POD ray path at each ht and used the simulated hTEC(ht) as the input to the OE/V6p retrieval. Figure 18 illustrates the integrating calculation of hTEC(ht) along the limb ray paths. Each hTEC(ht) in-

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.
NmF2 errors, which are induced by horizontal gradients across latitude, are generally larger during daytime, when the Ne in the F-region ionosphere is higher. As expected, the stronger the horizontal Ne 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 Ne gradient with the 1D algorithm.  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 Figures 19 and 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. Remote Sens. 2023, 15, x FOR PEER REVIEW 26 of 33 Figure 19. Differences in HmF2 (in km) and f0F2 (in MHz) between the IRI-2016 model truth and retrieved with V6p from the simulation. NmF2 = 1.24 × 10 10 · (foF2) 2 is used in the NmF2-f0F2 conversion, where NmF2 is el/m 3 and f0F2 is MHz. The mean and standard deviation of retrieval-model differences are shown in the scatter plots.
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 Ne 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 Ne 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 Figures 5-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

HmF2 Truth
HmF2 Diff (V6p -Truth) NmF2 Truth NmF2 Diff (V6p -Truth) Figure 19. Differences in HmF2 (in km) and f0F2 (in MHz) between the IRI-2016 model truth and retrieved with V6p from the simulation. NmF2 = 1.24 × 10 10 · (foF2) 2 is used in the NmF2-f0F2 conversion, where NmF2 is el/m 3 and f0F2 is MHz. The mean and standard deviation of retrievalmodel differences are shown in the scatter plots.
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.
inversion parameters in the OE algorithm.
In summary, impacts of ionospheric inhomogeneity on the 1D GNSS-POD Ne retrieval depend on the horizontal gradient of Ne 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. Figure 19, but for the simulations of cross-longitude inhomogeneity from the IRI-2016 model.

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 Figures 12-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(ht) 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 ht of 110 km.
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 Figures 5-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.

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 Figures 12-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 Ne 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 Ne 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.

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 Ne distribution. This a priori distribution would determine the partition between the near and far sides of Ne profiles, so that the retrieved 1D profile is a better representation of the Ne 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 Ne 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(ht) profile below the F2 peak (HmF2). As shown in Figure 22, the simulated hTEC(ht) profiles differ significantly if we flip the horizonal Ne gradient as used in Figure 18. In this flip, the Ne profile at the TP is identical, but 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.

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.
the tomographic technique. Unlike the conventional tomographic inversion problem, the proposed double-Ne-profile scheme (at the near and far sides of the TP) would make use of the hTEC(ht) 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(ht) are different between Ne profiles from the near and far sides of the TP. A fast-moving LEO satellite makes the nearside Ne sampling coarser than that of the far side, which is closer to a GNSS satellite. As a result, small differences in the hTEC(ht) profile can be used to infer the horizonal gradient of the Ne 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 Ne. Combining GNSS-POD and GNSS-RO measurements for retrieving a more complete ionospheric Ne profile would be very valuable in future algorithm development. Because the E-region Ne is prone to the residual errors from the F-region Ne result, a sequential OE method would be desired to first retrieve the E-region Ne and then the F-region Ne with a tight constraint on the E-region Ne.  Figure 18, with the same inhomogeneous Ne field but as observed from two opposite directions that yield the horizontally-flipped sampling with respect to the TP (a,b). The solid and dashed lines in (c) denote the hTEC(ht) profiles simulated for the two viewing directions. The Ne profile at the TP location [red line in (c)] is same in both cases. The dashed lines in (a,b) denote the location of two state vectors that may be retrieved from hTEC(ht) in the presence of horizonal inhomogeneity.  Figure 22. As in Figure 18, with the same inhomogeneous N e field but as observed from two opposite directions that yield the horizontally-flipped sampling with respect to the TP (a,b). The solid and dashed lines in (c) denote the hTEC(h t ) profiles simulated for the two viewing directions. The N e profile at the TP location [red line in (c)] is same in both cases. The dashed lines in (a,b) denote the location of two state vectors that may be retrieved from hTEC(h t ) in the presence of horizonal inhomogeneity.
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.

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 CDAACcalibrated 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 frequencywavenumber spectrum of the EIA was obtained to unambiguously resolve both eastwardand 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 tõ 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.