The Zenith Total Delay Combination of International GNSS Service Repro3 and the Analysis of Its Precision

Abstract

Currently, ground-based global navigation satellite system (GNSS) techniques have become widely recognized as a reliable and effective tool for atmospheric monitoring, enabling the retrieval of zenith total delay (ZTD) and precipitable water vapor (PWV) for meteorological and climate research. The International GNSS Service analysis centers (ACs) have initiated their third reprocessing campaign, known as IGS Repro3. In this campaign, six ACs conducted a homogeneous reprocessing of the ZTD time series spanning the period from 1994 to 2022. This paper primarily focuses on ZTD products. First, the data processing strategies and station conditions of six ACs were compared and analyzed. Then, formal errors within the data were examined, followed by the implementation of quality control processes. Second, a combination method is proposed and applied to generate the final ZTD products. The resulting combined series was compared with the time series submitted by the six ACs, revealing a mean bias of 0.03 mm and a mean root mean square value of 3.02 mm. Finally, the time series submitted by the six ACs and the combined series were compared with VLBI data, radiosonde data, and ERA5 data. In comparison, the combined solution performs better than most individual analysis centers, demonstrating higher quality. Therefore, the advanced method proposed in this study and the generated high-quality dataset have considerable implications for further advancing GNSS atmospheric sensing and offer valuable insights for climate modeling and prediction.

1. Introduction

Currently, the evolution of satellite Earth observing technologies has significantly advanced their applications in acquiring atmospheric parameters for weather and climate applications [1]. Originally designed for position, navigation, and surveying engineering, global navigation satellite systems (GNSSs) have already become valuable tools for atmospheric monitoring over several decades of groundbreaking advancements [2,3,4,5]. The ground-based GNSS atmospheric monitoring technique takes GNSS receivers as atmospheric sensors. These receivers track changes in satellite signals as they traverse Earth’s atmosphere, providing accurate, broad-coverage, and densely sampled atmospheric parameters of zenith total delay (ZTD) and precipitable water vapor (PWV) [6,7,8,9,10]. Compared to heritage techniques of measuring water vapor, such as radiosonde and water vapor radiometer, GNSS atmospheric sounding techniques offer distinct benefits like long-term stability, superior spatiotemporal resolution, and all-weather capability, making them well suited for weather and climate studies [11,12,13]. In recent years, with the innovative utilization of GNSS-derived ZTD and PWV, various statistical, numerical, and artificial intelligence-enhanced models have been developed for monitoring weather and climate extremes, especially heavy precipitation, tropical cyclones, and droughts [14,15,16,17,18], as well as analyzing climate change fingerprints [19,20,21,22]. For example, Ding et al. [23] studied the spatial–temporal variations in GNSS-derived PWV globally and examined how these variations reflect and impact global climate change. Wang et al. [24] used GNSS-derived PWV data from 56 stations near the ocean over more than 10 years to study the relationship between PWV and sea surface temperature. Based on this, they investigated the relationship between PWV and the El Niño–Southern Oscillation (ENSO). Many studies have shown that ZTD can serve as a substitute for PWV in studying climate change and meteorological forecasting events. For example, Li et al. [25] developed an improved model for detecting heavy precipitation using GNSS-derived ZTD, demonstrating the potential of these measurements for more accurate precipitation monitoring. Zhao et al. [26] showed that real-time precise point positioning (PPP)-based ZTD can be employed effectively for forecasting precipitation, highlighting the utility of GNSS-derived products in operational meteorology. Li et al. [27] introduced a new cumulative anomaly-based model that utilizes ZTD data, enhancing the detection of heavy precipitation through GNSS-derived tropospheric products, thereby improving early warning systems for severe weather. Therefore, these measurements hold substantial potential for effectively monitoring climate change and unraveling the intricate dynamics of weather and climate extremes.

As per previous research and findings, high-quality atmospheric parameters are of paramount importance in supporting weather and climate research [28,29,30]. However, these data, especially when covering a long period, often experience inconsistencies, i.e., temporal inhomogeneities, due to updates in international terrestrial reference frames (ITRF) and applied models, the use of varying elevation cut-off angles, the implementation of different mapping functions, and other changes in processing strategies [20,31]. Therefore, the GNSS data analysis strategy and applied models must remain consistent throughout the entire processing period by homogenously reprocessing to ensure the reliability and quality of the results. Many existing studies have proven that reprocessed atmospheric parameters are more suitable for weather and climate studies [28,32,33]. For example, Steigenberger et al. [34] compared a consistent time series of ZTD and tropospheric gradients from homogenously reprocessed GNSS and very long baseline interferometry (VLBI) solutions. They found that maintaining the homogeneity of these reprocessed time series is crucial for avoiding misunderstandings that may arise from changes in individual models. Thomas et al. [35] conducted a homogenous reprocessing of global GNSS data, focusing on 12 Antarctic stations. They also found that reprocessed GNSS-derived tropospheric estimates, using advanced models, now show significant potential for integration into both regional and global numerical weather models. Consequently, with the accumulation of nearly 30 years of GNSS observations since the early 1990s, this juncture presents an opportune moment to leverage the full potential of GNSS atmospheric monitoring techniques in climate applications. At the current stage, the first and foremost step forward is to generate long-term, reliable, and homogeneous GNSS climate records.

Over the past few years, international organizations have made considerable strides in enhancing the accuracy and consistency of long-term GNSS atmospheric parameters for climate applications. For example, the EUREF Permanent Network (EPN) facilitates high-quality GNSS data reprocessing for geodesy and climate applications. The second reprocessing campaign, known as EPN Repro2, covered all EPN stations from January 1996 to December 2013, with an expansion at the end of 2014 for atmospheric parameters. This effort involved about 280 stations, each processed by at least three analysis centers (ACs) to ensure the quality of their outputs [30,36]. Furthermore, the European Cooperation in Science and Technology (COST) Action is an intergovernmental framework aimed at fostering the coordination of nationally funded research activities across Europe. Typically, the EU COST Action ES 1206 project, particularly “Working Group 3: GNSS for climate monitoring”, aims to promote the use of reprocessed GNSS data for climate research, standardizing algorithms and methods to ensure the long-term stability and reliability of data [37,38]. More importantly, the International GNSS Service (IGS) has always been at the forefront, dedicating efforts to reprocessing activities that involve reanalyzing the entire dataset collected by the IGS network since 1994. Its primary goal is to maintain internal consistency by using the most current models and methods to reanalyze the GNSS data coherently [39,40,41]. Since late 2020, the IGS has completed the third reprocessing activity, known as IGS Repro3. In contrast to previous campaigns, i.e., Repro1 and Repro2 [42], Repro3 expands its dataset to include observations not only from the GPS and GLONASS constellations but also from the Galileo system. In addition, it was found that the Repro2 combination products were not suitable for long time span processing compared to the Repro3 [43]. Therefore, building upon all these endeavors, this study mainly focuses on the atmospheric parameters derived from the IGS Repro3 reprocessing initiatives.

It is reported that there is a total of eleven IGS ACs contributing to the Repro3 initiative (https://cddis.nasa.gov/archive/gnss/products/ accessed on 10 September 2024). However, data from five ACs, i.e., Natural Resources Canada (NRCan, EMR), the Massachusetts Institute of Technology (MIT), the National Geodetic Survey (NGS), Université de la Rochelle (ULR), and Wuhan University (WHU), were not included in this study due to the absence of tropospheric data. Therefore, this study focuses on examining the atmospheric parameters obtained from the remaining six ACs, i.e., the Center for Orbit Determination in Europe (CODE) [44,45], the European Space Agency (ESA) [46], GeoForschungsZentrum (GFZ) [47,48,49], Groupe de Recherche en Géodésie Spatiale (GRG) [50], the Jet Propulsion Laboratory (JPL), and the Graz University of Technology (TUG) [51].

The motivation of this study stems from the recognition that the ZTD data provided by each AC originate from different solutions, with some organizations generating official combined products [30,40,52]. Currently, ZTD combined products for IGS Repro3 have not been provided. Compared to ZTD estimates from a single AC, combined products can mitigate or even eliminate systematic errors from individual models or algorithms by integrating data from multiple ACs, thereby offering higher reliability and precision [40]. The combination process facilitates the evaluation of the consistency among different ACs. By comparing and integrating results from multiple ACs, any inconsistencies can be identified and corrected, thus improving the overall data quality [52]. Additionally, while the six ACs provide reprocessed ZTD estimates, relying on data from a single AC often results in time series for certain stations that are either short-term or sparsely populated with observations. Such limitations can affect their utility in climate research, where continuous and dense observations are crucial for ensuring the reliability and accuracy of the results. Furthermore, integrating data from multiple ACs also allows for the extension of the duration of the time series, ensuring continuous and dense observations necessary for climate research. This study not only focuses on combining data from multiple ACs to enhance overall data quality but also seeks to assess the results of each AC through rigorous quality control. This ensures that the data have been thoroughly validated before being used in data combination and climate change research.

Consequently, the main contribution of this study lies in the development of an advanced method that amalgamates data from all six ACs and the implementation of a rigorous data quality control process. This approach results in a more comprehensive and consistent GNSS climate dataset being generated, which not only improves the reliability and quality of the GNSS atmospheric parameters but also significantly expands their potential and uptake for climate applications.

The structure of the rest of the paper is outlined as follows: Section 2 and Section 3 detail the data and methodologies utilized in this study. Section 4 showcases the quality control outcomes and the precision evaluation of the combined results. Then, Section 5 offers discussions and conclusions.

2. Data

In this study, a total of ten sets of ZTD time series were used. Six of these were generated by different IGS ACs, i.e., the CODE, GFZ, GRG, JPL, TUW, and TUG, as part of the IGS Repro3 campaign. Additionally, a combined ZTD solution was generated by merging a ZTD series obtained from each individual AC. The remaining three sets comprise ZTD data from the ERA5 reanalysis dataset, VLBI, and radiosonde measurements.

2.1. IGS Data and Products

As stated in Section 1, ZTD data obtained from six ACs were included in this study. The data are provided in the form of daily files and stored in SINEX-TRO format, covering the period from the beginning of 1994 to the end of 2022. GNSS processing strategies vary among ACs, resulting in six distinct contributing solutions. Each solution uses different software and settings and covers different stations.

Table 1. IGS Repro3 processing options for each contributing solution at the six ACs.

	CODE	ESA	GFZ	GRG	JPL	TUG
Software	Bernese 5.3	NAPEOS 4.6	EPOS.P8	GINS	GIPSYX	GROOPS
System	GPS (2000)	GPS (1995)	GPS (1994)	GPS (2000)	GPS (1994)	GPS (1994)
	GLO (2002)	GLO (2009)	GLO (2012)	GLO (2008)		GLO (2009)
	GAL (2013)	GAL (2015)	GAL (2014)	GAL (2016)		GAL (2013)
Solution type	Network	Network	PPP	PPP	PPP	Network
Station number	493	718	331	236	559	1182
Secular pole	IERS2010	IERS2010	IERS2010	IERS2010	IERS2010	IERS2010
Ocean tidal	FES2014b	EOT11A	FES2014b	FES2014b	GOT4.8	FES2014b
TROPO mapping function	VMF1	GMF	VMF1	VMF2	VMF1	VMF3
GRADS mapping function	CHEN-HERRING	CHEN-HERRING	MacMillan	CHEN-HERRING		CHEN-HERRING
Elevation cut-off	3°	10°	7°	8°	7°	5°
Output rate	1 h	1 h	1 h	2 h	5 min/30 s	5 min

outlines the detailed processing characteristics of each AC’s solution. Each AC employs its own in-house developed software packages for processing, including Bernese GNSS Software package V5.3 [53], NAPEOS 4.6 [54], EPOS.P8 [55], GINS [56], GIPSYX [57], and GROOPS [58]. Regarding the employed constellations, five ACs, i.e., the CODE, ESA, GFZ, GRG, and TUG, incorporated GPS, GLONASS, and Galileo in their solutions, while JPL processed GPS alone. Additionally, it can be seen from this table that in the IGS Repro3 campaign, ACs have adopted advanced ocean tidal loading models. These models include EOT11A [59] developed by the ESA, GOT4.8 [60] from the JPL, and FES2014b [61] utilized by other ACs. For tropospheric modeling, two primary mapping functions were employed, the Global Mapping Function (GMF) [62] and the Vienna Mapping Function (VMF) [63,64]. All ACs used the VMF, with the exception of the ESA, which opted for the GMF. Furthermore, each AC contributes a solution at different output rates: the CODE, ESA, and GFZ provide hourly updates; GRG outputs every 2 h; the JPL offers solutions at 30 s or 5 min; and the TUG provides updates every 5 min.

Regarding the distribution of stations for these ACs, as depicted in Figure 1, the TUG operates the most with 1182 stations, followed by the ESA with 718, the JPL with 559, the CODE with 493, GFZ with 331, and GRG with the fewest, with 236 stations. Figure 1 employs varying dot colors to indicate a station is monitored by a different number of ACs. It can be found that most stations were intermittently monitored by a single AC, leading to sporadic data availability. Consequently, the time series for these sites exhibit short durations and sparse observations. Statistical analysis shows that 41.32% of the stations had three or more ZTD solutions, 15.33% had only two solutions, and 43.35% had a single solution. Furthermore, it can also be discovered from this figure that while the 1331 stations are globally distributed, there are significant differences in their density. Typically, North America has a high density of stations, particularly in the major cities and coastal areas of the United States and Canada. In contrast, South America has fewer stations, most concentrated in southeastern Brazil. Europe features a dense network, especially in Western and Central countries such as Germany, France, and the United Kingdom, whereas Eastern Europe shows a sparse distribution. In Africa, stations are mainly located in South Africa, albeit with a lower density. Asia shows dense station coverage in Japan and South Korea, while the Middle East and South Asia have more sparse distributions. Australia and New Zealand have a dense and moderate station distribution, respectively, particularly along southeastern coastal areas. Other scattered oceanic stations, while less dense, cover extensive areas.

In addition, Figure 2 shows the temporal evolution of the daily station count processed by each AC, highlighting the significant variation in the size of the daily station networks. For example, the JPL typically processes about 80 stations per day, while the TUG can accommodate over 800 stations on its peak days.

2.2. ERA5, VLBI, and Radiosonde Data

In this study, ZTD estimates derived from independent sources including ERA5, VLBI, and radiosonde observations were used to validate the performance of GNSS ZTD.

Typically, the ERA5, as the most recent atmospheric reanalysis from the European Centre for Medium-Range Weather Forecasts, offers significant improvements over its predecessor, i.e., the ERA-Interim, which was in use from 1979 to 2019. A major improvement in ERA5 is its enhanced spatial resolution. Additionally, ERA5 employs a four-dimensional variational data assimilation technique using an advanced model cycle (41r2), allowing for more accurate and detailed atmospheric, land surface, and ocean wave data [65]. In this study, global ZTD data derived from the ERA5 reanalysis over the period of 1994–2019 were taken as the reference to evaluate the performance of GNSS-derived ZTD estimates. For the evaluation of GNSS ZTD, we used hourly atmospheric estimates from ERA5 at 37 pressure levels, with a horizontal resolution of 0.25° × 0.25°. The ZTD was calculated following the methods outlined in [66,67]. For GNSS stations located above or below the lowest pressure level of ERA5, atmospheric parameters at the altitudes of the stations were estimated through interpolation or extrapolation from four adjacent grids. Subsequently, these parameters were horizontally interpolated to the location of GNSS sites. Detailed methodologies for these processes are described in [68,69].

Since 1984, VLBI has continuously provided high-quality observations, capturing accurate information about long-term ZTD estimates. Due to their high accuracy, the atmospheric parameters derived from VLBI are also widely used for the validation and calibration of parameters obtained using GNSS atmospheric monitoring techniques. In this study, the accuracy of GNSS-derived ZTD (at 13 GNSS stations) was evaluated by comparing it with those obtained at their respective co-located VLBI stations. The horizontal distances between each pair of GNSS and VLBI stations were maintained within 500 m.

Radiosonde observations are known for providing the most extensive historical record of ZTD estimates, with over 900 radiosonde stations deployed worldwide. In this study, in addition to the above-mentioned ERA5 and VLBI datasets, ZTD estimates derived from sounding balloons were also used in evaluating the accuracy of GNSS-derived ZTD at 30 pairs of co-located GNSS and radiosonde stations. The horizontal distances between each pair of stations were within 50 km. The radiosonde data utilized in this study, sourced from Integrated Global Radiosonde Archive (IGRA) Version 2, comprised quality-controlled observations of temperature, humidity, and pressure from stations worldwide. Detailed methodologies for the calculation of ZTD using atmospheric variables from radiosonde observations are provided in [67]. Similar to the ZTD calculation using ERA5, vertical interpolation or extrapolation was employed to derive atmospheric variables at the GNSS station altitude based on whether the station was positioned above or below the lowest radiosonde level.

3. Methodologies

As stated in Section 1, the main novelty of this study is the development of a robust method that amalgamates GNSS ZTD from all six IGS ACs and the implementation of a rigorous data quality control process. Therefore, the methodologies used in this study can be categorized into two parts, i.e., the combination of long-term GNSS-derived ZTD time series from the six ACs and the validation of the obtained datasets.

3.1. An Advanced Method for GNSS Data Combination

The combination process initiates with the reading and verification of the original SINEX-TRO files provided by the ACs. In this initial stage, all daily time series are merged into a unified long-term time series. Then, a preliminary quality control process is performed to scrutinize the data with the aim of identifying and eliminating any gross errors. Only stations that meet the criterion of having at least three distinct solutions are selected for further analysis, i.e., the subsequent determination of weight values in the combination process. In the second stage, a combination is conducted to derive optimal weighting factors for each contributing solution. These calculated weights are then incorporated in the final stage, in which the combined ZTD time series is finally determined. Figure 3 illustrates the flowchart for the detailed data combination methodology, providing a clear view of the sequential steps and their connections.

3.1.1. Preprocessing of the Individual AC Data

First, all daily series are merged into a unified long-term time series.

At each station, multiple time series exist, each representing the ZTD estimated provided by different ACs. A critical requirement for the combination process is that the time series from different ACs must align in terms of length and sampling interval. To address this, a set of evenly spaced time labels at 1 h intervals is established. Moreover, linear interpolation is also performed for the ZTD time series where consecutive time labels do not align with the predefined evenly spaced time labels [70]. This interpolation facilitates the alignment of ZTD estimates, ensuring that they are synchronized and contemporaneous across all contributing ACs. Following this, a rigorous data validation stage is performed to identify and rectify any gross errors present in the data.

3.1.2. Methodology of Combination

In the next phase, an average of the ZTD from various ACs is calculated. The formula for determining the combined time series of ZTD typically involves the integration of individual ZTD measurements from different sources, which is as follows:

$${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{j}}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left[{\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}\right]}{\mathrm{n}}}$$

(1)

where ${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{j}}$ denotes the mean ZTD value at station $\mathrm{s}$ for epoch $\mathrm{j}$. Here, $\mathrm{s}$ represents the station index, $\mathrm{i}$ refers to the specific AC, and $\mathrm{j}$ indicates each time epoch. $\mathrm{n}$ represents the number of ACs at a given epoch. ${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}$ is the ZTD value for station $\mathrm{s}$ from AC $\mathrm{i}$ at epoch $\mathrm{j}$.

Following the generation of a preliminary time series of the ZTD, the next step involves calculating the standard deviations (STDs) for each station over weekly intervals. These calculations compare the preliminary time series with those derived from different ACs.

$${\mathrm{S}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}=\sqrt{{\displaystyle \frac{{\sum }_{\mathrm{j}=1}^{\mathrm{W}}{\left({\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}-\overline{\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}\right)}^2}{\mathrm{W}}}}$$

(2)

$${\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}={\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}-{\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{j}}$$

(3)

$$\overline{\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}={\displaystyle \frac{{\sum }_{\mathrm{j}=1}^{\mathrm{W}}{\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}}{\mathrm{W}}}$$

(4)

where ${\mathrm{S}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}$ denotes the STD for station $\mathrm{s}$ over a weekly interval. ${\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}$ and $\overline{\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}$ can be found in Equations (3) and (4). $\mathrm{W}$ is the number of data points in the week.

Moreover, the mean value of the standard deviations is determined for all stations of this AC over a week.

$${\mathsf{\sigma }}_{\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}={\displaystyle \frac{{\sum }_{\mathrm{s}=1}^{\mathrm{M}}{\mathrm{S}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}}{\mathrm{M}}}$$

(5)

where ${\mathsf{\sigma }}_{\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}$ denotes the mean ${\mathrm{S}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}$ for all stations. $\mathrm{M}$ is the total number of stations for AC $\mathrm{i}$.

These average standard deviations are employed to assign weights to the individual AC solutions during the final combination process. The formula can be represented as

$${\mathrm{P}}_{\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}={\displaystyle \frac{1/{\mathsf{\sigma }}_{\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}1/{\mathsf{\sigma }}_{\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}^2}}$$

(6)

where ${\mathrm{P}}_{\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}$ is the weight of AC $\mathrm{i}$ in week $\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}$. $\mathrm{n}$ represents the number of ACs.

The final formula for the ZTD combination is

$${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{c}\mathrm{o}\mathrm{m},\mathrm{j}}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left[{\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}\mathrm{j}}\times {\mathrm{P}}_{\mathrm{i},\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{k}}\right]$$

(7)

where ${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{s},\mathrm{c}\mathrm{o}\mathrm{m},\mathrm{j}}$ demotes the final ZTD value of station $\mathrm{s}$ for epoch $\mathrm{j}$.

With the use of all the formulas, a combined ZTD dataset can be finally generated.

3.2. Methods for Validating the Quality of GNSS Climate Records

To validate the quality of the GNSS climate records, we compared the time series from the six ACs and the combined series with VLBI, radiosonde, and ERA5 data, calculating the bias, RMS, and STD between the GNSS data and the reference data. The bias, RMS, and STD are calculated as

$$\mathrm{b}\mathrm{i}\mathrm{a}\mathrm{s}={\displaystyle \frac{{\sum }_{\mathrm{s}=1}^{\mathrm{M}}\left(\frac{{\sum }_{\mathrm{j}=1}^{\mathrm{t}}\left({\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{G}\mathrm{j}}-{\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{R}\mathrm{j}}\right)}{\mathrm{t}}\right)}{\mathrm{M}}}$$

(8)

$$\mathrm{R}\mathrm{M}\mathrm{S}={\displaystyle \frac{{\sum }_{\mathrm{s}=1}^{\mathrm{M}}\left(\sqrt{\frac{{\sum }_{\mathrm{j}=1}^{\mathrm{t}}{\left({\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{G}\mathrm{j}}-{\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{R}\mathrm{j}}\right)}^2}{\mathrm{t}}}\right)}{\mathrm{M}}}$$

(9)

$$\mathrm{S}\mathrm{T}\mathrm{D}={\displaystyle \frac{{\sum }_{\mathrm{s}=1}^{\mathrm{M}}\left(\sqrt{\frac{{\sum }_{\mathrm{j}=1}^{\mathrm{t}}{\left({\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{j}}-\overline{\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}\right)}^2}{\mathrm{t}-1}}\right)}{\mathrm{M}}}$$

(10)

$${\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{j}}={\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{G}\mathrm{j}}-{\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{R}\mathrm{j}}$$

(11)

$$\overline{\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}={\displaystyle \frac{{\sum }_{\mathrm{j}=1}^{\mathrm{t}}{\mathrm{d}\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{j}}}{\mathrm{t}}}$$

(12)

where $\mathrm{M}$ is the number of the stations. $\mathrm{t}$ is the number of ZTD observations in the series. ${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{G}\mathrm{j}}$ and ${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{R}\mathrm{j}}$ represent the GNSS and reference ZTD value in the time series for epoch $\mathrm{j}$.

4. Results

4.1. Quality Control

Formal error is one of the important indicators for evaluating data quality [30,52,70,71,72]. By analyzing formal error, outliers or data points that significantly deviate from expected values can be identified and discarded. This approach has been previously utilized for quality control in related studies. Earlier research employed a formal error to detect gross errors in GNSS ZTD, with threshold values ranging from 10 mm to 100 mm. However, variations in software, observational data sources, and processing workflows across different ACs can influence the resulting formal error values. To establish a more precise threshold, this study conducts a comprehensive analysis of the formal error.

We computed the daily median formal errors for each AC. An analysis was carried out on the ZTD formal error values for the six ACs. Box plots were generated as shown in Figure 4, which encompassed the lower whisker, Q1 (first quartile), Q2 (second quartile), Q3 (third quartile), and the upper whisker. In addition, the interquartile ranges (IQRs) and upper atypical limits (Q3 + 3 × IQR) were also calculated, and their specific values can be found in

Table 2. Quartiles calculated for the formal error values [mm] of the ZTD estimated by each AC.

	CODE	ESA	GFZ	GRG	JPL	TUG	Mean
Lower whisker	0.0	1.7	2.1	0.07	1.2	0.0	0.8
Q1	1.0	3.1	2.7	0.53	1.7	0.2	1.5
Q2	1.3	3.9	3.1	0.63	1.9	0.2	1.8
Q3	1.8	4.7	3.4	0.84	2.3	0.4	2.2
Upper whisker	3.0	7.1	4.4	1.30	3.1	0.7	3.3
Upper outlier limit	4.2	9.5	5.5	1.77	4.1	1.0	4.4

. We investigated the distribution and variability of formal errors across these ACs, providing insights into the quality and consistency of data products generated by each center. The results indicate that the ESA exhibits a wider range of formal error values, suggesting greater variability in their data processing. Conversely, the CODE, GRG, GFZ, and JPL demonstrate relatively concentrated formal error distributions, while the TUG exhibits the narrowest range. These findings contribute to a better understanding of the performance and reliability of each AC’s data. From the box plots, the ESA shows the lowest quality, while the TUG exhibits the highest quality.

In general, the upper outlier limit is commonly used as the threshold for identifying outliers [52], and it was determined to be the appropriate tolerance in this case. Based on the information in Table 2, the maximum upper outlier limit is found to be 9.5 mm. Therefore, we established that the individual parameter formal error must be less than 10 mm as the acceptable tolerance. The tolerance value is similar to those used in previous applications by others [30,70]. Values that do not meet these tolerance criteria will be excluded or removed from the analysis. The percentages of data removal for each AC are shown in

Table 3. Data removal rate for each AC after quality control.

AC	Data Removal Rate (%)
CODE	0.61
ESA	1.62
GFZ	0.00
GRG	2.07
JPL	0.05
TUG	0.17

. GRG and the ESA exhibit relatively high removal rates, at 2.08% and 1.62% respectively, while GFZ has the lowest rate at 0.00%. The data removal rates for the other ACs range from 0% to 1%. From the perspective of data removal rates, GFZ and the JPL exhibit the highest quality.

Figure 5 shows the variations in formal errors before and after quality control. Figure 5a displays the time series plot of daily median formal errors, while Figure 5b shows the daily median formal errors of the ZTD after removing values with significant errors. Additionally, the figure also presents the average of all daily median formal errors. As time progresses, the STD values gradually decrease, indicating a gradual improvement in the data processing quality.

In Figure 5, the ZTD estimates of the TUG clearly stand out in terms of their precision compared to all other ACs. Preliminary speculation suggests that the higher precision of the TUG solution may be attributed to its utilization of the latest generation of geodetic analysis software packages, GROOPS [41]. On the other hand, the TUG dataset stands out as the only one that incorporates the refined VMF3 a priori zenith delays and mapping functions [73] and utilizes raw (undifferenced and uncombined) observables from all available signal frequencies [51].

4.2. Combined Results

The weight of each AC for each week is shown in Figure 6. The average weights are presented in

Table 4. Table of average weights for each AC.

AC	Average Weight
CODE	0.11
ESA	0.14
GFZ	0.25
GRG	0.05
JPL	0.27
TUG	0.18

. The JPL and GFZ have the highest average weights among the analysis centers, with values of 0.27 and 0.25, respectively. This indicates that they play the most significant role in the overall analysis process and contribute the most. The TUG and ESA also have relatively high weights, at 0.18 and 0.14, respectively. These centers perform well in terms of data quality and stability, although they are slightly less influential than the JPL and GFZ. The CODE’s average weight is 0.11. Although the CODE’s contribution is not as substantial as the top four centers, it still has a notable impact on the overall analysis process. GRG has the lowest average weight, at only 0.05. This suggests that GRG’s data quality is relatively lower or its data show greater variability throughout the analysis process, leading to it being assigned a smaller weight in the combination process.

Using the method described above, it is possible to determine a combined time series of the ZTD, as shown in Figure 7. The red line in the figures represents the combined ZTD time series. Compared to the individual submissions from each AC (CODE, ESA, GFZ, GRG, JPL, TUG), the combined ZTD shows greater smoothness. This indicates that by aggregating data from multiple ACs, the combined time series effectively reduces random errors and outliers. The combined ZTD time series aligns with the overall trend in the individual AC time series but with less variability, demonstrating higher consistency and reliability. This consistency makes the combined ZTD more reliable for applications requiring accurate and consistent ZTD estimates. Therefore, by integrating data from multiple analysis centers, the combined ZTD time series significantly enhances smoothness, consistency, and reliability, making it an effective method for error reduction and data optimization.

4.3. Time Span and Data Completeness Analysis

After the combination, the time span and data completeness of each station’s data are altered compared to the data submitted by individual ACs. The breakdown of station time spans is as follows: 22.84% of the stations have a time span exceeding 20 years, 16.83% have a time span between 15 and 20 years, 23.29% have a time span between 10 and 15 years, 21.26% have a time span between 5 and 10 years, and 15.78% have a time span of less than five years. Furthermore, it was found that 81% of the stations had a data missing rate within 25%, indicating that the data from these stations is relatively complete. The specific details are illustrated in Figure 8.

In the final analysis, it has been determined that 212 sites meet the criteria of having a data completion rate higher than 50%, a time span longer than 15 years, and the involvement of at least three different ACs for data processing. These 212 sites represent 16% of the total number of sites. These stations are geographically diverse, covering different regions and environmental conditions. Their distribution allows us to obtain reliable data from various geographic areas. Their geographical distribution is illustrated in Figure 9. Furthermore, the data from these stations span over 15 years, enabling us to access long-term trends and patterns. This is crucial for studying long-term changes and predicting future trends. Lastly, these stations have relatively low data missing rates, ensuring the integrity and reliability of the data. This allows us to confidently analyze and apply these data in various applications. In summary, the data from these 212 stations possess characteristics of high reliability and accessibility, offering significant value to various application domains.

4.4. Comparison of Combined Data with Respect to Single-AC Data

The ZTD reprocessing results were contributed to by the six ACs, covering a total of 1331 stations. The analysis process did not consider all Repro3 stations or the entire Repro3 period. For station selection, we ultimately selected the 212 stations mentioned in the previous section for analysis.

To assess the compatibility between the combined ZTD and the ZTD provided by each AC, we calculated the RMS values of the differences between the combined ZTD series and each AC’s series. The daily RMS values for all stations are shown in Figure 10. The figure includes points representing the daily RMS values and smoothed curves based on a 14-day window Gaussian filter. During the early phase of the Repro3 period (before 2005), the differences between individual ACs and the combined solution were significant. However, after 2005, these differences became more stable. The AC with the highest RMS value stabilized below 6 mm, while most ACs stabilized around 3 mm. The GFZ solution performed best, with an RMS value of approximately 1.7 mm, indicating the highest stability. Additionally, some ACs showed clear seasonal fluctuations in their RMS differences, particularly GRG and the TUG. GRG exhibited the greatest amplitude of fluctuations, which may be related to its ZTD output frequency of once every 2 h compared to other ACs that output once every hour or every 5 min. The daily RMS for individual stations also showed similar trends, as illustrated in Figure 11.

Next, we conducted a station-specific analysis. The bias, RMS, and STD were calculated between the ZTD sequences and the combined sequences for each AC. The bias calculation results between the ZTD of each station and the combined ZTD for different ACs are shown in Figure 12. It can be seen that the stations in the ESA and GRG tend to have positive biases, indicating that the ZTD estimates from these two ACs are generally higher than the combined solution. Conversely, the stations in the TUG tend to have negative biases, suggesting that the ZTD estimates from the TUG are generally lower than the combined solution. For the other three ACs (CODE, GFZ, JPL), the stations’ biases are roughly balanced between positive and negative. The RMS and STD calculation results of the differences are shown in Figure 13 and Figure 14, respectively. In these figures, GRG also exhibits high RMS and STD values. Some stations exhibited high RMS and STD values across multiple ACs and were removed during processing. From the perspective of station distribution, RMS and STD values decrease as latitude increases. This may be related to the greater variability of ZTD in low-latitude regions, which increases measurement uncertainty.

Table 5. Bias, RMS, and STD statistical table of the differences between the individual AC ZTD solutions and the combined solution in mm.

Mean over All Stations	CODE	ESA	GFZ	GRG	JPL	TUG	Mean
Bias	−0.07	0.13	0.04	0.16	0.04	−0.10	0.03
RMS	3.61	3.12	1.70	5.50	2.00	2.19	3.02
STD	3.58	3.07	1.66	5.47	1.96	2.16	2.98

presents the specific values. The bias values for each AC relative to the combined solution are all close to zero. The RMS and STD values provide a measure of the overall error and variability, respectively. Among the ACs, the TUG exhibits the lowest RMS and STD values, indicating the highest consistency with the combined ZTD solution. Conversely, GRG shows the highest RMS and STD values, suggesting greater variability in its ZTD estimates.

4.5. Evaluation of the ZTD Combined Products with Respect to Independent Datasets

4.5.1. ERA5

When comparing GNSS ZTD with ERA5 ZTD, significant differences were observed among various analysis centers. By calculating the bias, RMS, and STD between the ZTD series of each AC and the ERA5 ZTD series, these differences were quantitatively evaluated. The mean values for all stations of each AC are summarized in

Table 6. Statistical table of the bias, RMS, and STD between GNSS ZTD and ERA5 ZTD.

Mean over All Stations	COD-ERA5	ESA-ERA5	GFZ-ERA5	GRG-ERA5	JPL-ERA5	TUG-ERA5	Com-ERA5	Mean
Bias	−2.00	−1.80	−1.86	−1.75	−1.88	−2.07	−1.95	−1.90
RMS	11.22	11.39	9.98	11.43	9.60	10.41	10.14	10.60
STD	10.75	10.89	9.49	11.00	9.10	9.87	9.63	10.10

. The results of the bias and STD are shown in Figure 15. From the figures, it can be observed that the biases of different ACs exhibit different distribution trends across various latitudes. However, overall, the GNSS ZTD values are generally slightly lower than the ERA5 ZTD values. All ACs show the highest STD values near the equator (0°), at approximately 15–20 mm, while these values decrease to around 5–10 mm at higher latitudes (around ±60° and beyond). This indicates that errors and random errors are larger in equatorial regions, while they are smaller and more stable in high-latitude regions. This pattern may reflect the impact of complex meteorological conditions and high-humidity environments on ZTD estimation accuracy in equatorial regions. The high-humidity conditions in these regions indicate an increased water vapor content, which would subsequently leads to greater variability in ZTD measurements. Moreover, rapidly evolving atmospheric patterns and the frequent occurrences of weather extremes, like convective storms and heavy precipitation, further complicate the accurate estimation of ZTD in equatorial regions.

Analyzing the statistical values in the table, the bias analysis shows that the ZTD values of all ACs and the combined series are slightly lower than those of ERA5. Among them, GRG has the smallest bias at −1.75 mm, while the TUG has the largest bias at −2.07 mm. Further analysis using RMS and STD reveals the magnitudes of total and random errors. The JPL performs best across all metrics, with an RMS of 9.60 mm and an STD of 9.10 mm, indicating the smallest total and random errors. In contrast, although GRG has a smaller bias, its RMS and STD values are 11.43 mm and 11.00 mm, respectively, indicating larger random errors. Overall, the quality of the combined series is superior to that of most individual analysis centers.

4.5.2. VLBI

The study selected a total of 13 VLBI stations co-located with IGS stations, as shown in Figure 16. The distances between the IGS and VLBI stations are shown in

Table 7. List of IGS and VLBI stations with approximate distance between them and height differences (dh is the altitude of the GNSS station minus the altitude of the VLBI station).

GNSS Station	VLBI Station	Distance (m)	Dh (m)	GNSS Station	VLBI Station	Distance (m)	Dh (m)
BADG	7382	96.7	−10.2	MATE	7243	58	−7.7
BRFT	7297	57.9	−1.4	MEDI	7230	62.8	−17.1
FAIR	7225	94.3	−13.1	NYA1	7331	106.4	−3.1
HOB2	7374	107	0.1	SHAO	7227	92	−7.4
HRAO	7232	163.8	−1.5	TSKB	7345	303	−17.5
HRAO	7378	274.8	5.1	ZECK	7381	65.8	−8.8
KOKB	7298	46.5	−9.2

. The bias, RMS, and STD between the AC/combined data and the VLBI data were calculated, as shown in

Table 8. Statistical table of the bias, RMS, and STD between GNSS ZTD and VLBI ZTD.

Mean over All Stations	CODE-VLBI	ESA-VLBI	GFZ-VLBI	GRG-VLBI	JPL-VLBI	TUG-VLBI	Com-VLBI	Mean
Bias	0.74	0.83	1.06	0.94	0.97	0.92	0.94	0.91
RMS	8.03	7.50	6.95	9.18	7.01	7.25	6.96	7.55
STD	7.79	7.24	6.63	8.94	6.74	6.92	6.65	7.27

. Based on the comparative analysis between the 13 VLBI stations co-located with IGS stations, the following conclusions can be drawn: all the ACs’ ZTD values were overestimated compared to the VLBI observations, exhibiting a certain degree of systematic bias. Among them, GFZ had the largest bias at 1.06 mm, while the CODE had the smallest bias at 0.74 mm, reflecting the differences between the different analysis centers. In terms of the RMS and STD metrics, GRG had the largest total error and random error, while GFZ and combined data had results closest to the VLBI observations.

4.5.3. Radiosonde

This study selected a total of 30 radiosonde stations collocated with IGS stations, as shown in Figure 16. The distances between the GNSS and radiosonde (RS) stations are shown in

Table 9. List of GPS and RS stations with approximate distance between them and height differences (dh is the altitude of the GNSS station minus the altitude of the RS station).

GNSS Station	RS Station	Distance (km)	Dh (m)	GNSS Station	RS Station	Distance (km)	Dh (m)
ALRT	CAM00071082	0.6	272.1	MAC1	ASM00094998	0.5	229.8
ANKR	TUM00017130	12.7	94.4	MAW1	AYM00089564	0.5	10.4
BRAZ	BRM00083378	10.8	56.4	MCIL	JAM00047991	0.1	30.2
BRMU	BDM00078016	1.8	39.4	NICO	CYM00017607	4.1	221.9
CAS1	AYM00089611	0.2	200.5	NRIL	RSM00023078	1.4	300.1
CHAT	NZM00093986	0.7	69.0	NYA1	SVM00001004	15.1	90.3
COCO	CKM00096996	0.1	0.4	PARC	CIM00085934	35.8	161.9
GLSV	UPM00033345	6.3	134.3	REYK	ICM00004018	0.2	227.9
HERT	UKM00003882	3.7	88.6	SCOR	GLM00004339	42.5	21.7
INVK	CAM00071957	1.2	159.8	SHAO	CHM00058362	10.4	537.6
IRKJ	RSM00030715	43.1	215.7	SOFI	BUM00015614	2.2	17.4
KIRI	KRM00091610	0.9	2.4	TEHN	IRM00040754	6.4	35.6
KOSG	NLM00006260	43.9	161.7	TSKB	JAM00047646	3.1	67.2
LHAZ	CHM00055591	3.0	17.7	URUM	CHM00051463	1.4	−19.8
ALRT	CAM00071082	0.1	130.9	WIND	WAM00068110	8.9	124.9

. The bias, RMS, and STD between the AC/combined data and the RS data were calculated, as shown in

Table 10. Statistical table of the bias, RMS, and STD between GNSS ZTD and radiosonde ZTD.

Mean over All Stations	CODE-Radio	ESA-Radio	GFZ-Radio	GRG-Radio	JPL-Radio	TUG-Radio	Com-Radio	Mean
Bias	−1.11	−0.82	−0.84	−1.22	−0.94	−1.13	−1.01	−1.01
RMS	12.29	12.85	11.26	11.54	11.16	12.00	11.37	11.78
STD	11.44	11.70	10.32	10.50	10.29	11.10	10.44	10.83

. From the bias analysis, the ZTD values of all ACs were lower than the RS measurements, indicating that GNSS tends to underestimate compared to RS. Among them, GRG had the largest bias at −1.22 mm, while the ESA had the smallest bias at −0.82 mm, reflecting systematic differences between the different ACs. The RMS and standard deviation analyses revealed the magnitudes of the total errors and random errors, respectively. In the RMS metric, the ESA had the largest total error at 12.85 mm, while the JPL had the smallest at 11.16 mm. In the standard deviation analysis, the ESA had the largest random error at 11.70 mm and the JPL had the smallest at 10.29 mm. Considering all the indicators, the JPL’s measurements were closest to the radiosonde data. Moreover, the quality of the combined series remains superior to that of most individual analysis centers.

5. Discussion

The advanced method developed in this study has great potential to facilitate the generation of a robust, long-term, and homogenous GNSS climate record to serve the climate community. This would not only foster a more holistic understanding of atmospheric phenomena captured through GNSS atmospheric monitoring techniques but also significantly improve the quality and uptake of GNSS climate records for extensive climatic analysis. Despite the progress made in this study, there are some limitations that should be acknowledged. For instance, there is a lack of in-depth research on the trends in ZTD at individual stations, and the homogenization of long-term ZTD time series needs to be addressed. Additionally, while data quality has improved, there is still room for further enhancement in terms of precision and accuracy. Considering these limitations, we will conduct further research in the future.

6. Conclusions

In this study, an advanced method has been developed to effectively amalgamate six sets of GNSS-derived ZTD time series obtained from different IGS ACs. This method enabled the creation of a comprehensive and integrated GNSS climate record, encapsulating a robust dataset that reflects the combined expertise and data contributions of the involved ACs.

Specifically, during the use of this method, the data processing strategies and station conditions of the six ACs were initially compared and analyzed. Subsequently, the formal errors in the data were systematically analyzed, accompanied by the implementation of a rigorous quality control process. An analysis was conducted on the ZTD formal error values for six ACs. The Q1, Q2, Q3, IQR, and upper atypical limits (Q3 + 3 × IQR) were calculated for each AC. Therefore, it was determined that the formal error of an individual parameter must be less than 10 mm as an acceptable tolerance, and data exceeding this limit were removed.

The data were than combined using the method proposed in this paper; the ZTD values from the stations were combined, resulting in a smoother, more consistent, and reliable ZTD combined time series. After the generation of the final combined ZTD time series, its quality was evaluated by firstly comparing with the individual time series contributed by the six ACs, offering a comparative perspective on its consistency and accuracy. The mean bias of the ACs’ time series with respect to the combined time series is 0.03 mm and the mean root mean square is 3.02 mm. Among the six ACs, the TUG had the lowest RMS and STD values, indicating the highest consistency with the combined ZTD solution. Conversely, GRG exhibited the highest RMS and STD values, suggesting greater variability in its ZTD estimates, which may be related to its ZTD output frequency of once every 2 h. After the combination, each station’s data show an increase in the time span and an improvement in data integrity compared to the data submitted by each AC.

Finally, other external references, such as VLBI, radiosonde, and ERA5 data were employed to validate the quality and reliability of the data. When compared against the different reference data sources, the combined solution performs better than most individual analysis centers, indicating that the combined solution has high reliability. As a result, the advanced method proposed in this study, along with the generated high-quality dataset, holds significant potential for advancing GNSS atmospheric sensing.