Quality control of SIRGAS ZTD products

: The SIRGAS - CON network currently has more than 450 continuous GNSS stations, and it is used for geodetic purposes. In atmospheric studies, it is used for ionospheric monitoring and for the estimation of zenith tropospheric delays ( ZTDs ) . From the Neutral Atmosphere Analysis Center of SIRGAS, Centro de Ingeniería Mendoza Argentina, the ﬁ nal tropospheric products of this network are generated after several stages of quality controls and ﬁ ltering, in order to be published on a daily basis in the o ﬃ cial website of SIRGAS, since 2014 ( https://sirgas.ipgh. org/en/products/tropospheric - delays ) . These products arise from adjusting the solutions estimated by di ﬀ erent SIRGAS analysis centers. Prior to the combination, a quality control of the individual solutions is carried out, based on the precision estimator of each parameter and an internal control of each solution with respect to the combined value. In this work, we show the quality control process of the inputs, the selected tolerance and its justi ﬁ cation. The internal consistency analysis of tropospheric parameters for a period of 7 years was carried out. We also exposed the improvements in the estimation of tropospheric para meters implemented during 2021 and its impact in the generation of the ﬁ nal ZTD products ( in 99% of the stations the mean standard deviation of ZTD is less than 1 mm ) .


Introduction
The network of GNSS continuous stations in the Americas, called SIRGAS-CON (Figure 1), currently has more than 450 stations. From them, 86 stations have been removed, some of them have observations for more than 25 years and some have remained inactive for some periods. The majority of the stations have received maintenance and an increase in the regional station densification is observed in different American countries.
Since its origin (approximately around the year 2000), the SIRGAS-CON network has been used for geodetic purposes (Drewes et al. 2000, Seemueller et al. 2002, Sánchez and Drewes 2020, and in atmospheric studies, it has been used for the determination of total electron content on each site and in the modeling of this variable in the region (Brunini et al. 2004a, 2004b, 2008, Azpilicueta et al. 2005, 2012. In the last 10 years, efforts have been made in the estimation of tropospheric delay parameters, as a contribution to the knowledge of the neutral atmosphere. The zenith tropospheric delays (ZTDs) estimated in GNSS processing strongly depend on the instantaneous integrated water vapor (IWV) content of the atmospheric column at a site. Therefore, it is possible to use the ZTD to retrieve the IWV at each GNSS station, increasing the amount of data and improving its coverage. The accuracy results are equivalent to that expected from direct observational techniques, such as radiosondes and microwave radiometers (Bonafoni et al. 2013, Calori et al. 2016. Camisay et al. (2019) and Mackern et al. (2020) have exposed some contributions of the SIRGAS-CON network to the knowledge of this atmospheric variable.
The first ZTD estimations in the SIRGAS network were carried out, every 2 h, in the GNSS stations of South America, processed by the CIMA Analysis Center (Centro de Ingeniería, Mendoza, Argentina) (Calori et al. 2013(Calori et al. , 2015. Since 2014, SIRGAS implemented a strategy to combine the tropospheric parameters estimated by different analysis centers (ACs) with the Bernese GNSS Software V.5.2, in their weekly routine. The ZTDs obtained from the combination are considered final SIRGAS products and are published weekly on the SIRGAS website (https:// sirgas.ipgh.org/en/products/tropospheric-delays/). These values have been validated with respect to final International GNSS Service (IGS) ZTD products (at 15 IGS stations) and with respect to computed ZTD at radiosonde stations (10 sites in South America and the Caribbean), between 2014 and 2018 (Mackern et al. 2020).
During the 8 years since the implementation of the ZTD combination strategy, the number of stations in the SIRGAS-CON network has increased with a rate of 21% (considering only the active stations from 2014 to 2021). The number of ACs that participate in the weekly processing has also increased ( Table 1), and since 2021 the tropospheric parameters estimated with the GAMIT software have been incorporated by the GNA (Instituto Geográfico Nacional de Argentina) AC. The later ZTD data sets had not been included in this article.
The aim of this work is to present the results of the consistency analysis of the tropospheric parameters, prior to the incorporation of the GNA ZTDs (estimated with the GAMIT software). Section 2 describes the data used and the tolerance adopted to consider or reject the tropospheric parameters in the final combination and its justification. The methodology adopted to carry out the internal quality control of the individual tropospheric parameters is also explained.
Section 3 shows, as first results, the statistical indicators of applying quality control over 7 years of ZTD data (2014-2020) and the issues following its implementation. In Section 4, an analysis of the excluded parameters and the identification of the causes are described. The actions applied, during the year 2021, in three SIRGAS ACs, are mentioned.
Finally, in Section 5 the results obtained after applying the changes in three ACs and the precisions achieved during the last 6 months of 2021 are shown. As statistical indicators of improvement, we used the mean and standard deviation (SD) of ZTD parameters and the amount of parameters excluded by each AC.

Data used
ZTDs provided by seven SIRGAS ACs (Table 1), for a period of 7 years (2014-2020) were used.
Under the coordination of Working Group I (WG I) of SIRGAS, each AC carries out weekly processing of a group of assigned SIRGAS stations. In this process, a tropospheric product is estimated as a sum of three components, a priori model, estimated parameter and horizontal gradient parameters (equation (1)) with their SD (Dach et al. 2015), every hour (in some cases every 2 h).
The distribution of stations is done in such a way that each station is processed by at least three ACs and its distribution is coordinated by the WG I chief. Based on the AC guidelines, each station should be considered in the processing if, at least, 12 h of observation are available and if its observation file (RINEX file) is available in the corresponding Data Center at the starting time of the processing. This occurs with a delay of 14 days after the end of the observation week. Availability of final IGS orbits and earth orientation parameters is required. Therefore, each AC performs its individual processing with a delay between 15 and 21 days.
The Neutral Atmosphere Analysis Center, CIMA, analyzes the combination of the tropospheric parameters of different ACs, on 24th day after the end of the week of interest. All individual solutions delivered, which pass in the corresponding controls, participate in the combination. These individual parameters are the basic input for the calculation of the final SIRGAS tropospheric products, for which a rigorous quality control is carried out on them, prior to the combination, called "admission" (Figure 2).

Tolerance on the ZTD i precision estimator
Each AC delivers the ZTD ij (i corresponds to the AC and j to each epoch) calculated with its SD, σ ZTDij , for every station assigned to it. The calculation of the final tropospheric parameters of SIRGAS is performed as a combination per epoch, within a process that is carried out weekly. The first step in this process is "admission," for which a threshold was established on the SD corresponding to each parameter.
The mean upper atypical limit value (0.0023 m) was considered as the limit of outliers and it was decided to take as a tolerance ten times of this upper limit. In this way, we established as an admission tolerance that the SD of the individual parameter (σ ZTDij ) has to be less than 0.02 m (Figure 2), similar to that applied in EUREF (Pacione et al. 2011(Pacione et al. , 2017. Those parameters that do not satisfy this tolerance are excluded from the combination process.

Data redundancy
Each station is processed by at least three ACs, according to the guidelines of SIRGAS WG I. Therefore, it has been established as a condition for the combination, to have at  least two individual solutions per tropospheric parameter with their SD minor to the tolerance (Figure 2).

ZTD quality control, period 2014-2020
To verify the contribution of the individual solutions (ZTD ij ) in the combined final SIRGAS products, an analysis, from 2014 to 2020, was carried out on the total of parameters that passed the "admission" tolerance ( Figure 3). It was identified that on average for the 7 years analyzed, 31.3% of the stations had 3 or more tropospheric parameter solutions, 50.3% had only 2 solutions and 18.3% only had 1 solution. It was also found that in the last year (2020) the number of stations with only 1 and 2 parameter solutions analyzed had increased, which led us to suppose that the number of parameters rejected had also increased.
As part of the evaluation, two main causes of these shortages are identified. One cause is the delay in uploading the observation RINEX files from some stations within the maximum time established in the guidelines. Because of this, in the processing of some ACs (those who start processing earlier) the estimation of the ZTD of such stations is omitted. The second cause is the number of excluded tropospheric parameters with SD greater than the tolerance.
Regarding the first cause, from the SIRGAS WG I, it is emphasized in the Data Center guides on the importance of publicizing the observation files within the established times. Regarding the second cause, a specific analysis was carried out on the rejected parameters, which is detailed in the following section.

Analysis of the number of parameters out of tolerance
A statistical analysis was implemented on the parameter rejections made during a period of 6 months, from September 27, 2020 (GPS week 2,125) to March 27, 2021 (GPS week 2,150), trying to identify systematic effects on some SIRGAS ACs. The number of excluded parameters was counted weekly (Figure 4), for each AC (7 in the analyzed period). It was identified that the largest amount corresponded to the solutions of USC, CHL and URY, compared to what was expected (no more than 30 parameters excluded per week).
From WG I, a new redistribution of stations was carried out, establishing the subnets to be processed by each AC. This was implemented from GPS week 2,151. The calculation routines in the three mentioned ACs (USC, CHL and URY) were also checked over.
Prior to the calculation of week 2,151, some changes ( Table 3) were implemented in the USC processing routine, and it was observed that the magnitudes of the SD of the ZTDs were significantly reduced. The number of parameters with sigma greater than 0.02 m was also considerably reduced. Based on the improvements in the estimation of tropospheric parameters by USC for week 2,151, it was decided to keep on the new processing configuration. From SIRGAS WG I, CHL and URY were   In a daily processing, the parameters are estimated, with all coordinates loosely constrained (sigma a priori 1 m), in a 1-h interval and horizontal gradient parameters Tropospheric parameters are estimated through a 3-stage sequence: 1) Coordinates and tropospheric parameters are preliminary estimated (datum is defined from the IGS14 frame for selected IGS stations in each cluster). These are necessary in ambiguity resolution strategies 2) Tropospheric parameters are estimated with an interval of 1 h, after ambiguities resolution, per cluster 3) Final tropospheric parameters are estimated in the cluster combination per day, with all coordinates constrained, in a 1-h interval and horizontal gradient parameters advised to implement the same changes applied in the USC center into their corresponding routines. This was made by URY from GPS week 2,161 and by CHL from GPS week 2,166.
In order to evaluate the results achieved, the analyzed period was divided into four stages for better interpretation: Pre-improvement (GPS weeks 2,125-2,150), USC improvement (GPS weeks 2,151-2,160), URY improvement (GPS weeks 2,161-2,165) and CHL improvement (GPS weeks 2,166-2,181). Figure 5 shows the average number of excluded parameters for the period, for each AC. It can be seen how the implementation of the new strategy optimizes the contribution of ZTD from USC (greater number of parameters that are included with lower SD) from week 2,151. Likewise, the improvement in URY from week 2,161 and in CHL since week 2,166. In each station and for each epoch, after calculating the weighted estimator, the residual of each individual estimator is computed: A quality control is performed and those residuals greater than three times their SD are investigated as potential outliers. If that is verified, the final combination is excluded.
In order to evaluate the external precision of the individual AC solutions, the residuals of each parameter were analyzed with respect to the final products obtained from the combination, for each epoch. This analysis was carried out for a period of 10 months, from January 1 to October 30, 2021. This period was selected in such a way that it included 3 months prior to the change implemented in the USC strategy, and also involves 3 months during which the changes were applied in the three ACs mentioned in the previous item and four subsequent months of monitoring.
It was possible to observe that the residuals of the three ACs were greater than (in absolute value) the general average, before implementing the changes and reduced significantly after the improvement was applied. The residuals calculated at four stations (AMCO, BAPA, CORD and AGGO) are shown as an example.  Likewise, Figure 7 shows the ZTD residuals from the BAPA station (processed by GNA, IBG, URY and USC). Residuals of USC solutions between 0.03 and −0.03 m (before week 2,151) are identified and then the values are reduced between ±0.015 m. The residuals of URY solutions between 0.07 and -0.06 m (prior to week 2,151) are identified. Since week 2,151, this station does not have URY residuals because BAPA was not assigned to this AC after the aforementioned network redistribution. Figure 8 shows the individual AC residuals, corresponding to the CORD station (processed by DGF, GNA, IBG and URY), during the analyzed period. It can be observed how the residuals of the individual URY solutions appear from week 2,151, when its processing was assigned to this AC. The residuals of URY solutions between 0.08 and −0.07 m (from week 2,151 to week 2,161) are identified and then the values decrease notably (between ±0.01 m), from week 2,161 when the changes in the URY processing strategy were implemented.
Finally, the ZTD residuals in the AGGO station (processed by DGF, USC, GNA and CHL) are shown in Figure 9. Note the reduction in magnitude of residuals from week 2,151 when USC implemented the improvement and the reduction in magnitude of CHL residuals from week 2,166, since this AC applied the changes in its strategy.
6 Results and final conclusions 6.1 Increase in the number of individual ZTDs contributing to the final solution A monitoring was carried out for 3 months after the changes, during which it was verified that the exclusions due to SDs greater than 0.02 m were reduced in the three ACs. In average values, the number of rejected parameters in USC were reduced to 99.9%, in URY and in CHL were reduced to 99.6% with respect to the amounts that were excluded before the implemented routine ( Figure 5). Once the changes have been implemented in the three ACs (from week 2,166), a general agreement is  observed, compared to the exclusions of the other ACs ( Figure 10). The number of parameters excluded by each AC kept below 30 per week (sum of the 7 days, and the total number of corresponding assigned stations).
This directly influences the increase in the amount of ZTD itj added to the combination. Figure 11 shows, on average for the year 2021, that 75% of the stations participating in the combination with at least three individual ZTDs (optimal requirement) and for 98% of the stations it is possible to obtain the final ZTD using the combination of at least two individual ZTDs (minimum requirement).
A significant improvement is observed if we compare this result with the corresponding results obtained in the

ZTD precision estimation
The precision estimation of the final SIRGAS tropospheric parameters (ZTD SIRGAS_tj ) calculated in 2021 is reliable considering the previous controls. It is important to consider that the ZTD SIRGAS_tj have been combined from three or more individual solutions for 75% of the stations and have had at least one control from the combination of two solutions for 23.5% more stations. The ZTD precision is calculated from the mean annual SD obtained for each station. Such mean value was estimated, considering the SD of each parameter (1 parameter per hour) for every day of the year. Based on the total number of stations, the distribution of mean SD was carried out in five ranges ( Figure 12).
For a total of 440 stations with ZTD calculated in the year 2021 (after the improvement in the ZTD estimation), 48% (212) of the stations have a mean SD less than 0.1 mm, 48% (213) between 0.1 and 0.59 mm, 2.5% (11) between 0.6 and 1 mm and the remaining 1% (4 stations) have a mean SD less than 3 mm. This agrees more reliably with the precision estimate of the final SIRGAS tropospheric parameters (the mean SD is less than 1 mm in more than 84% of the estimated parameters) published by Mackern et al. (2020).  As a final conclusion, it can be said that the mean SD of the final products of the SIRGAS ZTD combination is less than 1 mm in 98.5% of the stations, according to the EUREF and IGS network solutions (Pacione et al. 2011, Brockmann et al. 2006. This precision confirms that these final SIRGAS tropospheric products can be used as a reference for scientific applications (e.g., validation of near real-time tropospheric products in the Latin American region) and their extended time series could be used to monitor trends and variability in the neutral atmosphere.