Modelling Global Ionosphere Based on Multi-Frequency, Multi-Constellation GNSS Observations and IRI Model

With the emergence of BeiDou and Galileo as well as the modernization of GPS and GLONASS, more available satellites and signals enhance the capability of Global Navigation Satellite Systems (GNSS) to monitor the ionosphere. However, currently the International GNSS Service (IGS) Ionosphere Associate Analysis Centers (IAACs) just use GPS and GLONASS dual-frequency observations in ionosphere estimation. To better determine the global ionosphere, we used multi-frequency, multi-constellation GNSS observations and a priori International Reference Ionosphere (IRI) to model the ionosphere. The newly estimated ionosphere was represented by a spherical harmonic expansion function with degree and order of 15 in a solar-geomagnetic frame. By collecting more than 300 stations with a global distribution, we processed and analysed two years of data. The estimated ionospheric results were compared with those of IAACs, and the averaged Root Mean Squares (RMS) of Total Electron Content (TEC) differences for different solutions did not exceed 3 TEC Unit (TECU). Through validation by satellite altimetry, it was suggested that the newly established ionosphere had a higher precision than the IGS products. Moreover, compared with IGS ionospheric products, the newly established ionosphere showed a more accurate response to the ionosphere disturbances during the geomagnetic storms.


Introduction
The ionosphere constitutes the upper part of the atmosphere, extending from approximately 60 to 1500 km above the Earth's surface, enriched with free electrons and ions. It is a particularly important region with regard to radio signal propagation and radio communications in general. The ionosphere is dispersive, which means that the time delay attributed to the ionosphere depends on the frequency of the signal. This dispersive property also makes the multi-frequency Global Navigation Satellite Systems (GNSS) efficiently sense the ionospheric activities and ionisation [1][2][3]. In the meantime, with the modernization of GPS and GLONASS as well as the rapid deployments of BeiDou System

Ionospheric Observation Equation
The unambiguous GNSS code measurements are widely used in ionosphere estimation. To significantly reduce the noise of code measurements, the so-called "phase-smoothed code" algorithm is used based on a continuous time series of code and phase measurements [5]. Assuming a common "cycle-slip-free" period of n epochs, the smoothed code measurements from receiver a to satellite i at epoch t are written as: where m denotes another frequency index except for the reference frequency index 1, f is the carrier frequency (Hz), P(t) and L(t) are the raw code and phase measurements (m) at epoch t, P and L denotes the averaged code and phase measurements (m) over the n epochs.
The undifferenced observation equations of smoothed code measurements at epoch t can be reformulated as follows: where ρ means the geometric distance between the antennas of receiver and satellite (m); c is the speed of light in vacuum (m/s); dt i and dt a represents the satellite and receiver clock bias (s); T depicts tropospheric delay (m); E is the Slant TEC (STEC) in the signal propagation path in unit of TEC Unit (1 TECU = 10 16 el/m 2 ); δh and δb are the hardware biases of satellite and receiver (m). A priori noise of the unsmoothed raw code measurement is set as 0.6 m, ε is the observation noise of the smoothed code measurement (m), which is reduced by about a factor of √ n with respect to the a priori noise. The ionospheric observation can be obtained through differencing the smoothed measurements at different frequencies [5]: P i a,m1 (t) = P i a,m (t) − P i a,1 (t) = 40. 3( where ε ion = √ 2 × ε, δh i m1 and δb a,m1 denote the satellite and receiver DCB between frequencies m and 1.

Initial Ionospheric Information from the International Reference Ionosphere (IRI)
The IRI (http://irimodel.org/) is an international project sponsored by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI), the latest version is IRI-2016 [15,16]. As an advanced empirical ionospheric model, the IRI-2016 can integrate the TEC from the ground to a selected height for a given location, time and date. A priori ionospheric information helps us to strengthen the stability and reliability of the solution, especially for the regions with less available GNSS data. Based on IRI-2016, the hourly GIM was pre-generated with a spatial resolution of 5 • × 2.5 • in longitude and latitude. The a priori ionospheric value for each grid is seen as a pseudo-observation and inserted into the normal equation. The pseudo-observation at a grid is formatted as: E IRI (B, L) = VTEC IRI (B, L) + ε(E IRI ) (5) Remote Sens. 2020, 12, 439 4 of 19 in which B, L denote the latitude and longitude of the grid, VTEC IRI is the corresponding Vertical TEC (VTEC) derived from IRI-2016, ε(E IRI ) represents the a priori precision of the VTEC. The precision of VTEC derived from the IRI model is around 2~5 TECU in medium and high latitudes and 5~10 TECU in low latitudes, which was deeply analyzed and well validated by various authors [19][20][21]. To avoid the absolute VTEC estimation, a relatively loose constraint of 10 TECU was employed for the region of latitudes lower than 25 • and 5 TECU for the region of latitudes higher than 25 • .

Global Ionospheric Representation
The single-layer model at an altitude of 450 km and SH expansion function are well suited to represent the global VTEC [21,22]. The VTEC was parameterised as: N nm P nm (sin β)( C nm cos(ms) + S nm sin(ms)) (6) where E v is the VTEC (TECU); β and s denote the sun-fixed geomagnetic latitude and longitude of the ionospheric pierce point (IPP); n max is the maximum degree of the SH expansion; N nm is the normalization function [23]; P nm is the classical, unnormalized Legendre function; C nm and S nm are the unknown SH coefficients. For absolute TEC estimation using ground-based GNSS data, the TEC along the vertical direction is of main interest. To convert the STEC to VTEC, the Modified Single-Layer Model (MSLM) mapping function was employed: where z and z are the zenith distances (i.e., the angle measured from zenith) of the satellite at the receiver and IPP, respectively. The best MSLM mapping function was achieved at the altitude H = 506.7 km and the scale factor α = 0.9782 [24], when using the mean earth radius R = 6371 km and assuming a maximum zenith distance of 80 degrees. Substituting (6) and (7) to (4), the multi-constellation GNSS ionospheric observations for global ionosphere determination was described as: with the initial ionospheric constraints from IRI-2016 for each grid of the hourly GIM: where B is from the latitude −87.5 • to 87.5 • with an interval of 2.5 • ; L is from the longitude −180 • to 180 • with an interval of 5 • ; i, j, k, l represents GPS (G), GLONASS (R), BDS (C) and Galileo (E) satellites, respectively. m, n, p, q are frequency indexes except for the reference frequency index 1, such as m = L2, L5 for GPS, n = G2 for GLONASS, p = B2, B3 for BDS and q = E5a, E5b, E5, E6 for Galileo. Considering the GLONSS inter-frequency bias [25] and BDS satellite-induced code bias [26], we adopted a scale factor of 1.5 for the noise of GLONASS and BDS ionospheric observations, which means ε R ion = ε C ion = 1.5ε ion and ε G ion = ε E ion = ε ion . In addition, the DCB parameters of the satellite and receiver for each observation in (8) are strongly correlated, causing singularity of the normal equation. Additional constraints need to be applied to separate the DCB of the receiver and satellites. According to the processing strategy of the IAACs, we assumed that the sum of satellite DCB parameters in one day for each constellation was equal to zero: The multi-constellation GNSS ionospheric observations are lumped to a normal equation epoch by epoch. The global ionospheric parameters are modelled as piece-wise linear function with an interval of 1 h and solved by the least square adjustment. Based on the estimated SH parameters and variance-covariance matrix, the hourly GIM and RMS maps with a resolution of 5 • × 2.5 • in longitude and latitude were generated through (6) and stored in IONEX format.

Data Sources
The Multi-GNSS Experiment (MGEX) was initialized in 2013, aiming to track, collate and analyse all available GNSS signals [27,28]. By January 2020, the MGEX network had grown to more than 300 stations and a continued expansion is expected in the future. Especially in the past two years, with the rapid deployment of BDS and Galileo, the MGEX stations have significantly increased. To achieve a meaningful conclusion, we collected data from 348 IGS and MGEX stations in 2017 and 2018, as illustrated in Figure 1. All the stations were able to track GPS signals, whilst 295, 208 and 141 stations could normally receive GLONSS, Galileo and BDS signals, respectively. Although we used 348 stations in total to model the global ionosphere, the number of stations we applied for each day was around 300, which was close to the number of stations used by the IAACs.
The Positioning and Navigation Data Analysis (PANDA) platform at Wuhan University in China [29] was extended to model the global ionosphere based on multi-frequency and multi-constellation GNSS, and the data of two years were processed with this software. Moreover, to deeply analyse, comprehensively assess and independently validate our global ionosphere products, the determined GIM was compared with the VTEC from satellite altimetry mission Jason-2/3. The final data products of satellite altimetry, issued by the French Space Agency-"Centre National d'Etudes Spatiales" (CNES), were collected from the database of the Archiving, Validation and Interpretation of Satellite Oceanographic (AVISO).
To have a clear sense of the increase of IPP brought by multi-constellation GNSS, as an example, Figure 2 plots the IPP at an epoch on the day of year (DOY) 001 in 2018. The IPP are noticeably increased since the BDS and Galileo satellites are introduced. Actually, the increase in observations is far more than the increase of IPP, because BDS and Galileo can transmit triple or more frequency signals. Each IPP of BDS or Galileo satellite in Figure 2 denotes at least two available ionospheric observations. Therefore, the observations for ionosphere modelling are obviously increased through introducing multi-frequency signals and multi-constellation Galileo and BDS satellites. Unfortunately, the IPP are not evenly distributed in space. The IPP are dense over the regions where there are many GNSS stations, and are scattered for the regions with limited GNSS stations. The VTEC values derived by IRI-2016 can somewhat make up for this deficiency and provide some effective pseudo-observations especially for the regions with less or no available GNSS observations.

Results and Analysis
For convenience, we defined the GIM derived by PANDA, IGS, CODE, ESA and UPC as PANG, IGSG, CODG, ESAG and UPCG, respectively. The summary of processing strategies for different types of GIMs are listed in Table 1. The determined PANG and RMS maps are first showed in this section, followed by the comparison analysis with IGSG, CODG, ESAG and UPCG products. Then, satellite altimetry was used to assess the precision of the GIM over the two years. Finally, the

Results and Analysis
For convenience, we defined the GIM derived by PANDA, IGS, CODE, ESA and UPC as PANG, IGSG, CODG, ESAG and UPCG, respectively. The summary of processing strategies for different types of GIMs are listed in Table 1. The determined PANG and RMS maps are first showed in this section, followed by the comparison analysis with IGSG, CODG, ESAG and UPCG products. Then, satellite altimetry was used to assess the precision of the GIM over the two years. Finally, the

Results and Analysis
For convenience, we defined the GIM derived by PANDA, IGS, CODE, ESA and UPC as PANG, IGSG, CODG, ESAG and UPCG, respectively. The summary of processing strategies for different types of GIMs are listed in Table 1. The determined PANG and RMS maps are first showed in this section, followed by the comparison analysis with IGSG, CODG, ESAG and UPCG products. Then, satellite altimetry was used to assess the precision of the GIM over the two years. Finally, the performance of the estimated ionosphere during two geomagnetic storms are examined and analysed.

Comparison with IAACs' Solutions
As an example, Figure 3 shows PANG, IGSG, CODG and their RMS snapshots at 20:00 UT for DOY 001, 2018. The VTEC results from different solutions had a similar distribution. The ionosphere conditions are in low and high ionisations at night and day time, respectively. The VTEC results ranged from 0 to 30 TECU and exhibited an apparent equator anomaly. The "peak-to-peak" equator anomaly was well aligned with the geomagnetic equator. Slight TEC differences among different solutions can be seen around the equator region. As for the RMS maps, we can see an apparently higher RMS for CODG over the regions with less GNSS data, and its RMS even exceeded 3 TECU. PANG introduced more available satellites, signals and a priori ionosphere model, and its RMS values were lower and achieved a better global consistency. IGSG is a combined solution from different IAACs, and its RMS map shows an irregular distribution. In brief, the solar-geomagnetic reference frame is well adopted to represent the global ionosphere, and the RMS is remarkably reduced after introducing additional BDS, Galileo signals and IRI. PANG, IGSG, CODG, ESAG, JPLG and UPCG have 25 1-hourly, 13 2-hourly, 25 1-hourly, 13 2hourly, 13 2-hourly and 13 2-houtly GIM snapshots for each day. One snapshot has 73 71=5183 × grids. To quantify the VTEC differences among different solutions, we compared the VTEC between the IAAC and PANG. The daily RMS of the VTEC differences for each IAAC's snapshot with respect to PANG are plotted in Figure 4, and the corresponding statistical values of the daily RMS are listed in Table 2. CODG had the lowest RMS values. This is because PANG and CODG adopted the same mathematical model and temporal resolution (1 hr). The mean values of daily RMS for IGSG and CODG were no more than 1.5 TECU, and the daily RMS for different solutions varied from 0.6 to 6.6 TECU in 2017 and 2018. The ionosphere in UPCG is independently modelled for each station with a 3D voxel model, not the SH expansion function. The different ionosphere representation models are possibly responsible for the larger RMS values between UPCG and PANG.  PANG, IGSG, CODG, ESAG, JPLG and UPCG have 25 1-hourly, 13 2-hourly, 25 1-hourly, 13 2-hourly, 13 2-hourly and 13 2-houtly GIM snapshots for each day. One snapshot has 73 × 71 = 5183 grids. To quantify the VTEC differences among different solutions, we compared the VTEC between the IAAC and PANG. The daily RMS of the VTEC differences for each IAAC's snapshot with respect to PANG are plotted in Figure 4, and the corresponding statistical values of the daily RMS are listed in Table 2. CODG had the lowest RMS values. This is because PANG and CODG adopted the same mathematical model and temporal resolution (1 h). The mean values of daily RMS for IGSG and CODG were no more than 1.5 TECU, and the daily RMS for different solutions varied from 0.6 to 6.6 TECU in 2017 and 2018. The ionosphere in UPCG is independently modelled for each station with a 3D voxel model, not the SH expansion function. The different ionosphere representation models are possibly responsible for the larger RMS values between UPCG and PANG. In addition, the solar and geomagnetic activity mainly modulates and controls the ionospheric evolution. The solar radio flux at 10.7 cm (2800 MHz), often called the F10.7 index, is an excellent indicator of solar activity. The Ap-index provides a daily average level for geomagnetic activity. The RMS for different solutions fluctuates with the trend of F10.7/AP index, which can be clearly seen from January to April and from September to November in 2017. The maximum RMS between UPCG and PANG was 6.6 TECU, which occurred on DOY 311 of 2017, in which day the maximum Ap-index could reach up to 94, reflecting a relative high level of geomagnetic activity. During 2018, the F10.7 index was always lower than 85, and the sun was in a low activity condition, inducing a relatively stable ionosphere environment. The VTEC among different solutions showed a good consistency with different solar and geomagnetic conditions, and the mean values of the RMS were all lower than 3 TECU.
Remote Sens. 2020, 12, 439 9 of 20 In addition, the solar and geomagnetic activity mainly modulates and controls the ionospheric evolution. The solar radio flux at 10.7 cm (2800 MHz), often called the F10.7 index, is an excellent indicator of solar activity. The Ap-index provides a daily average level for geomagnetic activity. The RMS for different solutions fluctuates with the trend of F10.7/AP index, which can be clearly seen from January to April and from September to November in 2017. The maximum RMS between UPCG and PANG was 6.6 TECU, which occurred on DOY 311 of 2017, in which day the maximum Ap-index could reach up to 94, reflecting a relative high level of geomagnetic activity. During 2018, the F10.7 index was always lower than 85, and the sun was in a low activity condition, inducing a relatively stable ionosphere environment. The VTEC among different solutions showed a good consistency with different solar and geomagnetic conditions, and the mean values of the RMS were all lower than 3 TECU.

Satellite Altimetry Validation
The IAACs' GIMs are generated from GNSS technology. The comparison analysis in the previous section just reflects an internal consistency or relative precision within the same technology, which cannot validate the absolute accuracy of the GIM. Satellite altimetry measures the vertical distance between a satellite and the ocean surface [30]. The range is derived from the traveling time of the radar impulse transmitted by the on-board radar altimeter and reflected from the ocean surface. The two widely separated frequencies on Jason3 allow TEC to be detected directly and precisely from nadir altimetry sampling data along the satellite track [31]. The satellite altimetry technology is one of the best ways to validate the ionospheric accuracy.
Considering IGSG is a weighted combination from different IAACs' solutions and has a relatively higher precision, we compared the VTEC from PANG and IGSG with that of Jason3. As an example, Figure 5 plots the VTEC of Jason3, PANG and IGSG from 11:00 to 14:00 in DOY 180, 360 of 2017 and 2018. The sampling rate of Jason3's data is 1 Hz, the PANG or IGSG VTEC at each track point of Jason3 was interpolated between two consecutively rotated GIM. The gaps existing in Figure  5 mean Jason3 was crossing the land or the satellite altimeter data were flagged as bad. The PANG and IGSG VTEC follow the VTEC variation of Jaons3. The VTEC could reach 30 TECU for DOY 360

Satellite Altimetry Validation
The IAACs' GIMs are generated from GNSS technology. The comparison analysis in the previous section just reflects an internal consistency or relative precision within the same technology, which cannot validate the absolute accuracy of the GIM. Satellite altimetry measures the vertical distance between a satellite and the ocean surface [30]. The range is derived from the traveling time of the radar impulse transmitted by the on-board radar altimeter and reflected from the ocean surface. The two widely separated frequencies on Jason3 allow TEC to be detected directly and precisely from nadir altimetry sampling data along the satellite track [31]. The satellite altimetry technology is one of the best ways to validate the ionospheric accuracy.
Considering IGSG is a weighted combination from different IAACs' solutions and has a relatively higher precision, we compared the VTEC from PANG and IGSG with that of Jason3. As an example, Figure 5 plots the VTEC of Jason3, PANG and IGSG from 11:00 to 14:00 in DOY 180, 360 of 2017 and 2018. The sampling rate of Jason3's data is 1 Hz, the PANG or IGSG VTEC at each track point of Jason3 was interpolated between two consecutively rotated GIM. The gaps existing in Figure 5 mean Jason3 was crossing the land or the satellite altimeter data were flagged as bad. The PANG and IGSG VTEC follow the VTEC variation of Jaons3. The VTEC could reach 30 TECU for DOY 360 in 2017, and did not exceed 20 TECU for DOY 180 in 2018 because of a low annual solar activity condition. Compared with IGSG, the VTEC of PANG appeared to have a better consistency with that of Jason3. This could be seen through the calculated VTEC Standard Deviation (STD) of IGSG and PANG. PANG showed a lower STD, which had improved by 0.2~0.6 TECU compared with the STD of IGSG. The daily system offset between satellite altimetry and GNSS mainly includes both the instrumental bias and the plasmaspheric component (i.e., VTEC above the height of the altimetry mission) [32]. Figure 6 clearly exhibits the daily offset and STD of VTEC for PANG-Jason3 (PANG minus Jason3) and IGSG-Jason3 (IGSG minus Jason3). The mean values of the VTEC offset and STD for PANG-Jaons3 and IGSG-Jason3 are listed in Table 3. Through Figure 6 and Table 3, it can be seen that the daily offsets were rather stable over 2017 and 2018, and ranged from 1 to 3 TECU. The daily offsets of PANG were slightly lower than those of IGSG. The mean VTEC offset for PANG-Jason3 was 1.1 TECU in 2017 and 1.2 TECU in 2018; For IGSG-Jason3, the mean VTEC offsets were 2.4 and 2.1 TECU in 2017 and 2018, respectively. The mean offset just reflect the system errors between the technologies of GNSS and satellite altimetry. Moreover, the mean STD of VTEC for PANG-Jason3 and IGSG-Jason3 were 2.45 and 2.65 TECU, respectively. In 2017, the mean STD were 2.5 and 2.8 TECU for PANG-Jason3 and IGSG-Jaon3, while in 2018 they were 2.4 and 2.5 TECU. All of these prove that PANG has a better consistency with Jason3 and the accuracy of PANG is more or less improved compared with IGSG. The daily system offset between satellite altimetry and GNSS mainly includes both the instrumental bias and the plasmaspheric component (i.e., VTEC above the height of the altimetry mission) [32]. Figure 6 clearly exhibits the daily offset and STD of VTEC for PANG-Jason3 (PANG minus Jason3) and IGSG-Jason3 (IGSG minus Jason3). The mean values of the VTEC offset and STD for PANG-Jaons3 and IGSG-Jason3 are listed in Table 3. Through Figure 6 and Table 3, it can be seen that the daily offsets were rather stable over 2017 and 2018, and ranged from 1 to 3 TECU. The daily offsets of PANG were slightly lower than those of IGSG. The mean VTEC offset for PANG-Jason3 was 1.1 TECU in 2017 and 1.2 TECU in 2018; For IGSG-Jason3, the mean VTEC offsets were 2.4 and 2.1 TECU in 2017 and 2018, respectively. The mean offset just reflect the system errors between the technologies of GNSS and satellite altimetry. Moreover, the mean STD of VTEC for PANG-Jason3 and IGSG-Jason3 were 2.45 and 2.65 TECU, respectively. In 2017, the mean STD were 2.5 and 2.8 TECU for PANG-Jason3 and IGSG-Jaon3, while in 2018 they were 2.4 and 2.5 TECU. All of these prove that PANG has a better consistency with Jason3 and the accuracy of PANG is more or less improved compared with IGSG. was 1.1 TECU in 2017 and 1.2 TECU in 2018; For IGSG-Jason3, the mean VTEC offsets were 2.4 and 2.1 TECU in 2017 and 2018, respectively. The mean offset just reflect the system errors between the technologies of GNSS and satellite altimetry. Moreover, the mean STD of VTEC for PANG-Jason3 and IGSG-Jason3 were 2.45 and 2.65 TECU, respectively. In 2017, the mean STD were 2.5 and 2.8 TECU for PANG-Jason3 and IGSG-Jaon3, while in 2018 they were 2.4 and 2.5 TECU. All of these prove that PANG has a better consistency with Jason3 and the accuracy of PANG is more or less improved compared with IGSG.   After eliminating the daily offsets between GNSS and satellite altimetry technologies, the VTEC density distribution of PANG and IGSG with respect to Jason3 in 2017 and 2018 are plotted in Figure 7. The point located in the cyan line denotes that the VTEC between GIM and Jason3 are equal. All the VTEC in 2017 and 2018 fall in the region of 0~50 TECU, with most of the points located in the region of 0~20 TECU. The density around the cyan line for PANG vs. Jason3 was higher than that of IGSG vs. Jason3, which meant that the VTEC of PANG VTEC was closer to those of Jason3.  After eliminating the daily offsets between GNSS and satellite altimetry technologies, the VTEC density distribution of PANG and IGSG with respect to Jason3 in 2017 and 2018 are plotted in Figure  7. The point located in the cyan line denotes that the VTEC between GIM and Jason3 are equal. All the VTEC in 2017 and 2018 fall in the region of 0~50 TECU, with most of the points located in the region of 0~20 TECU. The density around the cyan line for PANG vs. Jason3 was higher than that of IGSG vs. Jason3, which meant that the VTEC of PANG VTEC was closer to those of Jason3. To clearly see the VTEC discrepancies with respect to Jason3, Figure 8 plots the histogram of VTEC differences for PANG-Jason3 and IGSG-Jason3. All the discrepancies for PANG and IGSG were lower than 8 TECU. 79% of the VTEC differences were within ±3 TECU for IGSG-Jason3, which increased to 82.5% for PANG-Jason3. With regarding to RMS, the VTEC RMS for IGSG-Jason3 was 2.7 TECU, which declined to 2.5 TECU for PANG-Jason3. In summary, validated by Jason3, PANG has a higher accuracy compared with IGSG. To clearly see the VTEC discrepancies with respect to Jason3, Figure 8 plots the histogram of VTEC differences for PANG-Jason3 and IGSG-Jason3. All the discrepancies for PANG and IGSG were lower than 8 TECU. 79% of the VTEC differences were within ±3 TECU for IGSG-Jason3, which increased to 82.5% for PANG-Jason3. With regarding to RMS, the VTEC RMS for IGSG-Jason3 was 2.7 TECU, which declined to 2.5 TECU for PANG-Jason3. In summary, validated by Jason3, PANG has a higher accuracy compared with IGSG. Figure 7. VTEC density distribution of PANG and IGSG with respect to Jason3 in 2017 and 2018. The cyan line shows that the GIM VTEC is equal to Jason3 VTEC.
To clearly see the VTEC discrepancies with respect to Jason3, Figure 8 plots the histogram of VTEC differences for PANG-Jason3 and IGSG-Jason3. All the discrepancies for PANG and IGSG were lower than 8 TECU. 79% of the VTEC differences were within ±3 TECU for IGSG-Jason3, which increased to 82.5% for PANG-Jason3. With regarding to RMS, the VTEC RMS for IGSG-Jason3 was 2.7 TECU, which declined to 2.5 TECU for PANG-Jason3. In summary, validated by Jason3, PANG has a higher accuracy compared with IGSG.

Monitor Ionospheric Disturbance
Through the aforementioned analysis, we can see a slight overall improvement for PANG compared with IGSG. The improvement was relatively obvious during a period of ionospheric disturbance, which could be trigged by the geomagnetic storms or enhanced solar activities. The solar condition was quiet in 2017 and 2018, and we tried to find the ionospheric disturbance caused by geomagnetic storms. To analyse ionospheric features during storm-time, the storm onset time needed to be clearly defined. The Dst (Disturbance Storm Time) index is an index of magnetic activity derived

Monitor Ionospheric Disturbance
Through the aforementioned analysis, we can see a slight overall improvement for PANG compared with IGSG. The improvement was relatively obvious during a period of ionospheric disturbance, which could be trigged by the geomagnetic storms or enhanced solar activities. The solar condition was quiet in 2017 and 2018, and we tried to find the ionospheric disturbance caused by geomagnetic storms. To analyse ionospheric features during storm-time, the storm onset time needed to be clearly defined. The Dst (Disturbance Storm Time) index is an index of magnetic activity derived from a network of near-equatorial geomagnetic observatories that measures the intensity of the globally symmetrical equatorial electrojet. The geomagnetic storm in our manuscript is defined as the period that the Dst index startes decreasing until lower than −100 nT. Based on this definition, two geomagnetic storms were detected on 8 September 2017 and 26 August 2018, respectively. The responses of the ionosphere to the two geomagnetic storms were examined.

Ionospheric Responses for the Geomagnetic Storm that Happened on 8 September 2017
The geomagnetic storm on 8 September 2017 was the 18th strongest geomagnetic storms since January 1994 (https://www.spaceweatherlive.com/en/auroral-activity/top-50-geomagnetic-storms). The Dst index started to decrease from −1 nt at UT 22:00 on 7 September to −124 nT at UT 2:00 on 8 September. We selected the region of latitude from −30 • to 30 • and longitude from 120 • to 180 • . This region was chosen because its ionosphere was in a high level of ionisation during the onset time of the storm, and 15 points with a distance interval of 15 • in latitude and 30 • in longitude were selected. The VTEC variation of the 15 points from 6 September 2017 to 10 September 2017 are illustrated in Figure 9. The predicted VTEC values are calculated by the median values from the previous 15 days, and the gray shadow denotes the predicted VTEC ± 2*sigma. We are able to clearly notice the variations of VTEC in the normal days (6 September~7 September and 9 September~10 September) had a similar trend, and the predicted VTEC values had a better consistency with the estimated VTEC. Caused by the geomagnetic storm on 8 September, the estimated VTEC values showed a sudden increase compared with the predicted VTEC, especially for the points at B30 • L150 • , B-15 • L150 • and B-15 • L180 • . The ionosphere in low latitude is more active than that in medium and high latitudes, the maximum VTEC for the points at latitude 0 • and 15 • in a normal day usually exceeded 40 TECU. When the geomagnetic storm happened, the highest VTEC for the points at longitude 150 • could even reach up to 60 TECU. The VTEC differences between IGSG and PANG were close to zero or usually no more than 2 TECU for the normal days. During storm-time, IGSG and PANG still achieved a better consistency for the points at latitudes ± 15 • and ± 30 • , and the VTEC differences were within 5 TECU. However, for the three points at the equator (i.e., B0 • L120 • , B0 • L150 • and B0 • L180 • ), the VTEC differences were relatively larger, marked as black arrows in Figure 9, and the maximum VTEC differences even reached up to 10 TECU. This suggests PANG normally agrees well with IGSG, however, the VTEC discrepancies between them at the equator even exceed 10 TECU during a geomagnetic storm. The VTEC values derived from satellite altimetry were used to prove which solution was more reliable during the geomagnetic storm. The ionosphere was relatively active in low latitude, and the VTEC discrepancies between PANG and IGSG near the equator (from latitude -30 o to 30 o ) showed larger values than those in other regions. Therefore, we mainly analysed the VTEC around the equator from PANG, IGSG and Jason-3 on 08 September, where the VTEC differences of PANG minus Jason3 and IGSG minus Jason3 are illustrated in Figure 10. The number of GNSS stations used by PANG and IGSG to model the global ionosphere were 339 and 333, respectively. The distribution of the stations selected by PANG and IGSG were similar and without significant differences. The GNSS stations are sparsely distributed across the ocean area. This means that the available GNSS observations for modelling the ionosphere over the oceans are limited, as a consequence, the accuracy of the estimated VTEC above the ocean is relatively lower compared with that on the land area. In addition, the IGSG modelled the global ionosphere without an a priori model, just using GNSS observations. Therefore, it showed large VTEC differences for IGSG-Jason3 over the ocean area. While PANG introduced not only an a priori IRI-2016 model but also more available BDS and Galileo satellites, as a result, the VTEC differences of PANG-Jason3 were smaller than those of IGSG-Jason3. For PANG-Jason3, 83.42% of the TEC differences were within ± 5 TECU, while the percent for IGSG-Jason3 was just 52.38%, and 10% of the VTEC differences for IGSG-Jason3 exceeded 10 TECU. The VTEC values derived from satellite altimetry were used to prove which solution was more reliable during the geomagnetic storm. The ionosphere was relatively active in low latitude, and the VTEC discrepancies between PANG and IGSG near the equator (from latitude −30 • to 30 • ) showed larger values than those in other regions. Therefore, we mainly analysed the VTEC around the equator from PANG, IGSG and Jason-3 on 08 September, where the VTEC differences of PANG minus Jason3 and IGSG minus Jason3 are illustrated in Figure 10. The number of GNSS stations used by PANG and IGSG to model the global ionosphere were 339 and 333, respectively. The distribution of the stations selected by PANG and IGSG were similar and without significant differences. The GNSS stations are sparsely distributed across the ocean area. This means that the available GNSS observations for modelling the ionosphere over the oceans are limited, as a consequence, the accuracy of the estimated VTEC above the ocean is relatively lower compared with that on the land area. In addition, the IGSG modelled the global ionosphere without an a priori model, just using GNSS observations. Therefore, it showed large VTEC differences for IGSG-Jason3 over the ocean area. While PANG introduced not only an a priori IRI-2016 model but also more available BDS and Galileo satellites, as a result, the VTEC differences of PANG-Jason3 were smaller than those of IGSG-Jason3. For PANG-Jason3, 83.42% of the TEC differences were within ± 5 TECU, while the percent for IGSG-Jason3 was just 52.38%, and 10% of the VTEC differences for IGSG-Jason3 exceeded 10 TECU.
The VTEC comparisons of PANG-Jason2 and IGSG-Jason2 on 8 September 2017 are plotted in Figure 11. In the comparison with Jason2, because there are no GNSS observations in the Eastern Pacific Ocean, both of IGSG and PANG showed larger VTEC differences compared with Jason2. In the Western Pacific Ocean, Indian Ocean and Atlantic Ocean, PANG and IGSG have scattered GNSS stations, but PANG introduces BDS and Galileo signals in the global ionosphere estimation, which adds more IPP and somewhat improves the precision of the estimated ionosphere. Additionally, with the help of IRI-2016, the percent of VTEC differences within ± 5 TECU for PANG-Jason2 achieved 85.65% during the geomagnetic storm. The percent for IGSG-Jason2 was just 68.94%. Therefore, we can see that the PANG had a better performance during the geomagnetic storm. The VTEC comparisons of PANG-Jason2 and IGSG-Jason2 on 8 September 2017 are plotted in Figure 11. In the comparison with Jason2, because there are no GNSS observations in the Eastern Pacific Ocean, both of IGSG and PANG showed larger VTEC differences compared with Jason2. In the Western Pacific Ocean, Indian Ocean and Atlantic Ocean, PANG and IGSG have scattered GNSS stations, but PANG introduces BDS and Galileo signals in the global ionosphere estimation, which adds more IPP and somewhat improves the precision of the estimated ionosphere. Additionally, with the help of IRI-2016, the percent of VTEC differences within ± 5 TECU for PANG-Jason2 achieved 85.65% during the geomagnetic storm. The percent for IGSG-Jason2 was just 68.94%. Therefore, we can see that the PANG had a better performance during the geomagnetic storm. The ionospheric response for the geomagnetic storm that occurred on 26 August 2018 is shown in Figure 12. The level of the geomagnetic storm was not severe enough to be ranked in the top 50 geomagnetic storms since 1994. The Dst index started decreasing from 5 nT at UT 19:00 on 25 August to the lowest -156 nT at UT 10:00 on 26 August, then back to -78 nT at UT 15:00 and remained stable afterwards. We selected 15 points in the region of low latitude and the ionosphere over the 15 points

Ionospheric Responses for the Geomagnetic Storm Happened on 26 August 2018
The ionospheric response for the geomagnetic storm that occurred on 26 August 2018 is shown in Figure 12. The level of the geomagnetic storm was not severe enough to be ranked in the top 50 geomagnetic storms since 1994. The Dst index started decreasing from 5 nT at UT 19:00 on 25 August to the lowest −156 nT at UT 10:00 on 26 August, then back to −78 nT at UT 15:00 and remained stable afterwards. We selected 15 points in the region of low latitude and the ionosphere over the 15 points were relatively active during the sudden decrease of the Dst index. The solar activity was at a quiet level in 2018 and the geomagnetic storm was not as severe as that in 8 September 2017. Therefore, the VTEC values were normally lower than 30 TECU, even during the onset time of the geomagnetic storm, and the VTEC values at the equator were no more than 40 TECU. Because the ionosphere is quiet during this period, the VTEC differences between IGSG and PANG were very small, usually not larger than 2 TECU. The estimated VTEC values were consistent with the predicted VTEC until the onset of the geomagnetic storm at the day boundary between 25 August and 26 August. Then, during the geomagnetic storm, the estimated VTEC values were obviously larger than the predicted VTEC, and the VTEC differences between IGSG and PANG were also becoming larger. The large VTEC differences are marked by black arrows in Figure 12, and the maximum VTEC difference exceeded 5 TECU. After the storm, the consistency between the estimated and predicted VTEC values were recovered. Through comparing Figure 12 with Figure 9, it can be found that the VTEC differences between IGSG and PANG varied with the magnitude of the ionospheric disturbance. The greater the disturbance, the larger the difference. Similarly, the VTEC differences of PANG-Jason3 and IGSG-Jason3 during this storm time (from UT 18:00 on 25 August to UT 18:00 on 26 August ) are illustrated in Figure 13. We can observe that the VTEC discrepancies for PANG-Jason3 were smaller compared with those of IGSG-Jason3. 87.92% of the VTEC differences for PANG-Jason3 were within ± 5 TECU, and the percent for IGSG-Jason3 was 69.12%. Through comparing Figure 13 with Figure 10, it was noted that the percent of the VTEC Similarly, the VTEC differences of PANG-Jason3 and IGSG-Jason3 during this storm time (from UT 18:00 on 25 August to UT 18:00 on 26 August) are illustrated in Figure 13. We can observe that the VTEC discrepancies for PANG-Jason3 were smaller compared with those of IGSG-Jason3. 87.92% of the VTEC differences for PANG-Jason3 were within ± 5 TECU, and the percent for IGSG-Jason3 was 69.12%. Through comparing Figure 13 with Figure 10, it was noted that the percent of the VTEC differences within ± 5 TECU for PANG-Jason3 during the two storms were almost the same. Whilst the percent increased by 16.74% for IGSG-Jason3 during the geomagnetic storm in 2018. This implies that the accuracy of IGSG is lower during a geomagnetic storm, and declined even further during a severe ionospheric disturbance, while in a high or a low level of ionospheric disturbance, PANG showed a stable performance with a higher accuracy. The VTEC comparisons of PANG-Jason2 and IGSG-Jason2 during this storm-time are shown in Figure 14. PANG agreed better with Jason2 compared with IGSG. The large VTEC differences of IGSG-Jason2 mainly existed in the Western Pacific Ocean, where there were limited GNSS stations. Benefiting from the IRI-2016 model, PANG-Jason2 did not show obvious larger VTEC differences for the places with less available GNSS data. For the region with sparsely distributed GNSS stations, such as the Western Pacific Ocean, because PANG introduced more Galileo and BDS signals, the accuracy of PANG was slightly higher than that of IGSG. Compared with Figure 11, the improvement brought by PANG in Figure 14 was not as obvious as that in Figure 11. This was because the ionospheric disturbance was relative lower for the storm on 26 August 2018 compared with that on 8 September 2017, and the more intense the ionospheric disturbance was, the more obvious advantages PANG had. The VTEC comparisons of PANG-Jason2 and IGSG-Jason2 during this storm-time are shown in Figure 14. PANG agreed better with Jason2 compared with IGSG. The large VTEC differences of IGSG-Jason2 mainly existed in the Western Pacific Ocean, where there were limited GNSS stations. Benefiting from the IRI-2016 model, PANG-Jason2 did not show obvious larger VTEC differences for the places with less available GNSS data. For the region with sparsely distributed GNSS stations, such as the Western Pacific Ocean, because PANG introduced more Galileo and BDS signals, the accuracy of PANG was slightly higher than that of IGSG. Compared with Figure 11, the improvement brought by PANG in Figure 14 was not as obvious as that in Figure 11. This was because the ionospheric disturbance was relative lower for the storm on 26 August 2018 compared with that on 8 September 2017, and the more intense the ionospheric disturbance was, the more obvious advantages PANG had.

Summary and Conclusions
We applied all the available GNSS satellites and frequencies to model the global ionosphere. Meanwhile, the IRI-2016 model was introduced to enhance the strength of the ionosphere model derived by GNSS, especially for the region with sparsely distributed GNSS stations. We processed two years of data to analyse the performance of the estimated ionosphere. Through comparing the PANG with IAACs' ionosphere products, we concluded that the newly developed PANG had a higher precision. The mean RMS in 2017 and 2018 were always no more than 3 TECU no matter when compared with IGSG, CODG, ESAG or UPCG. Then, satellite altimetry was used to objectively evaluate the precision of the ionosphere. The results showed that the RMS of PANG was just 2.5 TECU compared with satellite altimetry data in 2017 and 2018. Additionally, it was better than the RMS of IGSG, which was 2.7 TECU.
The advantage of PANG was more obvious during the ionospheric disturbance. We selected two geomagnetic storms that occurred in 2017 and 2018 to examine the performance of PANG during storm-time. Both the geomagnetic storm caused ionosphere disturbance, showing a sudden VTEC increase. During the disturbance period, the VTEC differences were 5 ~ 10 TECU in the equatorial region, and the VTEC differences were larger with the increase in magnitude of the ionospheric disturbance. Through validation by the altimeter satellite Jason2/3, it was proven, compared with IGSG, PANG had a higher consistency with the VTEC from Jason2/3, and PANG showed a good performance over the region with limited GNSS stations, where the precision of IGSG obviously decreased. Therefore, through adding more available GNSS observations and introducing the advanced IRI model, the capability of GNSS in ionosphere modeling is enhanced. This can help us to better understand the ionospheric conditions and model the ionosphere in real-time with a higher accuracy in the future.

Summary and Conclusions
We applied all the available GNSS satellites and frequencies to model the global ionosphere. Meanwhile, the IRI-2016 model was introduced to enhance the strength of the ionosphere model derived by GNSS, especially for the region with sparsely distributed GNSS stations. We processed two years of data to analyse the performance of the estimated ionosphere. Through comparing the PANG with IAACs' ionosphere products, we concluded that the newly developed PANG had a higher precision. The mean RMS in 2017 and 2018 were always no more than 3 TECU no matter when compared with IGSG, CODG, ESAG or UPCG. Then, satellite altimetry was used to objectively evaluate the precision of the ionosphere. The results showed that the RMS of PANG was just 2.5 TECU compared with satellite altimetry data in 2017 and 2018. Additionally, it was better than the RMS of IGSG, which was 2.7 TECU.
The advantage of PANG was more obvious during the ionospheric disturbance. We selected two geomagnetic storms that occurred in 2017 and 2018 to examine the performance of PANG during storm-time. Both the geomagnetic storm caused ionosphere disturbance, showing a sudden VTEC increase. During the disturbance period, the VTEC differences were 5~10 TECU in the equatorial region, and the VTEC differences were larger with the increase in magnitude of the ionospheric disturbance. Through validation by the altimeter satellite Jason2/3, it was proven, compared with IGSG, PANG had a higher consistency with the VTEC from Jason2/3, and PANG showed a good performance over the region with limited GNSS stations, where the precision of IGSG obviously decreased. Therefore, through adding more available GNSS observations and introducing the advanced IRI model, the capability of GNSS in ionosphere modeling is enhanced. This can help us to better understand the ionospheric conditions and model the ionosphere in real-time with a higher accuracy in the future.