A Comprehensive Evaluation of Key Tropospheric Parameters from ERA5 and MERRA-2 Reanalysis Products Using Radiosonde Data and GNSS Measurements

: Tropospheric delay is a major error source in the Global Navigation Satellite System (GNSS), and the weighted mean temperature (T m ) is a key parameter in precipitable water vapor (PWV) retrieval. Although reanalysis products like the National Centers for Environmental Prediction (NCEP) and the European Center for Medium-Range Weather Forecasts (ECMWF) Re-Analysis-Interim (ERA-Interim) data have been used to calculate and model the tropospheric delay, T m , and PWV, the limitations of the temporal and spatial resolutions of the reanalysis data have affected their performance. The release of the ﬁfth-generation accurate global atmospheric reanalysis (ERA5) and the second Modern-Era Retrospective analysis for Research and Applications (MERRA-2) provide the opportunity to overcome these limitations. The performances of the zenith tropospheric delay (ZTD), zenith wet delay (ZWD), T m , and zenith hydrostatic delay (ZHD) of ERA5 and MERRA-2 data from 2016 to 2017 were evaluated in this work using GNSS ZTD and radiosonde data over the globe. Taking GNSS ZTD as a reference, the ZTD calculated from MERRA-2 and ERA5 pressure-level data were evaluated in temporal and spatial scales, with an annual mean bias and root mean square (RMS) of 2.3 and 10.9 mm for ERA5 and 4.5 and 13.1 mm for MERRA-2, respectively. Compared to radiosonde data, the ZHD, ZWD, and T m derived from ERA5 and MERRA-2 data were also evaluated on temporal and spatial scales, with annual mean bias values of 1.1, 1.7 mm, and 0.14 K for ERA5 and 0.5, 4.8 mm, and –0.08 K for MERRA-2, respectively. Meanwhile, the annual mean RMS was 4.5, 10.5 mm, and 1.03 K for ERA5 and 4.4, 13.6 mm, and 1.17 K for MERRA-2, respectively. Tropospheric parameters derived from MERRA-2 and ERA5, with improved temporal and spatial resolutions, can provide a reference for GNSS positioning and PWV retrieval.


Introduction
Tropospheric delay and weighted mean temperature (T m ) are key tropospheric parameters for high-accuracy positioning/navigation and precipitable water vapor (PWV) retrieval in space geodetic techniques [1][2][3][4][5][6]. Hence, this should be considered carefully. The Saastamoinen model can calculate with high accuracy the zenith hydrostatic delay (ZHD) based on the condition of obtaining accurate meteorological parameters. However, not all Global Navigation Satellite System (GNSS) stations are equipped with meteorological sensors. Meanwhile, the distribution of GNSS stations is uneven; most of them are distributed on land and few are distributed in the ocean. The high-accuracy prior zenith wet delay (ZWD) can effectively reduce the convergence time and improve the positioning precision in precise point positioning (PPP) [7][8][9]. Meanwhile, tropospheric parameter information of the ZTD, ZWD, ZHD, and T m derived from MERRA-2 and ERA5 reanalysis data. Finally, we summarize the conclusions of this work.

Data Description
The main data used in this study were the ERA5 reanalysis data from ECMWF; MERRA-2 reanalysis data from the National Aeronautics and Space Administration (NASA); GNSS observations collected from the International GNSS Service (IGS); radiosonde data from the websites of the University of Wyoming, both from 2016 to 2017. The performances of key tropospheric parameters derived from MERRA-2 and ERA5 reanalysis data were analyzed by comparing the GNSS observations and radiosonde data on a global scale.

Reanalysis Products
MERRA-2 and ERA5 are the latest reanalysis products from NASA and ECMWF, respectively. The surface-level and pressure-level data of MERRA-2 (https://goldsmr4 .gesdisc.eosdis.nasa.gov/data/MERRA2, accessed on 1 June 2021) and ERA5 (https:// www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5, accessed on 1 June 2021) were used in the current study. The MERRA-2 hourly surface-level data have a horizontal resolution of 0.625 • × 0.5 • including temperature, geopotential height, surface pressure, and specific humidity. The ERA5 hourly surface-level data have a horizontal resolution of 0.25 • × 0.25 • , containing temperature, geopotential, surface pressure, and surface specific humidity. The geopotential data are divided by gravity acceleration to obtain geopotential heights. The MERRA-2 6 hourly pressure-level data have a horizontal resolution of 0.625 • × 0.5 • and a vertical resolution of 42 levels from 1000 to 0.1 hPa and comprise temperature, geopotential height, and specific humidity. Compared to the MERRA-2 pressure-level data, the ERA5 pressure-level data have a temporal resolution of 1 h, the horizontal resolution of 0.25 • × 0.25 • , and vertical resolution of 37 levels from 1000 to 1 hPa involving temperature, geopotential, and specific humidity.

GNSS Data
IGS can offer ZTD products with a temporal resolution of 5 min and accuracy of 4 mm [29] (https://cddis.nasa.gov/archive/gnss/products/troposphere/zpd/, accessed on 1 June 2021). We obtained 6 hourly and hourly ZTD products from IGS to evaluate the ZTD from reanalysis products. The GNSS ZTD were collected at 290 GNSS stations in the IGS as shown in Figure 1.

Radiosonde Data
Radiosonde data can be used as a reference to evaluate the accuracies of ZWD and T m derived from reanalysis data [30][31][32]. In addition, radiosonde data can also provide measured meteorological parameters such as water vapor pressure, atmospheric pressure, and temperature. We can use these meteorological parameters to calculate ZHD and ZTD using the Saastamoinen model, and which are used as the reference to validate the performances of ZHD and ZTD from reanalysis products. The radiosonde data from the websites of the University of Wyoming can be accessed freely (http://weather.uwyo.edu/ upperair/sounding.html, accessed on 1 June 2021). Figure 1 shows the distribution of 545 radiosonde stations over the globe. The Saastamoinen model method was used to calculate the ZHD, ZTD, and ZWD from ERA5 and MERRA-2 surface-level data, and the integration method was used to calculate ZTD, ZWD, ZHD, and T m from ERA5 and MERRA-2 pressure-level data. After deriving the ZTD by the integration method from pressure-level data, the Saastamoinen model was used to calculate the ZTD above the pressure level and the meteorological data of the top level was used as input values in the model.
The ZTD, ZWD, T m , and ZHD values from ERA5 and MERRA-2 pressure-level data and the formulas are as follows [33]: T m = s (e/T)ds s (e/T 2 )ds where N w the wet refraction, N the total refraction, e denotes the water vapor pressure, and T denotes the atmospheric temperature. Then, the N w , N, and e were calculated using the following formulas [34,35]: e ≈ qP/0.622 (5) where P denotes the atmospheric pressure, q represents the specific humidity, k 3 = 377600 K 2 /hPa, k 2 = 64.79 K/hPa, k 2 = 22.97 K/hPa, and k 1 = 77.604 K/hPa [21,34]. The Saastamoinen model was used to calculate the ZTD, ZWD, and ZHD based on the ERA5 and MERRA-2 surface-level data, and the formulas are as follows [36]: where H denotes the height (m) of the station, T 0 denotes the surface temperature (K), ϕ represents the latitude (radian), and P 0 denotes the surface pressure (hPa).

Calculating ZTD, ZHD, ZWD, and T m Using ERA5 and MERRA-2 Data at Radiosonde and GNSS Stations
We can calculate the ZTD, ZWD, ZHD, and T m values of grid points by the above method. However, the radiosonde or GNSS stations and the ERA5 and MERRA-2 grid points are usually not collocated and have different elevations. The altitude difference between radiosonde or GNSS and the grid point can lead to large errors in evaluating the performances of ZTD, ZWD, ZHD, and T m derived from reanalysis data [37]. It is necessary to eliminate or reduce the influence of this altitude difference. Therefore, the ZTD, ZWD, ZHD, and T m values of grid points in this study were at the height of radiosonde and GNSS stations. It can be ensured that the height of stations was consistent with the height of the grid points, which significantly eliminates or reduces the influence of altitude differences in evaluating the accuracies of these tropospheric parameters.
Meteorological parameters at the same height need to be obtained before calculating the ZTD, ZWD, ZHD, and T m values of grid points at the height of radiosonde and GNSS stations. Notably, the different reference planes of the grid height, radiosonde height, and GNSS height should be unified before interpolating meteorological parameters [38]. A detailed introduction to calculating the temperature, pressure, and specific humidity of grid points at a station's height using pressure-level data can be found in Huang's paper [39].
The method for calculating the meteorological parameters at the height of radiosonde and GNSS stations using reanalysis surface-level data is different from that of using reanalysis pressure-level data. The surface-level data from ERA5 and MERRA-2 data do not directly provide water vapor pressure values. Therefore, the grid points of water vapor pressure values need to be calculated, firstly, by using Formula (5). The formulas for calculating the pressure, temperature, and water vapor pressure values of the grid points at a station's height is as follows [20,28]: e = e 0 (P/P 0 ) λ+1 (12) where λ denotes the water vapor pressure decrease factor, which can be obtained by using the GPT2w model; T 0 represents the temperature at the grid point height H 0 ; e 0 denotes the water vapor pressure at the same height; P 0 is the pressure at the grid height; H denotes the height of the radiosonde or GNSS site. The ZTD, ZHD, ZWD, and T m values of the grid points at the height of the radiosonde and GNSS stations can be calculated by the above method, whereas the corresponding values of these parameters at the radiosonde and GNSS stations were obtained by the inverse distance weighting (IDW) [40].
To access the performance of the ZTD, ZWD, ZHD, and T m derived from MERRA-2 and ERA5 reanalysis data, we compared the tropospheric parameters calculated by reanalysis data with those derived from the radiosonde data and GNSS observation. Then, the bias and RMS were selected to measure their performance, which could be computed via the following equations: where V i m and V i R are the values from the reanalysis data and reference data, respectively, and N is the number of samples.

Comparisons between ERA5/MERRA-2 ZTD and GNSS ZTD
The hourly ZTDs from 290 GNSS stations from 2016 to 2017 as a reference were employed to evaluate the accuracy of ERA5 and MERRA-2 surface-and pressure-level ZTD. For each GNSS station, we computed the annual mean bias and RMS values of ZTD. The results are shown in Figures 2 and 3 and Table 1. Table 1. Statistics on the bias and RMS for ERA5 and MERRA-2 ZTD from surface-and pressure-level data compared to GNSS ZTD from 2016 to 2017.

MERRA-2 ERA5
Pressure-ZTD Surface-ZTD Pressure-ZTD Surface-ZTD bias (mm)   The RMS values for MERRA-2 ZTD from pressure-level data compared to GNSS ZTD were less than 16 mm for 82% of the 290 GNSS stations. We also found that the RMS values between ERA5 pressure-level ZTD and GNSS ZTD were less than 16 mm for 90% of the 290 GNSS stations. Meanwhile, the annual mean bias and RMS were 4.5 and 13.1 mm for MERRA-2 pressure-level ZTD and -4.8 and 27.9 mm for MERRA-2 surfacelevel ZTD, respectively. The annual mean bias and RMS of ERA5 ZTD from pressure-level data with respect to GNSS ZTD were 2.3 and 10.9 mm, respectively, and those of ERA5 ZTD from surface-level data were 3.8 and 31.6 mm, respectively. These illustrate that the ZTD calculated from pressure-level data had superior consistencies with GNSS ZTD than surface-level ZTD. The reason for this phenomenon was that the error of the Saastamoinen model method was 10-30 mm greater than that of the integration method [37]. Furthermore, the performance of ZTD derived from ERA5 pressure-level data was slightly improved over MERRA-2 pressure-level data, with an RMS of 10.9 and 13.1 mm for ERA5 and MERRA-2, respectively. This small RMS difference between ERA5 and MERRA-2 may be caused by differences in data assimilation and physical models [41]. Thus, the following comparison and analysis were based on ERA5 and MERRA-2 pressure-level ZTD.

RMS
To analyze the seasonal variation features of ZTD derived from ERA5 and MERRA-2 data, the monthly average bias and RMS values of ERA5 ZTD and MERRA-2 ZTD compared with GNSS ZTD at all GNSS stations from 2016 to 2017 were calculated. There are 222 GNSS sites in the Northern Hemisphere and 68 GNSS sites in the Southern Hemisphere. The bias and RMS means of all stations in four seasons, namely, spring (i.e., March, April, and May), summer (i.e., June, July, and August), autumn (i.e., September, October, and November), and winter (December, January, and February) were also summarized, and the results are shown in Table 2. In the Northern Hemisphere, the biases of MERRA-2 ZTD were positive in different seasons and show apparent seasonal trends that have a larger bias in the summer and a smaller bias in the winter. Meanwhile, the seasonal variation in RMS of ZTD in the Northern Hemisphere was obvious. The RMS of MERRA-2 ZTD in the summer and winter were 16.4 and 9.6 mm, respectively, and those of ERA5 ZTD in summer and winter were 13.1 and 8.5 mm, respectively. These features are related to the complicated and changeable climate during summer when the content of atmospheric water vapor is abundant and changes violently, which leads to larger errors of MERRA-2 and ERA5 ZTD in the summer [21,23]. Seasonal variations in RMS of ZTD in the Southern Hemisphere are opposite to that in the Northern Hemisphere owing to the contrary seasons and climatic conditions [42]. In the Southern Hemisphere, the RMS of MERRA-2 ZTD in summer and winter were 12.2 and 16.1 mm, respectively, and those of ERA5 ZTD in the summer and winter were 10.8 and 13.1 mm, respectively. Although seasonal variations clearly affected the performance of ZTD derived from reanalysis data, the ERA5 ZTD and MERRA-2 ZTD still had high precision, with a smaller bias and RMS in all seasons. To analyze the daily variations of MERRA-2 and ERA5 ZTD, the times series of MERRA-2 and ERA5 ZTD bias and RMS in eight GNSS stations are shown in Figures 4 and 5. The details of eight GNSS stations are shown in Table 3. Overall, the ZTD derived from MERRA-2 and ERA5 data in the thu2 and maw1 stations in high latitudes showed better precision than that in the shao, lhaz, cnmr, wgtn, mgue, and samo stations in middle and low latitudes. This is attributed to the stable variation of meteorological parameters in high latitudes, which reduces the uncertainty of ZTD calculated from MERRA-2 and ERA5 data. It is evident that the RMS of ZTD in the shao and mgue stations show a distinct seasonal trend. The main reason is that the shao station in the Northern Hemisphere is located in a temperate monsoon climate and close to the ocean, and the abundant and changeable water vapor in summer will lead to a large error in the calculation of ZTD. While the mgue station in the Southern Hemisphere is located in the temperate grassland climate, the water vapor content is lower in summer and increases in winter, which causes fewer errors in the summer and large errors in the winter. The altitude of the lhaz station is above 3000 m and the climate is dry, which leads to small errors in ZTD calculation.  Figure 6 shows the variations in the latitude of the annual mean bias and RMS values of MERRA-2 ZTD and ERA5 ZTD compared with GNSS ZTD for all stations from 2016 to 2017. The variations in the biases of ZTD with latitude were not evident. While the RMS values of ZTD had clear trends with latitude, decreasing from the equator to both poles. Further, the distribution of the annual average bias and RMS values of ZTD over the globe from 2016 to 2017 are shown in Figure 7. We can find that both MERRA-2 and ERA5 ZTD RMS values were larger at low and middle latitudes and smaller at high latitudes, mainly influenced by the complex climate and drastic changes in water vapor at low latitudes [43]. The RMSs of ERA5 ZTD were generally smaller than MERRA-2 at low latitudes, especially in South Asia and Southeast Asia. Meanwhile, the accuracy of ERA5 ZTD was slightly improved over MERRA-2 at middle and high latitudes, especially in North America, Europe, Australia, and the eastern part of Asia. This indicates that ERA5 ZTD agrees well with GNSS ZTD more than MERRA-2 ZTD in these regions.

Evaluations of the ZHD, ZWD, and T m Derived from MERRA-2 and ERA5 by Radiosonde Data
Through Section 3.1, we know that the accuracy of ZTD calculations from pressurelevel data using the integration method is higher than that derived from surface-level data using the Saastamoinen model method. Thus, the following comparison and analysis were also mainly based on ERA5 and MERRA-2 pressure-level ZHD, ZWD, and T m . ZHD, ZWD, and T m calculated from radiosonde data from 545 stations from 2016 to 2017 were used to evaluate the accuracies of corresponding tropospheric parameters derived from ERA5 and MERRA-2 pressure-level data. The annual mean bias and RMS of ZHD, ZWD, and T m values at each station were calculated. The results are shown in Table 4 and Figure 9.  The RMSs of ZHD derived from MERRA-2 data were almost the same as those of the ERA5 data, with an annual mean RMS of 4.4 mm for MERRA-2 and 4.5 mm for ERA5. Meanwhile, the RMSs of ZHD from MERRA-2 and ERA5 data compared to radiosonde ZHD were less than 10 mm for 94% of the 545 radiosonde stations. The accuracy of T m calculated from MERRA-2 data was almost similar to that from the ERA5 data, with an annual mean RMS of 1.17 K for MERRA-2 and 1.03 K for ERA5. In addition, the RMSs of T m derived from MERRA-2 and ERA5 data, with respect to radiosonde T m , were less than 2.2 K for 93% of the 545 stations. Furthermore, the performance of ZWD derived from ERA5 data slightly improved over that of the MERRA-2 data, with an RMS of 10.5 mm for ERA5 and 13.6 mm for MERRA-2. The reasons for this are described in Section 3.1. In a word, the tropospheric parameters derived from MERRA-2 and ERA5 data had high precision.
To analyze the seasonal variation features of ZHD, ZWD, and T m derived from ERA5 and MERRA-2 data, the monthly average bias and RMS values of ZHD, ZWD, and T m at all radiosonde stations in four seasons are also summarized in Table 5. There are 470 radiosonde stations in the Northern Hemisphere and 75 radiosonde stations in the Southern Hemisphere. Seasonal variations did not evidently impact the accuracy of ZHD calculated from reanalysis data, with smaller differences in bias and RMS for different seasons over the globe, while a seasonal trend in bias and RMS for ZWD was evident in the Northern Hemisphere. The maximum and minimum RMS appeared in the summer and winter months, with an RMS of 19.1 and 8.0 mm for MERRA-2 and 14.8 and 6.6 mm for ERA5, respectively. This is caused by the abundant and changeable water vapor in summer [44]. Meanwhile, seasonal variations in the ZWD RMS in the Southern Hemisphere were opposite to that in the Northern Hemisphere, with an RMS of 14.6 and 18.8 mm in the summer and winter for MERRA-2 and 9.7 and 12.3 mm in summer and winter for ERA5, respectively. Furthermore, the biases of ERA5 T m in winter were higher than those in summer in the Northern Hemisphere. The major reason for this is that the variation in T m during winter is stronger than that in summer [21,45]. At the same time, the ERA5 T m has a larger bias in summer and a smaller bias in winter in the Southern Hemisphere. However, the seasonal variation in RMS of T m in the Southern Hemisphere and Northern Hemisphere was not obvious.  Figure 10 shows the variations in annual mean bias and RMS of ZHD, ZWD, and T m with latitude. The variations of biases of ZHD, ZWD, and T m with latitude were not distinct. However, the trend in RMSs of ZWD with latitude was clear, decreasing from the equator to the poles. Figures 11-13 display the distribution of the annual mean bias and RMS values of ZHD, ZWD, and T m over the globe. Across the board, the RMSs of ZHD were mostly within 1 cm, and those of T m were mostly within 2 K over the globe. Moreover, the RMSs of ZHD and T m represent slightly higher values in the north of Asia. This may be due to the seasonal variations of ZHD and T m at high latitudes that are stronger than at low latitudes [21]. There was less difference between the accuracy of ZHD calculated from MERRA-2 and that of ERA5 over the globe. While the RMSs of ERA5 T m were generally smaller than MERRA-2 at low latitudes, especially in South Asia, Southeast Asia, and the north of South America. In the meantime, the performance of ERA5 ZWD improved over MERRA-2 in these regions. In addition, the RMSs of ZWD were relatively large at low latitudes. This was due to the influence of the tropical climate, such as tropical rain forest and tropical monsoon, in the low latitude area, where the water vapor content is rich and changes violently. At the same time, the changes in the tropospheric meteorological parameters were more active as a result of the high temperature, which leads to an increase in the uncertainty in the calculation of tropospheric parameters from reanalysis data, and the accuracy of ZWD is reduced [46,47].     Figures 14 and 15. The trends in the bias and RMS of ZHD, ZWD, and T m with altitude were not clear. We can find that the radiosonde stations below 500 m showed relatively large bias and RMS values, mainly influenced by active meteorological parameters in low elevation areas, resulting in large errors of interpolation meteorological parameters. All the same, the RMSs of ZHD, ZWD, and T m were mostly within 15, 27 mm, and 2.7 K. This further shows that taking the station height as the integral starting height can effectively reduce the error of tropospheric parameters interpolation in the elevation direction.

Conclusions
The ZTD, ZWD, ZHD, and T m derived from MERRA-2 and ERA5 data were validated using radiosonde data and GNSS ZTD. The spatial and temporal distribution characteristics of bias and RMS values were also analyzed.
Taking hourly GNSS ZTD as reference values, the annual mean RMS of MERRA-2 ZTD from pressure-level data and surface-level data were 13.1 and 27.9 mm, respectively, and that of ERA5 ZTD from pressure-level data and surface-level data were 10.9 and 31.6 mm, respectively. This indicates that the performance of ZTD derived from pressure-level data was much better than surface-level data. Meanwhile, the accuracy of ERA5 pressure-level ZTD was slighter higher than MERRA-2 pressure-level ZTD, which may be caused by differences in data assimilation and physical models. Additionally, the temporal and spatial variation characteristics of the bias and RMS of ZTD derived from pressure-level data were analyzed over the globe. The seasonal variations in ZTD from ERA5 and MERRA-2 were clear over the globe, with larger errors in the summer and smaller errors in the winter in the Northern Hemisphere, while there were smaller errors in the summer and larger errors in the winter in the Southern Hemisphere. Furthermore, it was found that the trend in RMS of ZTD with latitude was obvious, decreasing from the equator to both poles. However, the trend in the bias and RMS of ZTD with altitude was not evident.
Compared with ZHD, ZWD, and T m calculated from radiosonde data, the annual average RMS of ZHD, ZWD, and T m derived from MERRA-2 data were 4.4, 13.6 mm, and 1.17 K, respectively, and those of ERA5 were 4.5, 10.5 mm, and 1.03 K, respectively. The bias and RMS of ZWD had evident seasonal variation which was similar to that of ZTD. Meanwhile, the trend for RMS of ZWD with latitude was also similar to that of ZTD. However, the RMS values for ZHD and T m were relatively larger at high latitudes in the Northern Hemisphere. This may be due to the seasonal variations of ZHD and T m at high latitudes that were stronger than at low latitudes. In addition, the variations in ZHD, ZWD, and T m with altitude were not clear.