Assessment of SIRGAS-CON tropospheric products using ERA5 and IGS

2.1 ZTD derived from GNSS observations

In order to achieve accurate GNSS measurements, it is necessary to account for the tropospheric delays inherent to them. To account for these effects, over the last few years, the analysis centers from GNSS networks have been computing and providing the ZTD for their GNSS stations’ observations systematically and operationally. Processing methods may be different for each network and center and are well established, reaching an accuracy of a few millimeters (Niell et al. 2001; Pacione et al. 2011).

The ZTD is the sum of the hydrostatic and wet delays in the zenith direction, namely the Zenith Hydrostatic Delay (ZHD) and the Zenith Wet Delay (ZWD), respectively. It is modeled according to equation (1), transforming the Slant Total Delay (STD), which is the tropospheric delay suffered by the GNSS signal when crossing the atmosphere at different angles of elevation (E), using the mapping functions mf h and mf w , corresponding to the hydrostatic and wet components, respectively.

SIRGAS-CON tropospheric products are estimated for all its stations in highly reliable combined weekly solutions (Mackern et al. 2020). The processing methodology adopted in the estimation of GNSS-derived ZTD has a significant impact on the accuracy of the estimated time series (Stępniak et al. 2022). SIRGAS-CON analysis centers use Bernese GNSS Software v5.2 (BSW52, Dach et al. 2015), and ZTDs are calculated with the final IGS products (orbits and earth rotation parameters, ERP). The operational SIRGAS-CON products are estimated according to the methodology and parameters shown in Table 1 (Mackern et al. 2020).

SIRGAS-CON tropospheric products are available with a latency of approximately 30 days on the SIRGAS webpage, in SINEX TRO files with 1 h temporal resolution since January 2014. The products are freely available for download on the webpage ftp://ftp.sirgas.org/pub/gps/SIRGAS-ZPD (SIRGAS 2022).

Nowadays, the SIRGAS-CON network has a total of more than 560 stations spread all over Latin America and also some of them in the USA, Antarctica, Europe and Africa, as shown in Figure 1, with a historical series of 7 years of data. For this work, all ZTDs available from January 2014 to December 2020 have been considered for analysis, forming a subset of 467 stations, of which 63 are common to the IGS network.

Figure 1

SIRGAS-CON stations represented in blue points and SIRGAS-CON stations common to IGS in red points.

2.2 ERA5 ZTD computation

The NWMs are physical and mathematical models that represent the atmospheric variables, computed from equations that describe their dynamics, using a combination of data from different sources.

ECMWF provides different NWMs, such as the Operational model (ECMWF-Op), with 6 h temporal resolution and 0.125° × 0.125° (16 km × 16 km) spatial resolution (Miller et al. 2010); ERA-Interim (Dee et al. 2011), which has a 6 h temporal resolution and 0.75° × 0.75° (80 km × 80 km) spatial resolution; and ERA5 (Hersbach et al. 2018), currently the best and most recent reanalysis model, with 1 h temporal resolution and 0.25° × 0.25° (30 km × 30 km) spatial resolution. ERA5 replaced ERA-Interim, which has been discontinued since August 31, 2019 (ECMWF 2021). ERA5 is available in single-level variables at the surface, some of which are representative of the total atmospheric column (integrated variables), or in 137 vertical levels (3D) for some variables, from the surface up to 0.01 hPa (around an altitude of 80 km), and also interpolated to standard levels, such as pressure levels, at 37 levels with constant pressure, from the surface (1,000 hPa) to 1 hPa (around an altitude of 45–50 km; Vieira et al. 2019a).

ZTD can be obtained from NWM as the sum of ZHD and ZWD. Since the NWM has a different spatial and temporal resolution from the GNSS observations, they need to be interpolated in space and time to derive values for the same instant and location as the GNSS measurements, making them directly comparable. The NWM parameters used to calculate the ZHD and the ZWD are given at the model orography level, making it necessary to reduce the values first computed at the orography level to the height of each GNSS station.

The ZHD corresponds to approximately 90% of the ZTD (Fernandes et al. 2014) and basically depends on the pressure of atmospheric dry gases, which has an almost constant distribution throughout the atmosphere, enabling its computation with an error of less than 1 cm. The ZHD can be computed from ERA5 sea level pressure fields using the modified Saastamoinen model (Davis et al. 1985), according to equation (2):

In equation (2), the ZHD results in meters, φ is the geodetic latitude, h _s is the surface height over the geoid, p _s is the surface pressure, calculated from the surface pressure at sea level p ₀, by equation (3), which denotes the pressure variation with altitude (Hopfield 1969).

In equation (3), T _m is the mean temperature value in Kelvin (K) of the layer between heights h ₀ and h _s and can be estimated as the mean value between T ₀ and T _s temperatures at heights h ₀ and h _s, respectively, from T ₀ values at mean sea level from Global Pressure and Temperature models (Böhm et al. 2007). R ₀ is the specific dry air constant and g is the mean gravity, resulting from equation (4) as follows:

The ZWD corresponds to approximately 10% of the ZTD (Fernandes et al. 2014, 2021), and it can be computed from NWM single-layer fields of total column water vapor (TCWV), given in mm or the equivalent, kg/m², and 2-m temperature (T ₀), which is the air temperature close to the surface, given in Kelvin (K), as presented in equation (5) (Bevis et al. 1992; Bevis et al. 1994):

T _m is the mean temperature of the troposphere, which can be modeled from T ₀, also in K, according to equation (6) (Mendes et al. 2000; Mendes 1999):

The ZWD obtained by equations (5) and (6) is given at the model orography level. The altitude dependence of ZWD is much more complex than ZHD because it depends essentially on the water vapor, which has non-constant and much more variable distribution throughout the atmosphere, unlike the dry gases. Nevertheless, ZWD height reduction can be modeled according to equation (7) (Kouba 2008) as follows:

Even though it is widely used, Kouba reduction is not recommended for height differences greater than 1,000 m (Kouba 2008) and needs some improvements in its modeling (Vieira et al. 2019b).

3.1 Comparison of SIRGAS-CON ZTD with ERA5 ZTD

For each of the 467 SIRGAS-CON stations with available ZTD between January 2014 and December 2020, the mean and standard deviation of the SIRGAS-CON ZTD and the ERA5 ZTD differences were calculated.

The means of the ZTD differences for the various stations are in a range from 0 to 130 mm in absolute value, from which 97 stations have an absolute mean value greater than 1 cm, 44 stations greater than 2 cm and 19 stations greater than 3 cm, representing 21, 9 and 4% of the total number of analyzed stations, respectively. This indicates that 79% of the stations have a mean difference with respect to ERA5 with an absolute value of less than 1 cm. The mean value of the set of mean differences is 8 mm.

Figure 2 shows the boxplot of the means of the ZTD differences, with the median very close to zero and a large number of outliers over the lower and upper whiskers, with values of −13.6 and 15.1 mm, respectively.

Figure 2

Boxplot of the mean of differences between SIRGAS-CON ZTD and ERA5 ZTD for each station.

Selecting only the outliers from the boxplot for spatial visualization, as shown in Figure 3, it is noted that the occurrence zone is mainly in the Andes Mountains and other mountainous areas of Central America and the USA. It is possible that the model is not suitable for these zones, as the model’s orography is a smoothed version of a Digital Elevation Model (DEM) and may not fit well in mountainous areas because its resolution may not be able to map highly variable topography gradients (Bento et al. 2017). Therefore, only stations that are not in these zones, stations with orthometric height less than 1,000 m and stations for which the difference between orthometric height and the model orography level is also less than 1,000 m are considered for the subsequent analysis since ZWD reduction for height differences greater than 1,000 m can cause gross errors.

Figure 3

107 stations located in zones where the NWM may not fit the relief variations well due to its resolution.

After eliminating the stations located in areas where the comparison with the model may be inadequate (107 in total), 360 stations are left for analysis. For this subset, the means of the differences are in a range of 0–14 mm in absolute value, with 15 stations having values greater than 1 cm, representing 4% of the total analyzed stations. This shows that about 96% of the stations have sub-centimeter precision, regarding the mean of differences with the ERA5.

Regarding the standard deviations of the ZTD differences, considering the original set of 467 stations, the range is from 5 to 50 mm and the mean of standard deviations is 16 mm, with 55 stations greater than 2 cm and eight stations greater than 3 cm, representing 12% and 2% of the total, respectively.

For the subset of 360 stations selected to analyze the means of the ZTD differences, the standard deviations are in a range of 5–40 mm (only one station has a standard deviation of 50 mm in the original set, so the intervals are practically the same), and the mean of standard deviations remains 16 mm, with 36 stations greater than 2 cm and three stations greater than 3 cm, representing 10% and 1% of the subset, respectively. This analysis shows that the problem verified for the means does not happen for the standard deviation. Figure 4 shows the behavior of the differences between the ABRA (panel a) and COEC (panel b) stations, which are examples of stations that have a high bias, but a stable time series with a low standard deviation of 11 and 9 mm, respectively.

Figure 4

SIRGAS-CON ZTD and ERA5 ZTD differences for stations (a) ABRA, with a mean of 59 mm and standard deviation of 11 mm, and (b) COEC, with a mean of 37 mm and standard deviation of 9 mm.

The root mean square (RMS) of the ZTD differences considering the original set of 467 stations have a range from 5 to 140 mm and the mean value of the set of RMS of 20 mm, with 130 stations greater than 2 cm and 35 stations greater than 3 cm, representing 28 and 8% of the total, respectively.

For the subset of 360 stations, the RMS is in a range from 5 to 40 mm and the mean value of the set of RMS decreases to 17 mm, with 55 stations greater than 2 cm and five stations greater than 3 cm, representing 15% and 1% of the subset, respectively. It is noted that the percentage of stations with RMS greater than 2 cm on the selected set decreases to almost half that on the original set.

These statistics are partly in agreement with the results of similar studies, such as Fernandes et al. (2013a,b), which present for a set of 11 IGS/EUREF Permanent Network (EPN) common stations the mean and standard deviation of the ZTD differences between GNSS and ERA-Interim, the preceding ECMWF reanalysis model, of 4 and 12 mm, respectively. In the same study, the standard deviation of the ZTD differences has a range from 8 to 14 mm. It is important to highlight that the set of stations analyzed in this work is much larger and that the GNSS stations of the previous work are located in the temperate zone, while most of the SIRGAS-CON stations are located in the tropical zone, where the ZWD is much more variable and can reach much higher values, which may cause the much larger range of the standard deviation.

The variability of the standard deviation of the ZTD differences between SIRGAS-CON and ERA5 is shown on the map in Figure 5, for all 467 stations with available ZTD. The green points represent the stations with a standard deviation between 10 and 20 mm, the interval that contains most of the stations.

Figure 5

SIRGAS-CON stations with colors representing the standard deviation of the differences between SIRGAS-CON ZTD and ERA5 ZTD for each station.

The blue points represent stations with a standard deviation lower than 10 mm. The regions that have most of the stations with standard deviations in the lowest range, represented by the blue points, are the extreme south of South America, in the Patagonia region and the Antarctic continent. The yellow and red points represent stations with the highest standard deviations. There is a greater concentration of these stations in the equatorial zone, and some are in other regions, such as Uruguay. The pattern of the standard deviation shown on the map is in accordance with the natural ZWD variability, which is greater in the tropics and lower in areas closer to the poles. However, those stations with a high standard deviation in Uruguay are not supposed to have those values, since they are located in the temperate zone, which indicates that these stations may have some problems with their GNSS ZTD solutions.

Although the ZHD corresponds to 90% of the ZTD, it has a standard deviation of up to 3 cm, and it is much more constant throughout the atmosphere and easier to model than the ZWD, which corresponds only to 10% of the ZTD and has a standard deviation value of up to 15 cm (about five times larger than that of ZHD), evidencing that it is much more variable and difficult to model, with values ranging from 0 to 50 cm (Fernandes et al. 2021). Therefore, most of the ZTD variability comes from ZWD.

The stations located in latitudes inside or near the polar zones have a lower and less variable ZTD, since the water vapor in the atmosphere in this region is also lower and less variable, while the opposite happens in the tropics and close to the equator, where the water vapor becomes larger and more variable. Figure 6 shows the OHI2 (panel a) and LIMN (panel b) stations. The OHI2 station is located in Antarctica (63° 19′ 15.960″ S, 57° 54′ 4.680 W″) and has a continuous time series over the years, and it is a good example of a stable station that has ZTD differences between both sources agreeing very well, with a mean very close to zero and low dispersion (standard deviation of 6 mm). In contrast, the LIMN station, located in Costa Rica, is in the tropical zone (9° 59′ 35.125″ N, 83° 01′ 34.930″ W) and has also a continuous and stable time series, with a mean that is also close to zero but high dispersion (standard deviation of 20 mm). This exemplifies the variation of water vapor as the latitude changes.

Figure 6

SIRGAS-CON ZTD and ERA5 ZTD differences for stations (a) OHI2 and (b) LIMN.

The regions closer to the equator are warmer and wetter and have the highest ZWD absolute values (Fernandes et al. 2013a,b). Southwards of the tropics, as the latitudes increase, the amount of water vapor in the atmosphere decreases as well as its variability, so it is expected that the standard deviation of the corresponding ZTD differences also decreases.

On the map in Figure 5, it is possible to observe some stations with a high standard deviation greater than 20 mm in Uruguay and along the west of South America. This can indicate the presence of instrumental errors, processing errors or even errors in station coordinates (the coordinates used are from SIRGAS-CON log files, also available on the SIRGAS webpage [SIRGAS 2022]), since they are located in the temperate zone and the standard deviations are expected to be lower.

The stations UYPT and UYSO, as shown in Figure 7a and b, respectively, have a stable series for most of the period, with an expected dispersion caused by the seasonal signal, characterized by an increase in the ZTD variability in the summer and a decrease in the winter, but there is a period in the time series when the dispersion deviates from their normal pattern and remains with very large ZTD differences until the end of the time series. These stations have a high standard deviation due to these periods of high dispersion, but the mean of differences is not affected, having values close to zero, since the dispersion of the differences, even in the periods when it is highly dispersed, is often almost symmetrical.

Figure 7

SIRGAS-CON ZTD and ERA5 ZTD differences for stations (a) UYPT and (b) UYSO.

Some stations have short periods of high dispersion of differences, as for BCH1 and CXEC stations, shown in Figure 8a and b, respectively. Despite this inconsistency, the standard deviation is 13 mm for the BCH1 station and 15 mm for the CXEC station, below the mean of standard deviations of the whole set. This happens because there is a low percentage of GNSS ZTD observations in specific short periods of the time series causing the high dispersion. In these cases, graphical analysis proves to be helpful to detect these types of inconsistencies. Another problem detected by graphical analysis which is sometimes not detected by general statistics is a discontinuity, detected by a tendency change in the means of the ZTD differences. These discontinuities can be caused by various reasons, such as changes in data acquisition equipment (antenna, receiver, etc.), changes in the processing method or errors in station coordinates. An example of this case is the AUTF station, which has for its differences with ERA5 a mean of −1 mm and a standard deviation of 11 mm. These statistics can be considered a good indicator but it is possible to see in Figure 9a that at the beginning of 2017 there is a sudden decrease in the values of the differences when compared to the period just before this occurrence. There was also an increase at the beginning of 2018, but the mean did not return to the value it had before the first decrease. The BEJA station, as shown in Figure 9b, is another example of a discontinuity. For this analysis, a LOESS (locally weighted non-parametric regression fitting using a second-order polynomial) smoothing filter (Marsh 2022) has been used to highlight the discontinuities.

Figure 8

SIRGAS-CON ZTD and ERA5 ZTD differences for stations (a) BCH1 and (b) CXEC.

Figure 9

SIRGAS-CON ZTD and ERA5 ZTD differences for stations (a) AUTF and (b) BEJA.

3.2 Intercomparison between SIRGAS-CON, ERA5 and IGS ZTD

There are 63 stations in SIRGAS-CON that are also part of the IGS global network. To support the understanding of the differences between SIRGAS-CON ZTD and ERA5 ZTD analyzed in the previous section, a comparison between SIRGAS-CON ZTD and the corresponding IGS ZTD was performed.

Regarding these differences, considering the whole set of 63 stations common to both networks, the means are in the range from −1 to 6 mm. The standard deviations are in the range from 3 to 16 mm, and the RMS are at most 1 mm larger than the standard deviation, except for MDO1, PALM and PALV stations, which are 2 mm larger. The differences between both GNSS solutions are mostly small, but six stations have standard deviations and RMS of ZTD differences greater than 1 cm, shown in Table 2, which may indicate the presence of errors in one of the solutions. The mean values of the means, standard deviations and RMS for this set of stations are 1, 7 and 7 mm, respectively.

These statistics show a good agreement between both GNSS solutions. Results are in agreement with other studies, as in Fernandes et al. (2013a,b), which present 3 mm in mean and 5 mm in the standard deviation of the ZTD differences in a comparison between different GNSS solutions using 51 IGS and EPN common stations spanning 2002–2009. The same work presents absolute means between 2 and 4 mm and standard deviations between 2 and 7 mm, considering each station individually. Mackern et al. (2020) present 2 mm in mean and 7 mm in RMS using 15 IGS and SIRGAS-CON common stations, spanning 2014 to 2019.

The differences between GNSS and model ZTD are usually greater than those between two different GNSS solutions (Fernandes et al. 2013a,b). One of the reasons for the first differences being greater than the second is the need to reduce the ZHD and ZWD derived from the model orography level to the station level, which has an associated error, whereas there is no height reduction in the comparison between different solutions for the same GNSS station from different networks. Figure 10 shows the ZTD differences for station PIE1 when comparing data from the three different sources, showing that these differences have different biases. The bias is associated with the height reduction, while the range of the differences and the seasonal signal are associated with the low capacity of the NWM to describe the peaks of the seasonal signal since the NWM is a smoothed representation of the ZTD field.

Figure 10

SIRGAS-CON ZTD and ERA5 ZTD (blue) and SIRGAS-CON ZTD and IGS ZTD (cyan) differences for station PIE1.

A case that presents an inconsistency between both GNSS solutions is the station MGUE, with a period of high dispersion in the differences between SIRGAS-CON ZTD and IGS ZTD from the end of 2016 to the first half of 2018, shown in Figure 11a. An intercomparison of both solutions with NWM helps to find out if these instabilities come from SIRGAS-CON or IGS data. The differences for this station between IGS ZTD and ERA5 ZTD and between SIRGAS-CON ZTD and ERA5 ZTD are shown in Figure 11b and Figure 11c, respectively.

Figure 11

Differences between (a) SIRGAS-CON ZTD and IGS ZTD; (b) IGS ZTD and ERA5 ZTD; and (c) SIRGAS-CON ZTD and ERA5 ZTD for station MGUE.

When comparing the three plots, it is clear that the high dispersion in the differences for the aforementioned period is present in the SIRGAS-CON ZTD differences with ERA5 ZTD. It is concluded that the data dispersion comes from SIRGAS-CON ZTD. Another case of instability is the SALU station, as shown in Figure 12. The three plots analyzed together suggest that the discontinuity comes from IGS ZTD.

Figure 12

Differences between (a) SIRGAS-CON ZTD and IGS ZTD; (b) IGS ZTD and ERA5 ZTD; and (c) SIRGAS-CON ZTD and ERA5 ZTD for station SALU.

This analysis shows that the comparison of ZTD from different types of sources may have different characteristics and magnitudes. The intercomparison helps to understand these characteristics and come to the conclusion that a comparison just between two different sources may be insufficient.