**1. Introduction**

According to the three-step development strategy [1,2], the Beidou navigation satellite system (BDS) has experienced BDS-1 (1990–2003), BDS-2 (2003–2012) and BDS-3 (2016–2020). BDS-3 has started to run formally for the whole world on 31 July 2020, providing the Position, Navigation and Timing service (PNT), Radio Determination Satellite Service (RDSS), Global Short Message Communication (GSMC), Satellite-Based Augmentation Service (SBAS), Precise Point Positioning (PPP), Search and Rescue (SAR), and so on. Among them, RDSS is an important component and a distinctive feature of BDS that has advantages and a competitiveness which are different from Global Positioning Systems (GPS), Global Navigation Satellite Systems (GLONASS) and the Galileo system [3–5].

BDS-2 is composed of fourteen satellites and has provided services since 2012. Among the satellites are five geosynchronous earth orbit (GEO) satellites, five inclined geosynchronous orbit (IGSO) satellites and four medium earth orbit (MEO) satellites. Five GEO satellites can broadcast RDSS signals, and each satellite has two wide beams. BDS-2 RDSS can provide a timing service. The official claimed accuracy of BDS-2 RDSS is 50 ns for the one-way timing service, and 20 ns for the two-way timing service [6,7].

BDS-3 is composed of thirty satellites and has provided services since 2020. Among the satellites are three GEO satellites, three IGSO satellites and twenty-four MEO satellites.

**Citation:** Guo, R.; Wang, D.; Xing, N.; Liu, Z.; Zhang, T.; Ren, H.; Liu, S. Preliminary Analysis and Evaluation of BDS-3 RDSS Timing Performance. *Remote Sens.* **2022**, *14*, 352. https://doi.org/10.3390/rs14020352

Academic Editor: Jose Moreno

Received: 20 November 2021 Accepted: 8 January 2022 Published: 13 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Three GEO satellites can broadcast RDSS signals. In order to improve the capacity of the system, BDS-3 RDSS adopts narrow beams, and each satellite has six narrow beams to provide inbound and outbound links. At the same time, BDS-3 RDSS separates civil signals and military signals, which are modulated on different carriers to improve their isolation, and to realize the security and compatibility. The official claimed accuracy of BDS-3 RDSS is 50 ns for the one-way timing service, and 10 ns for the two-way timing service [6,7].

As a spatial information infrastructure, the availability and continuity of satellite navigation system services are important. Up to now, there are rare existing studies on BDS-3 RDSS timing performances. Based on BDS-3 RDSS signal measurements, this paper analyzed and evaluated the performance of the one-way timing service and two-way timing service for the first time, which provides technical support for RDSS timing application.

#### **2. Materials and Methods**

#### *2.1. One-Way Timing Principle*

This section describes the basic principle of the RDSS one-way timing service. In the one-way timing process, the user starts a timer at the beginning of a second (the local 1-pps). At the same time, it receives the *n*-th frame of the outbound signal (which carries a time stamp) transmitted by the central control system (CCS) through the satellite, and uses the time stamp of the *n*-th frame to stop the timer. Then, the timer measures the time interval Δ<sup>1</sup> [8,9]. The schematic diagram of the one-way timing principle is illustrated in Figure 1.

**Figure 1.** The one-way timing principle.

The clock difference between the user and CCS, denoted by Δ*ε*, can be obtained as:

$$
\Delta \varepsilon = \Delta\_1 - \tau - n\Delta t \tag{1}
$$

where Δ*t* = 125 ms is the frame period according to the interface control document (ICD) of the BDS-3 RDSS outbound signal, and *τ* is the forward propagation delay of the signal, which can be calculated as:

$$
\pi = t\_{equ} + t\_{R01} + t\_{R\tau1} \tag{2}
$$

where *tequ* is the one-way time delay of the equipment, which is stored in the user. *tR*<sup>01</sup> = *R*01/*c* + *tR*<sup>01</sup>−*pro* is the propagation delay from CCS to the satellite; *R*01 is the geometrical distance between CCS and the satellite; *<sup>c</sup>* is the light speed; *tR*<sup>01</sup> = *R*01/*c* + *tR*<sup>01</sup>−*pro* is the ionosphere and troposphere delay, which can be accurately calculated using the parameters broadcast by CCS. *tR<sup>τ</sup>*<sup>1</sup> = *Rτ*1/*c* + *tR<sup>τ</sup>*1−*pro* is the propagation delay from the satellite to the user; *Rτ*<sup>1</sup> is the geometrical distance between the satellite and the user; *tR<sup>τ</sup>*<sup>1</sup> = *Rτ*1/*c* + *tR<sup>τ</sup>*1−*pro* is the ionosphere and troposphere delay, which is calculated by the user according to the satellite position, user position and broadcast correction parameters [8,9].

From the above Equations (1) and (2), it can be seen that the main factors affecting the accuracy of RDSS one-way timing include GEO satellites' ephemeris errors, ionosphere and troposphere delay correction errors, various equipment delay errors, user position errors, and so on.

#### *2.2. Two-Way Timing Principle*

This section describes the basic principle of the RDSS two-way timing service. In a two-way timing process, the user also starts and stops a timer at the start of the local second and at the time stamp of the *n*-th frame received. What differs is that the forward propagation delay is calculated by CCS and is transmitted to the user. As shown in Figure 2, when the user receives the *n*-th frame, it sends an inbound signal to CCS immediately. CCS receives the inbound signal, measures the interval Δ2, and calculates the forward propagation delay *τ* [8–10], based on Equation (3).

**Figure 2.** The two-way timing principle.

Assuming that the signal drift of satellite in the process is ignored, the forward propagation delay can be obtained as:

$$
\tau = t\_{equ\\_one-way} + \left(\Delta\_2 - t\_{equ\\_two-way} - t\_+ - t\_-\right)/2\tag{3}
$$

where *tequ*\_*one*−*way* and *tequ*\_*two*−*way* refer to the equipment time delay of the one-way and two-way (including the time delay of CCS, satellite transmitters and users), Δ<sup>2</sup> refers to the signal round-trip time measured by CCS, *t*+ refers to the ionosphere and troposphere delay from CCS to the user through the i-th satellite and *t*<sup>−</sup> refers to the ionosphere and troposphere delay from the user to CCS [8–10].

From the above Equation (3), it can be seen that the main factors affecting the accuracy of RDSS two-way timing include ionosphere and troposphere delay correction errors, various equipment delay errors, and so on. Therefore, the two-way timing accuracy is higher than the one-way timing accuracy. It is worth noting that the signal drift of the satellite in the process is ignored when the forward propagation delay is calculated in Equation (3). However, the drift is related to the speed of motion, which means that periodic fluctuations can also appear in the two-way timing.

### **3. Results**

#### *3.1. Analysis and Evaluation Method*

The performance index of BDS-3 RDSS timing service includes timing accuracy, availability and continuity. The analysis methods of each index are described briefly below.

#### 3.1.1. Timing Accuracy

An RDSS receiver, of which its position is known precisely, is used to receive the outbound RDSS signal and then calculate the one-way timing accuracy. Meanwhile, the same receiver is used to send inbound two-way timing signal. CCS will calculate the timing accuracy and compare the results with the known reference 1 pps, which provide a time-frequency reference for the receivers, and measure the timing error of each epoch. The mean value, standard deviation and root mean square (RMS) can be obtained as [8]:

$$\overline{\tau} = \frac{1}{n} \sum\_{i=1}^{n} \tau\_i \tag{4}$$

$$\text{STD} = \sqrt{\frac{1}{n-1} \sum\_{i=1}^{n} |\tau\_i - \overline{\tau}|^2} \tag{5}$$

$$\tau\_{RMS} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} |\tau\_i|^2} \tag{6}$$

where *τ<sup>i</sup>* is the discrete point of timing accuracy, *n* is the number of sample points, *τ* is the mean value which reflects the constant deviation of timing results, STD is the standard deviation, which reflects the stability and variation range of timing results, *τRMS* is RMS, which is the comprehensive result of mean value and standard deviation, and reflects comprehensively the deviation and stochastic error of timing results.

#### 3.1.2. Availability

The availability of the timing service is the time percentage of when the timing accuracy meets the requirements, during the specified period, conditions and service area [6]. For BDS-3 RDSS one-way timing and two-way timing, the accuracy requirements are 50 ns and 10 ns respectively.

Taking the timing error sequence as the evaluation object, the available time percentage of timing service is calculated as:

$$Ava\_l = \frac{\sum\_{t=t\_{start}, inc=T}^{t\_{end}}bool\{EPk \le f\_{Acc}\}}{1 + \frac{t\_{end} - t\_{start}}{T}} \tag{7}$$

where *tstart* and *tend* is the start time and end time respectively, *T* is the sampling interval, *EPEk* is the timing error at a certain time, *fAcc* is the timing accuracy threshold and *bool*{∗} is the Boolean function.

#### 3.1.3. Continuity

RDSS messages are broadcast very second, and the continuity is evaluated without considering the message type. If the receiver receives a message at the current time but it is lost at the next time, it would be recorded as an interrupt, and the interruption duration is the time that receivers cannot receive the message [6].

In order to ensure that the statistical results reflected the messages continuity strictly, three different receivers were used in the experiment. This indicates that an interruption happened only when all three receivers lost the message at the time.

It should be known that the continuity shoud exclude the planned interrupted exception time, which will be published with the warning announcement on the BDS official website, subsequently.

#### *3.2. One-Way Timing Performance Analysis*

In order to evaluate the performance of BDS-3 RDSS one-way timing, the following analyses were conducted by the real measured data from the Beijing RDSS receiver. The assessment scenarios included two different states: normal state and satellite orbit maneuver.

#### 3.2.1. Normal State

From 14 February 2021 to 16 February 2021, the receiver collected data from beam two of the BDS-3 C59 satellite. During this period, the satellite was in the normal state, without satellite orbit maneuver or other abnormal status. In this evaluation, the receiver collected timing results for 3 days successively with a 1 s sampling interval, and the one-way timing accuracy curve can be seen in Figure 3.

**Figure 3.** RDSS one-way time service error (normal state).

It can be seen from Figure 3 that: (1) the one-way timing error of all epoch points is less than 30 ns, and the RMS is 5.45 ns, which is 20% smaller than that of the BDS-2 one-way timing error of 6.81 ns [11]; (2) The mean value of one-way timing error is 0.39 ns, which means that there is no significant deviation between the reference 1 pps and the timing results, the equipment delay was calibrated accurately and both ionosphere and troposphere delays were corrected properly; (3) The standard deviation of one-way timing error is 5.43 ns, which actually reflects the white noise in the signal. Moreover, the orbital errors, ionosphere and troposphere delay residual errors always appear as bias.

Furthermore, Figure 3 shows that one-way timing results had a period of 1 day, of which the characteristic was related to the GEO satellite orbit characteristic. In order to further analyze the characteristics of one-way timing error, the spectrum of the one-way timing is shown in Figure 4 and the error removing orbit period is shown in Figure 5.

**Figure 4.** Spectrum analysis results of RDSS time service error.

**Figure 5.** RDSS one-way timing error (remove orbit period).

It can be seen from Figure 4 that the maximum spectrum is 1.01 days (about 24 h 14 m 24 s), which is basically consistent with the orbital period of the GEO Satellite, of 23 h 56 m, and the reason for the inconsistency may be the observation data size and the measurement error. When compared with Figure 3, It can be seen from Figure 5 that the characteristics of the GEO satellite orbital period have been eliminated. However, there is still a periodic jump of the timing error, of which its magnitude is about 10 ns. This phenomenon is related to the RDSS timing principle and the receiver's algorithm. According to the Interface Control Document (ICD), GEO satellite updates the ephemeris every 6 s. Receivers need

to interpolate the ephemeris to obtain the satellite position at the timing epoch, and then calculate the timing results based on the fixed point. The interpolation algorithm is the main factor causing the timing error jump.

In order to further evaluate the performance under normal state, the availability and continuity results are given in Table 1. The availability and continuity of the RDSS one-way timing service are 100%.

**Table 1.** The availability and continuity of one-way timing service (normal state).


#### 3.2.2. Orbit Maneuver State

The ephemeris accuracy of GEO satellites is lower than that of MEO and IGSO satellites due to the orbital characteristics. The orbit of the GEO satellite needs to be adjusted and maintained regularly [12,13]. One-way timing accuracy is mainly limited by orbit determination accuracy during GEO satellite orbit maneuver, because it is related to satellite orbital error. Therefore, the data of one-way timing services were used for analysis and evaluated the orbit maneuver state and recovery of the C59 satellite from 8 February 2021 to 10 February 2021, where the maneuver period was from 11:30 to 13:30 on 8 February 2021.

It can be seen from Figure 6 that: (1) during the orbit maneuver state, the accuracy recovery time of one-way timing is 7.68 h. Specifically, from the start of orbit maneuver, the timing error is within 100 ns about 2.18 h, then the timing error increases to 1000 ns at about 5.08 h and then the timing error gradually decreases and finally returns to the level of 50 ns about 0.52 h. (2) Affected by orbit maneuver, the RMS of the RDSS one-way timing error is 15.33 ns, and its standard deviation is 15.31 ns. This is because in the process of RDSS one-way timing, receivers extrapolate the orbit data according to the ephemeris broadcast by satellites, which is greatly affected by the ephemeris errors of GEO satellites.

**Figure 6.** BDS-3 RDSS one-way timing accuracy (orbit maneuver state).

In order to further evaluate the orbit maneuver, the timing accuracy results are given in Table 2. The RMS of normal 8 February, and the recovery on 9 and 10 February are 5.52 ns, 5.49 ns and 5.46 ns, respectively.

**Table 2.** The accuracy of the RDSS one-way timing service (orbit maneuver state, RMS, unit: ns).


#### *3.3. Two-Way Timing Performance Analysis*

Similar to the one-way timing, the performance analysis of two-way timing also used the signal measurements from the Beijing RDSS receiver [13]. The assessment scenarios included two different states: normal state and satellite orbit maneuver.

#### 3.3.1. Normal State

From February 14 to 16, 2021, the receiver collected data from beam two of the BDS-3 C59 satellite. During this period, the satellite was in the normal state, without satellite orbit maneuver or other abnormal status. In this evaluation, the receiver adopted timing results for 3 days successively, with z 1 min sampling interval, and the one-way timing accuracy curve can be seen in Figure 7.

**Figure 7.** RDSS two-way timing error (normal state).

It can be seen from Figure 7 that: (1) the two-way timing error of all epochs is less than 8 ns, and the RMS is 3.59 ns, which is 34% smaller than that of the one-way timing error of 5.45 ns; (2) the mean value of the two-way timing error is 0.92 ns, which is not much larger than that of the one-way timing of 0.39 ns. This indicates that the equipment delay was calibrated accurately, and both the ionosphere and troposphere delays were corrected properly; (3) The standard deviation is 2.09 ns, which is reduced by 62% when compared with the one-way timing, indicating that GEO satellites orbital errors, ionosphere and troposphere delay errors are greatly reduced in the two-way timing process.

Furthermore, Figure 7 shows that two-way timing results also have a period which is connected to ionosphere and troposphere parameters, which are related to the sun zenith angle, since the GEO satellite is located in geostationary orbit. The spectrum of the two-way timing is shown in Figure 8 and the error removing the orbit period is shown in Figure 9.

**Figure 8.** Spectrum analysis results of RDSS timing error.

**Figure 9.** RDSS two-way timing error (remove orbit period).

It can be seen from Figure 8 the maximum spectrum is 1.01 days (about 24 h 14 m 24 s), which is consistent with the analysis results of the one-way timing error series, indicating that the two-way timing is still affected by the GEO satellite orbital error. It is worth noting that the drift is related to the speed of motion, which means that periodic fluctuations can

also appear in the two-way timing. When compared with Figure 7, it can be seen from Figure 9 that the characteristics of the GEO satellite orbit period have been eliminated, and there is a 5 ns timing error, which is a comprehensive reflection of the ionosphere and troposphere correction residual errors.

In order to further evaluate the performance under normal state, the availability and continuity results are given in Table 3. The availability and continuity of the RDSS two-way timing service are 100%.

**Table 3.** The availability and continuity of two-way timing service (normal state).


#### 3.3.2. Orbit Maneuver State

In order to comprehensively analyze the performance of BDS-3 two-way timing, we also analyzed and evaluated the two-way timing results under the orbit maneuvering state. The data was also from 8 February 2021 to 10 February 2021, where the maneuver period was from 11:30 to 13:30 on 8 February 2021.

It can be seen from Figure 10 that the two-way timing under orbit maneuver is consistent with that of the normal state, which means the GEO satellite ephemeris errors were eliminated by the two-way subtraction. Therefore, RDSS two-way timing accuracy has a certain impact on orbit residuals, but it is not obvious. Specifically, the two-way timing error of all epochs during orbit maneuver is smaller than 8 ns, and the RMS, mean value and standard deviation are 3.60 ns, 0.95 ns and 2.07 ns, respectively, which are very close with the results of the normal state. It is worth noting that whether it is maneuvering or not, the GEO satellite moves in the geostationary orbit, and the two-way timing still appears in the periodic fluctuation.

**Figure 10.** RDSS two-way timing accuracy (orbit maneuver state).

In order to further evaluate the performance of the orbit maneuver state, the availability and continuity results are given in Table 4. The availability and continui ty of theRDSS two-way timing service are 100%.

**Table 4.** The availability and continuity of two-way timing service (orbit maneuver state).


#### **4. Preliminary Conclusions**

This paper discusses the principles of one-way timing and two-way timing, and provides analysis and evaluation methods of timing accuracy, availability and continuity indexes. Using RDSS signal measurements, the performance of BDS-3 RDSS one-way timing and two-way timing were evaluated for the first time under both the normal state and orbit maneuver state. The preliminary conclusions are obtained as follows:


In conclusion, the preliminary evaluation results show that the timing performance of BDS-3 is better than the officially claimed one, which can provide technical reference for BDS-3 RDSS users in the normal state and orbit maneuver state.

**Author Contributions:** Conceptualization, R.G. and D.W.; methodology, Z.L.; software, N.X.; validation, T.Z. and H.R.; formal analysis, R.G.; resources, S.L.; data curation, Z.L.; writing—original draft preparation, R.G. and D.W.; writing—review and editing, D.W.; funding acquisition, R.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 41874043, 61603397, 41704037.

**Data Availability Statement:** The data sources are supported by Beijing Satellite Navigation Center.

**Acknowledgments:** Thanks for the data sources support of Beijing Satellite Navigation Center, and thanks for the fund support of the National Natural Science Foundation of China (Nos. 41874043, 61603397, 41704037).

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

#### **References**


**Mengjie Wu 1,2, Peng Guo 3,\*, Wei Zhou 4, Junchen Xue 1, Xingyuan Han 5, Yansong Meng <sup>5</sup> and Xiaogong Hu 1,2**


**Abstract:** 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 defined a new mapping function assuming the vertical ionospheric distribution as an exponential profiler 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.

**Keywords:** radio occultation; TEC; mapping function; plasmaspheric scale height; GCPM

## **1. 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–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–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-

**Citation:** Wu, M.; Guo, P.; Zhou, W.; Xue, J.; Han, X.; Meng, Y.; Hu, X. A New Mapping Function for Spaceborne TEC Conversion Based on the Plasmaspheric Scale Height. *Remote Sens.* **2021**, *13*, 4758. https:// doi.org/10.3390/rs13234758

Academic Editor: Chung-yen Kuo

Received: 27 October 2021 Accepted: 19 November 2021 Published: 24 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

gast (called F&K hereafter) [13], and some extended investigations [14–18]. The vertical structure of the ionosphere is usually assumed concentrating on one single thin layer. The height of the layer is called the ionospheric effective height (IEH) or shell height, which is calculated by the centroid method or integral median method or simplified as a constant value. The performance of mapping functions differs when receivers are located at different orbit heights or applying different IEH selections. How to choose the optimized IEH for the LEO-based TEC conversion needs systematical investigation. According to Zhong et al. [19] and Huang & Yuan [20], the IEH selection becomes more significant for the mapping function with an increasing zenith angle. Xiang & Gao [16] reviewed several mapping functions and proposed a mapping function that utilizes the ionospheric varying height assisted by the International Reference Ionosphere (IRI) model [21]. The study by Gulyaeva et al. [22] explored the center-of-mass of the ionosphere as the effective varying shell height produced by the ionospheric equivalent slab thickness. Zhong et al. [19] examined the applicability of three mapping functions for LEO-based GNSS observations: the thin layer model, the F&K model, and the Lear model, and found that the F&K geometric mapping function together with the IEH from the centroid method is more suitable to convert the LEO-based TEC measurements. The optimized IEH for the F&K can be approximately expressed as (2.5*hLEO* + 110) km under medium solar activity (MSA) conditions. As a consequence, the COSMIC-based TEC conversion (with satellites running at ~800 km) will get minimum mapping errors when the IEH is set to be around 2,100 km during MSA years. The current data processing in the COSMIC data analysis and archive center (CDAAC) applied the F&K mapping model, and the IEH is fixed at several hundreds or thousands of kilometers to simplify the uses [3,23]. However, one common drawback of these thin-layer mapping functions is that they ignore the vertical structure of the ionosphere and plasmasphere, and the IEH without spatial and temporal variations will affect the accuracy of the TEC conversion. One option to improve the TEC mapping is to estimate the 3D ionospheric electron density distribution by tomography or data assimilation based on the dense measuring network of sufficient spatial and temporal resolution and deduce realistic and persistent estimates of the ionospheric effective height [24]. Another approach is to assume multiple layers in the ionosphere implemented with specific mapping functions [17,20]. To further consider the vertical distribution of the ionosphere, Hoque & Jakowski [15] proposed a multi-layer mapping function according to the typical structure described by the Chapman layer. Nevertheless, the above methods are dependent on either extensive data coverage or accurate ionospheric models and parameters that are difficult to obtain in the global scope. Therefore, more advanced ionospheric mapping functions with simple free parameters deserve continuous attention and potential development.

The scale height in the plasmasphere (*HP*) is an important factor to present the dynamic nature and variations of the plasmasphere [25]. In Wu et al. [26], we have developed a scale height model depending on parameters, including the month of the year, local time, geographic latitude, and the solar radio flux index F10.7. This study is aimed to propose a new mapping function taking advantage of the former scale height model (named as the scale-height-based mapping function) to give an accurate obliquity factor to convert LEO-based slant TEC to vertical TEC and vice versa. The ionosphere is no longer collapsed into any shell but assumes an exponential vertical structure. The scale height introduces realistic spatial and temporal variations in the ionospheric and plasmaspheric electron density into the TEC conversion. With the proposed method, we expect to solve the problem of optimized IEH specification and neglect of ionospheric physical variations in the current single layer mapping models and, finally, facilitate the TEC and DCB estimations of LEO satellites. Section 2 introduces the mathematical solutions of the new mapping function, and Section 3 gives the distribution of plasmaspheric scale height and the associated mapping function values; in Section 4, global assessments of the scale-height-based mapping function and the F&K model are performed according to different seasons, locations, solar activity levels, and orbit altitudes, and the conclusion is drawn in Section 5.

#### **2. Materials and Methods**

The scale height is one of the most important characteristics in the ionosphere– plasmasphere system due to the role in shaping the electron density profile [25,27–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:

$$N\_{\varepsilon}(h) = N\_{\varepsilon}(h\_0) \cdot \exp\left(-\frac{h - h\_0}{H\_P}\right) \tag{1}$$

where *Ne*(*h*0) is the base electron density of plasmasphere (here set to be the electron density at the receiver height *h*0), and *HP* 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:

$$M(z) = \frac{STEC}{VTEC} = \frac{\int\_{\vec{r\_0}}^{\vec{r\_c}} N\_c \left(\vec{r}\right) d\vec{r}}{\int\_{h\_0}^{h\_c} N\_c(h) dh} \tag{2}$$

where <sup>→</sup> *r*<sup>0</sup> and <sup>→</sup> *rc* represent the locations of LEO and GNSS satellites from the Earth center, respectively; *h*<sup>0</sup> and *hc* 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

$$M(z) = \frac{\int\_{r\_0}^{r\_c} \exp\left(-\frac{r - r\_0}{H\_P}\right) \cdot \frac{r}{\sqrt{r^2 - r\_v^2}} dr}{\int\_{h\_0}^{h\_c} \exp\left(-\frac{h - h\_0}{H\_P}\right) dh} \tag{3}$$

where *rv* = *r*0·*sin*(*z*); *r*<sup>0</sup> 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 *HP* 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:

$$\begin{split} M(z) &= \frac{\sqrt{2r\_0/H\_P}}{\sin(z)} \exp\left(I^2\right) \cdot \frac{\sqrt{\pi}}{2} \operatorname{erfc}(I) \\ \operatorname{erfc}(I) &= \frac{2}{\sqrt{\pi}} \int\_I^\infty \exp\left(-y^2\right) dy \end{split} \tag{4}$$

where *I* = - *r*0 <sup>2</sup>*HP 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

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.

**Figure 1.** Geometric schematic of LEO-based GNSS TEC mapping model; LEO—Low Earth orbit.

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

$$M(z) = \frac{1 + R\_{shell} / R\_{orbit}}{\sqrt{(R\_{shell} / R\_{orbit})^2 - \sin^2(z)} + \cos(z)}\tag{5}$$

where *Rshell* and *Rorbit* 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.

#### **3. Results**

We calculate a mapping factor grid of zenith angle *z* varying between 5◦ and 80◦ with 5◦ resolution, and *HP* 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 *HP*. The grid map along with zenith angle and *HP* 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 *HP* 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 *HP* values (shown in panel (d)).

**Figure 2.** The mapping factor grid of zenith angle and the plasmaspheric scale height *HP*; panel (**a**,**c**) are the numerical and analytical results of scale-height-based mapping function, respectively; (**b**) is those for F&K model; (**d**) is the absolute relative difference between (**a**,**c**).

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 *HP* model is suitable to provide the essential parameter in the scale-heightbased mapping function. Besides the COSMIC-derived *HP*, 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*<sup>+</sup> *<sup>e</sup>* , 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 *HP* 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) could be approximately written as *TEC* ≈ *Ne*(*h*0)·*HP*. The *HP* retrieved from the observations and empirical model (represented by *HP*\_*COSMIC* and *HP*\_*GCPM*, respectively) are shown according to four seasons, local time, and dipole geomagnetic latitude in the low solar activity (LSA) year 2008 (Figure 3) and high solar activity (HSA) year 2013 (Figure 4), respectively. The four seasons are represented by March equinox (March, April, May), June solstice (June, July, August), September equinox (September, October, November), and December solstice (January, February, December) in this work. The *HP*\_*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

when applying *HP*\_*COSMIC* and *HP*\_*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 *HP* derived from satellite observations is assumed more reliable when dealing with the realistic mapping issues, but, in the following study, the *HP*\_*GCPM* is discussed because we are aiming at the comprehensive validation of the scale-height-based mapping function in a simulation way.

**Figure 3.** The distribution of the COSMIC and GCPM retrieved *HP* (represented by '*HP*\_*COSMIC*' and '*HP*\_*GCPM*') along with the local time, geomagnetic latitude, and season in the LSA year; panels (**a**,**c**,**e**,**g**) are for *HP*\_*COSMIC*, panels (**b**,**d**,**f**,**h**) are for *HP*\_*GCPM*.

**Figure 4.** Same as Figure 3 but for the HSA year.

#### **4. Discussion**

#### *4.1. 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 *VTECmod*). The relative root mean square (RMS) of each zenith angle *z* is calculated as the assessment criteria:

$$RMS(z) = \sqrt{\frac{\sum\_{i=1}^{N} \left[ \left( \frac{STEC\_i(z)}{M\_i(z)} - VTE\_{mod} \right) / VTE\_{mod} \right]^2}{N}} \times 100\% \tag{6}$$

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 '*HP*\_*N*' and '*HP*\_*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]:

$$IEH = (0.0027F\_{107} + 1.79)h\_{LEO} - 5.52F\_{107} + 1350\tag{7}$$

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 *HP* 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 *HP*-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.

**Figure 5.** The relative RMS statistics of the scale-height-based mapping function with numerical and analytical solutions (denoted as '*HP*\_*N*' and '*HP*\_*A*', respectively), and the F&K model ('F&K') in four seasons under low and high solar activity conditions (represented by 'LSA' and 'HSA'); panels (**a**–**d**) are for the LSA year, and (**e**–**h**) are for HSA year.

#### *4.2. 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 *HP*, 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◦.

**Figure 6.** The local time and latitude-dependent variations of the retrieved VTEC mean deviation for the mapping experiments in different seasons in LSA years. The zenith angle is fixed at 40◦. The '*HP*\_*N*', '*HP*\_*A*', and 'F&K' represent the numerical and analytical mapping functions and the F&K model, respectively; panels (**a**,**d**,**g**,**j**) are *HP*\_*N*-based mapping errors; (**b**,**e**,**h**,**k**) are *HP*\_*A*-based mapping errors; (**c**,**f**,**i**,**l**) are F&K mapping errors.

**Figure 7.** Same as Figure 6 but for HSA year.

Figures 8 and 9 demonstrate the mean deviation of the VTEC along with the latitude and longitude of three chosen zenith angles, 20◦, 40◦, and 60◦, in the December solstice. 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.

**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**–**c**)), 40◦ (panels (**d**–**f**)), and 60◦ (panels (**g**–**i**)), respectively. The '*HP*\_*N*', '*HP*\_*A*', and 'F&K' represent the numerical and analytical mapping functions and the F&K model, respectively.

#### *4.3. 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–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 11.** Same as Figure 5 but for LEO satellite at 1400 km.

To understand the distinctive performance of the scale-height-based mapping function refering to different LEO orbit height, it should be noted that the *HP* 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*<sup>+</sup> 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*<sup>+</sup> or *H*<sup>+</sup> *<sup>e</sup>* in the plasmasphere [26]. Therefore, the averaged *HP* 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 *HP*. 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 *HP* 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 *HP* 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.

#### **5. 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, *HP*, 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.

**Author Contributions:** Conceptualization, M.W. and J.X.; methodology, W.Z.; software, M.W. and J.X.; validation, X.H.(Xingyuan Han), Y.M.; formal analysis, P.G. and M.W. ; investigation, M.W.; resources, X.H.(Xiaogong Hu); data curation, X.H.(Xiaogong Hu); writing—original draft preparation, M.W.; writing—review and editing, P.G., W.Z. and X.H.(Xingyuan Han); supervision, P.G.; project administration, X.H.(Xiaogong Hu); funding acquisition, M.W. and Y.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China, grant number 2020YFA0713501; the National Natural Science Foundation of China, grant number 11903064, U1831116; the stability support fund project of national key laboratory, 2020SSFNKLSMT-06; the stability support fund project of Science and Industry Bureau, grant number HTKJ2020KL504006; the opening project of Shanghai Key Laboratory of Space Navigation and Positioning Techniques, NO. KFKT201906.

**Data Availability Statement:** The data presented in this study are openly available in ResearchGate at 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).

**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.

#### **Abbreviation List**


#### **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 \approx X \cos(z) + \frac{X^2 \sin^2(z)}{2r\_0} \tag{A1}$$

where *r*<sup>0</sup> 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,

$$\begin{split} STEC &= \int\_{0}^{\infty} N\_{0} \exp\left(-\frac{H}{H\_{\text{P}}}\right) dX = N\_{0} \cdot \int\_{0}^{\infty} \exp\left[-\left(\frac{X\cos(z)}{H\_{\text{P}}} + \frac{X^{2}\sin^{2}(z)}{2r\_{0}H\_{\text{P}}}\right)\right] dX \quad \text{(A2)}\\ \text{If } a &= \frac{\sin^{2}(z)}{2r\_{0}H\_{\text{P}}} \text{ and } b = \frac{\cos(z)}{2H\_{\text{P}}}, \text{ then} \end{split} \tag{A2}$$

$$\begin{split} STEC &= N\_{0} \cdot \int\_{0}^{\infty} \exp\left[-\left(2bX + aX^{2}\right)\right] dX \\ &= N\_{0} \exp\left(\frac{b^{2}}{a}\right) \cdot \int\_{0}^{\infty} \exp\left[-a\left(X + \frac{b}{a}\right)^{2}\right] dX \\ &= \frac{N\_{0}}{\sqrt{a}} \exp\left(\frac{b^{2}}{a}\right) \cdot \int\_{\sqrt{a}}^{\infty} \exp\left(-y^{2}\right) dy \end{split} \tag{A3}$$

*N*<sup>0</sup> is the electron density at the receiver. The vertical TEC is easy to obtain with integration as *VTEC* ≈ *N*0·*HP*; therefore, the analytic solution of *M*(*z*) relates to the scale height *HP*, zenith angle *z*, and the complementary error function (*er f c*(*I*)):

$$\begin{split} M(z) &= \frac{\sqrt{2r\_0/H\_P}}{\sin(z)} \exp\left(l^2\right) \cdot \frac{\sqrt{\pi}}{2} \text{erfc}(I) \\ &\text{erfc}(I) = \frac{2}{\sqrt{\pi}} \int\_I^\infty \exp\left(-y^2\right) dy \end{split} \tag{A4}$$

$$\text{where } I = \sqrt{\frac{r\_0}{2H\_P} \cot(z)}.$$

#### **References**

