Global GNSS-RO Electron Density in the Lower Ionosphere

:


Introduction
The D-(60-90 km) and E-region (90-150 km) ionosphere plays important roles in radio wave propagation because of their effective absorption of short-wave signals (e.g., high-frequency, HF radio) [1] and reflection of long-wave communication and navigation radios (e.g., very low frequency, or VLF) [2].It overlaps with the neutral atmosphere, known as the mesosphere, and is associated with large variability.The atmospheric temperature decreases with height in the mesosphere (or the D-region ionosphere), and the atmospheric thermal and density structures can be significantly perturbed by the atmospheric forcings from below (e.g., tides, planetary and gravity waves) [3].
The D-region is arguably one of the most critical interface layers between the lower atmosphere and the ionosphere.It is the layer where the electrical coupling takes place between the ionosphere and the phenomena such as sprites, elves and blue jets from tropospheric thunderstorms.The D/E-region conductivity is important to complete a so-called global electric circuit (GEC), an electrical coupling between the ionosphere and the planetary surface with a DC current driven by tropospheric thunderstorms [4,5].Because the D/E-layer conductivity can vary substantially between day and night, as well as between solar minimums and maximums over the 11-years cycle, it can short-circuit the F-region dynamo during the day [6,7,8] and modulate the GEC with the atmosphere throughout the day.However, it remains unclear if GEC has a significant feedback effect on thunderstorms, electrified clouds, or even on global cloud cover.If this feedback exists, the coupling between the ionosphere and the troposphere would act an important process to change Earth's radiation and climate system [9,10,11].
The D-region plasma is weakly ionized during daytime, mainly from photoionization of nitric oxide (NO) by the solar Lyman-α (Ly-a) radiation [12].When the photoionization is absent at night, the D-region electron density (Ne) decreases substantially but is maintained at a low level from complex sources [13].During the quiet-time, the daytime D-region Ne is balanced largely by positive molecular and cluster ions, where the upper D-layer mostly consists of molecular ions and the lower is dominated by cluster ions [14].A strong solar activity can bring hard X-rays (wavelength < 1 nm) to ionize mesospheric N2 and O2, but even in this situation the recombination rates are so high in the D-region that the environment tends to relax back to the neutral condition quickly once the ionization sources are gone.Hence, the number concentration in this altitude region is dominant by neutral molecules.
The D/E-region ionosphere can be often perturbed by irregular small-scale dynamic and electromagnetic variability and instability.In the E-region, for instance, sporadic E (Es) can occur as a transient disturbance to the background Ne profile at 90-120 km altitudes [15].The wind shear is thought to play a critical role in Es formation by converging the plasma density via the v×B drift in this region [16,17].Global distribution and seasonal variation of Es have been obtained and studied extensively with global navigation satellite system radio occultation (GNSS-RO) [18,19,20] and ground-based observations [21].These sporadic layers can even occur in the D-layer [22].At high latitudes, influences on the D/E-region come from magnetospheric disturbances through polar electron precipitation with enhanced ionization.This process is thought to play a role in increasing HOx and NOx productions through atmospheric chemical chain reactions involving electrons, positive and negative ions [23,24].The short-lived HOx tends to produce direct effects on the O3 loss that are localized in space (i.e., the mesosphere) and time (i.e., precipitation events) [25,26], whereas the long-lived NOx can be transported down to the stratosphere to produce indirect effects on the O3 loss [27,28].
The D/E-region Ne data are rare and limited, largely because of poor instrument sensitivity to low Ne concentration and incapability of remote sensing from space to profile a sharp Ne vertical gradient (i.e., four orders of magnitude over 30 km).Inconvenient and costly access to this altitude region makes in-situ measurements infrequent and sparse.Thus, the majority of lower-ionospheric Ne measurements were obtained from sounding rockets (D/E-region) and ionosondes (E/F-region), which constitute the basis of empirical models such as FIRI (Faraday-International Reference Ionosphere) [29] and IRI-2016 (International Reference Ionosphere-2016) [30].
The objective of this study is to develop a robust algorithm to retrieve global D/Eregion Ne, from the high-rate GNSS-RO data that have been increasing in spatiotemporal sampling from a large number of SmallSat/CubeSat constellations.An improved spatiotemporal coverage is critical for studying the lower ionosphere because of substantial and complex diurnal variations and influences from both the lower atmosphere and upper ionosphere.As in the algorithm developed previously for the E-region by Wu [31], the new retrieval method removes F-region contributions in the excess phase profile with a linear function fit to the data, before they are input for the D/E-region Ne retrieval.The GNSS-RO Ne retrievals tend to produce a lower Ne in the E-region compared to the IRI-2016 but in a general agreement with FIRI [32].Similar results are found in the GNSS-RO Ne retrieval from this study through a brief comparison with ionosonde observations at the Hermanus site.The paper is organized to describe the improved algorithm of Ne retrieval with examples in section 2, followed by the retrieval results and climatology of the GNSS-RO D/E-region Ne in section 3. Brief evaluations of the GNSS-RO Ne with ionosonde observations at the Hermanus site and a novel retrieval technique for total electron content (TEC) are discussed in section 4, before it is concluded in section 5.

Data and Methods
The algorithm developed in this study is an extension of the earlier attempt to retrieve the Ne profile in the lower ionosphere from the GNSS-RO excess phase ( !" ) profile [31].The RO excess phase is the difference between the measured and the free-space straight-line distance of GNSS (transmitter) and LEO (receiver) satellites.In the so-called "bottom-up" technique, a linear function is fitted to the  !" profile over a narrow ht range below the D-region to remove the F-region effects, before the data are input for the D/Eregion Ne retrieval.The F-region effects need to be removed first as cleanly as possible on a profile-by-profile basis, because a small  !" residual from the F-region variability can induce large errors in the D/E-region retrieval.The bottom-up method, self-sufficient in removing F-region effects, is specially designed for retrieving low Ne concentrations in the lower ionosphere [Appendix A].It has several important differences from the conventional onion-peeling or top-down approach.
First, the top-down method assumes spherical homogeneity for all Ne layers along the RO path and requires absolute calibration of the excess phase ( !" ) measurements, whereas the bottom-up method assumes only homogeneity in the D/E-region.The homogeneity assumption on the F-region Ne is not required, because it can become problematic for retrieving the Ne below it, as a small error in the F-region Ne retrieval can induce a large effect on the D/E-region Ne.This error propagation is caused by a long tail of the Abel weighting function, which does not vanish even at the RO tangent height (ht) far below the F-region [31].As a result, small residual errors of the F-region Ne retrieval can propagate down to the D/E-region Ne in the data inversion process, and generate a large error there [33,34].A manifestation of this error propagation is clearly seen in the Ne profiles retrieved from the onion-peeling approach, which often exhibit large oscillatory, sometimes negative values in the lower ionosphere [35,36,37].On the other hand, the bottomup method shows little contribution from the F-region in its weighting function in the D/E-region, and therefore its Ne retrieval is less impacted by the F-region error.The differences between the Abel and bottom-up weighting functions are detailed in section 2.3.
Second, the bottom-up method is resilient to complex F-region variability with capability to cope with most of the F-region effects (e.g., inhomogeneity, structural variations) on a profile-by-profile basis.It is a self-sufficient approach, without relying on any auxiliary data.In other words, the algorithm contains a stand-alone data process using the same  !" measurement profile for both E-region Ne retrieval and F-region effect removal.Previous studies found that the E-region Ne retrieval using an auxiliary model or constrained to an assimilated 3D Ne distribution had limited improvements for mitigating unrealistic F-region Ne features [38] and gradient effects on the derived E-region Ne [39].These model-aided approaches are prone to ionospheric variability and tend to bias towards the quiet-time condition, not working reliably for a disturbed ionosphere.
Third, the bottom-up algorithm does not require absolute  !" calibration that is needed to derive a line-of-sight (LOS) or horizontal TEC (hTEC) profile prior to the Ne inversion.It uses a differential  !" profile, or the hTEC difference (Δℎ), with respect to a linear function fitted to the  !" or hTEC data below the D/E-region.Thus, the linear fit serves a reference of each individual measurement profile, and the derived Δℎ is used for the Ne retrieval.The algorithm assumes that the F-region contributions to the hTEC in the D/E-region vary linearly over a narrow ht range, and this linear variation can be estimated and removed using the hTEC measurements below the D-region.As a result, the new algorithm is perhaps more sensitive to residual ionospheric errors (RIEs) and phase scintillations that are uncorrectable with the dual-frequency GNSS RO receivers.
But these residual errors are relatively smaller than the systematic error from the spherical homogeneity assumption about the F-region.
Fourth, the bottom-up retrieval should not be used to replace the Abel inversion for the whole ionospheric Ne.Because the linear function used to remove F-region effects is a fit to the lower-ionospheric data, its uncertainty tends to increase with height and causes a large error in the Ne retrieval in the upper ionosphere.Thus, the Abel inversion from hTEC, which requires absolute calibration of the  !" profile, still works better for the whole-ionospheric Ne retrieval.

RO Excess Phase (𝜙 !" ) Measurements
As a key variable for the Ne retrieval, hTEC(ht) is the integral of Ne along the RO path at each ht.It is related to excess phase  !" profile by Errors in hTEC can come from higher-order dependence (e.g.,  $% ) and multi-path propagation.In the E-region ht, the hTEC error contains contributions from both the F-and Eregion, of which the former is much larger because its residual effects cannot be fully removed.Thus, any small error from the F-region Ne retrieval would cause a large error in the E-region Ne retrieval.The hTEC profile can be derived from the L1 ( # =1.57542 GHz) and L2 ( $ =1.22760 GHz) excess phase ( !"%# and  !"%$ in meters) by ℎ = ( !"*( −  !"*+ ) In the previous algorithm, Wu [31] used only the  !"%# measurements to derive hTEC.Therefore, the linear fit for F-region effect removal had to be carried out at ht = 50-80 km, to avoid the atmospheric bending ( &'( ) at ht <60km.Because the  &'( contribution is same in the L1 ( !"%# ) and L2 ( !"%$ ) profiles, Eq.2 uses both L1 and L2 measurements to cancel the  &'( contribution, such that the linear fit to be carried out from lower ht.In this study we fit the hTEC data at ht =30-60 km to a linear function, and extrapolate it up and subtract it from hTEC to yield Δℎ for the D/E-region Ne retrieval.Hence, the use of hTEC data at ht =30-60 km allows the Ne retrieval extended to the D-region. To better understand and evaluate error sources and their contributions to hTEC, we broke down  !" into the ionospheric and atmospheric components.First, the neutral atmospheric component ( )*+*_-.!! ) in the  !" profile is defined by Here  )*+*_-.!! includes the atmospheric bending  &'( and the residual ionospheric effect (RIE), or  /01 , that cannot be accounted by the assumption from Eq.1 Because  /01 is not from the neutral atmosphere and is not negligible at ht >80km, it should be considered as a measurement error of hTEC in Eq.2.With the dual-frequency approach, one would hope to correct the  2$ dependent ionospheric effects in radio propagation.However, multi-path propagation and higher-order frequency dependence may induce the RIEs at these altitudes [40,41].In Fig. 1a, the  &'( contribution is obtained from fitting an exponential function to the  )*+*_-.!! data at ht =20-50 km, and extrapolated to the heights above.
Furthermore, it is important to understand the contributions to  !"_)*+* induced by ionospheric bending ( 3!+4_)*+* ) and phase advance ( 5_)*+* ) effects as radio wave propagates in the ionosphere.The contributions from these processes are significant in the  !"_)*+* profile shown in Fig. 1b, and these components have the opposite sign to each other.
This is because bending delay (due to a longer pathlength) and phase advance (due to faster phase speed) in the ionosphere produce the opposite effects in the measured excess phase  !" .These components are defined by where RO and SL in the integral represent the RO and straight-line (SL) path of radio propagation, and n is refractive index of the ionosphere.Although the RO  !"_)*+* is dominated by − 5_)*+* , ionospheric bending effects are not negligible because of the long pathlength in RO wave propagation.On the other hand, phase advance results from a faster phase speed of wave propagation in the ionospheric plasma, which will shorten  !"_)*+* .The two competing effects both manifest themselves in the GNSS-RO  !" profile shown in Fig. 1b.
The number of Ne retrievals may decrease with the dual-frequency approach, to some extent, compared to the  !"%# -only approach.Data quality control is obviously more demanding from the joint use of  !"%# and  !"%$ profiles, because large spikes in either measurement can cause a bad Ne retrieval.Because GNSS-RO  !" measurements generally have better signal-to-noise ratio (SNR) for  !"%# than for  !"%$ from GPS (Global Positioning System), the required  &'( calculation may subject to the additional error induced from the  !"%$ profile, which can lead to an unsuccessful Ne retrieval.In addition, the relative signal strengths between L1 and L2 transmitters may vary with the GNSS systems.For example, GLONASS (Global Navigation Satellite System) satellites sometimes transmit stronger signals at the L2 frequency than L1.
The linear fitting for the F-region effect removal is similar between this study and one developed by Wu [31], except for slightly different height ranges.In this study a linear function is fitted to the hTEC profile at ht =30-60 km on a profile-by-profile basis for ht =30-60 km (8) The derived linear fit is extrapolated to the heights above at ht > 60km, as an estimated Fregion contributions.By subtracting it from ℎ, we can obtain the D/E-region Δℎ profile (in TECu) where 1 TECu = 10 16 m -3 .Note that the conversion in Eq.9 neglects higher-order frequency dependence and ionospheric bending effects.
In essence, the bottom-up approach uses the difference Δℎ(ℎ ' ) with respect to the empirical reference (i.e., a linear function derived from ht=30-60 km), instead of absolute ℎ(ℎ ' ) values, for the D/E-region Ne retrieval.As shown in the next section, the Δℎ(ℎ ' ) at the D/E-region ht have much smaller contributions from the F-region Ne compared to the Abel weighting functions, which allows to yield a more reliable inversion for the D/E-region Ne.

Low-Rate (1 Hz) POD and High-Rate (50-100 Hz) RO Data
Both precise orbit determination (POD) and RO antennas can track GNSS signals at negative elevation angles or, in other words, at ℎ ' below the LEO receiver altitude.Generally speaking, the POD antennas have a wider field-of-view (FOV), lower SNR and a coarser sampling resolution, compared to the RO antennas.Depending on the POD receiver design, receiver satellite altitude and data acquisition schedule, the POD 1-Hz occultations usually have a fewer number of profiles (10%-30% in COSMIC-1 podTec) reaching the D/E-region.There are slightly more POD occultation profiles reaching the D/Eregion when the satellite is in a lower orbital altitude.Thus, the higher RO cutoff tangent height and noisier hTEC measurements make the POD data less attractive for the D/Eregion study.
In this study we focus mainly on utilizing the high-rate (50-100 Hz) RO data, of which the sampling and spatiotemporal coverage are much higher than the POD data.The dense spatiotemporal sampling is critically needed to characterize the lower ionosphere variability.In addition, the high-rate data have advantages to detect and correct  !" measurement errors such as clock time and cycle slips in the presence of ionospheric scintillations [42].Smaller  !" measurement errors help to retrieve a low Ne concentration in the lower ionosphere.
In recognition of the high-rate RO measurements for ionospheric observations, the MetOp-A operation conducted a short experiment with the high-rate RO acquisition up to ht ~290 km during 2020d161-2020d254 (June 9-September 10, 2020).The success of the experiment allows the extended acquisition mode to continue from 2021d201 to the decommission in December 2021.As shown later in section 3, the high-rate provides a significant benefit to the Ne retrievals in the E-to-F region transition and the retrieval of TEC.In addition to the extended high-rate profile from MetOp-A, the Spire GNSS-RO constellation, which NOAA uses to fill the polar latitudes outside the COSMIC-2 coverage, has been routinely acquiring the high-rate RO profile up to ht~175 km.Among the earlier missions, the C/NOFS (Communications/Navigation Outage Forecasting System) had a frequent high-rate RO data reaching up to ht~165 km from a low-inclination (13°) elliptical orbit,.
The high-rate data in the D/E/F-region also provide a valuable diagnostic for better understanding and quantifying the RIEs of  !" measurements in the ionosphere.The high-rate data help to reduce some of the scintillation-induced  !" measurement noise, but the high-order frequency dependence and multi-path effects remain as a source of the Δℎ errors in Eq.9 [31].As shown in Fig. 2, Es-induced phase scintillation can become a large random error source of Δℎ at ht below the Es layer.A RIE from the F-region can sometimes induce a bias.If the bias is confined in a narrow height region, their impact on the Ne retrieval is generally small.But if the bias is present at all heights, the induced bias can be significant in the retrieved Ne.
RIEs can also have a non-negligible impact on the neutral atmospheric refractivity retrieval.Bending angle (), which is used to derive the atmospheric refractivity, is essentially a vertical derivative of the  )*+*_-.!! profile, i.e.,  ≅ −  )*+*_-.!! ℎ ' ⁄ .Therefore, a bias from RIEs would not be a problem for the  retrieval, but a slope or vertical gradient of RIEs would.The RIE-induced  error has been a major issue in the climate-quality data analysis of the upper stratosphere [43,44].As shown in the second case in Fig. 2, a vertical gradient (~1.25´10 -7 radians) and bias (~0.03 m) are evident in the RIE profile.For this study, therefore, we simply allow a 0.3 m uncertainty for these RIEs in the D/E-region retrieval algorithm.Full characterization and correction for RIEs is an ongoing research and will be reported in a separate study.m/km, or ~1.25´10 -7 in radian, along with a bias of 0.03 m.

Inversion from Δℎ𝑇𝐸𝐶
The inversion algorithm for the Ne retrieval in this study uses the optimal estimation inversion method [45].The measurement vector y is the Δℎ derived from Eqs.3-8 as a function of ht, i.e., y ≡ {Δℎ(ℎ When the measurement vector is changed from ℎ to Δℎ, the weighting functions of the bottom-up algorithm is also transformed from the Abel functions to a new set of functions as shown in red in Fig. 3d.Unlike the Abel function, the new weighting functions for Δℎ have a very small contribution from higher altitudes, which provide the advantage of minimizing F-region impacts. The linearized forward model for the measurement vector () is related to the state vector () through a weighting function K, or the Jacobians as in  = /, as follows where the linearization is based on the state vector   (Ne in m -3 ) and the measurement vector   (Δℎ in TECu).A single Ne profile is used for the   in Eq.10, which is a merged profile from the annual mean of IRI-2016 from 2008 for the F-region [30] and the mean FIRI E-region profile for the solar minimum condition [29].We assume that the weighting function K is invariant with the coefficient  # (i.e., the slope of the linear function in Eq.8).The linearization model   ,   , and K are depicted in Fig. 3.The measurement uncertainty of y, or   , is 0.3 m for the  !"%# , or 1.9 TECu for Δℎ, to account for multipath, RIE, scintillation, neglection of ionospheric bending, and other measurement errors.The multipath error from reflections off the spacecraft structure is estimated to <0.2 m for COSMIC-2 [46] and <1 m for COSMIC-1 [38].The RIE error, as shown in Figure 1, is on the order of ~0.1 m but can be larger in a strong scintillation condition.The  inversion in the optimal estimation method [45] can be expressed by where  K is an optimal solution to x in Eq.10 and a is the a priori of x, which is set to be same as   in this study.  = 0.02 • (5 • 10 > +   ) •  and   =    •  are covariance matrices for a and y respectively.Note that   is a height-dependent parameter, which is empirically determined for stable Ne inversion over a large dynamic range and variability in the D/E-region altitudes.The detailed data flow diagram of the new algorithm and changes can be found in Appendix A.
For computation efficiency of the inversion from Δℎ to Ne, we employ a variable height range for the state vector such that  (&" ≤ ℎ (&" +10 km.Because Δℎ(ℎ ' ) does not have much information on Ne at altitudes above ℎ (&" , it would waste resources to carry out the Ne retrieval at extra levels above ℎ (&" .This adaptive height range in the retrieval helps to reduce the number of Ne levels in the state vector, and therefore save CPU time for the inversion calculation in Eq.11.

D/E-Region
The bottom-up algorithm is developed to enable and optimize the D/E-region Ne retrieval in the situations where the F-region hTEC measurements are unavailable for correcting their effects on the lower-ionospheric retrieval.It is difficult to make this correction with the onion-peeling method when the entire hTEC profile is not available.An example from Eqs (9-10) is shown in Fig. 4   The D/E-region Ne. retrievals show a good sensitivity to the diurnal variations [Fig.5], since a constant linearization profile is used in Eq.11 without any a priori information on solar zenith angle (c).As expected for the low Ne.concentration, the retrieval result is noisier at lower altitudes.More scatters are evident in the early evening hours, especially in the E-region, suggesting the potential impacts from ionospheric scintillations [47].The local time coverage from the COSMIC-2 is fairly uniform on a single day.By the time of this day (i.e., January 1, 2021), the 6-satellite constellation had spread out in local time sampling from their initial launch orbits on June 25, 2019.The low inclination (24°) orbit also helps to cover the tropical diurnal cycle in a relatively shorter period compared to those higher-inclination orbiters.
The diurnal variation of the retrieved D/E-region Ne is generally consistent with FIRI, although it exhibits a large deviation at 60 km [Fig.5].Unlike FIRI, which is symmetric about the local noon, the RO Ne retrievals are skewed slightly toward the afternoon hours.It is clear that the morning-hour onset of the daytime RO Ne lags the FIRI curve, while the retreat from the RO in the evening hours appears lagging behind the noon-symmetric FIRI model.The daytime RO Ne retrievals are slightly higher than the FIRI at 110 km.Given that F10.7 is at 77.7 solar flux unit (sfu) on Jan 1, 2021, the RO Ne retrievals are generally higher, compared to the FIRI prediction for F10.7=75 sfu.

Lower F-Region
The bottom-up algorithm can be extended to retrieve the Ne in the lower F-region if the RO hTEC profile reaches a higher altitude.However, most of the high-rate RO measurements do not go above ht~120 km in the normal operation.On the other hand, MetOp-A has acquired some high-rate RO data with a top up to ~290 km, allowing to test the new algorithm for the lower F-region retrieval.Fig. 7 shows the extended Ne retrievals from the ascending (nighttime) and descending (daytime) orbits in October 2021, revealing the Eto-F-region transition across the meridional plane at two local times (8-11h and 20-23h).
Because the E-and F-region dynamos are electrically connected, this global coverage in the lower ionosphere with a high spatiotemporal sampling is critically needed for studying and understanding their coupling processes.
In addition to MetOp-A, the Spire constellation also produces ROs with a top routinely up to ~175 km, higher than most of the GNSS-RO operations.The Spire's high-rate data occasionally go up to an even higher altitude for large scintillation cases.As shown in Fig. 7, these high profiles provide only sparse sampling at altitudes above 175 km, not enough to produce a reliable climatology.Below that altitude, the Spire sampling seems to be uniform.Different from MetOp-A, the Spire constellation covers all local times.In other words, the Spire zonal mean is truly daytime and nighttime from all local times, whereas the MetOp-A zonal means are from two sun-synchronous local times [Appendix B].Nevertheless, the Spire daytime and nighttime climatologies below 175 km exhibit the distributions similar to those from MetOp-A.

TEC and Vertical Gradient of ℎ𝑇𝐸𝐶(ℎ ' )
The vertical gradient of ℎ(ℎ ' ) profiles at the lower ionospheric ht contains quantitative information on TEC and can be used to retrieve TEC without requiring absolute calibration of ℎ(ℎ ' ).In this section we further analyzed the TEC information using the synthetic data from the January 2008 IRI-2016 Ne profiles for different geolocation and solar local time conditions.We compute TEC and ℎ(ℎ ' ) from these Ne profiles and applied the same retrieval algorithm as in Eq.8 to the modeled ℎ(ℎ ' ) and coefficient  # .As illustrated in Fig. 8a, the coefficient  # from the synthetic data and TEC are highly correlated.There is a clear difference between the daytime and nighttime  # -TEC relationships, showing a smaller slope from the daytime profiles.On the other hand, TEC appears to correlate well the hTEC difference or DhTEC between ht =200 and ht =60 km.These sensitivities have an important application for the TEC remote sensing, because the hTEC differencing (DhTEC) approach does not require absolute calibration about hTEC measurements.It can be simply achieved by using the Δℎ with respect to the value at a reference ht.In a broader sense, it is a similar technique to the measurement vector as formulated in Eqs.7-8.
The model simulation shows that TEC can be retrieved from the lower-ionospheric Δℎ measurements with the reference at ht =60 km.The technique appears to work best for the Δℎ measurements with the top reaching ht =200 km, the altitude slightly below the F2 peak.Assuming a measurement uncertainty of 1 TECu, we test this idea with the synthetic IRI-2016 data and derive a set of regression coefficients for the TEC retrieval.The synthetic Δℎ and TEC data represent all variabilities from 2008, and the linear regression is used to develop the retrieval coefficients for Δℎ.Fig. 8c and Fig. 8d show two TEC retrieval methods from different 5-level Δℎ data.The regressed coefficients for the five-level Δℎ data are given by Eq.12 and Eq.13, for the RO top reaching 220 km and 180 km respectively, where TEC' is the TEC retrieval from the 5-level Δℎ measurements.In summary, the TEC retrieval method demonstrated in Fig. 8c is practically useful for the high-rate RO data with the top above 220 km.On the other hand, the method in Fig. 8d is applicable to the Spire near-real-time (NRT) data when the top of RO profiles reaches ht ~175 km.For accurate TEC retrievals, these high-ht RO observations (close to the F2 peak) are more desirable.Thus, it is recommended in future GNSS-RO operations to extend the RO acquisition up to at least ht =220 km.

GNSS-RO Data
Global Navigation Satellite System (GNSS) provides continuous broadcast of positioning and time data from space to the ground and other satellite receivers that can be used to determine their location.There are six major GNSS constellations either in operation or development: U.S. Global Positioning System (GPS), Russian Global Navigation Satellite System (GLONASS), European Navigation Satellite System Galileo (Galileo), Japanese Quasi Zenith Satellite System (QZSS), Chinese BeiDou System (BDS), and Indian Regional Navigation Satellite System (IRNSS).In addition, civil Satellite-based Augmentation System (SBAS) for aviation safety uses geostationary (GEO) satellites to improve the accuracy and reliability of GNSS information by broadcasting the augmentation data with signal measurement error corrections.Each of these GNSS constellations is assigned with satellite system code: G: GPS, R: GLONASS, E: Galileo, J: QZSS, C: BDS, I: IRNSS, and S: SBAS, which is used in naming every GNSS-LEO link for the RO data.
The GNSS-RO data used in this study are summarized in Appendix B (Table B1), which lists different RO missions with information such as operation period, LEO (low Earth orbit) satellite initial and final altitude, orbital precession, latitude coverage, RO top height, and number of daily RO profiles.The majority of GNSS-RO data is obtained from the COSMIC Data Analysis and Archive Center (CDAAC) distribution, which processed the data to Level-2.The Level-1b data are processed with different versions of software, namely, NRT, postprocessing (postProc), and reprocessing (e.g., with subfix '2016').The satellites that contributed to the global GNSS-RO acquisition include CHAMP (Challenging Minisatellite Payload), SAC-C (Satellite de Aplicaciones Cientifico-C), COSMIC-1/2 (Constellation Observing System for Meteorology, Ionosphere, and Climate 1 or 2), METOP-A/B/C (Meteorological Operational Polar Satellite A/B/C), GRACE (Gravity Recovery and Climate Experiment), C/NOFS, TSX (TerraSAR-X), TDX (TanDEM-X), Korean KOMPSAT-5, Spanish PAZ, and Chinese Fengyun-3 (FY3).

Diurnal Variations
The D/E-region Ne is a strong function of c as the lower ionosphere is largely influenced by the Ly-a ionization of mesospheric species such as NO.To characterize the Ne morphology, we averaged all the COSMIC-1 RO retrievals from 2006-2020 to generate a monthly climatology for each 2-hourly local time, 4° latitude and 2-km altitude bin.As shown in Fig. 9 for four selected months and four selected altitudes, the new Ne data reveal several interesting features and variations in the lower ionosphere.A complete set of monthly climatology can be found in Fig. S1.
Although the D/E-region Ne shows a clear c-dependence, it exhibits strong latitudinal variations for the same c.For example, almost at all altitudes there is an enhanced daytime peak at the summer midlatitude during the solstice months (e.g., June and December).This peak occurs at a late morning hour slightly before noon.For the months near equinox (e.g.March and September), a dip at the magnetic equator is evident in the daytime Ne, especially at 100 and 120 km.This low Ne band, known as the "fountain effect" in the Fregion, is likely an imprint of an upward plasma draft induced by the equatorial electrojet (EEJ).As expected for an eastward electric field build-up from EEJ, the E´B force could create an upward plasma drift [48].Because of the strong vertical gradient in the D/E-region Ne, this vertical drift would manifest itself as a low Ne value compared to those at nearby latitudes without the drift.
Another noticeable feature is the asymmetry of Ne distribution about the noon.The D/E-region Ne climatology clearly shows that Ne arises slightly late in the morning with a lagged decay into the evening.The D-region Ne asymmetry about noon was also reported in ground-based observations [49].In all monthly climatologies, the extended Ne enhancement from evening to midnight is observed.During the months near the equinox, such a prolonged Ne enhancement in the tropics is punched by a short decrease around LST=19h.

D-Region H' and 𝛽 Parameters
The D-region (60-90 km) ionosphere is poorly characterized due to lack of Ne observations.In-situ (e.g.sounding rockets) [50,51], remote sensing from incoherence scattering radar (ISR) [52] and radio wave propagation [53] are the common methods for measuring the D-region Ne [54], but their locations are sparse and the operation cost can be high.In addition, low electron concentrations in the D-region requires high power, large-size antenna radars to obtain useful Ne measurements, which makes such systems unaffordable for wide deployment.On the other hand, high-frequency radio waves, although useful for the Ne sounding in the E-region and above, tend to pass through the D-region.At very low frequency (VLF, 3-30 kHz) and extremely low frequency (ELF, 0.3-3 kHz), radio waves are completely reflected by the D-layer and the ground surface, of which the two conducting boundaries act an Earth-ionosphere waveguide (EIWG) to allow a low-loss long-range radio propagation.As a result, the VLF/ELF technique has been increasingly used for remote sensing of the D-region properties including Ne [55,56,57].
To characterize the D-region Ne, the VLF/ELF observation community often uses a two-parameter D-region Ne model.As described by Wait and Spies [1964], the formula has characteristic reference height H′ (km) and exponential growth rate of Ne with height or  (km −1 ), i.e., The model parameters from fitting to VLF/ELF observations yield a typical H' value of 65-90 km and  of 0.3-0.5 km -1 with a lower H′ value and slightly higher  value during daytime [58][59][60][61][62][63][64].Recently, an update was proposed to characterize a height split in the Dregion Ne profiles with a four-parameter scheme, which appears to represent the daytime profiles better [65].
We used the same 2-parameter model (Eq.10) to fit the D-region Ne retrieved from COSMIC-1 [Fig.10].The H′ parameter derived from COSMIC-1 are slightly lower than those derived from the VLF observations [59,60], showing a daytime variation of 70-73 km, compared to 71-77 km from VLF.The COSMIC-1  parameters exhibit a little variation around 0.31 km -1 during day, compared to 0.25-0.4km -1 reported from the VLF observations.For nighttime, the COSMIC-1 H′ and  values are also lower, with H′ varying between 71 and 76 km compared to 82-86 km from VLF, and  varying between 0.29 and 0.30 km -1 compared to 0.4-0.6 km -1 from VLF [61,62].

Magnetic-Field Dependence
In addition to the diurnal variation, the D/E-region Ne can be significantly modulated by the magnetic field.To illustrate the magnetic-field effect, we averaged the COSMIC-1 Ne data from all LSTs and produced a monthly map for each altitude as a function of geographical longitude and latitude.Because the Ne at low and midlatitudes is dominated by the daytime climatology, the all-time averages look similar to the daytime maps except for high latitudes.Fig. 11 shows the 70-km maps , but the Ne distributions at the E-region altitudes are similar.A complete monthly climatology of all-LST maps can be found in Fig. S2.The features observed below are generally applicable to the COSMIC-1 data at 70-120 km.
As shown in Fig. 11, magnetic-field modulation of the D/E-region Ne is present in all months where the longitudinal variations follow the field line at most latitudes.There is an important distinction between the high-latitude and low-and-midlatitude variations on what magnetic fields are in control.At high latitudes, as expected for the domination of auroral electron precipitation, the Ne distribution is correlated better with the magnetic field in the magnetosphere (e.g., at one Earth radius Re).As seen Fig. 11 (left panels), the correlation looks worse if it is compared to the magnetic pitch angles from the field at 100 km.But the distribution of auroral electron precipitation correlates better with the pitch angle derived from the magnetic field in the magnetosphere.In addition to the direct contribution from electron precipitation, the winter polar vortex is capable of transporting the aurora-produced NO down to the mesosphere and further enhance the nighttime D-region Ne [66].The fact that the polar Ne enhancements are closely associated with the magnetic pitch angle implies that such vertical transport from the polar vortex only plays a secondary role in the polar winter anomaly.As shown in the RO observations, the auroral electron precipitation does not produce a very strong polar Ne distribution as simulated under different hypothesized scenarios for the D-region ion chemistry [67].
It is worth noting that the summer-midlatitude Ne distributions in Fig. 11 (June and December) bear a striking resemblance to the global Es morphology reported in the earlier studies with GNSS-RO [18,19,20].As stressed by Wu [31], the new D/E-region Ne data will allow a more comprehensive study of the Es formation and its connection to the background E-region Ne concentration.The similar distribution between Es and Ne suggests that a common process is modulating their longitudinal variations in the summertime midlatitudes.
It remains challenging to fully understand what processes drive large longitudinal variations along the same field pitch angle [Fig.11].Planetary waves, such as the migrating and non-migrating tidal waves, can perturb the D/E-region Ne distribution by modulating the NO density and temperature in the upper mesosphere.In addition, some observations show that the lightning from tropospheric thunderstorms can induce a direct decrease in D-region Ne [68] or sharpening in the Ne vertical gradient [69,70].The lightning-induced Ne reduction may explain the relatively lower Ne values (along the same pitch angle) in June over the American and Indian Monsoon regions where deep convection and lightning activity are likely higher than elsewhere.It remains challenging with the GNSS-RO Ne observations to separate the direct lightning effect from others, as the D/E-region Ne is strongly coupled with other processes (e.g., solar radiation, planetary waves).

Seasonal Variations
To characterize the Ne seasonal variations, we divided the monthly data into daytime, post-sunset (LST=19h), and nighttime bins, and plotted these binned data as a function of time and magnetic latitude.Fig. 12 shows the time series of Ne at 60, 80, 100, and 120 km.To the first order, the seasonal variation of daytime Ne at these altitudes is driven by solar radiation as the c maximum moves around the equator with season.The midlatitude maximums near the solstice months are clearly seen in the time series, whereas the seasonal minimums at the midlatitudes appear to be quite different.The minimum in the SH winter is much lower than that in the northern winter.The seasonal asymmetry in Ne is extended to the nighttime time series where the c effect can still be seen.The nighttime time series provides a better view of the seasonal variation of polar Ne.Similar to the post-sunset result, the polar Ne has a larger concentration in the summer month, as expected for a stronger field-aligned current from the sunlit hemisphere [71,72].
A strong semiannual variation is observed in the nighttime D-region Ne with peaks near the equinox months.Several studies have suggested that the semiannual variation in the D-region could be induced by the NO concentration modulated by dynamic transport processes in the mesosphere and lower thermosphere [73,74].In this case the modulated NO concentration and the enhanced nighttime ionization in the D-region are correlated with a peak around the equinox.Consistent with the semiannual variation reported in these earlier studies, the COSMIC-1 observations seem to have two peaks in the tropics, but with a slightly higher concentration in October than that in March.

E-to-F-Region Transition
The onion-peeling approach often produce negative or large oscillatory Ne values in the lower ionosphere [38].These large errors make the onion-peeling retrieval unsuitable for studying the E to F-region transition.The new Ne data, such as those with a higher RO top (e.g., Spire), can be used to describe the plasma transition from the E to F region with continuous height coverage and high spatiotemporal sampling.Fig. 13 presents an example of monthly Ne statistics from such observations where a high RO-top (e.g., Spire) RO data are used.The retrieval covers an altitude range between 60 and 170 km and the data are averaged to 2-hourly bins in local time.More displays of Spire monthly climatology can be found in Fig. S3 for the zonal mean and in Fig. S4 for diurnal variations at the magnetic equator.One of the interesting features in the 2-hourly zonal mean [Fig.13] is the shift of the daytime Ne maximum from the SH E-region to the NH lower F-region.This shift implies a height-varying c dependence of the Ne morphology in the lower F-region.Again, the equatorial dip in Ne is obvious in the late afternoon and early evening hours (i.e., LST=15h, 17h and 19h), continuing through all altitudes in the E-to lower F-region.The polar features from auroral electron precipitation are more pronounced in the nighttime zonal mean in Fig. 13.It is also known as the polar winter anomaly where ionospheric absorption of radio waves becomes extraordinarily strong [75].Relatively speaking, they have a nearly uniform concentration in the lower F-region, but with a greater extension to the E-region compared to the low-to-midlatitude Ne.The penetration to the D/E-region seems to be deeper around the midnight hours in the NH.Compared to the D/E-region, the nighttime Ne distribution in the lower F-region is less symmetric about the equator.In the early evening hours, the lower F-region appears to have a higher Ne in the winter subtropics.But the distribution peak shifts to the summer subtropics during the predawn hours.
The results presented in Fig. 13 are based on the sampling from a subset of the Spire constellation data provided to the NOAA's CWDP program, which have ~4000 daily RO profiles.Processed by the same COSMIC software to Level-1 and Level-2, these data only become available on CDAAC for the research community since September 2021.As the constellation continues to expand, it is anticipated that more RO profiles will yield a further improved diurnal coverage on a monthly basis, perhaps to enable one-hour intervals in local time and a spatial resolution better than the 8°´4° longitude-latitude grid used in this study.One of the advantages with dense spatial sampling is to allow a tomographic retrieval of Ne, which can help to mitigate the horizontal inhomogeneity problem associated with the poor horizontal resolution of the GNSS-RO technique.Characterizing the long-term D/E-region Ne variations requires stable satellite operation and uniform sampling/quality in GNSS-RO observations.Because of the strong diurnal variability in the D/E-region Ne, a non-uniform or variable sampling of local times would make the observed Ne variation difficult to interpret and prone to the sampling error from aliasing to the diurnal variation.Although COSMIC-1 had a full diurnal coverage during the early period of mission, its sampling degraded over time as several of its operating satellites stopped working.As a result, the COSMIC-1 data are not ideal for studying the solar-cycle variation.On the other hand, the MetOp observations have been stable with the longest record at two sun-synchronous local times around 8-10AM and 8-10PM.Since the routine GNSS-RO acquisition from MetOp goes up to ~90 km, this long record is only available for the D-region Ne.To fill in the early period (2007)(2008) of the MetOp record, we added the SACC data when they had the sampling at the similar equator-crossing time [see Appendix B].

Solar-Cycle Variations
The daytime and nighttime solar-cycle variations at 70 km from the MetOp/SACC observations are highlighted in Fig. 14 as a function of magnetic latitude.As expected for the sunlight-driven ionization in the D-region, the D-region Ne correlates well with the Ly-a flux in terms of the solar-cycle variation.In addition, there is a strong semiannual modulation in the magnetic tropics from both daytime and nighttime time series.Because of the presence of the dominating solar-cycle and semiannual variations, it is not straightforward to discern a quasi-biennial oscillation (QBO) modulation in the D-region Ne, which was reported in some of the early studies [76,77].

Comparisons with Ionosonde Data at Hermanus
As a preliminary validation effort, the GNSS-RO E-region Ne profiles are compared against ionosonde Ne.The Hermanus (HE13N), South Africa (34.4°S, 19.2°E) ionosonde site was selected for the comparison because of the low fmin (minimum return frequency observed by the ionosondes) values over the period of interest: July 2008 -Aug 2010.A low fmin value in the ionogram measurement is vital to derive the E-region Ne profiles.To collocate the GNSS-RO observations, RO tangent points were required to be within 2° of the Hermanus site and within 30 minutes of an ionosonde measurement that has an fmin below 1.2 MHz.The ionosonde data were obtained from the Digital Ionogram Database (DIDBASE, accessed 17 January 2022, https://giro.uml.edu/didbase/), and all ionograms were hand-scaled using the SAO Explorer application for ionosondes (https://ulcar.uml.edu/SAO-X/SAO-X.html) to eliminate uncertainties caused by automatic scaling and ensure the quality of ionograms used in the comparison.The ionogram with sporadic-E or an uncertain E-region profile was excluded in the comparison.Over this twoyear period, a total of 60 ionograms were hand-scaled for the comparison against the RO profiles.
Two examples of the ionosonde-vs-RO Ne comparisons are displayed in Fig. 15.The vertical line corresponds to the minimum return frequency measured by the ionosonde, fmin.The ionograms show E-region returns from fmin up to foE, but no measurements are obtained for the E-to F-region valley.Therefore, this comparison is limited to the altitude range between the fmin altitude and hmE, which varies with the time of day, location, and ionospheric conditions.As shown in Fig. 15, the typical starting height for ionosonde measurements is 90 km while GNSS-RO is capable of measuring densities at D-region altitudes (see section 3).The NmE estimated from the RO and ionosonde techniques are in close agreement for these two examples, and the bottomside E-layer Ne gradients with respect to altitude (dNe/dz) are also in agreement.
While many of the profiles show agreement in dNe/dz and NmE, the profiles are shifted in altitude with RO profiles being slightly higher by 2 km on average above the ionosonde profiles.This relatively small altitude difference may be due to the "internal uncertainty of the starting height of the profile" for ionosonde measurements [78].The ionogram inversion assumptions and uncertainties were described in detail in several studies [79,80,81], revealing challenges in converting measured virtual heights to real heights for Ne.On the other hand, errors from ionospheric bending and derivation of the Earth ellipsoid model may also induce an uncertainty in the GNSS-RO HSL, but it is expected to be ≤ 1 km.
For this comparison, we show the results obtained for both unshifted (original) and shifted profiles.The shifted profiles are obtained by adjusting the ionosonde Ne altitude to minimize the altitude difference between profiles at fmin.GNSS-RO and ionosonde derived Ne at each height within the fmin altitude and hmE for the 60 profiles analyzed in this comparison are displayed in Fig. 16.Here, the ionosonde Ne profiles are interpolated at each GNSS-RO altitude for a direct comparison.While the R 2 and slope values for the linear regressions are nearly equal for the shifted and unshifted profiles, the y-intercept is closer to zero for the shifted profiles, closer to the 1:1 line that represents a perfect match between the two techniques.Overall, the densities show reasonable agreement with a larger variation for elevated densities.Histograms of the differences in density (Ne,GNSS-RO -Ne,ionosonde) as a function of altitude are displayed in Fig. 17.The altitude bias in the unshifted profiles results in an underestimation of the electron densities by the GNSS-RO technique for the altitude range of roughly 90-110 km.However, when the ionosonde profiles are shifted in altitude, this underestimation bias is mostly removed, and the density differences are slightly less than zero for most altitudes.The mean absolute error (MAE) for the unshifted profiles is 1.9x10 10 m -3 while the shifted profiles MAE is reduced to 1.4x10 10 m -3 .Finally, the averaged daytime Ne for both techniques are displayed in Fig. 18 for the shifted and unshifted profiles.Daytime profiles are constrained between 0800-1600 in LST, providing a total of 22 profiles for comparison.The altitude bias is obvious in the unshifted profiles and shifting the ionosonde profiles up in altitude helps to improve agreement.While the shifted ionosonde densities are larger than the RO derived densities on average, the mean Ne profiles agree within their estimated uncertainty at E-region altitudes.Encouraged by the preliminary comparison, we plan to carry out a larger further validation study to include more ionosonde sites from different latitudes.The new GNSS-RO Ne data at the altitudes near the ionosonde fmin may provide an alternative pathway for pinning the starting heights during ionogram inversion.

TEC Comparisons and Maps
The TEC retrieval methods described in Fig. 8 are applied to the real RO data and compared with the TEC measurements from ground-based GPS receivers.To validate the TEC retrieval methods, we use the MetOp-A data from June 9-18, 2020 when it acquired the high-rate RO data with the top up to ~290 km.To find collocated RO profiles with ground-based TEC observations, the RO TEC measurements are required within 1°´1° latitude-longitude and 5-min windows.Fig. 19 shows the RO TEC retrievals for the same four methods as illustrated in Fig. 8. Similar to the conclusion in Fig. 8, the method from Eq.12 works best, showing the lowest bias and standard deviation.The RO-TEC and ground-TEC are in generally good agreement along the 1:1 line.The standard deviations from the MetOp-A retrieval are higher than the simulation from the synthetic data, which may be caused by additional variability and differences in their observing times and spatial averaging where horizontal inhomogeneity can play a role.To develop a generalized method beyond Eqs.11-12, we obtained a set of regression coefficients for the TEC retrieval from the RO profiles with different top cutoffs, using the synthetic IRI-2016 data.This method uses all Δℎ measurements below the RO top and the retrieval coefficients for different RO top cutoffs are stored in a lookup table.For the RO Δℎ profiles with the top below 180 km, the TEC retrieval error tends to increase.We found that the RO profiles with a top above 120 km can still produce a scientifically useful TEC retrieval, but the retrieval error becomes too large if the top reaches only to 90 km such as the normal operation with MetOp-A/B/C.One of the advantages with the high-rate RO data is to provide a better spatiotemporal coverage for the TEC measurement.Spatiotemporal sampling is critical for characterizing the ionosphere with large diurnal and geographical variabilities.Because of the recent high number of high-rate RO profiles [see Appendix B], the 2-hourly sampling of global TEC maps can be obtained on a monthly basis from GNSS-RO observations.Fig. 20 demonstrates the latest capability of dense sampling from the combined Spire and COS-MIC-2 constellations plus a few additional GNSS-RO satellites in October 2021.More examples of monthly TEC maps can be found in Fig. S6.While the daily sampling from the GNSS-RO observations is still sparse, the monthly sampling from these RO observations can provide nearly global coverage over a 4°´8° latitude-longitude grid every 2 hours.Thanks to the COSMIC-2 constellation, there is little data gap in the tropical/subtropical latitudes (< 45°N/S).Also, a small number of gaps are found in extratropical and polar latitudes largely because of the coverage from the Spire constellation.FY3D helps to fill in a small sampling gap in the Spire constellation around noontime at mid-to-high latitudes  Compared to the recent large number of RO profiles, the sampling from the COS-MIC-1 constellation plus a few additions still had many regional data gaps in the 2-hourly maps for August 2008 when they were in a mission prime time [Fig.21].The unsampled regions are clearly more than in the case for October 2021.The regions with missing data appear to be randomly distributed, as the COSMIC-1 constellation had been fully deployed and spread in local time coverage.More 'apparent' data gaps in the polar region are partly due to the reduced equivalent area from the 4°´8° latitude-longitude grid.
A few salient features emerge from the 2-hourly TEC maps from October 2021 [Fig.20].First, there appears to be a trough in the NH mid-latitude that evolves with time.October 2021 was generally a quiet month, but this trough feature is absent in the background TEC model [82].Second, an auroral contribution to TEC is evident in the nighttime SH during UTC=11Z-17Z.The auroral contribution is generally too weak to observe in TEC except from the enhancement during substorms [83].Third, there is UTC-dependent symmetry in the subtropical TEC about the magnetic equator.Asymmetric distributions are more pronounced for the UTC time before 07Z.
Current global TEC nowcast and forecast rely primarily on the data analyses from ground-based measurements at the International GNSS Service for Geodynamics (IGS) stations.Because the ground stations are mostly over lands, the IGS-associated analysis centers, including GFZ (GeoForschungsZentrum), JPL (Jet Propulsion Laboratory), and USNO (U.S. Naval Observatory), use various interpolation-filtering methods to fill the data-sparse regions in producing a global TEC map [82,84,85].Thus, a large TEC uncertainty is expected over the data-sparse regions.Over oceans satellite altimeters may fill sampling voids with a better sensitivity to TEC measurements [86].Nevertheless, the dense spatiotemporal sampling is needed to provide globally uniform quality for the TEC measurement.As demonstrated in this study, the high-rate RO data can provide such sampling and contribute to the Global Ionosphere Map (GIM) analysis for improved accuracy over data-sparse regions.

Further Improvements for the Bottom-Up Inversion
As described section 2, in the Ne inversion the modified measurement vector is a differential hTEC profile relative to its linear variation determined from ht=20-50km.This change yields two advantages: 1) significantly reducing the impact of F-region residuals on the D/E-region Ne retrieval; and 2) making the D/E-region Ne retrieval less dependent on absolute error or accuracy of the hTEC measurement.Apart from these changes, the core inversion formulation (i.e., Eq.11) is same in both approaches.
The same principle is applied to the TEC retrieval using hTEC differences between the measurement at ht=60 km and those.at ht < 220 km [Section 2.3.4].This DhTEC-based TEC retrieval will allow a full use of the high-rate RO data for both neutral atmospheric and ionospheric applications, and the retrieved TECs agree reasonably well with the ground-based observations in the case where the RO profile top reaches 220 km.Because the number of RO profiles is generally greater and the RO data quality (i.e., SNR and hTEC) is better than POD measurements, the high-rate RO profiles with a higher top will attract more attention in the dual applications for atmospheric and ionospheric sciences.
The bottom-up inversion algorithm has several similarities to the Abel inversion in formulation.A further development with this algorithm will be to apply the inversion technique to the 1-Hz POD measurements, some of which have the RO profile reaching to ht = 100km or lower.Although these hTEC profiles may not have the same quality as with the high-rate RO data, they can serve as a surrogate to verify the D/E-region Ne and TEC retrievals from the bottom-up method.If both POD and RO data can be used independently for the TEC retrieval, it will increase the sampling and coverage by a significant amount.Previous Ne retrievals from the POD hTEC have been struggling in the lower ionosphere, showing large oscillations in the Ne profile [38].Preliminary results show that the sequential optimization inversion performs better than the onion-peeling approach for retrieving Ne retrieval without artificial oscillations in the lower ionosphere.
Although it is feasible to jointly retrieve the Ne profile from the 1-Hz POD and highrate RO hTEC measurements, the LEO-PRN pair in both observations needs to be consistent.For COSMIC-2 there are ~50% of the POD profiles that can be paired to the highrate RO profile.Out of these matches, approximate 50% are overlapped in tangent heights.As a result, there are only ~25% of the RO profiles that may be used for jointly retrieving Ne from the D-region to F-region.

Conclusions and Future Work
We have extended the bottom-up Ne retrieval to the D-region (60-90 km) using the ΔhTEC measurement derived from high-rate GNSS-RO data.The new algorithm can produce a Ne profile up to 180 km if the high-rate RO top reaches this height.The optimal estimation inversion has been updated to handle both measurement vectors hTEC with the Abel weighting function and ΔhTEC with the "bottom-up" weighting function.This flexibility allows the algorithm to process the podTEC data, if needed, for the Ne retrieval.In addition to the D/E-region Ne, a new TEC retrieval technique was developed using the high-rate RO data with a top above 120 km.For better characterizing the E-to-F Ne transition and the novel TEC retrieval, it is recommended to have the routine GNSS-RO acquisition up to ht =220 km.
The RO Ne retrieval was compared briefly with the ionosonde observations at the Hermanus, South Africa site.A generally good agreement was found in the E-region, with the RO Ne profile being lower in density and shifted upwards slightly in altitude.Encouraged by the preliminary comparison, a greater validation effort is planned to include more ionosonde sites from different latitudes.
The morphology of new GNSS-RO Ne data revealed many interesting features in the diurnal, seasonal, solar-cycle and magnetic-field-dependent variations.While the D/E-region Ne is a strong function of solar zenith angle (c), it exhibits strong latitudinal variations for the same c and its distribution is asymmetric about noon.The magnetic-field modulation of Ne varies from month to month.There are large longitudinal variations along the same pitch angle of magnetic field.The summer-midlatitude Ne distributions in June and December resemble closely the global Es morphology, suggesting that the same processes likely drive these geophysical variations.In the polar region the Ne distribution correlates better with the pitch angle derived from the magnetic field in the magnetosphere than from the 100 km, as expected for the magnetosphere-driven auroral electron precipitation.There is a clear semiannual variation in the nighttime zonal mean Ne in the D-region with peaks near the equinox months.We were able to derive the daytime and nighttime solarcycle variations in the D-region Ne from the MetOp RO data that had stable operation and uniform sampling at ~9AM and ~9PM local times.There is a good correlation between the long-term Ne variations and the F10.7 flux during 2007-2021.
Spatiotemporal sampling is critical to characterize large diurnal and geographical variabilities in the ionosphere.Because there are more high-rate RO measurements than podTEC, it provides a clear advantage to sample the D/E-region ionosphere better with the high-rate RO observations.The combined spatiotemporal observations from COSMIC-2 and Spire constellations provide an unprecedented coverage for studying the ionosphere since 2021.As the constellation continues to expand, it is feasible to improve the diurnal sampling to one-hour intervals and a spatial sampling better than the 8°´4° longitude-latitude grid as used in this study.The dense spatial sampling will perhaps allow a tomographic retrieval of Ne, to improve the horizontal resolution of the GNSS-RO technique.
For future algorithm developments we plan to apply the retrieval to the podTEC and ionPhs data to obtain a full Ne profile in the ionosphere.The podTEC data from CDAAC contain the hTEC profile with absolute calibration that can be readily processed with the algorithm developed in this study for an Abel inversion of Ne.The Abel inversion can be further improved to include the podTEC data from both positive and negative elevation angles, making it less dependent on the orbital altitude of the RO receiver.It is also feasible to develop an algorithm for the joint retrievals of TEC and Ne using the Δℎ profile, which should work in the cases where the absolute hTEC calibration is challenging.Finally, we will conduct a retrieval exercise for the tomographic inversion of TEC and Ne once the GNSS-RO constellations begin to provide a sufficient sampling density.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1,algorithm is extended to use of both  !"%# and  !"%$ profiles.Slightly different thresholds for quality control are used in the new algorithm.Because of the need in quality control for both  !"%# and  !"%$ profiles, the total number of Ne retrievals is slightly lower than the algorithm that uses the  !"%# profile only.
Although the inversion module bears the essentially same architecture, the weighting functions and linearization vectors are somewhat different.The new algorithm is linearized on a climatological profile derived from the IRI and FIRI models, instead of zero values.The weighting functions in the new algorithm are extended to a higher altitude such that it can handle the measurement and state vector up to 800 km.Both 'bottom-up' and Abel weighting functions are computed and stored in the new algorithm, making it flexible to retrieve Ne from different input measurement vectors, namely Δℎ (without absolute  !" calibration) and hTEC (with absolute  !" calibration) profiles.The capability of retrieving the entire ionospheric Ne from hTEC allows the new algorithm to process the podTEC data produced at CDAAC for COSMIC-1 and COSMIC-2.Initial testing of the retrieval from podTEC data shows that the optimal estimation inversion with the Abel weighting function produces better-quality (much less oscillations) Ne in the lower ionosphere than the inversion from the onion-peeling approach [38,87,88].The algorithm architecture in this study can readily accommodate further extensions of the Ne retrieval from the podTEC data, such as the use of podTEC from both positive and negative elevation angles.
This study added a new module to retrieve TEC from the high-rate RO data with a top up to 120 km.Preliminary results from this module were presented in section 4.2, revealing scientifically-useful TEC measurements on a global basis.Unlike the TEC retrieval derived from the Ne profile retrieved from the podTEC data, this TEC retrieval method does not require absolute calibration in the hTEC or  !" profile, as in most of the high-rate RO measurements.However, it does require the RO top to reach at least 120 km in order to observe a sufficient vertical gradient in the hTEC profile.This added TEC retrieval capability will significantly improve the spatiotemporal sampling of global ionosphere on a daily basis.Thus, further development and refinement are expected for this TEC module as the GNSS-RO constellations produce more the high-quality high-rate RO profiles.

Appendix B: GNSS-RO Data Availability and Statistics
The global GNSS-RO observations can be perhaps divided into three periods in terms of the total number of daily RO profiles [

D! /_0
Fit linear F-region contributions to % !"_#$_&'(' 30km< h t <60km F-region: Fit linear F-region contributions to ! !"#$ 50km< h t <75km F-region: space.The BlackJack RO receiver on CHAMP, developed by NASA Jet Propulsion Laboratory (JPL), was able to track a dual-frequency GPS signals for precise (cm accuracy range) orbit determination and continuous coverage [89].In its occultation mode, the receiver software was also able to schedule every 50 Hz tracking of setting occultations of up to four GPS satellites.The BlackJack on CHAMP has only aft antennas for the RO sounding, which yielded ~250 profiles per day.The sampling was further improved by tracking rising occultations with the open-loop (OL) tracking successfully demonstrated on the SAC-C and later COMSIC-1 satellites.This advance led to the dual-antenna (aft and fore) and OL tracking IGOR (Integrated GPS and Occultation Receiver) implemented by the COMSIC-1 constellation [90], which produced ~700 daily ROs per satellite, or an average total of ~4000 daily ROs from the six-satellite constellation.A recent boost in RO acquisitions comes from the combination of a new four-antenna TriG (Tri-band GNSS) receiver [91] implemented by the COSMIC-2 constellation and commercial CubeSat constellations (e.g., Spire) [92].  (6 Constellation 2018-varying varying 90°S/N 170-400 G,R,E,J ~4000 (7)  ~12,000 (8)  Other (9)  Constellations 2020-varying varying varying varying varying varying    The full diurnal coverage with GNSS-RO was not achieved on a monthly basis, until the COSMIC-1 constellation was launched.Prior to COSMIC-1, CHAMP made a diurnal sampling through a single orbital precession from multi-month or seasonal composites.The lengthy averaging for ionospheric observations can create aliasing problems between diurnal and seasonal variations, whereas the averaging over the same month from different years may bring interannual variability to the diurnal characteristics.Hence, a dense spatiotemporal sampling on a monthly basis is needed and critical for characterizing ionospheric variability.
While COSMIC-1 achieved the needed monthly diurnal coverage initially, the 6-satellite constellation degraded after 2016 with a fewer satellite in operation.As seen in Fig. B2, large gaps appeared in the GNSS-RO LST sampling.This creates a challenge to use the COSMIC-1 data for studying solar-cycle variations of the ionosphere.As a result, the long-term monitoring of ionospheric variability is limited to the observations at a few LSTs such as the 9AM and 21AM at which the combined operation from MetOp-A/B/C has provided more than 14 years of RO data.
A great boost in diurnal coverages comes from the recent COSMIC-2 and Spire constellation.Both constellations are capable of sampling a full LST range within a month, but COSMIC-2 covers only low latitudes (<44°S/N) with a good longitudinal sampling.On the other hand, the Spire constellation has a good pole-to-pole coverage with a fewer number of longitudinal samples at low latitudes, which was a similar situation from COSMIC-1.The combined COSMIC-2/Spire constellation, therefore, will provide an unprecedented spatiotemporal coverage for studying the ionosphere in the years to come.In ionospheric limb sounding the orbit altitude of RO receivers can have a significant effect on the retrieval, because of the extensive vertical range of ionosphere.Satellites are often orbiting inside the ionosphere.Depending on the orbital altitude, a significant portion of the ionosphere can reside above the satellite.Although the orbit altitude may affect little on the bending angle measurement for the atmosphere, the retrieval algorithms for TEC and Ne need to take into account the ionosphere portion that is not in the RO line-ofsight, which is dependent on orbit altitude.
As shown in Fig. B3, the GNSS-RO satellite altitude varies widely between 300 and 800 km.The orbit of satellites such as CHAMP, GRACE and C/NOFS decreased gradually during their mission lifetime, whereas others (e.g., operational weather satellites) were able to keep the orbital altitude nearly constant.To form their constellation with a uniform LST coverage, COSMIC-1 and COSMIC-2 had to raise and lower the altitude of its member satellites, respectively.After more than one year of the constellation formation, the COSMIC-2 satellites are now all operating in 525-540 km altitudes.Its companion polarorbiting constellation from Spire, under contract with NOAA CWDP, is also operating at the similar orbital altitude.Other operational features, such as the raised orbit by FY3C after 2018 and the MetOp-A decommission toward the end of 2021, are also evident.

Figure 2 .
Figure 2. Two examples of the RIE in the MetOp-A RO profiles from 2020d162: SNR profiles (left panels) and RIE profiles (right panels).The red line in the right panels is  -./ .The F-region and Es-induced RIEs are evident in the first case (top panels), whereas the sloped RIE profile is shown in the second case (bottom panels).The vertical gradient  &'('_*+,, profile from RIEs is ~1.25´10 -4

Figure 3 .
Figure 3. (a) State vector   and (b) measurement vector   in the linearization model, and (d) the associated weighting functions K (red) and their comparison with the Abel functions (black) for selected tangent heights between 50 and 500 km.The weighting functions K are flipped in sign purposely, to help the comparison.The intermediate step for the calculated RO hTEC and its linear fit (red line) from ht = 30-60 km is illustrated in (b), from which   is derived as their difference.See text for more discussions.
, where the Ne profile is retrieved up to the RO top.Since most of the COSMIC 50-Hz RO data have a top height around ht~125 km [see Appendix B], the D-region and a good portion of the E-region can be resolved.But the number of COSMIC-1 ROs started to decrease after 2017.The sampling reduction creates a problem for uniform local time coverage during the later years of COSMIC-1 mission.On the other hand, the operation from sun-synchronous MetOp satellites has been stable throughout the mission, acquiring 100-Hz ROs up to ht~90 km for the D-region.Therefore, the MetOp-A/B/C data are quite useful for studying long-term D-region Ne variations at two sunsynchronous local times (8-11h and 20-23h).

Figure 4 .
Figure 4. Illustration of a COSMIC-1 RO Ne profile retrieved from January 1, 2008.Differences between the Δℎ measurement profile and the linearized model   are inverted to produce a new state vector  for Ne.

Figure 5 .
Figure 5. COSMIC-2 Ne retrievals in the tropics (5°S-5°N) as a function of local time from Jan 1, 2021.The Ne value in the linearization state vector   is the red line, and the FIRI models for three tabulated F10.7 (75, 130, 200 sfu) levels are the blue curves.

Figure 7 .
Figure 7. Zonal mean Ne profile from MetOp-A (top panels) and Spire (bottom panels) as a function of magnetic latitude for October 2021.
Fig.8c produces a mean standard deviation of 1.03 TECu of the retrieved TEC from five ht levels at 140, 160, 180, 200, 220 km, whereas Fig.8d the retrieved TEC error increases to 1.75 TECu if the lower levels (100, 120, 140, 160, 180 km) are used.The Δℎ measurements near the F2 peak (ht =180-220 km) are generally weighted more than those from other heights, as indicated by the regression coefficients in the simulated data.The weight coefficient decreases for the Δℎ measurements at ht above the F2 peak (not shown), as expected for the measurements not going through the bulk of the ionosphere.

Figure 9 .
Figure 9. COMSIC-1 Ne at 60, 80, 100 and 120 km as a function of LST and latitude for March, June, September and December.The scaling factor for the color bar is indicated by the number in the parenthesis.The solid line denotes solar zenith angle c = 90°.

Figure 10 .
Figure 10.The model parameters (H' and ) for September and December as a function of magnetic latitude and local time from the fitting of COSMIC-1 D-region (60-80 km) Ne to Eq.10.The solid line is the terminator where the solar zenith angle c at 100km is equal to 90°.

Figure 11 .
Figure 11.COSMIC-1 monthly Ne maps at 70 km from all local times for March, June, September and December.Magnetic field-line pitch angles (a.k.a., dip angles) are contoured at an interval of 20°.The pitch angles on the cylindrical-projection map (left) are derived from the magnetic field at a 100-km altitude whereas the ones on the azimuthal polar maps are from the field at one Earth radius above the surface.The longitude-latitude grid size is 8°´4° for these maps.

Figure 13 .
Figure 13.E to F-region transition in Ne for December 2021 from the Spire constellation.The Spire receivers acquire GNSS-RO routinely with a top up to ~175 km, allowing the Ne retrieval to be extended to the F1-region.

Figure 14 .
Figure 14.Solar cycle variations of daytime (top) and nighttime (bottom) Ne at 70 km as a function of magnetic latitude from the combined MetOp-A/B/C and SAC-C (before 2008) monthly averages.The long-term variations at other altitudes (60 and 80 km) can be found in Fig.S5.The time series of daily Ly-a is the overplotted white line.The Ly-a data are obtained from the GSFC's Space Physics Data Facility (SPDF) OMNIweb service.

Figure 15 :
Figure 15: Example electron density profiles for the GNSS-RO and ionosonde measurements.The vertical line corresponds to the minimum return frequency observed by ionosondes, fmin.(a) corresponds to a daytime profile while (b) shows a nighttime profile.

Figure 16 :
Figure 16: Scatter plot of the electron densities predicted for each GNSS-RO altitude within the ionosonde fmin altitude and hmE.The ionosonde profiles are interpolated at each RO altitude for comparison.Unshifted profiles are compared in (a) while the ionosonde profiles shifted in altitude are displayed in (b).The linear regression parameters are also shown for both panels.

Figure 17 :
Figure 17: Histograms of the Ne difference (Ne,GNSS-RO -Ne,ionosonde) at each GNSS-RO altitude between the fmin altitude and hmE.(a) shows the unshifted Ne differences while (b) shows the shifted profile differences.The colorbar corresponds to the number of measurements within each bin.

Figure 18 :
Figure 18: Average daytime (LST=0800-1600) Ne for the GPS-RO and ionosonde techniques for (a) unshifted and (b) shifted profiles.The average fmin is displayed as a vertical line while the average hmE is displayed as a horizontal line.

Figure 19 .
Figure 19.RO-TEC retrievals using the modeled coefficients in Figure 8 and Eqs.(11-12), and their comparison with Ground-TEC.The MetOp-A RO data from June 9-18, 2020 are used because its high-rate ROs had a top up to ht ~290 km.The 1:1 line is plotted to guide the comparisons where daytime and nighttime measurements are in red and blue respectively.

Figure 20 .
Figure 20.TEC maps from multiple GNSS-RO satellites/constellations (Spire, Metop-A, COSMIC-2, FY3-C/D, Kompsat-5, TSX, PAZ) for October 2021.Satellite data with the high-rate RO acquisition top < 100 km are excluded in this TEC climatology.The monthly data are aggregated into a 4°´8° latitude-longitude and 2-hourly local time bins.The gridboxes without any measurements are colored in white.

Figure 21 .
Figure 21.As in Figure 20 but for August 2008 observed by COSMIC-1, GRACE, TSX, and SAC-C.

Figure S1 :
Figure S1: As in Fig.10 but for COMSIC-1 monthly Ne in the D/E-region as a function of local solar time (LST) and latitude from 2006-2020.

Figure S2 :
Figure S2: As in Fig.12 but for COSMIC-1 Ne maps at 70 and 110 km from all months in 2006-2020.

Figure S3 :
Figure S3: As in Fig.13 but for the 2-hourly zonal mean from Spire data in Sep, Oct, and Nov 2021, and in Jan 2022.

Figure S4 .
Figure S4.Spire Ne monthly climatology at the magnetic equator (6°S-6°N) as a function of local time and height from September, October, November, December of 2021, January of 2022.The NOAA CWDP NRT data delivered to CDAAC for September 2021 have 17 days, and the other months are complete.

Figure S5 .
Figure S5.As in Fig.14 but for the solar-cycle variations at 60 and 80 km.

Figure S6 .
Figure S6.As in Fig.20 but for TEC maps from Nov and Dec 2021 and Jan 2022.Author Contributions: Conceptualization, D.Wu; methodology, D.Wu; software, D.Wu; validation, D.Emmons and N.Swarnlingam; formal analysis, D.Wu; investigation, D.Wu; resources, D.Wu; writing, D.Wu and D.Emmons; funding acquisition, D.Wu.All authors have read and agreed to the published version of the manuscript.Funding: The work is supported by NASA's Living With Star and Sun-Climate research funds to Goddard Space Flight Center (GSFC).

Figure A1 .
Figure A1.Algorithm flow diagram in this study in comparison with the algorithm developed earlier by Wu [2018].

( 7 )
NOAA's Commercial Weather Data Pilot (CWDP) NRT data published by CDAAC from September 2021 to the present.The daily RO number is from November 2021 statistics.

( 8 )
NASA's Commercial Smallsat Data Acquisition (CSDA) data from November 2019 to the present.The daily RO number is from November 2021 statistics, of which approximately 880 are overlapped with the NOAA CWDP profile.

( 8 )
Other commercial GNSS-RO data include the constellations such as GeoOptics and PlanetIQ.

Figure B1 .
Figure B1.Number of daily GNSS-RO profiles since the CHAMP mission.Recent rising numbers are mainly from the COSMIC-2 and Spire constellations.The Spire RO numbers are based on the data archived under NASA's CSDA program.

Figure B2 .
Figure B2.Time series of daily mean equator-crossing LSTs of the RO measurement from different missions.The LSTs from satellite's ascending and descending orbits are separated roughly by 12h.Each orbit node has two close LSTs separated by ~1.5h, which comes from the RO measurements made by the fore and aft viewing antennas.Most of the modern RO instruments have at least two receivers to take measurements from the fore and aft directions.

Figure B3 .
Figure B3.Daily average orbital altitude of the past and current GNSS-RO missions.The altitudes are estimated from the LEO satellite position in the level-1b data.The Spire data contain only satellites (cubesats) from NOAA CWDP.

Table A1 .
Summary of GNSS-RO data in this study