Status of CAS global ionospheric maps after the maximum of solar cycle 24

As a new Ionosphere Associate Analysis Center (IAAC) of the International GNSS Service (IGS), Chinese Academy of Sciences (CAS) started the routine computation of the real‐time, rapid, and final Global Ionospheric Maps (GIMs) in 2015. The method for the generation of CAS rapid and final GIMs and recent updates are presented in the paper. The quality of CAS post‐processed GIMs is assessed during 2015–2018 after the maximum of solar cycle 24. To perform an independent and fair assessment, Jason‐2/3 Vertical Total Electron


Introduction
The Global Navigation Satellite Systems (GNSSs), especially the Global Positioning System (GPS), have provided the opportunities to continuously monitor the variability of the global ionosphere. To provide a public service for the reliable generation of ionospheric Total Electron Contents (TECs), the Ionosphere Working Group (IWG, Feltens 2003) of the International GNSS Service (IGS, Dow et al. 2009) was created in 1998. The initial analysis centers of the IWG include the Center for Orbit Determination in Europe (CODE, Dach et al. 2009), European Space Agency/European Space Operation Center (ESA/ ESOC), Jet Propulsion Laboratory (JPL), and Polytechnic University of Catalonia (UPC-IonSAT). The Energy, Mines and Resources/Natural Resource Canada (EMR/ NRCan) also contributed as an Ionosphere Associate Analysis Center (IAAC) during 1998-2002. Different from the computation of regional ionospheric models with a dense network of permanent receivers (Bergeot et al. 2014), the distribution of GNSS receivers used in the generation of Global Ionospheric Maps (GIMs) is notably inhomogeneous given the sparse ground receivers over

Open Access
Satellite Navigation https://satellite-navigation.springeropen.com/ *Correspondence: wangningbo@aoe.ac.cn 1 Aerospace Information Research Institute (AIR), Chinese Academy of Sciences (CAS), Beijing, China Full list of author information is available at the end of the article oceanic regions and continents of the southern hemisphere (Hernández-Pajares et al. 2010). With different mapping techniques e.g., the use of solar-geomagnetic reference frame (Schaer 1999) or Spline/Kriging interpolations (Mannucci et al. 1998;Orús et al. 2005), rapid and final GIMs are generated by individual IAAC and provided in the IONosphere EXchange (IONEX, Schaer et al. 1998) format v1.0 with a temporal resolution of 15 min to 2 h and a spatial resolution of 2.5° × 5° in geographic latitude and longitude. The IGS GIM is a weighted combination of the TEC maps from two IAACs (for the final GIM: CODE and JPL, and for the rapid GIM: ESA and JPL), which is initially generated by UPC and presently by the University of Warmia-Mazury (UWM) in Olsztyn, by applying an independent assessment of different GIMs as described in Hernández-Pajares et al. (2017) and Krankowski et al. (2017). The latency of the file accessibility is commonly 1-2 days and 1-2 weeks for the rapid and final GIMs, respectively. Predicted GIMs with 1, 2 or 5 d leading time are also provided by CODE, ESA/ESOC, and UPC since 2009 Li et al. 2018;Schaer 1999).
Different global ionosphere mapping techniques are employed by different IAACs. For instance, the Spherical Harmonic (SH) expansion is used by CODE (Schaer 1999) and ESA/ESOC (Feltens 2007), and a three-shell model and tomographic approach are used by JPL (Komjathy et al. 2005;Mannucci et al. 1998) and UPC-IonSAT (Hernández-Pajares et al. 1999). With the inclusion of GLONASS data (Dach et al. 2009;Vergados et al. 2016;Yasyukevich et al. 2020) and the transition from low (2 h) to high (15 min to 1 h) temporal resolutions in the generation of individual IAAC GIMs, the accuracy of the IGS combined GIMs has been continually improved (Hernández-Pajares et al. 2009). To support more reliable global ionospheric service, EMR/NRCan resumed GIM computation since 2016, and at the same time, the Chinese Academy of Sciences (CAS) and Wuhan University (WHU) were included as new IAACs of the IGS (Hernández-Pajares et al. 2016). The GIMs from the three new and the first four IAACs were compared during the past solar cycle 23 with the TOPEX/Jason Vertical TECs (VTECs), and a good consistency among the seven global TEC maps was reported in Roma-Dollase et al. (2018). Additionally, the German Geodetic Research Institute of the Technical University of Munich (DGFI-TUM) was recommended as a new contributor after the quality assessment of their GIMs (Schmidt 2018). The contribution from more than seven analysis centers with multiple ionospheric mapping techniques shall benefit the combination and application of the IGS GIM. It is worth mentioning that Massachusetts Institute of Technology (MIT) Haystack Observatory (Rideout and Coster 2006) and the German Aerospace Center (DLR) in Neustrelitz ) also produce global VTEC maps for scientific applications.
As discussed in Li et al. (2015), CAS uses the SH expansion plus Generalized Trigonometric Series (GTS, Yuan and Ou 2004) function (SHPTS) to generate its global TEC map. The generation of the real-time, rapid, and final GIMs started in 2015. Some modifications, e.g., the inclusion of BeiDou Navigation Satellite System (BDS) data and the use of a Modified GTS (MGTS) function for local ionospheric modeling , have been made for its GIM computation since then. The routine delivery of CAS GIMs to the Crustal Dynamics Data Information System (CDDIS, Noll 2010) started from January 2017, with a latency of 1 day and 3 days for the rapid and final products, respectively. CAS rapid and final GIMs cover the time spans 2015-now and 1998-now, respectively, which are accessible from the IONEX archive of CDDIS. Aside from the generation of predicted and Real-Time (RT) GIMs , the routine assessment of different GIMs with respect to the GPS differential Slant TECs (dSTEC, Feltens et al. 2011;Hernández-Pajares et al. 2017;Krankowski et al. 2017) and the Jason VTECs are also implemented since 2017. While CAS predicted and RT GIMs have not yet been officially provided to the IGS, the products are downloadable from CAS repository itself (ftp:// ftp. gipp. org. cn/ produ ct/ ionex/). The purpose of the paper is to update the GIM computation implemented at CAS and its GIM products (as of January 2020). Since the rapid and final GIMs have been routinely provided to the IGS, these two products are the focus of the work. The description of CAS global ionosphere mapping technique is first presented, followed by the data sets and processing strategies used for the assessment of different GIMs. Here, GPS dSTECs derived from the receivers of the Multi-GNSS Experiment network (MGEX, Montenbruck et al. 2017) of the IGS as well as the Jason VTECs are selected as the references following Hernández-Pajares et al. (2017). The quality assessment of the CAS rapid and final GIMs is performed during a 4-year period after the maximum of solar cycle 24 with proper interpretation and discussion. Finally, the summary and conclusions are given.

CAS global ionosphere mapping technique
The SHPTS method is employed by CAS to generate global TEC maps, which constitutes the extraction of raw ionospheric observables, the global VTEC modeling, the local VTEC modeling, as well as the generation of global TEC and Root-Mean-Square (RMS) maps. The ionospheric observables can be directly extracted from the geometry-free combination of GNSS dual-frequency pseudorange and carrier phase measurements (Afraimovich et al. 2013;Hernández-Pajares et al. 2011). The pseudorange measurements are significantly affected by the multipath and observation noises, whereas carrier phase measurements suffer from the unknown integer phase ambiguities. To improve the quality of the derived ionospheric observables, the Carrier-to-Code Leveling (CCL) approach is used (Mannucci et al. 1998;Schaer 1999). The global and local ionospheric modeling techniques as well as the generation of global TEC/RMS maps are discussed in the following subsections, followed by the overview of CAS ionospheric products.

Global ionospheric modeling
For global ionospheric modeling using the SH expansion, the VTEC is modeled in a solar-geomagnetic reference frame where the electron distribution is assumed to be mostly stationary (Schaer 1999). The Slant TEC (STEC)S ϕ ′ , ′ , z between a satellite and a receiver is mapped into the VTEC V SH ϕ ′ , ′ of the Ionospheric Pierce Point (IPP) under the thin-layer assumption with a mapping function M(z) of satellite zenith angle z as follows: with the brief notation c n,m for SH coefficients to be estimated,Y n.m ϕ ′ , ′ for the spherical harmonic expansion,m and n for the order and degree of the SH expansion, ϕ ′ and ′ for the latitude and longitude of the IPP in the solar-geomagnetic reference frame,P n,|m| (·) for the normalized associated Legendre function, R = 6378 km for the mean radius of the Earth, and H ion = 450.0 km for the altitude of the assumed single-layer ionosphere. For the conversion between earth-geographic and solar-geomagnetic reference frames, we refer to Schaer (1999).
The spherical harmonic expansion up to degree 15 is used. The time spacing between model sessions is set to 1 h, and a piecewise linear interpolation function (1) between the SH coefficients of consecutive sessions is used (Schaer 1999). Satellite and receiver Differential Code Bias (DCB) parameters are assumed constant within one/three days. The intra-frequency biases, i.e., GPS and GLONASS P1-C1 DCBs, are directly estimated from the two associated pseudorange measurements simultaneously tracked by the receivers (Wang et al. 2016a). For the receivers only supporting the simultaneous tracking of GPS C/A and P2 or GLONASS C1 and P2 pseudoranges, the P1-C1 corrections of the satellite are applied for the GPS/GLONASS P1-P2 DCB estimation. While the estimation of GLONASS satellite-toreceiver pair biases is employed by CODE and JPL (Dach et al. 2009;Vergados et al. 2016), GPS and GLONASS satellite-plus-receiver DCBs are split into the sum of satellite-dependent and receiver-dependent biases in our case. To remove the rank deficiency in the estimation of satellite-and receiver-specific DCBs, a zero-constellation-mean constraint is applied in the least-squares solution. Considering that receiver P1-P2 DCBs listed in the header of IONEX files might be the C1-P2 DCBs, satellite and receiver DCB solutions are also provided in the Bias-SINEX (Schaer 2016) format v1.0 for unambiguous purpose. In our rapid and final GIM computation, satellite and receiver DCBs are estimated together with global ionospheric parameters (Li et al. 2015), whereas in our real-time GIM computation, satellite DCBs are fixed to CAS MGEX DCB solutions (Wang et al. 2019a and receiver DCBs are estimated together with local ionospheric parameters, i.e., the MGTS function . The classification of rapid and final GIMs is based on the input observation length as well as the output product latency. The rapid solution is estimated in 24 h observation arcs and a latency of 1 day. To ensure the smoothness at day boundaries, the final GIM computation employs additional observations of 4 h lengths prior and after the day of the estimated ionospheric parameters with a latency of 3 days.

Local ionospheric modeling
The local VTEC is modeled as the sum of two-dimensional polynomials of latitude and local time and a finite Fourier series of local time in the original GTS function (Li et al. 2012;Yuan and Ou 2004). The MGTS function is formed by the inclusion of the spherical cap coordinate system and a local-time dependent weighting function, which has been employed for the local ionosphere modeling in the CAS GIM computation since early 2018. Using the same mapping function and altitude of the assumed single-layer ionosphere as given in Eq.
(1), the variability of the local VTEC in MGTS V MGTS (ϕ d , d , h) is described as where n,m and k are the degrees of the polynomial and Fourier series functions, E n,m ,C k and S k are the model coefficients to be estimated,t denotes the local time, and ϕ d and d are the spherical distance between IPPs and the station along the latitudinal and longitudinal directions, which are computed as in which (ϕ, ) are the geographic latitude and longitude of the IPP,(ϕ 0 , 0 ) are the geographic latitude and longitude of the station, and (ϕ c , c ) are the latitude and longitude of IPPs in the spherical cap coordinate system (Haines 1988). A satellite-elevation (e) and local-time (t) dependent stochastic model, as shown in Eq. (4), is employed in MGTS for local VTEC modeling.
where the maximum degrees n max ,m max and k max in Eq.
(2) are set to 2, 2 and 4, respectively. With the use of a satellite-elevation and local-time dependent weighting function, the effects of observation noises and day-night differences of the ionospheric variability shall be reduced.
The MGTS function has also been used in the generation of CAS MGEX DCB products since mid-2018 ).

Generation of global TEC and RMS grid maps
In the IGS convention, the two-dimensional GIMs are generated with a spatial resolution of 2.5° in geographic latitude and 5° in geographic longitude . As reported in Li et al. (2015), the method of Differential Areas for Differential Stations (DADS, Yuan and Ou 2002) is employed to generate CAS global TEC and RMS maps. (2) The VTEC V i and its corresponding RMS σ i at the i-th grid point are computed using global SH and local MGTS models as follows: with the notations V SH,i and σ SH,i for the VTEC and RMS computed from the global SH expansion for the i-th grid point,V MGTS,i,m and σ MGTS,i,m for the VTEC and RMS computed from the local MGTS model corresponding to the i-th grid point and the m-th station, M for the number of contributing stations, which is determined by the elevation (> 25°) between stations and line-of-sight grid point on the assumed single-layer sphere, P m for the weight in the computation of MGTS-based grid VTEC, and RMS values as follows: where σ 2 0,m is the posterior variance of the unit weight of the local MGTS least-squares solution, and e denotes the elevation between station and the corresponding ionospheric grid point.
The ionospheric VTEC at the grid point near GNSS stations is calculated from the local MGTS model. In the case that two or more stations contribute to the grid VTEC computation, an elevation and ionosphere-modeling-accuracy dependent weighting function is used to adjust the VTECs computed from local ionospheric models. For those grid points far from GNSS receivers, e.g., the global ocean, VTECs are directly calculated by the global SH expansion. Considering that local ionospheric models generally exhibit better performance than the global model, the accuracy of the corresponding grid VTEC estimates around GNSS stations shall be improved in SHPTS-based GIMs (Li et al. 2015).

Overview of CAS ionospheric products
The routine computation of the SHPTS-based rapid and final GIMs started in 2015, and the delivery of the two products to CDDIS started in 2017. Since satellite and receiver biases are simultaneously estimated in the GIM computation, the daily, weekly, and monthly GPS/GLONASS DCB solutions are generated in the Bias-SINEX format v1.0 (Schaer 2016) at the same time and provided in the DCB archive of CAS repository (ftp:// ftp. gipp. org. cn/ produ ct/ dcb/). The 1, 2 and 5 d predicted GIMs are generated at CAS by predicting each SH coefficient using the Fourier series expansion. When individual SH coefficients are generated, the predicted GIM can be also reconstructed by the SH expansion . A "predicting-plusmodeling" method is developed at CAS to generate the RT-GIM with multi-GNSS (including GPS and GLO-NASS L1 + L2, BDS B1 + B2 and Galileo E1 + E5a RT data streams provided by the IGS, MGEX and other regional GNSS tracking networks ). The VTEC is modeled in a solar-geographic reference frame (with 2 h shifted in the mean sun fixed longitude) by the SH expansion to degree15. To mitigate the impacts of the unstable RT data streams, e.g., the interruption of data streams, the predicted TEC information is also included in the RT-GIM computation of CAS . Triggered by the Klobuchar-style ionospheric coefficients produced by CODE, Klobuchar-like coefficients of GPS and NeQuickG coefficients of Galileo are also re-computed at CAS (Wang et al. 2019b). A detailed description of the method for the Klobuchar-like and NeQuickG coefficient estimation using dual-frequency GPS/GLONASS measurements is presented in Wang et al. (2016bWang et al. ( , 2017. As the third generation BeiDou Navigation Satellite System (BDS-3) employs a new BeiDou global ionospheric model (BDGIM, Yuan et al. 2019) for single-frequency ionospheric delay error mitigation, the estimation of BDGIM coefficients is also supported since early 2018. The re-estimated ionospheric coefficients of GPS, Galileo, and BeiDou-3 are routinely generated in the RINEX format (available at ftp:// ftp. gipp. org. cn/ produ ct/ brdion/), and the corresponding IONEX files, along with the predicted, real-time, rapid, and final GIMs are also provided at the ionex archive of CAS repository (ftp:// ftp. gipp. org. cn/ produ ct/ ionex/). The status of different GIM products generated by CAS and other IAACs of the IGS is summarized in Table 1.

Data sets and processing
As the rapid and final GIMs of CAS have been routinely delivered to the CDDIS, we mainly focus on the performance analysis of these two products during their routine computation period. The method for the independent and fair assessment of ionospheric electron content models was discussed and summarized in Hernández-Pajares et al. (2017). Here, Jason-2/3 altimeter VTEC and GPS dSTEC are used as the references to evaluate the quality of the rapid and final GIMs over the oceanic and continental regions, respectively.
Jason-series altimeters operate at a mean orbit height of around 1 330 km with a Ku-band primary frequency and a C-band auxiliary frequency (Fu and Haines 2013). Different from the GNSS technique, Jason VTEC observations are directly obtained in the vertical direction (1, 2, 5 d) n/a n/a n/a
, Schmidt (2018) between the ocean surface and the satellite orbit altitude, which inhibits the mapping errors in the conversion from slant to vertical TECs. Because of the lower orbit height of Jason altimeters compared with GNSS satellites, the Jason VTECs are expected to be smaller than the GNSS-derived VTECs as the latter ones also contain the contribution of plasmaspheric electron contents (PECs, Lee et al. 2013). To reduce the pronounced measurement noises in the Jason derived VTECs, a moving averaging procedure with a sliding window of 16 s is used to generate smoothed VTECs. As reported in Hernández-Pajares et al. (2017), the Jason VTEC presents a Standard Deviation (STD) of around 1.0 TECu after applying the smoothing procedure. As illustrated in Fig. 1, the daily mean VTECs from Jason-2 and -3 altimeters are computed along the satellite ground track during the overlap period of the two satellite missions. The global mean VTECs from Jason-2 and -3 are almost identical, with a mean value of 15.7 and 15.4 TECu, respectively, during the Day of Year (DOY) 077, 2016 and DOY 063, 2017. Compared to the UPC rapid GIM (UQRG), Jason-3 VTECs exhibit significant negative biases and a pronounced systematic offset of around 2.7 TECu is found between Jason-2 and -3 VTEC values. Since Jason-2 and -3 satellites operate on different orbits during the overlap period, the spatial and temporal discrepancies between the two satellites might result in the pronounced deviation between Jason-2 and -3 VTECs. The systematic bias is removed by subtracting 2.7 TECu from Jason-3 smoothed VTECs to generate the consistent results for comparison.
Although the CCL approach has been widely applied for GNSS STEC computation, the accuracy of the CCL derived STECs is largely related to the length of the smoothing arc (Schaer 1999), which is also significantly affected by the intra-day variation of receiver DCBs as well as the leveling errors introduced by the pseudorange noise and multipath effect through the leveling arc (Ciraolo et al. 2007). The dSTEC analysis is developed to weight the GIMs from different IAACs and used to evaluate the performance of different ionospheric models (Feltens et al. 2011;Hernández-Pajares et al. 2017;Roma-Dollase et al. 2018). As illustrated in Fig. 2, the dSTEC is defined as the difference between the STEC measured at a given epoch along a continuous arc and the STEC measured when the GNSS satellite reaches its highest elevation. After the pre-processing of raw carrier phase measurements, including the detection of gross errors and cycle slips, the derived dSTEC o (t i ) at epoch t i is formed as where I L i ,L j (t i ) denotes the geometry-free combination of dual-frequency carrier phase measurements on L i and L j frequencies at epoch t i, α L i ,L j is the frequencydependent factor, which equals , and t ref denotes the epoch at which the satellite zenith angle is the smallest through the observation arc. In the same manner, the GIMbased dSTEC at the same epoch and location is computed as with the brief notation V (ϕ i , i , t i ) for the VTEC value obtained from the corresponding GIM products, and M(z i ) for the ionospheric mapping function as explained in Eq.
(1). The dSTEC assessment is performed by analyzing the differences between dSTEC m (t i ) and dSTEC 0 (t i ) across the entire observation arc of all involved stations. The dSTEC is directly computed from the dual-frequency carrier phase measurements, which exhibits much smaller noises and multipath effects than pseudorange observations. To remove the unknown integer ambiguity within the continuous observation arc, a single-difference between the given STEC and the STEC at the highest satellite elevation is formed. As discussed in Hernández-Pajares et al. (2017), the dSTEC involves different geometries and different times of the ionosphere, which is reliable for the assessment of different electron content models. To perform a fair dSTEC assessment, only the GNSS stations which have not contributed to the GIM computation should be used. Since GIMs are commonly generated using IGS stations, 55 globally distributed MGEX stations are selected as the test sites for the GPS dSTEC assessment.

Results and discussion
As GIMs are produced with a temporal resolution of 15 min to 2 h whereas Jason and GPS measurements have different sampling intervals and geographic locations, the GIM TECs are interpolated to the same epoch and location of Jason VTEC and GPS dSTEC observations by a bivariate spatial interpolation and a linear temporal interpolation between two consecutive TEC maps . Aside from CAS rapid and final GIMs, GIM products of the other six IAACs as well as the IGS combined ones are also included in the assessment. The comparison with respect to the Jason-2 and -3 VTECs over global oceans is first presented, followed by the assessment with the GPS dSTEC at the selected MGEX stations.

Assessment with the Jason VTEC
The comparison of Jason-2 and GIM VTECs was first performed for the geomagnetically quiet days (Kp < 3) during a low solar activity period, i.e., the year 2016. As Jason-2 employs a non-sun-synchronous orbit advancing around 2° per day, it takes about 90 days to cover all local times. The data was then binned in three seasons, i.e., equinox (DOY 50-110 and 234-294), June solstice (DOY 111-233) and December solstice (DOY 1-50 and 295-366), to check the seasonal variation of the ionosphere. The global TEC maps generated with the VTEC data from Jason-2 data and CAS final GIMs (CASG) and their relative differences are depicted in the geographic latitude and local time coordinate system (with 2° × 0.25 h bins) as shown in Fig. 3. As a quick comparison between Jason-2 and CAS final GIM VTECs, CAS GIM can reproduce the temporal and latitudinal variations of the global ionosphere. The large VTECs appear in daytime and low-latitude and equatorial regions in all seasonal cases. The notably small VTECs are found at high latitudes during the southern and northern winters. As the orbit altitude of GNSS satellites is much higher than that of Jason-2 altimeter, it is suspected that GIM VTEC should be larger than the Jason-2 VTEC. However, the difference maps show both positive values (i.e., VTEC[GIM] > VTEC[Jason-2]) and negative values (i.e., VTEC[GIM] < VTEC[Jason-2]). The negative value reaches more than 30% of the Jason-2 VTEC for the latitudes higher than 50° during the nighttime of all seasons, and even larger (> 50%) at high latitudes of the southern winter. The negative discrepancy can be partly explained by the influence of the nighttime plasmaspheric electron density, especially for the southern hemisphere during low solar activity conditions (Jee et al. 2014). The positive differences are found at low and middle latitudes in most of daytime, which shows a close correlation with equatorial anomaly structures reproduced by the two TEC sources. During June and December solstices, the positive differences are around 10%, while during the equinox the corresponding values reach more than 40%. The difference can be largely attributed to the low accuracy of the GIMs in equatorial regions caused by the inadequacy of the single-layer assumption in the presence of large latitudinal gradients (Hernández-Pajares et al. 1999;Juan et al. 1997). The difference map also presents a notable hemispheric asymmetry in June solstice, which shows a significantly negative difference at middle and high latitudes of the southern hemisphere. As the number of the GPS/GLONASS receivers in the ocean and the southern hemisphere are limited for GIM computation, the worse performance is expected over the corresponding regions (Hernández-Pajares et al. 2009), which is one of the main limitations of the GIMs provided by the IGS.
Since seven IAACs have routinely provided their rapid and final GIMs to the IGS, we first select the GIMs generated with different techniques as representatives to analyze the long-term variation and latitudinal dependency of GIM errors. The selected GIMs are from CAS, CODE, JPL, UPC as well as the IGS. The long-term variation of the bias and STD of the selected GIMs in comparison with the Jason-2/3 VTECs is depicted in Fig. 4. Both positive and negative biases can be found in GIM-minus-Jason VTEC data series, and the positive difference is more pronounced during high solar activity conditions (2015) than low solar activity conditions (2018). JPLG exhibits significantly positive deviation with respect to the Jason-2/3 VTECs, which varies in the range of 2.0-5.5 TECu across the entire test period. The biases of CASG and CODG are at the same level, slightly smaller than that of UQRG, UPCG and IGSG. As the altimeterderived VTECs contain the uncalibrated bias of around few TECu (Azpilicueta and Brunini 2008;Jee et al. 2014), the STD of GIM-minus-Jason VTEC differences is also checked, which provides a measure of how well the GIM can reproduce changes in the Jason VTECs while neglecting the biases. As shown in the bottom panel of Fig. 4, UQRG presents the smallest STD and no notable differences are found among the remaining GIMs. The result indicates an overall better performance of UQRG (in terms of STD) compared to other GIMs over the ocean. As few GPS/GLONASS stations contribute to the GIM computation over oceanic regions, the interpolation or extrapolation techniques shall be of great importance. Considering that the Kriging interpolation is employed for UQRG computation, it suggests that Kriging interpolation exhibits better performance than Spline interpolation in the generation of UPC's global ionospheric TEC maps (Hernández-Pajares et al. 2017). Note that UQRG and UPCG are generated with a temporal resolution of 15 min and 2 h, respectively. The higher temporal resolution of UQRG also reduces interpolation errors. In addition to the pronounced correlation with the solar activity (i.e., reduced solar activity levels correspond to decreased STD values), the seasonal variation is also observed in the STD series. The largest and smallest STD values appear around March/September equinoxes and June solstice, respectively, which keeps in proper accord with the seasonal variability of the ionosphere.
The latitudinal variation of GIM errors compared to Jason VTECs is checked and plotted in Fig. 5. The results are computed from the daily mean values in the entire test period within individual 3° latitude bins. While JPLG shows notably positive biases with respect to the Jason-2/3 VTECs at different latitudes, CASG, CODG and IGSG present positive biases in low latitude and equatorial regions but negative biases in high latitude regions. Note that the biases of UQRG and UPCG show less latitude-dependent characteristics, which exhibit the smallest variation range compared to other GIMs. In terms of STDs, the largest discrepancy is found in low-latitude and equatorial regions for all GIMs, and the discrepancy in the southern hemisphere is more significant than that in the northern hemisphere. The results keep in proper accord with the previous findings generated during different test periods Li et al. 2015Li et al. , 2018Rovira-Garcia et al. 2016), which can be attributed to the latitudinal dependence of the plasmaspheric electron content, in addition to the pronounced ionosphere variation in low-latitude and equatorial regions and a small number of GPS/GLONASS stations in continents of the southern hemisphere.
The Jason-2/3 VTEC provides an independent way to assess the quality of GIMs over oceans where few GPS/GLONASS stations are used for GIM computation. In the comparison of the long-term and latitudinal variation of GIM-minus-Jason differences, it is found that JPLG shows significantly positive biases (i.e., VTEC[JPLG] > VTEC[Jason]) compared to the GIMs from CAS, CODE and UPC. While a three-layer model is used by JPL, a single-layer assumption and a tomographic method are employed by CAS/CODE and UPC, respectively. The pronounced biases of JPL GIMs are more likely related to their specific data processing strategies rather than the use of a multi-layer ionospheric model, although the reasons are not well explained. The quality of GIMs is low in equatorial and southern hemisphere regions. Although the multi-layer assumption or tomographic method have been employed for global ionosphere modeling, the current IONEX v1.0 is only used to generate the two-dimensional TEC maps in the IGS community, which inhibits the proper reproduction  of ionosphere variation in the presence of large ionospheric latitudinal gradients in equatorial regions. The small number of contributing stations in the ocean and southern hemisphere is the main reason for the poor performance of the GIMs in the corresponding regions. The use of advanced interpolation techniques, e.g., the Kriging interpolation employed by UPC, shall improve the performance of GIMs in the case of no additional data sources contributing to GIM computation (Hernández-Pajares et al. 2017). A comparison of the final and rapid GIMs from the seven IAACs with respect to the Jason VTECs is given in Tables 2 and 3, respectively. The short time span of EMR/NRCan (EMRG) and WHU (WHUG) final GIMs, especially for WHU rapid GIMs (WHRG), is due to the missing of files from CDDIS during the test period. CAS final GIM presents a comparable performance with CODE final GIM (21.3% versus 21.7% in relative RMS errors), which is slightly worse than UPC final GIM (21.1%) but better than the GIMs of other IAACs and the IGS combined one. JPL final GIM shows significantly positive biases (3.04 TECu), but performs at the same level in STD as other GIMs except for WHUG. The large STD of WHUG is caused by the poor performance of GIM during DOY 309, 2017 and DOY 075, 2018. The relative RMS error of WHUG can be reduced from 27.3% to 23.0% if the above time period is excluded in the analysis. For the rapid GIM assessment, CAS rapid GIM (23.0% in relative RMS error) still performs better than other GIMs except for UPC rapid GIMs (20.1-20.6%) during the test period.
As summarized in Tables 2 and 3, the quality of the final GIMs (see STD and relative RMS errors) is slightly better than the rapid GIMs from individual IAACs except for the UPC GIMs. This can be explained by the longer observation arcs used for the final GIM computation compared to the rapid ones (three-day versus one-day solutions). UQRG and UPCG are generated using tomography with Kriging and Spline interpolations (Hernández-Pajares et al. 1999;Orús et al. 2005), respectively, whereas UHRG and UPRG are generated by reducing the temporal resolution of UQRG to 1 and 2 h (Roma-Dollase et al. 2018). The use of Kriging interpolation and higher temporal resolution results in the overall better performance of UPC rapid GIMs than its final GIM (Hernández-Pajares et al. 2017). For the GIM generated with the same modeling technique, its accuracy depends on the temporal resolution, the input data (GPS versus GNSS), contributing stations (number and distribution), as well as other factors (Rovira-Garcia et al. 2016). This can be justified in the comparison of the GIMs computed by CODE, EMR/NRCan and ESA/ ESOC, which used the identical SH expansion method but with different networks of GNSS receivers and data processing strategies. The IGS combined final and rapid GIMs perform at the same level (23.7% versus 23.6%) during the test period. Considering the comparable performance of the GIMs from the three new and the first four IAACs, one can foresee the reliability of the IGS combined GIMs can be improved with the inclusion of the GIMs produced by different IAACs.

Assessment with GPS dSTEC
GPS dSTEC assessment provides an alternative way to evaluate the ability of ionospheric models to reproduce both spatial and temporal gradients in the ionosphere (Feltens et al. 2011). Compared to the conventional GPS STEC assessment derived by the CCL method, GPS dSTEC measurements are free of leveling errors and intra-day variation of receiver biases, which exhibit much better accuracy (< 0.1 TECu). GPS dSTEC assessment was performed during a 16-month period starting from September 2017. Six GIMs (same as the ones used in Jason-2/3 VTEC assessment) are selected as representatives to examine the temporal and latitudinal variations of GIM errors in comparison with the GPS dSTECs. The daily bias and STD of the involved GIMs during the fourth quarter of 2017 are depicted in Fig. 6. JPLG still presents notably larger biases (> 0.5 TECu) compared to the GIMs from other IAACs (− 0.1 to + 0.2 TECu) as well as the IGS (0.3 TECu). The notable jumps on DOY 264 and 275 are found in the bias series of CASG (see the top panel), which is caused by the fact that the computation of SHPTS-based GIM is failed and only the SH expansion is used for the GIM generation during the corresponding period. The bias of the UPC rapid and final GIMs also shows a significant jump on DOY 311 and DOY 314, which presents a proper correlation with the geomagnetic event on DOY 311 (see https:// omniw eb. gsfc. nasa. gov/ form/ dx1. html). The sudden increase of STDs is found for all GIMs on DOY 251, which is more likely related to the X-class solar flare on September 6, 2017 (Berdermann et al. 2018). The result indicates that the dSTEC assessment is sensitive to the large gradients of the ionosphere, especially in the presence of solar and/ or geomagnetic events. Since the dSTEC addresses how well spatial and temporal ionospheric gradients can be captured, the statistics are likely sensitive to the spatial and temporal resolutions of the ionospheric models. As a result, the dSTEC assessment might benefit for the GIMs with a higher temporal resolution, e.g., the UQRG.
The latitudinal variation of GIM errors compared to the GPS dSTEC is depicted in Fig. 7 during the fourth quarter of 2017. Note that the result is not presented at the individual station but for specific latitude bins considering that less test sites are located at high latitudes whereas more stations at middle latitudes, especially for the northern hemisphere. Significantly positive biases are found at low-latitude and equatorial stations, and slightly negative biases occur at mid-and high-latitude stations. The pronounced deviation in the GIM-minus-GPS dSTEC differences in equatorial and low-latitude regions can be largely attributed to the inferior reproduction of latitudinal gradients with the GIMs. As for the latitudinal variation of STDs, the low-and mid-latitude stations present the largest and smallest STDs, respectively, and the STDs at high-latitude stations of the southern hemisphere are comparable to those at low-latitude stations. The pronounced STDs further confirm the poor performance of the GIMs in high latitude regions of the southern hemisphere, which should be considered by individual IAACs in their future GIM computation, and by users when carrying out associated applications.
As discussed in Hernández-Pajares et al. (2017), the altimeter VTEC and GPS dSTEC assessments exhibit good consistency in the evaluation of UQRG with 26 GPS stations on the islands. While the altimeter VTEC provides a fair assessment of the GIMs over global oceans, the uncalibrated biases will induce a systematic offset of few TECu (Azpilicueta and Brunini 2008). The VTEC measurements from different altimetry missions may also present systematic biases, which must be removed before any applications. As a complementing reference dataset to altimeter VTEC assessment, the GNSS dSTEC is computed from dual-frequency carrier phase measurements in a continuous observation arc, which is free of pseudorange noises and satellite/receiver biases, thus exhibiting a high level of accuracy (< 0.1 TECu). Note that the external GNSS stations which are not used in the ionospheric model estimation should be employed, e.g., MGEX stations used for the GIM dSTEC assessment in the present study. Additionally, more detailed characteristics can be found in GNSS dSTEC assessment, e.g., the decreased performance of the GIMs during space weather events and the inferior accuracy of the GIMs in high latitude regions of the southern hemisphere as shown in Figs. 6 and 7, respectively. Considering its external and fair assessment of ionospheric correction models, the dSTEC assessment has been used for the routine ranking of different GIMs and the combination of IGS GIMs with an independent network of GPS receivers (Hernández-Pajares et al. 2017).
The comparisons of final and rapid GIMs from the seven IAACs with respect to the GPS dSTECs are summarized in Tables 4 and 5, respectively. JPL GIMs exhibit the best performance with a relative RMS error of 9.9% and 10.4% for the final and rapid products, respectively. CAS GIMs are in the middle with a relative error of 15.3% for CASG and 15.7% for CARG. The worst performance of the WHU GIMs relates to the pronounced errors during the period of late 2017 and early 2018. In the assessments compared to the Jason-2/3 VTECs and GPS dSTECs, JPL GIMs present the best performance over the continents, but worst performance over the oceans (except for WHU GIMs). The UPC rapid GIM UQRG shows the best performance over both continental and oceanic regions, which outperforms its final GIM UPCG by 1.0%-3.5%. The quality of the new GIMs from CAS is comparable to the CODE GIMs, which is slightly better than the GIMs from ESA, EMR/NRCan and WHU during the test period.
During 2015-2018, the IGS final GIM is produced by the combination of JPLG and CODG results, and the rapid GIM by JPRG and ESRG. The correlation coefficient between IGSG and JPLG/CODG is 99.52%/99.53%, and the corresponding value is 99.32%/99.27% between IGRG and JPRG/ESRG. It is understandable that the JPL GIMs are used for the generation of the IGS combined rapid and final GIMs as of their best performance in the GPS dSTEC assessment (see Tables 4 and 5). With the consideration of the significantly positive bias of the JPL GIMs over the global ocean (see Tables 2 and 3), it is suggested that the independent GNSS stations with a more homogeneous distribution should be used for the ranking the GIMs of individual IAAC as well as the potential combination of the IGS GIM.

Summary and conclusions
The SHPTS method is employed by CAS to generate the rapid and final GIMs. This paper updates the generation of the CAS rapid and final GIMs and summarizes the latest status (as of January 2020) of the CAS ionospheric products. In addition to the rapid and final GIMs, the predicted and real-time GIMs as well as the re-estimation of global broadcast ionospheric coefficients of GPS, Galileo and BDS-3 are also provided by CAS. Since CAS rapid and final GIMs have been routinely provided to the IGS, the quality of the two GIMs is assessed during a four-year period starting from January 2015.
Compared to the Jason-2/3 VTEC, the relative RMS error is 23.0% and 21.3% for the CAS rapid and final GIMs, respectively. The GPS dSTEC is free of pseudorange noises and satellite/receiver biases, and therefore has a high level of accuracy. Compared to the GPS dSTEC, the relative RMS error is 15.7% for the rapid GIM and 15.3% for the final GIM of CAS. The CAS GIM shows comparable performance with the CODE GIM,