*Article* **A New Method to Determine the Optimal Thin Layer Ionospheric Height and Its Application in the Polar Regions**

**Hu Jiang 1,2,3, Shuanggen Jin 1, Manuel Hernández-Pajares 4, Hui Xi 5,6, Jiachun An 2,\*, Zemin Wang 2, Xueyong Xu <sup>3</sup> and Houxuan Yan <sup>3</sup>**


**Abstract:** The conversion between the line-of-sight slant total electron content (STEC) and the vertical total electron content (VTEC) depends on the mapping function (MF) under the widely used thin layer ionospheric model. The thin layer ionospheric height (TLIH) is an essential parameter of the MF, which affects the accuracy of the conversion between the STEC and VTEC. Due to the influence of temporal and spatial variations of the ionosphere, the optimal TLIH is not constant over the globe, particularly in the polar regions. In this paper, a new method for determining the optimal TLIH is proposed, which compares the mapping function values (MFVs) from the MF at different given TLIHs with the "truth" mapping values from the UQRG global ionospheric maps (GIMs) and the differential TEC (dSTEC) method, namely the dSTEC- and GIM-based thin layer ionospheric height (dG-TLIH) techniques. The optimal TLIH is determined using the dG-TLIH method based on GNSS data over the Antarctic and Arctic. Furthermore, we analyze the relationship between the optimal TLIH derived from the dG-TLIH method and the height of maximum density of the F2 layer (hmF2) based on COSMIC data in the polar regions. According to the dG-TLIH method, the optimal TLIH is mainly distributed between 370 and 500 km over the Arctic and between 400 and 500 km over the Antarctic in a solar cycle. In the Arctic, the correlation coefficient between the hmF2 and optimal TLIH is 0.7, and the deviation between them is 162 km. Meanwhile, in the Antarctic, the correlation coefficient is 0.60, with a phase lag of ~3 months, with the hmF2 leading the optimal TLIH, and the deviation between them is 177 km.

**Keywords:** thin layer ionospheric height (TLIH); mapping function; dG-TLIH technique; global navigation satellite system (GNSS); height of maximum density of the F2 layer (hmF2)

#### **1. Introduction**

The ionosphere is an important part of the solar-terrestrial space environment, and it is closely related to human production and life. Due to the influence of free electrons in the ionosphere, global navigation satellite system (GNSS) signals are affected by the reflection, refraction, and delay of the ionosphere [1–3]. In order to improve the correction accuracy of ionospheric error for single-frequency users, it is necessary to obtain the temporal and spatial distribution of the ionospheric total electron content (TEC) [4].

**Citation:** Jiang, H.; Jin, S.; Hernández-Pajares, M.; Xi, H.; An, J.; Wang, Z.; Xu, X.; Yan, H. A New Method to Determine the Optimal Thin Layer Ionospheric Height and Its Application in the Polar Regions. *Remote Sens.* **2021**, *13*, 2458. https:// doi.org/10.3390/rs13132458

Academic Editor: Stefania Bonafoni

Received: 10 May 2021 Accepted: 21 June 2021 Published: 23 June 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/).

Most of the ionospheric correction information broadcasted by GNSS and their augmentation systems (such as WAAS and EGNOS) provides users with the vertical TEC (VTEC) based on two-dimensional ionospheric models [5–7]. GNSS users need to convert the VTEC to slant TEC (STEC) on the satellite-to-receiver line-of-sight (LOS) path, under the thin layer ionospheric model (TLIM) [7]. The conversion between VTEC and STEC is performed through a mapping function (MF) related to satellite elevation and thin layer ionospheric height (TLIH).

In order to realize the conversion between STEC and VTEC, the ionosphere is assumed to be an infinitely thin fixed-height layer. The intersection of the straight line from the satellite to the receiver and the thin layer is called the ionosphere piercing point (IPP) [7]. The STEC at the IPP represents the ionospheric TEC in the LOS path from the satellite to the receiver. Two-dimensional ionospheric models can be used to calculate the VTEC at the IPP. Then, the thin layer ionospheric model (TLIM) assumes a sample vertical gradient of ionospheric electron density, typically for ground-based receivers [8–10]. Nava et al. [11] detected the influence of electron density gradients for the thin layer model in the American sector using the coinciding pierce point (CPP) technique. In the calm period of the ionosphere, the TLIM can successfully simulate VTEC in the middle latitudes. However, in low-latitude regions or in the case of ionospheric storms, the TLIM can produce ionospheric mapping errors (IMEs) of up to dozens of TEC units (TECU) [12–14]. In order to reduce the impact of IMEs, many scholars have detected the optimal TLIH under different temporal and spatial conditions [11,15–17].

The TLIH is an important parameter of the mapping function, and it directly affects the conversion accuracy between STEC and VTEC [11]. It is usually selected at a ground height of between 350 km and 450 km [4,18]. Klobuchar [4] suggested that the global positioning system (GPS) broadcast ionospheric model uses a TLIH of 350 km. Widearea augmentation systems also establish ionospheric TEC grid models with a TLIH of 350 km [19]. The global ionosphere map (GIM) products provided by the International GNSS Service (IGS) Ionospheric Analysis Center use a TLIH of 450 km [18,20,21]. However, in comparison with the Chapman profile mapping function, the mapping function value based on a TLIH of 428.8 km is the closest to the true value on the global scale [7]. Li et al. [17] determined the optimized TLIH using combined IGS global ionospheric maps over China. Birch et al. [15] determined the TLIH by comparing the TEC from a pair of satellites observed simultaneously along slant and vertical paths over a ground station. However, due to the limitations of satellite orbit inclination and the accuracy of the GIM model [22], the aforementioned methods are not applicable in the polar regions. In addition, the height of maximum density of the F2 layer (hmF2) has also been used to model the ionospheric TEC, and it is superior to the constant TLIH [23,24].

The polar regions, namely the Arctic and the Antarctic, are the locations of the geographical and geomagnetic poles. The magnetic field lines over the polar caps penetrate the Earth or extend outward to connect with the interplanetary magnetic field, and the high-energy particles are mapped into the high-latitude ionosphere along magnetic field lines [25]. The polar ionosphere is controlled by the solar wind and the interplanetary magnetic field; meanwhile, it is coupled with the magnetosphere and thermosphere [26,27]. Many factors cause the complex characteristics of the ionosphere in the polar regions [28]. Moreover, the large-scale convection of the ionosphere in the polar regions also affects the ionospheric electron density [29]. The temporal and spatial variations of the ionosphere in the polar regions is obviously different from those at middle and low latitudes. However, there are still some problems that have not been systematically investigated over the Antarctic and Arctic, for example, whether the most widely used TLIH of 350–450 km is valid, and how to choose the range of optimal TLIH values for TEC conversion.

This paper presents a new method for detecting the optimal TLIH based on the dSTEC (differential STEC) measurements and UQRG (UPC quarter-of-an-hour rapid GIM) GIMs [6,30]. Meanwhile, we analyzed the reference range of the optimal TLIH values during almost one and a half solar cycles from 2003 to 2018 over the Antarctic and Arctic to mitigate ionospheric delay errors, using ionospheric modeling under the TLIM based on the new method, and the relationships between the optimal TLIH and the hmF2 were analyzed. To provide a detailed analysis of the variations of the optimal TLIH during a solar cycle, Constellation Observing System for Meteorology, ionosphere, climate (COSMIC), and GNSS data are used in this paper.

#### **2. Methods and Data**

#### *2.1. The CPP Method*

The conversion between VTEC and STEC is given by a mapping function related to satellite elevation and the TLIH, as shown in Equation (1).

$$\text{mf}(E) = \frac{\text{STEC}}{\text{VTEC}} \cong \frac{1}{\cos z} = \frac{1}{\sqrt{1 - \left(R \cos(E)/(R+H)\right)^2}}\tag{1}$$

where *m f*(•) denotes the mapping function, *R* is the geocentric distance of the receiver antenna, *H* is the thin layer ionospheric height, *z* is the satellite zenith distances at the IPP, and *E* is the satellite elevation at the observation point receiver antenna phase center (see Figure 1).

**Figure 1.** Geometry for estimating the ionospheric mapping errors (IMEs) [16].

The CPP method determines the optimal TLIH by analyzing the IMEs based on all pairs of CPP for a given epoch [11]. As illustrated in Figure 1, for practical purposes, if two pierce points, IPP1 and IPP2, satisfy Equation (2), they are considered to be a "coinciding" pierce point:

$$\begin{cases} |\varphi\_1 - \varphi\_2| < 0.2^\circ\\ |\lambda\_1/\cos(\varphi\_1) - \lambda\_2/\cos(\varphi\_2)| < 0.2^\circ \end{cases} \tag{2}$$

where (*ϕ*1, *λ*1) and (*ϕ*2, *λ*2) are the geomagnetic latitude and longitude of ionospheric IPP1 and IPP2, and the unit is degrees.

According to the thin layer approximation for the ionosphere, the converted VTEC values at IPP1 and IPP2 should be equivalent. However, due to the strong gradients of ionospheric electron density, conversion between the STEC and VTEC using the mapping function could result in TEC conversion errors related to the TLIH. The difference between the two converted VTEC values at IPP1 and IPP2 is defined as the IMEs, as shown in Equation (3). The optimalTLIH is defined as the one that minimizes the slant-to-vertical TEC conversion errors. In order to obtain enough CPP ionospheric TEC, dense GNSS monitoring stations are required to use this method.

$$
\Delta VTEC = \left| VTEC\_1 - VTEC\_2 \right| = \left| \frac{STEC\_1}{mf(Z\_1)} - \frac{STEC\_2}{mf(Z\_2)} \right| \tag{3}
$$

The STEC is derived from the GPS dual-frequency pseudo-range and carrier–phase measurements based on carrier-to-code leveling methods [7]. In addition, the differential code biases (DCB) products provided by the Center for Orbit Determination in Europe [31] were adopted to eliminate the GPS satellite DCBs. The receiver DCBs are estimated in conjunction with ionospheric model parameters using the generalized triangular series function [32].

#### *2.2. The dG-TLIH Method*

To overcome the limitations of the existing methods used in the polar regions, a new approach, based on the UQRG GIMs and dSTEC method [6], is proposed for optimal TLIH determination.

The GNSS dual-frequency carrier–phase observation equation can be expressed as follows:

$$\boldsymbol{\phi}\_{r,\boldsymbol{k}}^{s}(t) = \boldsymbol{\rho}\_{r}^{s}(t) + \boldsymbol{\varepsilon}(\delta\boldsymbol{t}\_{r}(t) - \delta\boldsymbol{t}^{s}(t)) - \kappa\_{k}\delta\boldsymbol{I}\_{r,1}^{s}(t) + \delta T\_{r}^{s}(t) + \boldsymbol{\varepsilon}(\boldsymbol{\gamma}\_{k}^{s} + \boldsymbol{\gamma}\_{r,k}) + \boldsymbol{\theta}\_{k}(t) + \lambda\_{k}N\_{r,k}^{s}(t) + \boldsymbol{\varepsilon}\_{\boldsymbol{\theta}} \tag{4}$$

where *s*, *r*, *k* denote the satellite PRN, receiver, and frequency, respectively; *φ<sup>s</sup> <sup>r</sup>*,*<sup>k</sup>* represents the carrier–phase measurement at frequency *fk*; *ρ<sup>s</sup> <sup>r</sup>* is the geometric range from satellite s to receiver r; *c* is the speed of light in a vacuum; *δtr* and *δt <sup>s</sup>* are the receiver and satellite clock errors; *<sup>κ</sup><sup>k</sup>* <sup>=</sup> *<sup>f</sup>* <sup>2</sup> 1 *f* 2 *k* is a constant factor; *δI<sup>s</sup> <sup>r</sup>*,1 <sup>=</sup> 40.3*STEC<sup>s</sup> r f* 2 1 denotes the ionospheric delays at frequency *f*1; *δT<sup>s</sup> <sup>r</sup>* denotes the tropospheric delays; *γ<sup>s</sup> <sup>k</sup>* and *γr*,*<sup>k</sup>* are the satellite and receiver– phase instrumental delays; *ϑ<sup>k</sup>* is the carrier–phase windup and is considered corrected; *λ<sup>k</sup>* denotes the wavelength at frequency *fk*; *N<sup>s</sup> <sup>r</sup>*,*<sup>k</sup>* is the carrier–phase ambiguity; and *εφ* is the noise and multipath effect of carrier–phase measurements.

According to Equation (4), the geometry-free combination of the carrier phase can be calculated as follows:

$$\begin{cases} \begin{array}{c} \phi\_{r, \rm GF}^{s}(t) = \phi\_{r,1}^{s}(t) - \phi\_{r,2}^{s}(t) = \boldsymbol{\upsilon} \cdot \boldsymbol{STEC\_{r}^{s}}(t) + \tilde{N}\_{r}^{s} + \varepsilon\_{\phi\_{\rm GF}}\\ \tilde{N}\_{r}^{s} = (\lambda\_{1}\tilde{N}\_{r,1}^{s} - \lambda\_{2}\tilde{N}\_{r,2}^{s}) + \boldsymbol{c}(\gamma\_{r,1} - \gamma\_{r,2}) + \boldsymbol{c}(\gamma\_{1}^{s} - \gamma\_{2}^{s}) \end{array} \tag{5}$$

where *<sup>υ</sup>* <sup>=</sup> 40.3 <sup>×</sup> <sup>10</sup><sup>16</sup> <sup>×</sup> (*<sup>f</sup>* <sup>−</sup><sup>2</sup> <sup>2</sup> <sup>−</sup> *<sup>f</sup>* <sup>−</sup><sup>2</sup> <sup>1</sup> ) is the conversion factor between the ionospheric delay (*δI<sup>s</sup> r*,1) and the ionospheric TEC (*STEC<sup>s</sup> <sup>r</sup>*). Along a phase-continuous satellite-receiver arc, as shown in Figure 2, there is an observation time (*tE*max ) with the maximum elevation. In theory, an observation value with higher elevation has lower errors, due to the much lower relevance of the mapping function. Based on Equation (5), in a continuous satellite– receiver arc of measurements, the difference in the given *φ<sup>s</sup> <sup>r</sup>*,*GF*(*t*) and the *<sup>φ</sup><sup>s</sup> <sup>r</sup>*,*GF*(*tE*max ) can be expressed as:

$$
\Delta S(t) = \phi^s\_{r, \zeta F}(t) - \phi^s\_{r, \zeta F}(t\_{E\_{\text{max}}}) = \upsilon (STEC^s\_r(t) - STEC^s\_r(t\_{E\_{\text{max}}})) \tag{6}
$$

where the definitions of the variables in Equation (6) are the same as in Equation (5).

According to Equation (6), the amount of STEC change between the given time *t* and the time *tE*max can be expressed as:

$$\Delta STEC(t) = \frac{1}{\upsilon} \Delta S(t) = STEC(t) - STEC(t\_{E\_{\text{max}}}) \tag{7}$$

where the definitions of the variables in Equation (7) are the same as that in Equation (6). Since the dSTEC is calculated from carrier–phase observations, its accuracy is less than 0.1 TECU.

According to Equation (7), the STEC at the given time *t* can be expressed as Equation (8) for each given pair of satellites and receivers, and for a common arc of measurements.

$$STEC(t) = STEC(t\_{E\_{max}}) + dSTEC(t) \tag{8}$$

Indeed, the "true" mapping value at the given time *t* can be expressed as:

$$MF\_T(t) = \frac{STEC(t)}{VTEC(t)} = \frac{VTEC(t\_{E\_{\text{max}}}) \cdot mf(E\_{\text{max}}) + dSTEC(t)}{VTEC(t)}\tag{9}$$

where *VTEC*(*tEmax* ) and *VTEC*(*t*) are the vertical TEC at the maximum elevation and time *t*, respectively, calculated using UQRG–GIM; and *m f*(*Emax*) is the mapping function value, calculated by Equation (1) at the maximum elevation. As can be seen from the equation, the "true" mapping value is independent of the TLIH.

**Figure 2.** Schematic diagram of dSTEC method.

The mapping function error indicator (*σH*) at a given TLIH can be expressed as

$$
\sigma\_H = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left( M F\_T(t\_i) - m f\_H(E\_{t\_i}) \right)^2} \tag{10}
$$

where *n* is the total number of samples; *MFT*(*ti*) denotes the mapping value calculated by Equation (9) at time *ti*; *Eti* is the satellite elevation; and *m fH*(·) denotes the mapping function with the TLIH of H.

According to Equation (10), the indicator *σ<sup>H</sup>* at different TLIHs can be calculated, and then the TLIH corresponding to the minimum *σ<sup>H</sup>* can be selected as the optimal TLIH. According to previous studies, it has been found that when the elevation angle is increased to approximately 40◦, the mapping errors caused by ionospheric gradients can be ignored. Then, the STEC calculated by Equation (8) is not affected by the ionospheric gradient. Therefore, for the present study, only a continuous satellite–receiver arc with a time length of more than 2 h and a maximum elevation angle greater than 50◦ were used.

In this paper, only the STEC at the maximum elevation angle in the continuous arc segment is calculated using the GIM model and the mapping function, to reduce the influence of the ionospheric gradients. Therefore, the reliability of this method depends on the accuracy of the GIM model used. In this study, the VTEC value was calculated by the UQRG–GIM model, which is provided by UPC, one of the IGS Ionosphere Working Group members. According to the comparison with the VTEC altimeter and globally distributed dSTEC GPS data, the RMS errors of the UQRG–GIM model are 2.0 TECU and 0.5 TECU, respectively [6]. In comparison with other GIM models, UQRG is the most accurate, and it has been used to detect ionospheric events in the polar regions [30].

#### *2.3. Observation Data*

In order to investigate the variation of optimal TLIH at high latitudes, datasets under different spatio-temporal conditions were used, including GPS data for the 2003–2018 period and COSMIC data for the 2007–2016 period. The time series of solar activity index F10.7 is shown in Figure 3.

**Figure 3.** Time series of the solar activity index F10.7 from 2003 to 2017.

GPS observations provided by the IGS and the Polar Earth Observing Network (POLENET) were used in this study. Figure 4 shows the geographical locations of the GPS ground stations over the Arctic and Antarctic, where the blue points are used for mapping ionospheric total electron content (TEC), the red circles are used for detecting the optimal thin layer ionospheric height (TLIH), and the green points are used for evaluating the ionospheric TEC model. For the present study, the sampling rate of the GPS measurements was set to 30 s. The maximum elevation of the GPS measurements was greater than 60◦ at latitude 75◦. Before 2010, only a few stations with non-uniform distribution were available in the Antarctic [28]. In addition, in order to show the relationship between the IMEs and the TLIH, global GPS ground observation data from the IGS (www.igs.org/station-resources/#site-guidelines, accessed date: 1 December 2019) were used.

**Figure 4.** Locations of GPS ground stations over the Arctic (left) and Antarctic (right), with the blue points mapping ionospheric total electron content (TEC), the green points detecting the optimal thin layer ionospheric height (TLIH), and the red circles evaluating the ionospheric TEC model.

To better understand how the hmF2 varies depending on solar activity and season in the polar regions, COSMIC data for the years 2007–2016, which were obtained from the COSMIC Data Analysis and Archive Center (CDAAC, http://cosmic-io.cosmic.ucar.edu/ cdaac/index.html, accessed date: 1 January 2018), are also considered.

#### **3. Results and Discussion**

#### *3.1. Optimal TLIH Determination Using CPP Method*

The CPP technique requires a dense GNSS monitoring network; otherwise, it is difficult to obtain GPS-derived TEC observations for the "coinciding" pierce points for a given epoch. To investigate the effects of the electron density gradients on the global ionosphere, GPS-derived STEC data for DOY 203, 300, and 355 in the year 2014 were used. The STEC data were obtained using a GPS satellite tracked globally by 460 IGS ground stations. The experiment was repeated for three different levels of solar activity: low (F10.7 = 92.6 sfu), medium (F10.7 = 187.8 sfu), and high (F10.7 = 205.8 sfu), corresponding to DOY 203, 300, and 355 in the year 2014, respectively. For the present study, the sampling interval of the data was set to 30 s, and only the STEC data corresponding to the lineof-sight of the receiver to the satellite with elevation angles between 15◦ and 40◦ at the observation points were used.

Values were computed for all epochs available during the analysis period. Examples of IMEs as a function of local time for different thin layer heights are illustrated in Figure 5. It can be seen that the solar activity had a significant impact on the IMEs: during low solar activity (DOY 203), most of the IMEs were between 0 and 8.0 TECU and there was no significant fluctuation; during medium solar activity (DOY 300), most of the IMEs were between 0 and 10 TECU, but the IMEs in the daytime were obviously greater than those in the nighttime, and some error values were higher than 20 TECU; during high solar activity (DOY 355), the IMEs increased significantly and large fluctuations appeared. In particular, during local times 00:00–04:00 and 13:00–16:00, the maximum IME reached up to 80 TECU, which was larger than at other times. Therefore, ignoring the ionospheric gradient when using TLIM will introduce significant errors caused by the conversion between STEC and VTEC. In addition, the standard deviations of the IMEs were smallest at the thin layer height of 400 km, during the selected days.

According to previous studies [15,16,33], it has been found that IMEs are related to latitude and TLIH. In the Northern Hemisphere, IMEs at high latitudes are greater than those at middle latitudes [16]. Although the TEC at high latitudes is lower than that at middle latitudes, due to the complex ionospheric variations in the polar regions, both the ionospheric gradient and the IMEs are larger than those at middle latitudes. In order to improve the conversion accuracy between the STEC and VTEC, it is necessary to deeply study the variation of the thin layer ionospheric height in the polar regions. All IGS tracking stations north of 55◦N were used in 30 s intervals. The mapping errors at the thin layer ionospheric height ranging from 300 to 700 km in steps of 50 km were calculated, and the thin layer height corresponding to the minimum error was selected as the optimal TLIH. Figure 6 presents the time series and histograms of the optimal TLIH for the Arctic region from 2003 to 2014, based on the CPP method. The optimal TLIH increased with the increase in solar activity. As illustrated in Figure 6, the highest proportion of the optimal TLIH was 400 km. The TLIH had a relatively low resolution of 50 km, mainly because this method was too computationally intensive and time-consuming. Meanwhile, the low resolution could introduce some errors, which could decrease the accuracy and reliability of the results.

Although the use of the CPP method to determine the TLIH has been investigated by many researchers in recent years, there are still some limitations with this method when it is used in the polar regions. First, this method depends on the dense GNSS monitoring network, but the available stations are few in the polar regions and their distributions are uneven. Therefore, it is difficult to get enough evenly distributed GNSS-derived TEC observations for the "coinciding" pierce points for a given epoch. Second, ionosphere

activity is intense, as ionospheric scintillation and patches often occur in the polar regions, both of which have greater effects on the statistical results of the CPP method.

**Figure 5.** IMEs as a function of local time for different layer heights. The data are related to the periods DOY 203, 300, and 355 in the year 2014.

**Figure 6.** Time series (top) and histograms (bottom) of daily optimal TLIH for the Arctic region from 2003 to 2014 based on the CPP method.

#### *3.2. Optimal TLIH Determination Using dG-TLIH Method*

Considering the limitations of the CPP method, a new approach for detecting the optimal TLIH is proposed. Figure 7 shows the time series and histograms of the daily optimal TLIH based on the dSTEC method and UQRG–GIM during a whole solar cycle: days 001, 2002 to 223, 2018 in the polar regions. In this experiment, the resolution of the TLIH was set to 1.0 km. As illustrated in Figure 7, the optimal TLIH typically varies from 300 to 600 km, and the difference between the maximum and the minimum is about 100 km in one year in the polar regions. In addition, solar activity has a considerable influence on the optimal TLIH, which is manifested as an upward optimal TLIH with the increase in solar activity, and which is similar to the optimal TLIH from the CPP method variation against solar activity shown in Figure 6. In the Arctic region, the optimal TLIH is mainly distributed between 370 and 500 km. During high solar activity years (2011–2014), the optimal TLIH mainly varies from 400 to 500 km, and during low solar activity years (2007–2009), the optimal TLIH mainly ranges from 300 to 500 km. In the Antarctic region, the distribution of the optimal TLIH has been more dispersed than that in the Arctic, especially before 2010. This can be attributed to the fact that the GNSS observation stations in the Antarctic are few and unevenly distributed, resulting in the low accuracy of the GIM model in this region, which affects the reliability of the mapping value. The time series of optimal TLIH can reflect its overall trend of change over the Antarctic. In addition, the optimal TLIHs show different seasonal characteristics: in the Arctic, the peak of optimal TLIH occurs from February to March and the valley appears from September to October, while in the Antarctic, the peaks appear from April to May and the valley appears from November to December. The seasonal variations of optimal TLIH are caused by many factors. In the summer, with the increase of local solar activity, the ionization of the middle and lower atmosphere increases, which leads to the downward movement of the optimal TLIH. In the winter, with the weakening of the local solar activity, the ionization of the middle and lower atmosphere weakens, resulting in the upward movement of the optimal TLIH. In addition, the polar ionosphere is affected by the magnetosphere and the interplanetary magnetic field, and the electron density also fluctuates.

**Figure 7.** Time series (left) and histograms (right) of daily optimal TLIH during a whole solar cycle: days 001, 2002 to 223, 2018 over the Arctic (top) and Antarctic (bottom).

To further analyze the variations of the optimal TLIH across different seasons and solar activity conditions, we computed the "true" mapping values based on Equation (9) and the mapping function values based on Equation (1) corresponding to the different TLIHs. Figures 8 and 9 show the variation in the mapping values computed by Equations (1) and (9) as a function of elevation angle for different TLIH under three different solar activity conditions: low (2009), medium (2016), and high (2014), for the months of March, June, September, and December over the Antarctic and Arctic, respectively. The "true" mapping values vary considerably across different seasons and years. During a low solar activity year (2009), there were significant differences in the mapping values in the different seasons, which is consistent with the characteristic of the optimal TLIH shown in Figure 7. When the elevation angles were greater than 50◦, the mapping values tended to be uniform. When the elevation angles are the same, the largest mapping value occurs in the spring and the smallest in the autumn in the Antarctic, as contrasted to the Arctic region. In addition, we also find that the mapping function values corresponding to different TLIHs were very different when the elevation angle was less than 15◦.

**Figure 8.** Variation of the mapping values computed by Equation (1) for different TLIHs and Equation (9) against elevation angles for three different solar activity conditions: low (2009), medium (2016), and high (2014), for the seasons of spring (September), summer (December), autumn (March), and winter (June) over the Antarctic.

Table 1 shows the optimal TLIH and the corresponding mapping function error indicator (*σH*) estimated by Equation (10) in March, June, September, and December during low (2009), medium (2016), and high (2014) solar activity years in the polar regions. The standard deviations at low solar activity are relatively larger than those at high solar activity, which indicates that the mapping values fluctuate greatly during low solar activity years. The difference between the maximum optimal TLIH in March and the minimum in September is more than 120 km. There are two main causes for the large fluctuations of the mapping values during low solar activity. First, the ionospheric TEC is smaller during low solar activity, and a slight change in the TEC will cause a great change in the mapping value. Second, there are few available GNSS tracking stations in Antarctica (especially before 2010) and the stations are unevenly distributed, resulting in low accuracy of the GIM model, which affects the reliability of the mapping value. At both medium and high solar activities, the variation of optimal TLIH in one year was only about 50 km. On the same date, the optimal TLIH in the Antarctic was higher than in the Arctic.

**Figure 9.** Variation of the mapping values computed by Equation (1) for different TLIHs and Equation (9) against elevation angles for three different solar activity conditions: low (2009), medium (2016), and high (2014), for the months of spring (March), summer (June), autumn (September), and winter (December) over the Arctic.

**Table 1.** The optimal TLIHs (km) and the corresponding standard deviations in March, June, September, and December of 2009, 2014, and 2016, respectively, over the Antarctic and Arctic.


The ionospheric gradient does not have a significant influence on the thin layer model with an elevation angle above 50◦. Therefore, only cases with an elevation angle lower than 56◦ were considered in the experiment. Figure 10 shows the variations of optimal TLIH against the day of year (DOY) and elevation angle in low (2009), medium (2016), and high (2014) solar activity years in the polar regions. It can be seen that the optimal TLIH increases with the increase in elevation angle. At low solar activity, the optimal TLIH with elevation between 25◦ and 55◦ in the Antarctic was significantly higher than that in the Arctic. The maximum value of the optimal TLIH in one year occurs in local winter, and the minimum occurs in local summer.

**Figure 10.** Variations of the optimal thin layer ionospheric height (TLIH) for three different solar activity conditions: low (2009), medium (2016), and high (2014), over the Antarctic (left) and Arctic (right), and as a function of the day of year (DOY) and elevation angle.

In order to verify the reliability of the dG-TLIH method, we used three kinds of TLIH (350 km, 450 km, and optimal TLIH from the dG-TLIH method) to build the ionosphere TEC model in the polar regions. The spherical cap harmonic (SCH) model was used to construct the ionospheric TEC model. The specific parameters of the SCH model are shown in previous research [28]. The numerical experiment was repeated for three different levels of solar activity: low (DOY 136, 2017), medium (DOY 127, 2014), and high (DOY 46, 2014), as detailed in Table 2.


**Table 2.** The information of solar activity and optimal TLIH for the 3 days selected to verify the reliability of the optimal TLIH over the Antarctic and Arctic.

In addition, we compared the performances of the three models based on different TLIHs for mapping the polar ionospheric TEC using dSTEC method [6]. Table 3 presents the mapping bias and RMS errors of the three models for the Arctic and Antarctic regions, respectively. The statistical results were based on the dSTEC measurements provided by the verification stations shown in Figure 3. When the elevation angle of the satellite was greater than 40◦, the influence of ionospheric gradient on the thin layer ionospheric model was weak. Therefore, statistical results only considered the data in the range of 10◦ to 40◦ for the elevation angle. It can be seen from the table that the performance of the model based on the optimal TLIH was better than the model based on 350 and 450 km, especially

in the case of high solar activity. By statistical analysis of the average of bias and RMS error on these three days, it can be found that in the Antarctic, the model bias error from the optimal TLIH was better than 1.21 and 0.20 TECU for 350 and 450 km, respectively, and the RMS error was better than 0.79 and 0.06 TECU, respectively; meanwhile in the Arctic, the bias error was better than 1.22 and 0.25 TECU for 350 and 450 km, respectively, and the RMS error was better than 0.90 and 0.20 TECU, respectively.

**Table 3.** Statistical results of ionospheric TEC models for different TLIHs using dSTEC measurements in the Antarctic and Arctic.


#### *3.3. Relationship between hmF2 and Optimal TLIH*

According to previous research [17,23,24], hmF2 is related to the TLIH, but it is typically lower due to the asymmetry in the electron density profiles at GNSS transmitter heights, and it can reflect the characteristics of the TLIH to a certain extent. In general, the centroid height of the ionospheric electron density was used as the TLIH. In practical applications, the centroid height of the ionospheric electron density is equal to the hmF2 plus 80 km [23]. Therefore, as an indicator of the TLIH, we analyzed the relationship of the hmF2 from COSMIC data and the optimal TLIH from the dG-TLIH method in the Arctic and Antarctic.

Figure 11 shows the time series of the daily optimal TLIH from the dG-TLIH method, the hmF2 from the COSMIC data, and the solar activity index F10.7, from 2007 to 2016 in the Arctic. In nearly one solar cycle, the daily mean hmF2 ranges from 200 to 330 km. It can be seen that the hmF2 is highly correlated with optimal TLIH, and increases with the increase in the solar activity. The correlation coefficient between these two variables is up to 0.70 in the Arctic.

Figure 12 shows the time series of the daily optimal TLIH from the dG-TLIH method, the hmF2 from the COSMIC data, and the solar activity index F10.7, from 2007 to 2016 in the Antarctic. It can be seen that the time series of the optimal TLIHs are relatively dispersed. In order to better analyze the change characteristics of the optimal TLIH, we use a Fourier series to fit it, as shown by the green line in Figure 12. Although both the optimal TLIH and the hmF2 have an upward trend with the increase in solar activity, the change trend of the two variables in one year is opposite, especially in the high solar activity year. Therefore, we measured the relationship between the optimal TLIH and the hmF2 by

the maximum correlation coefficient and corresponding phase discrepancy, as shown in Figure 13. The optimal TLIH positively relates to the hmF2, with a correlation coefficient of 0.60 and a phase lag of ~3 months, with the hmF2 leading optimal TLIH.

**Figure 11.** Time series of optimal TLIH (blue line), hmF2 (red line), and the levels of solar activity, F10.7 (bottom), from 2007 to 2016 over the Arctic.

**Figure 12.** Time series of optimal TLIH (blue dots), fitting TLIH (green line), hmF2 (red line), and levels of solar activity, F10.7 (bottom), from 2007 to 2016 over the Antarctic.

In the Arctic, the average values of optimal TLIH and hmF2 were 420.43 and 258.43 km during 2007–2016, respectively; meanwhile in the Antarctic, they were 441.87 and 264.60 km, respectively. The differences between the optimal TLIH and the hmF2 were 162 and 177.27 km over the Arctic and Antarctic, respectively, which is greater than the empirical

value of 80 km. The main reason for this lies in the existence of the plasmasphere outside the ionosphere, and the TEC observed by GNSS will inevitably reflect the contribution of the plasmasphere [23]. The contribution of the plasmasphere to the vertical TEC is about 10~20% during the day and up to 50% at night [34].

**Figure 13.** The correlation between the optimal TLIH and the hmF2; positive lag means hmF2 leads, and vice versa.

#### **4. Conclusions**

In this paper, a new method (dG-TLIH) to determine optimal TLIH was proposed, and the optimal TLIH was determined using the dG-TLIH method and the CPP technique over the Antarctic and Arctic based on GNSS data during a period of about one solar cycle. Relationships between the optimal TLIH from the dG-TLIH method and the hmF2 from the COSMIC data in the polar regions were analyzed.

According to the tests of the CPP technique, the results indicate that the optimal TLIH ranges from 300 to 600 km, and the height of 400 km ws the most frequent TLIH during 2003–2014 in the Arctic. Because the method requires dense and uniform ground tracking stations, it is not suitable for use in the Antarctic. When compared to the CPP method, the dG-TLIH method used in this study overcomes the limitation of requiring dense observation stations. According to the test of the dG-TLIH method, it is also confirmed that the optimal TLIH is related to solar activity, which mostly ranges from 320 to 500 km during one solar cycle. The optimal TLIH values present obviously seasonal variation characteristics, and most of the maximum and minimum values occurred in local winter and summer, respectively. Compared with the fixed TLIH (350 and 450 km), the bias errors of the ionospheric TEC model based on the optimal TLIH are decreased by 0.05 to 3.0 TECU, and the RMS errors are decreased by 0.02 to 2.11 TECU under different solar activity. According to the analysis of the relationship between the optimal TLIH and the hmF2, it was found that in the Arctic, the correlation coefficient between the hmF2 and optimal TLIH was 0.7, and the deviation between them was 162 km. In the Antarctic, they had a correlation coefficient of 0.60 with a phase lag of ~3 months, with the hmF2 leading the optimal TLIH, and the deviation between them was 177 km.

In order to facilitate ionospheric TEC modeling, we give an optimal TLIH each day over the Antarctic or Arctic. However, because the variation of ionospheric electron density is affected by solar activity, latitude, season, time, and other factors, it leads to the complex spatiotemporal variation of the optimal TLIH. In addition, the shape of the earth also needs to be taken into account when determining the optimal TLIH at the global scale. These problems will be the focus of future study.

**Author Contributions:** Conceptualization, H.J. and J.A.; methodology, M.H.-P.; data curation, H.J. and H.X.; validation, H.X. and H.J.; formal analysis, S.J. and M.H.-P.; writing—original draft, H.J. and H.X.; writing—review and editing, H.J., S.J., M.H.-P., J.A., Z.W., X.X. and H.Y. 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 41776195, 41941010, 41531069; and the State Key Laboratory of Geodesy and Earth's Dynamics, Innovation Academy for Precision Measurement Science and Technology, CAS, grant number SKLGED2021-2-3.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank the following organizations and projects for providing the data used in this work: International GNSS Services (IGS, accessed date: 1 December 2019), the Polar Earth Observing Network (POLENET, accessed date: 1 December 2019), and the COSMIC Data Analysis and Archive Center (CDAAC, accessed date: 1 January 2018).

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

#### **References**

