Characteristic differences in tropospheric delay between Nevada Geodetic Laboratory products and NWM ray-tracing

Numerical weather models (NWMs) are important data sources for space geodetic techniques. Additionally, the Global Navigation Satellite System (GNSS) provides many observations to continuously improve and enhance the NWM. Existing comparative analysis experiments on NWM tropospheric and GNSS tropospheric delays suffer from being conducted in highly specific regions with limited spatial coverage; furthermore, the length of time for the experiment is too short for analyzing seasonal characteristics, and the insufficient number of stations limits spatial density, making it difficult to obtain the equipment-dependent distribution characteristics. After strict quality control and data preprocessing, we have calculated and compared the bias and standard deviation of tropospheric delay for approximately 7000 selected Nevada Geodetic Laboratory (NGL) GNSS stations in 2020 with the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5) hourly ray-traced tropospheric delay for the same group of stations. Characterizations in time, space, and linkage to receivers and antennas reveal positive biases of approximately 4 mm in the NGL zenith tropospheric delay (ZTD) relative to the NWM ZTD over most of the globe; moreover, there is a seasonal amplitude reaching 6 mm in the bias, and an antenna-related mean bias of approximately 1.6 mm in the NGL tropospheric delay. The obtained results can be used to provide a priori tropospheric delays with appropriate uncertainties; additionally, they can be applied to assess the suitability of using NWMs for real-time positioning solutions.

Ray-tracing via NWMs is a highly accurate method for obtaining tropospheric delays (Landskron and Böhm 2018b); it is inseparable in both the direct calculation of the slant path delay (SPD) and the indirect mapping delay from the zenith direction down to the slant direction through mapping functions (Zhou et al. 2020). High-resolution NWM raytracing tropospheric delay has been used to evaluate and validate GNSS tropospheric delay (Andrei and Chen 2009;Li et al. 2015). In addition, comparing GNSS tropospheric and NWM tropospheric delays is seemingly useful for detecting possible increases in volcanic activity (Cegla et al. 2022). However, these advantages do not indicate that NWM tropospheric delays are more accurate than GNSS tropospheric delays. Douša et al. (2016) showed that GNSS products could provide more detailed structures in the atmosphere than the state-of-the-art numerical weather models, demonstrating that the two techniques have very high similarity and complementarity. GNSS tropospheric delays are also usually used to assess and compare the differences between different NWMs (Li et al. 2015;Zhou et al. 2020). In recent years, some scholars have compared and analyzed the two types of tropospheric delays; however, due to the limited amount of data, these experiments are usually performed on only a few to a few dozens of stations in the region, and the data duration covers only a few days to several months (Kačmařík et al. 2017;Hordyniec et al. 2018;Lasota et al. 2020). Other problems arise because the data qualities of the stations are not strictly controlled (Andrei and Chen 2009;Elsobeiey 2020;Zhou et al. 2020) (stations with large amounts of data missing are not excluded, the threshold of integrity rate is set too low, the mean value of all stations with uneven distribution is used to represent the global level, etc. However, these neglects cannot be ignored in millimeter-level accuracy comparisons). These problems lead to large differences between the results of these studies, and the spatiotemporal and equipment-dependent distribution characteristics of the differences between the two types of data cannot be reliably analyzed and discussed.
In this research, we calculate and analyze the bias and standard deviation (std) values from the tropospheric delay data of approximately 7000 GNSS stations that are strictly quality-controlled by Nevada Geodetic Laboratory (NGL) solutions and ray-traced tropospheric delay data of the same locations to compare the differences in the two most accurate types of tropospheric delays. The results can be used to provide a priori tropospheric delays with their uncertainties and to assess the suitability of using NWMs for realtime positioning solutions. We introduce the data sources and preprocessing methods in the following section. In the next section, we analyze and discuss the spatiotemporal and equipment-dependent distribution characteristics of bias and std in the following section. Finally, the analysis results are summarized, and some perspectives are provided.

Data and methods
This section introduces the NGL tropospheric products and their preprocessing methods; furthermore, the NWM used for ray-tracing is introduced. In addition, the bias and std values of the two tropospheric delays are calculated and analyzed.

ERA5 ray-traced tropospheric products
Ray-tracing by NWMs, which tracks as approximately as possible to the real travel path and considers the effects of the whole atmosphere, is one of the most accurate ways of obtaining tropospheric delays (Landskron and Böhm 2018b). European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5) is the latest generation of ECMWF reanalysis for the global climate and weather over the past 4 to 7 decades (https:// cds. clima te. coper nicus. eu/ cdsapp# !/ datas et/ reana lysis-era5-press ure-levels? tab= overv iew, 2022.4). ERA5 replaces the ERA-Interim reanalysis, and ERA5 zenith tropospheric delay (ZTD) performs better than ERA-Interim ZTD at all scales (Zhou et al. 2020). In this research, the ECMWF ERA5 hourly products with the horizontal resolutions of 1° × 1° are used for ray-tracing to derive zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD), and the ZTD is the sum of ZHD and ZWD.

NGL tropospheric products
The NGL has produced tropospheric products from more than 16 thousand sites for over 34 million station-days since November 5, 2017 (Blewitt et al. 2018). On November 28, 2019, the NGL reprocessed and updated the data products with improved models, including the VMF1 (Vienna Mapping Functions 1) mapping function and nominal troposphere, and improved JPL Repro 3 orbits, and the latest global reference frame IGS14. The number of updated stations exceeded 19 thousand, and the data exceeded 43 million station days in total (http:// geode sy. unr. edu/, 2022.4). In this research, NGL ZTD and NGL ZWD were derived from NGL products, and NGL ZHD was equal to NGL ZTD minus NGL ZWD. Figure 1 shows approximately 7000 NGL stations selected for this research. These stations included data from more than 80 receiver models and more than 200 antenna models in 2020. In addition, these stations had good data quality, and neither the receiver nor the antenna was replaced throughout 2020.

Data preprocessing and analysis methods
Unlike reanalysis data, such as ERA5, direct measurements, such as GNSS, always inevitably have missing epochs of data, with the times and numbers of missing epochs being unequal across separate stations. The error caused by this missing data cannot be ignored in highprecision data processing. To reduce the impacts of these missing data as much as possible, we selected stations with data integrity rates (ratios of the number of epochs with data to the total number of epochs) exceeding 95% (the reason 95% was chosen as the threshold is explained later). In addition, we used the annual and semi-annual fitting model (Chen et al. 2020) to complete the missing epochs at all the selected stations. Figure 2 shows the ZTD time series of station VIUA as an example. This station was chosen here because its data integrity rate just met the threshold we set, and almost all of the missing epochs were at one end (which was the worst of all the missing cases), suggesting that none of the other stations was worse than this station regarding data completeness.
The top left panel of Fig. 2 shows the fitting of the ZTD time series of the NGL, where the blue line is the fitted line, and the green points are the model values. The top right panel of Fig. 2 shows the comparison of the NGL ZTD and NWM ZTD, where the blue line is the linear fit line, and the green dots are the comparisons of the model value of the NGL ZTD and its corresponding epoch NWM ZTD. From the linear fitting results, the slope was 1.0142, and the coefficients of the two time series were 0.98, indicating that the two types of data had very high agreements and correlations.
The bottom panels of Fig. 2 present the residual series of the NGL-NWM and their absolute values. According to the residual series, the mean of residuals was very close to zeromean and the mean of the residual series is 4 mm, indicating a very small systematic deviation between NGL ZTD and NWM ZTD. The short-term variation in the ZTD difference series was disordered and it was almost impossible to find a pattern. However, the NGL-NWM residual series after the data preprocessing remained essentially approximately zeromean, and its distribution range was within the distribution of residuals throughout the year (setting a lower threshold would not guarantee this). Therefore, the thresholds for our integrity rates were chosen reasonably well, and the model we used outperforms methods that only use data from epochs near the missing epoch for interpolation or extrapolation.
The time series of the absolute values of the NGL-NWM biases showed the apparent seasonal signals. To eliminate the influences of systematic errors, the absolute value here was the absolute value of the NGL-NWM bias time series minus its mean. The fitting amplitude and mean absolute value of the fitting residual of each station data were used as indicators to analyze the temporal distribution characteristics. Fitting results for all stations are shown and analyzed in the following section on temporal distribution characteristics. Figure 3 shows the mean bias (left panels) and standard deviation (right panels) of the bias series between NWM and NGL after preprocessing by the above method. For better Fig. 1 Global distribution of the selected NGL stations and corresponding receiver types visualization, we use different color bars for biases exceeding zero and less than zero; additionally, different scales for the panels are used.

Global distribution of bias and std
The bias of the ZHD between the NGL and NWM is shown in the top left panel of Fig. 3. NGL tropospheric products do not include ZHD; instead, they use the VMF1/ NWM grid as model input, suggesting that the bias here is the difference between the different NWMs used by NGL and ERA5. Over 75% of the stations show negative ZHD bias values; however, in high-altitude regions, such as western South America and North America, northern and eastern Africa, northern India, and parts of Antarctica, the ZHD biases exhibit larger positive values. These results suggest that the NWM used by VMF1 may overestimate the ZHD in the high-altitude region relative to ERA5. To further confirm this phenomenon, we compared the VMF1 ZHD with the VMF3 (Vienna Mapping Functions 3, the latest version of the VMF model) ZHD. We find that there is indeed a centimeter-level difference between the two models in some high-altitude regions. This phenomenon, which has issues when applied to the above model over regions with highly variable topography, has also been found by the research of Zhang et al. (2021), and our calculations are much closer to the VMF3.
In GNSS data processing, the hydrostatic delay is typically obtained through models, and the wet delay (normally ZWD) is estimated as an unknown parameter. The errors of hydrostatic delay are absorbed by the wet delay. This absorption of errors is clearly displayed by the ZWD bias. In regions with large positive and negative values of the ZHD bias, the ZWD bias behaves as a correction to ZHD in the opposite direction approaching zero.
The biases of the ZTD between the NGL and NWM are approximately − 10 to 15 mm, and in most regions, the ZTD bias is positive; only in a few stations is the bias negative, suggesting that in most regions, the ZTD calculated by GNSS is overestimated relative to the NWM ZTD. In addition, the bias is generally larger in low-latitude regions than in high-latitude regions. In summary, due to the absorption of ZHD bias, applications requiring GNSS ZWD, such as GNSS meteorology, require high-accuracy ZHD models. However, the ZTD bias is less influenced by the ZHD bias, suggesting that applications based on GNSS ZTD, such as GNSS ZTD-based tropospheric delay models, are only slightly affected.
Std represents the degree of aggregation of the residual series toward the mean bias and the seasonal differences between NGL and NWM. The larger the std is, the greater the variation in the difference between the NGL tropospheric delay and the NWM tropospheric delay. The ZHD std is less than 5 mm for most regions, especially in the lower latitudes where there are few large values. In contrast, many more scattered large values of ZHD std exist in the mid-and high-latitude regions. In addition, large values of ZHD std occur in the high-altitude regions of Antarctica and in eastern South America.
Stations with larger std values have larger absolute bias values; stations with larger absolute bias values do not always have larger std values. ZWD std is approximately 0-25 mm over most of the region, exhibiting a significant negative correlation with latitude (decreasing with increasing latitude). The larger ZWD std occurs mostly in southeastern North America and in north-central South America. ZTD std is regarded as a concatenation of the ZWD std and ZHD std values. Since the ZHD std is relatively smaller than the ZWD std, the ZTD and ZWD std values are highly similar. Stations in the polar regions have smaller ZTD std values, and their ZWD std values are closer to those of ZHD std.

Results and analysis
In this section, the temporal distribution characteristics of the bias and std between the NWM tropospheric and NGL tropospheric delays are statistically analyzed. Then, the latitude and height dependences of bias and std are derived and analyzed. Finally, the tropospheric biases associated with the receiver/antenna category are extracted. Figure 4 shows the amplitude (left panels, hereafter referred to as the amplitude) and the mean of the absolute residuals (right panels, hereafter referred to as the mean) of the absolute values of bias between the NGL and NWM tropospheric delays under the annual and semi-annual fitting models of all the selected stations around the globe. For the fitted amplitude and the mean absolute values of the residuals, refer to the bottom right panel of Fig. 2. Note that different scales are used for different panels to better show the differences between regions.

Temporal distribution characteristics of bias and std
The ZHD amplitude is mostly less than 1 mm. However, due to the presence of many scattered stations with larger values at mid-to-high latitudes and high altitudes, the overall amplitude is 0-2 mm. The ZWD amplitude ranges from 0 to 6 mm, there is no significant latitude-related pattern, and almost all values are very small in the high-latitude Stations with larger values are mainly distributed in the continental centers of North and South America, the western sides of the Japanese Islands and the center of Australia; additionally, some stations near the Black Sea and the Caribbean Sea have larger amplitudes. We suggest that this regional dependence may be a combination of latitudinal dependence and the heterogeneity of the data assimilated by ERA5. The ZTD amplitude is broadly considered as a concatenation of the ZWD and ZHD amplitudes, combining all the characteristics of the numerical distribution of the two. Since the ZHD amplitude is much smaller than the ZWD amplitude, the ZTD bias is very close to the ZWD bias. By combining the amplitudes of the three delays, it is seen that only north-central Europe falls within the region where all three types of amplitudes are relatively small; this region has the flattest tropospheric variability during the year.
The ZHD mean range is the smallest (0-4 mm), and its distribution is very similar to the ZHD amplitude. However, in mid-and high-latitude regions, the ZHD mean differences between stations are more pronounced than those of the ZHD amplitude. The mean of the fitting residual absolute value of the ZWD residual absolute value ranges from 0 to 12 mm, showing a significant negative latitude correlation (decreases with increasing latitude), especially in the low-latitude regions of the American continent, where the value reaches its maximum. The ZTD mean is the same as the ZWD mean, both ranging from 0 to 12 mm. The ZTD mean is the union of ZWD and ZHD; however, since ZHD is much smaller than ZTD, ZTD only reflects the distribution characteristics of ZWD.

Spatial distribution characteristics of bias and std
To obtain more detailed information on the relationship between bias and std with station latitude and height, we grouped the bias and std values of the stations according to their latitude and height; the results are presented in Fig. 5. Note that GNSS stations are only distributed over land, and the distributions of GNSS stations are very heterogeneous (both horizontally and vertically). This phenomenon results in a relatively small number of stations in the high-latitude and high-altitude groups. If there are stations with outlying values in these groupings, the outliers may increase or decrease the values of the whole grouping. In addition, neither the height grouping nor the latitude grouping removes the influences of the other side; thus, the analysis here is conducted from the perspective of the overall trend.
The top panels of Fig. 5 demonstrate the latitude dependence of bias and std. The ZTD bias appears positive in all latitude groupings and decreases with increasing latitude. In low-latitude regions, the ZTD bias typically exceeds 5 mm; in the middle-and high-latitude regions, the ZTD bias is below 4 mm. The ZHD bias is negative in most latitude groups, and its magnitude and absolute value have no significant latitude correlation. However, in the middle-and low-latitude regions, the absolute value of the ZHD bias is usually smaller than the ZTD bias; in high-latitude regions, the absolute value of the ZHD bias becomes larger than the ZTD bias. Since ZWD absorbs the error of ZHD, the ZWD bias should be analyzed with the ZTD and ZHD biases. The figure shows that the ZWD bias is generally similar to the ZTD bias; in the region where the absolute value of the ZHD bias is large, to compensate for the ZHD error, the ZWD bias begins to increase.
Unlike bias, std exhibits a more significant latitude correlation. The figure shows that ZTD std and ZWD std are very similar, both of which decrease with increasing latitude. Additionally, the differences between the values in the same latitude group are very small; the difference is only larger in the high-latitude group. Numerically, the ZTD std and ZWD std values exceed 5 mm in all latitudinal regions and 15 mm in the lower-latitudes regions. The change in ZHD std with latitude is opposite to the previous two; that is, it increases with increasing latitude. However, the value of ZHD std is the smallest, and it is also less than 5 mm in high-latitude regions and usually only 2-3 mm in low-latitude regions.
The bottom panels of Fig. 5 show the height dependence characteristics of bias and std. The values of the ZTD bias in all height groups are approximately 4 mm; these values are positively correlated with height and increase slowly with increasing height. The ZHD bias has positive and negative values in regions with heights of less than 2 km, and these values are less than ± 2 mm; in regions with heights exceeding 2 km, all of the values are negative, and the values decrease as the height increases. These results indicate that in regions with heights exceeding 2 km, the agreement between ERA5 and the NWM used by NGL to calculate ZHD decreases rapidly with increasing height. This decrease in consistency results in rapid increases in ZWD biases in regions with heights exceeding 2 km.
When calculating the height dependence of std, we find that the tropospheric delay std decreases with increasing height; this phenomenon is obviously unreasonable because it contradicts the fact that the water vapor fluctuates more near the surface of the earth. Since the tropospheric delay decreases rapidly with height, we replace std with relative std (the ratio of std at each station to the annual mean of the tropospheric delay) to eliminate the impact. The results show that the relative std of the ZTD decreases as the height increases. Additionally, the relative std values of the ZHD and ZWD increase as the height increases, and the relative std of the ZTD changes more slowly than the other two. The ZTD, ZHD and ZWD relative std values are approximately 0.5%, 0.3% and 15%, respectively. The relative std of ZWD is more than an order of magnitude larger than the other two, indicating that the intra-annual variation in bias is mainly contributed to by ZWD. Ejigu et al. (2019) found an average deviation of 1.8 mm in ZWD in European experiments, which may be caused by different receiver equipment. The experiments of Ejigu et al. (2019) are based on double-difference (DD) network solutions; however, Stępniak et al. (2022) found that DD solutions contain more numerous and larger ZTD outliers. The NGL solution follows the precise point positioning (PPP) strategy, and the NGL gradients are more tightly constrained; these make NGL products more suitable for exploring equipment-dependent biases. In addition, the large number of high-density NGL stations allows us to analyze the correlation of tropospheric delays with receivers and antennas.

Linking of tropospheric delay bias to receivers and antennas
To minimize the effect of spatial distance on the analysis, we selected 1411 pairs of stations by finding the horizontal distance and vertical distance between two stations. Each Average bias and (relative) std of the NGL tropospheric delay and NWM tropospheric delay at different latitudes (top panels) and heights (bottom panels). The relative std is the percentage of std to the annual mean of the tropospheric delay at each station. Note that ZTD and ZHD correspond to the left Y-axis, and ZWD corresponds to the right Y-axis pair of stations has a horizontal distance of less than 20 km and a vertical distance of less than 1 m (in theory, under this strict constraint, the difference should be almost zero). These station pairs are shown in Fig. 6. Note that each pair of stations is represented as a single point on the map due to the close proximity.
These station pairs are divided into four categories according to the type of antenna and receiver: same antenna with the same receiver (SS; red dots, 435 pairs), same antenna with a different receiver (SD; green dots, 164 pairs), different antenna with the same receiver (DS; blue dots, 263 pairs) and different antenna with a different receiver (DD; magenta dots, 549 pairs). As seen from the graph, the first three types of station pairs are heavily distributed in the three regions with the largest number of stations-the USA, Europe and Japan-making the results of the analysis more fair and reliable.
The horizontal and vertical distances have been constrained to reduce the effects of station spatial distances, and the differences in tropospheric parameters between stations (defined as single-difference) of the NGL are free from the influences of satellite-related errors. However, the location difference of the station pair will still inevitably introduce tropospheric deviations. To further reduce these spatially related differences, we further calculate the differences between the single-differenced NGL and NWM tropospheric parameters, defined as double-difference below. Through the above processing, the double-differenced tropospheric parameters are actually affected mostly by the equipment.
We categorize the statistics of the double-differenced tropospheric parameters of each station pair by the receiver and antenna types. The results show that there are no direct patterns related to a specific type of receiver or a specific type of antenna. We further analyze the bias and std values by station pair according to the SS, SD, DS, and DD groups defined above; the results are shown in Fig. 7. Note that only ZTD and ZWD are considered here, as ZHD is not a direct output of GNSS estimates and is not related to the receiver or antenna. Figure 7 shows the error bar exceeds the corresponding column in all receiver and antenna combinations, which is expected because it is the result of the annual and global averages. However, the large error bar values indicate that the variations in the biases between combinations cannot be ignored. The figure shows the antenna is the main reason for the difference in tropospheric delay rather than the receiver; this phenomenon results in a ZTD/ZWD bias of approximately 1.6 mm, while the receiver only results in a bias of approximately 0.2 mm. The mean std behaves similarly to the mean bias, suggesting that the differences in antennas lead to greater tropospheric delay biases and to greater seasonal bias fluctuations.

Conclusions and perspectives
In this study, GNSS tropospheric delay data from approximately 7000 stations in the NGL for 2020 and tropospheric delay data calculated using ray-tracing (NWM from ECMWF ERA5 hourly data) at the same location, are compared and analyzed. Spatiotemporal and Fig. 6 Global distribution of station pairs of different types of receivers and antenna combinations. SS, SD, DS, and DD represent the same antenna for the same receiver, the same antenna for different receivers, different antennas for the same receiver, and different antennas for different receivers, respectively Fig. 7 Average bias (left panel) and std (right panel) of different types of receiver and antenna combinations. SS, SD, DS, and DD represent the same antenna for the same receiver, the same antenna for different receivers, different antennas for the same receiver, and different antennas for different receivers, respectively equipment-dependent distribution characteristics of bias and std values between these two types of tropospheric delay data are derived regarding time, spatial, receiver and antenna dependences.
Our results show that the ZTD biases of the two types of data exhibit positive values in most regions, with values of approximately 4 mm, showing a significant latitudinal correlation. The ZHD bias caused by the differences between the NWM adopted by VMF1 and ERA5 will rapidly increase (reaching the centimeter level) in high-altitude regions exceeding 2 km, further increasing the ZWD bias. This increased bias may lead to the precipitable water vapor (PWV) values obtained by the GNSS method in high-altitude regions being inaccurate while having little effect on the ZTD bias.
The annual amplitudes of the absolute values of the tropospheric delay biases are 0-6 mm. These results can be used to set the initial tropospheric variance values in the PPP to accelerate the positioning convergence time. In addition, we find that the receiver antenna contributes approximately 1.6 mm or so of tropospheric delay bias, which PCO/PCV accuracies may cause.
The results in this research help comprehensively understand the differences between the two types of tropospheric delays with the highest accuracies thus far; the results also determine the spatiotemporal and equipment-dependent distribution characteristics of the bias. In addition, the large amount of station data used in this research helps to provide a priori values with appropriate uncertainties for positioning; these values help to decorrelate other parameters, such as station coordinates and receiving clocks, to achieve better accuracy and reliability characteristics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.