A New Mapping Function for Spaceborne TEC Conversion Based on the Plasmaspheric Scale Height

: The mapping function is crucial for the conversion of slant total electron content (TEC) to vertical TEC for low Earth orbit (LEO) satellite-based observations. Instead of collapsing the ionosphere into one single shell in commonly used mapping models, we deﬁned a new mapping function assuming the vertical ionospheric distribution as an exponential proﬁler with one simple parameter: the plasmaspheric scale height in the zenith direction of LEO satellites. The scale height obtained by an empirical model introduces spatial and temporal variances into the mapping function. The performance of the new method is compared with the mapping function F&K by simulating experiments based on the global core plasma model (GCPM), and it is discussed along with the latitude, seasons, local time, as well as solar activity conditions and varying LEO orbit altitudes. The assessment indicates that the new mapping function has a comparable or better performance than the F&K mapping model, especially on the TEC conversion of low elevation angles.


Introduction
Benefitting from the accomplishment of the global navigation satellite system (GNSS) and many low Earth orbit (LEO) satellite missions, ionospheric exploration is greatly facilitated by various spaceborne measurements. Several successful LEO missions, such as Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC), contribute greatly to the ionospheric modeling and data assimilation system [1][2][3][4][5]. One significant product provided by these satellites is the sounding measurements that contain the total electron density content (TEC) along the signal paths. To obtain the absolute TEC from the raw GNSS/LEO observations, data analysis centers and scientific researchers have devoted great efforts into the main procedures, including the cycle slip detection and correction, carrier-phase to pseudorange leveling, multipath effect correction, and the differential code bias (DCB) estimation [6][7][8][9][10]. The conversion between the slant TEC (STEC) along the ray path and the vertical TEC (VTEC) in the zenith is an essential procedure during the data processing. An obliquity factor called mapping function (MF) is commonly used to do the conversion. MF is a crucial parameter in the estimation of receiver DCB by assuming the simultaneous VTECs transformed from two GNSS slant observations of one LEO antenna are equal [8]. Furthermore, the mapping function plays a significant role in ionospheric and plasmaspheric TEC modeling due to the fact that the accuracy of the MF is highly correlated with the estimated VTEC and DCB [11].
Several mapping functions are widely used in the TEC conversion, such as the thin layer model (TLM) [12], the simple "geometric" model proposed by Foelsche and Kichen-The scale height is one of the most important characteristics in the ionosphereplasmasphere system due to the role in shaping the electron density profile [25,[27][28][29]. There are various definitions of scale heights in published literature, including the theoretical plasma scale height, the vertical scale height, and the effective scale height. Alessio Pignalberi et al. [25] discussed the possible relationships among the different definitions of scale height. In this study, the plasmaspheric scale height adopts the definition of effective scale height, in which the analytical formulation of fitting the electron profile is chosen as exponential function. The new scale-height-based mapping function is not associated with specific IEH or thin-layer assumption but relates to the plasmaspheric scale height deduced from the exponential layer ionospheric model proposed by Stankov et al. [30]; i.e., the vertical electron density profile can be represented as: where N e (h 0 ) is the base electron density of plasmasphere (here set to be the electron density at the receiver height h 0 ), and H P is the plasmaspheric scale height. Thus, the mapping function can be written as the ratio of slant and vertical integral of electron density along the slant ray and the vertical path: where → r 0 and → r c represent the locations of LEO and GNSS satellites from the Earth center, respectively; h 0 and h c are the corresponding satellite heights above Earth's surface; z is the zenith angle of the ray at the LEO receiver. Figure 1 shows the geometric schematic of the mapping function. If the ionospheric spherical symmetry assumption is adopted, together with Equation (1), the numerical expression of the scale-height-based mapping function is where r v = r 0 ·sin(z); r 0 represents the radius of LEO satellite orbit from the Earth center; z is the zenith angle of the ray at the LEO receiver. M(z) can be obtained by numerical integral if H P is known. Since the plasmaspheric scale height varies with the solar activity, season, local time, and location [26], this mapping function is embedded with temporal and spatial variations. It should be mentioned that the scale-height-based mapping function implicitly assumes the exponential distribution of the ionosphere and plasmasphere. Considering the general orbit range of spaceborne satellites, the assumption is reasonable when dealing with LEO-based TEC conversion. Besides the numerical expression, an analytical solution for the mapping function can be written as: where I = r 0 2H P cot(z) and er f c( ) is the complementary error function. The process to obtain Equation (4) is shown in Appendix A. There are differences between the numerical and analytical expression of M(z) because the analytical solution adopts some approximations. We noticed that Hoque & Jakowski [15,31] presented a multi-layer mapping function approach for potential space-based augmentation system (SBAS) service under the assumption of Chapman layer distribution of the vertical electron density depending Remote Sens. 2021, 13, 4758 4 of 16 on the peak ionization height and atmospheric scale height. Here, we adopt the hypothesis of an exponential scale height to deal with the LEO-based slant to vertical TEC conversion.
the hydrostatic path delay with only one free parameter: "effective height". When applied to the LEO-based TEC conversion, the IEH is assumed at the height of a thick spherical layer above the receiver, hypothetically concentrating the whole impact of the ionosphere. The F&K mapping function is written as where ℎ and are the radii of the assumed layer and LEO satellite, respectively; is the zenith angle. The mapping factor may differ a lot with diverse IEH specifications.

Results
We calculate a mapping factor grid of zenith angle varying between 5° and 80° with 5° resolution, and varying in (100, 6000) km with a 50 km step according to equation (3) and (4). The obliquity factor of arbitrary TEC conversion can be obtained by interpolation with specified and . The grid map along with zenith angle and retrieved from the numerical and analytical solutions are demonstrated in Figure 2 (panel (a) and (c)). As an example, the spaceborne receiver is assumed at 800 km and the transmitter is assumed at 20,000 km. Moreover, the F&K model with the shell height varying from 100 km to 6000 km is displayed in panel (b). The blank in panel (b) is caused by the unreasonable result of the F&K model when the IEH is lower than the assumed LEO satellite height at high zenith angles. Similarly, the analytical resolution is not available when the ray path is approaching the zenith direction. We see the mapping factor increases when the zenith angle gets greater, and the difference brought by varying the value or shell height selection becomes significant when the zenith angle is larger than 40°. The mapping factor has similar variations along with the scale height or IEH. The numerical The scale-height-based mapping function is compared with the F&K method in this study. The F&K mapping function was originally proposed to describe the dependence of the hydrostatic path delay with only one free parameter: "effective height". When applied to the LEO-based TEC conversion, the IEH is assumed at the height of a thick spherical layer above the receiver, hypothetically concentrating the whole impact of the ionosphere. The F&K mapping function is written as where R shell and R orbit are the radii of the assumed layer and LEO satellite, respectively; z is the zenith angle. The mapping factor may differ a lot with diverse IEH specifications.

Results
We calculate a mapping factor grid of zenith angle z varying between 5 • and 80 • with 5 • resolution, and H P varying in (100, 6000) km with a 50 km step according to Equations (3) and (4). The obliquity factor of arbitrary TEC conversion can be obtained by interpolation with specified z and H P . The grid map along with zenith angle and H P retrieved from the numerical and analytical solutions are demonstrated in Figure 2 (panel (a) and (c)). As an example, the spaceborne receiver is assumed at 800 km and the transmitter is assumed at 20,000 km. Moreover, the F&K model with the shell height varying from 100 km to 6000 km is displayed in panel (b). The blank in panel (b) is caused by the unreasonable result of the F&K model when the IEH is lower than the assumed LEO satellite height at high zenith angles. Similarly, the analytical resolution is not available when the ray path is approaching the zenith direction. We see the mapping factor increases when the zenith angle gets greater, and the difference brought by varying the H P value or shell height selection becomes significant when the zenith angle is larger than 40 • . The mapping factor has similar variations along with the scale height or IEH. The numerical and analytical results differ more with a large zenith angle and higher H P values (shown in panel (d)). and analytical results differ more with a large zenith angle and higher values (shown in panel (d)). The only free parameter in the scale-height-based mapping function is the plasmaspheric scale height. In Wu et al. [26], we have calculated the plasmaspheric scale height from the COSMIC measurements from 2007 to 2014 and proposed an empirical monthly climate model that has reasonable accuracy validated by the scale height retrieved from the International Satellites for Ionospheric Studies (ISIS) observations during its mission period. This model is suitable to provide the essential parameter in the scale-heightbased mapping function. Besides the COSMIC-derived , the scale height obtained by the global core plasma model (GCPM) provides an alternative option [32]. The GCPM provides empirically derived core plasma density and ion composition ( + , + , and + ) as a function of the geomagnetic and solar conditions throughout the inner magnetosphere. The model merges with the IRI model at low altitudes and is composed of separate models for the plasmasphere, plasmapause, trough, and polar cap. The is calculated from the overhead vertical TEC simulated with the GCPM electron density field and the base electron density at the assumed receiver altitude since the integral in equation (1)   is generally greater than the GCPM simulations, but they have generally similar variations along with the latitude, local time, season, and solar activity. The nighttime scale height is larger than that in the daytime and usually achieves the maximum before the sunrise, The only free parameter in the scale-height-based mapping function is the plasmaspheric scale height. In Wu et al. [26], we have calculated the plasmaspheric scale height from the COSMIC measurements from 2007 to 2014 and proposed an empirical monthly climate model that has reasonable accuracy validated by the scale height retrieved from the International Satellites for Ionospheric Studies (ISIS) observations during its mission period. This H P model is suitable to provide the essential parameter in the scale-heightbased mapping function. Besides the COSMIC-derived H P , the scale height obtained by the global core plasma model (GCPM) provides an alternative option [32]. The GCPM provides empirically derived core plasma density and ion composition (H + , H + e , and O + ) as a function of the geomagnetic and solar conditions throughout the inner magnetosphere. The model merges with the IRI model at low altitudes and is composed of separate models for the plasmasphere, plasmapause, trough, and polar cap. The H P is calculated from the overhead vertical TEC simulated with the GCPM electron density field and the base electron density at the assumed receiver altitude since the integral in Equation (1)  December solstice (January, February, December) in this work. The H P_COSMIC is generally greater than the GCPM simulations, but they have generally similar variations along with the latitude, local time, season, and solar activity. The nighttime scale height is larger than that in the daytime and usually achieves the maximum before the sunrise, which is consistent with the conclusion drawn in Wu et al. [26]. The seasonal variations are more predominant for COSMIC scale heights, with higher values in the winter hemisphere and lower at the summer half. Under high solar activity conditions, the scale height decreases for both the observational and empirical models. The obliquity factor will differ a bit Remote Sens. 2021, 13, 4758 6 of 16 when applying H P_COSMIC and H P_GCPM in the mapping function. At the high latitudes (±60 •~± 90 • ), the electron density is quite small and rarely changed, so the scale height value can be regarded as a constant. The H P derived from satellite observations is assumed more reliable when dealing with the realistic mapping issues, but, in the following study, the H P_GCPM is discussed because we are aiming at the comprehensive validation of the scale-height-based mapping function in a simulation way.
which is consistent with the conclusion drawn in Wu et al. [26]. The seasonal variations are more predominant for COSMIC scale heights, with higher values in the winter hemisphere and lower at the summer half. Under high solar activity conditions, the scale height decreases for both the observational and empirical models. The obliquity factor will differ a bit when applying _ and _ in the mapping function. At the high latitudes (±60°~±90°), the electron density is quite small and rarely changed, so the scale height value can be regarded as a constant. The derived from satellite observations is assumed more reliable when dealing with the realistic mapping issues, but, in the following study, the _ is discussed because we are aiming at the comprehensive validation of the scale-height-based mapping function in a simulation way.

Assessment by Global Statistics
To assess the performance of the new mapping function, the GCPM is used to generate simulating LEO-based TEC observations and examine the retrieved errors of both mapping models. The LEO satellite orbit is chosen at 800 km, which is the typical orbit altitude of the COSMIC-1 constellation that made a great contribution to the ionosphere and plasmasphere exploration. The signal transmitter is at 20,000 km, referring to the GNSS satellites. It should be noted that both the scale-height-based and the F&K mapping function assume the spherical symmetric distribution of the ionosphere, which means the electron density on the same sphere is identical everywhere. Thus, the simulated STEC of a different zenith angle is actually equal to the integral of vertical electron density along the slant ray path with a 50 km step. The VTEC retrieved from the slant TEC and mapping function is compared to the VTEC 'truth' integrated directly from the GCPM background (regarded as ). The relative root mean square (RMS) of each zenith angle is calculated as the assessment criteria:

Assessment by Global Statistics
To assess the performance of the new mapping function, the GCPM is used to generate simulating LEO-based TEC observations and examine the retrieved errors of both mapping models. The LEO satellite orbit is chosen at 800 km, which is the typical orbit altitude of the COSMIC-1 constellation that made a great contribution to the ionosphere and plasmasphere exploration. The signal transmitter is at 20,000 km, referring to the GNSS satellites. It should be noted that both the scale-height-based and the F&K mapping function assume the spherical symmetric distribution of the ionosphere, which means the electron density on the same sphere is identical everywhere. Thus, the simulated STEC of a different zenith angle is actually equal to the integral of vertical electron density along the slant ray path with a 50 km step. The VTEC retrieved from the slant TEC and mapping function is compared to the VTEC 'truth' integrated directly from the GCPM background (regarded as VTEC mod ). The relative root mean square (RMS) of each zenith angle z is calculated as the assessment criteria: where N is the TEC observation number of global grids at the zenith angle z. Given the difference between the numerical and analytical representation of mapping function, we checked the performance of the scale-height-based mapping function of both solutions, denoted as 'H P_N ' and 'H P_A ', respectively, hereafter. To introduce independent reference into the validation, the F&K mapping function is also considered with adjustable IEH represented as a function of the LEO orbit height and solar flux proxy F10.7 (F 107 ) according to Zhong et al. [19]: Thus, the IEH specifications take different satellite altitudes and solar activity levels into consideration, and the F&K mapping model is optimized by choosing a constant IEH. The zenith angle of the simulated slant ray path varies from 0 • to 80 • in steps of 5 • , and the geophysical location of the receiver changes with horizontal resolution of latitudinal 5 • from −90 • to 90 • and longitudinal 10 • from −180 • to 180 • . The simulating experiments are executed on one day in four seasons with a temporal resolution of 2 h under both the low and high solar activity conditions. The mean retrieved RMS errors for each season are calculated according to the zenith angles and shown in Figure 5. The F&K mapping function depends on the fixed IEH globally; thus, the retrieved errors are closely associated with the ionospheric variations of the GCPM background. According to Figure 5, the behavior of the numerical and analytical mapping functions differs from each other while indicating better retrieved results than the F&K model, especially with increasing zenith angles greater than 40 • . In the LSA years, the retrieved errors of the numerical mapping function are slightly less than those of the analytical solution at lower zeniths. When the solar activity is more active, the mapping errors of all the methods decrease to some extent, and the analytical mapping function achieves the best performance throughout the year. The difference between the two scale-height-based models is associated with the variation in the H P value under different solar activity conditions and the zenith angles. From the statistics, we conclude that, with the LEO satellite located at the 800 km, the H P -based mapping function is a competitive approach in comparison with the F&K model, especially at low elevation angles; i.e., at high zenith angles. where is the TEC observation number of global grids at the zenith angle . Given the difference between the numerical and analytical representation of mapping function, we checked the performance of the scale-height-based mapping function of both solutions, denoted as ' _ ' and ' _ ', respectively, hereafter. To introduce independent reference into the validation, the F&K mapping function is also considered with adjustable IEH represented as a function of the LEO orbit height and solar flux proxy F10.7 ( 107 ) according to Zhong Thus, the IEH specifications take different satellite altitudes and solar activity levels into consideration, and the F&K mapping model is optimized by choosing a constant IEH. The zenith angle of the simulated slant ray path varies from 0° to 80° in steps of 5°, and the geophysical location of the receiver changes with horizontal resolution of latitudinal 5° from −90° to 90° and longitudinal 10° from −180° to 180°. The simulating experiments are executed on one day in four seasons with a temporal resolution of 2 hours under both the low and high solar activity conditions. The mean retrieved RMS errors for each season are calculated according to the zenith angles and shown in Figure 5. The F&K mapping function depends on the fixed IEH globally; thus, the retrieved errors are closely associated with the ionospheric variations of the GCPM background. According to Figure  5, the behavior of the numerical and analytical mapping functions differs from each other while indicating better retrieved results than the F&K model, especially with increasing zenith angles greater than 40°. In the LSA years, the retrieved errors of the numerical mapping function are slightly less than those of the analytical solution at lower zeniths. When the solar activity is more active, the mapping errors of all the methods decrease to some extent, and the analytical mapping function achieves the best performance throughout the year. The difference between the two scale-height-based models is associated with the variation in the value under different solar activity conditions and the zenith angles. From the statistics, we conclude that, with the LEO satellite located at the 800 km, the -based mapping function is a competitive approach in comparison with the F&K model, especially at low elevation angles; i.e., at high zenith angles.

Global Variations of Mapping Errors
To further analyze the time and spatial dependence of the mapping functions, the global maps of mean retrieved error with the zenith angle fixed at 40 • are displayed in Figures 6 and 7 of the LSA and HSA years, respectively. The scale-height-based mapping function of the numerical or analytical solutions and the F&K model exhibit quite different characteristics with the local time, geographical latitude, and seasons. Corresponding to the predominant seasonal asymmetric features of H P , the VTEC retrieved by the new method is overestimated in the nighttime winter hemisphere, especially when using the analytical mapping factor. Both the scale-height-based models underestimate the VTEC during the daytime at low and equatorial latitudes. Consistent with the statistics in Figure 5, the general performance of the numerical mapping function is fairly good, with the mapping error varying between −5% and 5%. The latitudinal and temporal variations of the mapping errors reflect the reliability of the exponential assumption in the GCPM model. It proves that the variations in the ionosphere and plasmasphere are more complicated and nonuniform at low latitudes. For the F&K model, the relative error maps indicate the inhomogeneous variations of the ionospheric effective height of the background. The VTEC is significantly underestimated during the whole day at low latitudinal areas, while it is overestimated at higher latitudes. The performance of all the methods is remarkably improved under higher solar activity, as shown in Figure 7. Generally, the scale-height-based mapping function achieves better conversion results than the F&K method, especially in the equatorial anomaly regions between −30 • and 30 • .   The VTEC at low latitudinal areas beside the geomagnetic equator rather underestimated when applying the numerical mapping function and the F&K model. At higher elevation angles, the TEC conversion errors brought by various mapping functions are negligible with minor differences. However, the residuals increase pronouncedly when the elevation angle is lower, especially for the F&K model. The performance of the analytical solution is generally the best at mid-and high latitudes, especially in a high solar activity year. The winter hemispheric overestimation along the geomagnetic equator is associated with the seasonal variation in scale height. As mentioned earlier, the geographical error distribution of the F&K model actually reflects the inhomogeneity of the IEH and the ionospheric background model. The results depend greatly on the IEH selection and the ionospheric model used to specify the IEH. The scale-height-based mapping model considers the vertical structure and variations in the ionosphere by involving the plasmaspheric scale height, and it is more convenient and reliable to provide realistic TEC conversion results. Moreover, the degradation in the performance is minor for scale-height-based mapping along with an increasing zenith, so it is especially suitable to be used on observations of low elevation angles. Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 18

Influence of LEO Orbit Altitudes
Given the diversity of the orbit altitudes of different LEO satellites, we set up another two scenarios for more detailed discussion: the receiver being located at 500 km and 1400 km, respectively. The performances of the scale-height-based mapping function and the F&K model with the IEH specified according to Equation (7) are investigated in the same way as in Section 4.1. In Figure 10, when the LEO orbit altitude is 500 km, the mapping errors increase. The reason is that, at 500 km of altitude, the reliability of considering an exponential vertical profile for the electron density (with the associated scale height) is quite low [33][34][35]. This hypothesis adopted in the scale-height-based mapping will produce higher errors than at 1400 km of altitude. In an LSA year, both the numerical and analytical mapping function achieve fewer RMS errors than the F&K model, and the analytical solution is slightly better than the numerical one. Under an HSA condition, the numerical mapping factor results in similar or even a bit worse results than the F&K model below a zenith angle of 60 • , while the analytical mapping function is promising to obtain the least mapping errors. The situation is quite different with the LEO satellite higher at 1400 km. According to Figure 11, the numerical mapping function remarkablely improves the mapping errors in both LSA and HSA years. The analytical solution has a relatively poor peformance at high elevation compared to the F&K method but still functions well with greater zenith angles. Figure 8. The longitude-and latitude-dependent variations of the retrieved VTEC mean deviation in the December solstice in LSA years. The zenith angle is chosen at 20° (panels (a), (b), and (c)), 40° (panels (d), (e), and (f)), and 60° (panels (g), (h), and (i)), respectively. The ' _ ', ' _ ', and 'F&K' represent the numerical and analytical mapping functions and the F&K model, respectively.

Conclusions
This paper proposed a new model of an LEO-based TEC mapping function based on a prior model of the plasmaspheric scale height. The numerical and analytical forms of the mapping function are illustrated, and the mapping factor grid is obtained as a function of the scale height and zenith angle. The scale height and mapping factor are accessible to the public at the database: https://www.researchgate.net/publication/353703384_map-ping_factor_grid(accessed on 08/2021) and https://www.researchgate.net/publication/353704297_plasmapsheric_scale_height(accessed on 08/2021).
The new mapping model is driven by the only free parameter, , obtained either from realistic observations or the empirical model. The performance of this mapping function is assessed by the simulated TEC conversion experiments based on the GCPM electron density field. This method is promising to be applied to LEO-based TEC conversion as a comparable or better alternative to the F&K model according to the assessments under various spatial and temporal specifications. Instead of collapsing the ionosphere into To understand the distinctive performance of the scale-height-based mapping function refering to different LEO orbit height, it should be noted that the H P obtained in this work is basically the averaged scale height under the assumption that the ionosphere and plasmasphere are varying roughly in an exponential trend above the satellite orbit. In fact, the actual scale height is varying with the altitude gradually below the O + -H + transition height, which peaks at 700 km and achieves an average value of 862 km [36], and approaching constant when the dominant ion becomes H + or H + e in the plasmasphere [26]. Therefore, the averaged H P has a larger deviation with the varying scale height, especially when the satellite is locating below the transition height, such as in the 500-km-scenario. The retrieved VTEC is generally less than the truth since the mapping factor is overestimated because of the averaged H P . It agrees with the underestimation of the VTEC of the numerical solution in Figures 6-9. Along with the increase in the orbit height, the averaged H P is closer to the realistic scale height and leads to fewer conversion errors. That accounts for the improvement in the performance of the numerical solution with the orbit varying from 500 km to 1400 km. Due to the assumption and approximation adopted in the analytical mapping function, a higher H P value leads to a smaller mapping factor ( Figure 2). Therefore, the analytical solution basically compensates the overestimation of M(z), sometimes presenting a better performance than the numerical one in an HSA year. Given the overall performances of the three methods, the scale-height-based mapping function with the numerical solution is the most stable approach under varying solar activity conditions and in changing LEO orbits, and it improves the F&K model significantly, especially at lower elevation angles.

Conclusions
This paper proposed a new model of an LEO-based TEC mapping function based on a prior model of the plasmaspheric scale height. The numerical and analytical forms of the mapping function are illustrated, and the mapping factor grid is obtained as a function of the scale height and zenith angle. The scale height and mapping factor are accessible to the public at the database: https://www.researchgate.net/publication/353703384_mapping_ factor_grid (accessed on 1 August 2021) and https://www.researchgate.net/publication/ 353704297_plasmapsheric_scale_height (accessed on 1 August 2021).
The new mapping model is driven by the only free parameter, H P , obtained either from realistic observations or the empirical model. The performance of this mapping function is assessed by the simulated TEC conversion experiments based on the GCPM electron density field. This method is promising to be applied to LEO-based TEC conversion as a comparable or better alternative to the F&K model according to the assessments under various spatial and temporal specifications. Instead of collapsing the ionosphere into a thin shell, the new formulation considers the exponential layer assumption of the electron density distribution and applies the plasmaspheric scale height to introduce seasonal, local time, and geolocation-dependent variations into the mapping model. The performance of the F&K model is highly correlated with the IEH specification at different latitudes during varying solar activity periods, as well as the orbit altitudes of LEO satellites. The scale-height-based mapping model has advantages in dealing with the spaceborne TEC conversion, especially with low-elevation-angle observations of higher orbit height.
We recognize that the major errors remaining in the TEC conversion are associated with the horizontal inhomogeneous distribution of the ionosphere. Therefore, a comprehensive mapping model involving the azimuth angle variation of the signal ray paths is now under development based on this investigation.

Acknowledgments:
The authors are very thankful for the open access provided by CDAAC for all the COSMIC data involved and the national aeronautics and space administration for GCPM sources codes. The processed data and figure datasets used for this paper are available at the site: https://doi.org/10.5281/zenodo.5205255 (accessed on 1 August 2021).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Assuming an arbitrary point P on the ray path in Figure 1, the relationship between the slant ray X and the vertical projection H to the receiver is approximately written as H ≈ Xcos(z) + X 2 sin 2 (z) 2r 0 (A1) where r 0 represents the radius of LEO satellite orbit from the Earth center; z is the zenith angle of the ray path. Under the assumption of ionospheric spherical symmetry, the slant TEC is obtained by integrating the electron density along the ray path, If a = sin 2 (z) 2r 0 H P and b = cos(z) 2H P , then N 0 is the electron density at the receiver. The vertical TEC is easy to obtain with integration as VTEC ≈ N 0 ·H P ; therefore, the analytic solution of M(z) relates to the scale height H P , zenith angle z, and the complementary error function where I = r 0 2H P cot(z).