Determination of Tropospheric Parameters from ERA Surface Data for Space Geodetic Techniques

: This study investigates methods of deriving meteorological parameters needed in space geodetic applications, from the surface data of the numerical weather model (NWM). It is more efﬁcient than pressure level data in terms of storage and transmission. Based on more realistic assumptions for the structure of the troposphere, formulas for accurate vertical reduction of pressure ( P ) and precipitable water vapor (PWV) are deduced, and they are applied with the gridded lapse rate data provided by the GPT2w model. The new method achieves better accuracy especially when a large height difference between the grid point and station exists. Validation with global radiosonde observations shows that the RMS errors of P , temperature ( T ), and water vapor pressure ( e ) derived from 2.5 ◦ × 2.5 ◦ ERA surface data are 1.16 hPa, 1.95 K, and 1.76 hPa respectively; zenith tropospheric delays (ZTDs) calculated from derived P , T , and e values have a mean RMS error of 3.26 cm, comparable to that obtained from in situ measurements; adding PWV will increase ZTD estimation accuracy to 1.52 cm, comparable to that obtained from NWM pressure level data. Validations with Global Navigation Satellite System estimated ZTDs from global and regional station networks display similar results on the globe, as well as features for localized regions. Using higher spatial resolution NWM seems to have little effect on the accuracy of ZTDs calculated from P , T , and e , while it apparently improves the accuracy of ZTDs calculated from P , T , e , and PWV.


Introduction
The neutral atmosphere effect, mainly the signal transmission time delay, exists in a variety of radio space geodetic techniques, such as Global Navigation Satellite System (GNSS), Very Long Baseline Interferometry (VLBI), as well as optical techniques like Satellite Lase Ranging (SLR) [1]. To compensate for such an effect, usually information about atmospheric parameters at the observation station should be known beforehand or sometimes be estimated in geodetic data analysis. Improper handling of the neutral atmosphere effect would lead to obvious degradation of measurement accuracy, for example precise positioning and precipitable water sensing. Those atmospheric parameters include pressure, temperature, humidity, tropospheric delay, weighted mean temperature, and so on. Besides utilization in the determination of zenith hydrostatic delays (ZHD), pressure is also related to atmosphere pressure loading [1], and temperature can be applied in modeling the thermal deformation of VLBI telescopes or geodetic monuments [1], and it also affects the GNSS receiver-dependent hardware delays between different frequencies [2].
Obviously, the most accurate and direct way to obtain useful atmospheric information is meteorological sensors at geodetic sites; however, local measurements are usually not available or accessible. Currently, most International GNSS Service (IGS) reference stations are not equipped with collocated meteorological sensors, and some meteorological sensors at IGS stations are not properly calibrated, leading to erroneous measurements [3,4]. In usage limitation, we will not obtain the vertical reduction method by statistical analysis on local climate for certain regions, as was done in some studies [13]. Instead, we deduce algorithms for parameters from basic atmospheric functions, and then apply them with gridded coefficients derived from the empirical model GPT2w, mainly the temperature and water vapor pressure lapse rates.
Section 2 describes the characteristics of different data sets used in this study. Section 3 introduces the vertical reduction approach for atmospheric parameters. Section 4 presents a comprehensive assessment of the results obtained using our method by a variety of validation data. Finally, the conclusions are summarized in Section 5.

Data
In this study, the estimation values of atmospheric parameters at the target location are retrieved from ERA-Interim surface data together with the lapse rate data from the GPT2w model. Then, the derived parameters are assessed through two worldwide meteorological data sets, namely the Integrated Global Radiosonde Archive (IGRA) and IGS GNSS-ZTD data, and two regional GNSS-ZTD data sets over China and Europe, respectively. All data cover a time period of one year from January 2015 to December 2015.

ERA-Interim Surface Data
ECMWF has been producing a variety of forecasts, analysis, climate reanalysis, and specific data sets since the 1970s by assimilating a wide range of meteorological observations collected by satellites and earth observation systems via atmospheric models. ERA-Interim is a global atmospheric reanalysis product of ECMWF [26], which contains gridded atmospheric and surface parameters. This study utilizes the ERA-Interim surface fields, including both the analysis fields and forecast fields, which are available every 6 h (at 00, 06, 12, 18 UTC) and every 3 h (at 00, 03, 06, 09, 12, 15, 18 UTC), respectively. These forecasts are produced twice daily (from 00 and 12 UTC) at forecast steps from 3 to 12 h. Surface pressure (P), dewpoint temperature at 2 m (d2m), temperature at 2 m (T), and precipitable water vapor are chosen to determine the meteorology of the target station with the aid of the lapse rates provided from the GPT2w model. This study is primarily based on data field on a global grid of 2.5 • × 2.5 • ; however, the 0.75 • × 0.75 • grid that represents the original horizontal resolution of ERA-Interim is also considered for investigating the effect of spatial resolution. Here, we do not employ the advanced ERA5 hourly products [27] provided by ECMWF, because the higher resolution will greatly increase the burden of data storage and transmission, and then weaken the advantage of this method.

GPT2w
The GPT2w model is an empirical tropospheric delay model, capable of providing estimations for air pressure (P), temperature (T) and its lapse rate (β), water vapor pressure (e) and its height decrease factor (λ, also called the lapse rate of water vapor pressure), weighted mean temperature (T m ), and mapping function coefficients [7]. These parameters are modeled using a harmonic function representing annual and semiannual variations on a global grid of 1.0 • × 1.0 • . The GPT2w model exhibits good performance in the vicinity of the earth surface, with generally low system offsets for most regions over the global and a mean RMS error of about 3.6 cm for ZTD estimations [7]. In this study, we mainly use the lapse rate values of temperature and water vapor pressure from GPT2w to adjust P, T, e, and PWV from the ERA-Interim surface grid height to the height of the target height.

IGRA
IGRA, which we obtained from the National Climatic Data Center of NOAA, is the most complete and comprehensively used radiosonde data source. It is comprised of 1-4 quality-assured radiosonde observations per day at more than 1500 globally distributed stations with varying periods of record [28]. Radiosonde observations contain profiles of pressure, temperature, and water vapor pressure from the surface to the top of the sounding Remote Sens. 2021, 13, 3813 4 of 24 (around 30 km). Those parameters are utilized to calculate the atmospheric refractivity, and then the tropospheric delays (ZHD, ZWD, and ZTD) at the surface level are obtained through integration of refractivity along the whole profile. The detailed algorithm is shown in Section 3.1. According to a previous analysis, the ZHD and ZWD calculated from the radiosonde data have mean accuracies of 2.4 and 4.4 mm, respectively [14]. By using observations at global IGRA stations, the atmospheric parameters retrieved from ERA-Interim surface fields can be comprehensively evaluated. The reference values include P, T, e, ZHD, ZWD, and ZWD at the surface level of radiosonde. IGRA data are available from ftp://ftp.ncdc.noaa.gov/pub/data/igra/v1/derived-v2/data-por/, accessed on 13 November 2019.
Although various quality assurance algorithms have already been applied to remove gross errors in the IGRA data, we still need further checks in order to obtain reliable ZHD, ZWD, and ZTD values. In this study, a radiosonde sounding is eliminated when any of the following criteria is satisfied: (1) the highest level of sounding is lower than 15 km or the number of vertical levels is no more than 4; (2) the pressure difference between two successive vertical levels is great than 300 hPa; (3) humidity observation or water vapor pressure is not available at the surface (i.e., the lowest level of sounding); (4) the highest level where valid water vapor pressure observation is available shows a value greater than 0.1 hPa; (5) the height difference between two successive vertical levels where valid water vapor pressure observation is available exceeds 10 km. In addition, the radiosonde data only provides geopotential height, which should be transformed to geometric height before numerical integration is performed to calculate the tropospheric delays. Only IGRA stations that have valid records on more than 70% of the days in 2015 are considered in the validation, and this leads to 530 stations in all, with their geographic locations being shown in Figure 1a. In this study, only the surface (the lowest level) meteorological observations of radiosonde soundings and the calculated surface tropospheric delays are used to assess the retrieved parameters.

GNSS-ZTD
Since this study aims to provide meteorological parameters for geodetic users, ZTD time series derived from GNSS observation data are also employed to test the performance at those GNSS stations. The IGS ZTD data, which is now routinely produced over hundreds of stations with a temporal resolution of 5 min, have a formal error of 4 mm [29]. In the statistics, GNSS-ZTD data with a formal error exceeding 20 mm were eliminated and the comparisons were conducted every 30 min. To avoid the seasonal bias, only stations with valid record lengths longer than 70% of the entire period are considered, and this leads to 302 stations in total. These criteria are also applied in the other two GNSS-ZTD data sets. The geographic locations of those IGS stations are shown in Figure 1c.
The Crustal Movement Observation Network of China (CMONOC), which includes more than 260 continues GNSS reference stations as well as other geodetic technologies such as VLBI and SLR over the Chinese mainland, have provided significant basic data for crustal movement and geodynamic research, earthquake prediction, wide-area differential GNSS services, atmospheric science, etc. [30]. The GNSS data product service platform of China Earthquake Administration routinely collects the CMONOC GNSS observation data and calculates and distributes the hourly tropospheric products (http://www.cgps.ac.cn, accessed on 10 December 2019). Similar criteria are applied to exclude improper data and this leads to 237 stations in all. The geographic locations of those stations are shown in Figure 1e.
The EUREF Permanent Network (EPN) is now made up of more than 280 continuously operating GNSS reference stations distributed over Europe, the North Atlantic, and Greenland. Since 2001, the EPN analysis centers routinely archive and distribute the tropospheric zenith delay products in time interval of an hour [31]. These ZTD data are available from https://igs.bkg.bund.de/root_ftp/EUREF/products/, accessed on 10 December 2019. Applying the above-mentioned criteria leads to 256 stations in all and Remote Sens. 2021, 13, 3813 5 of 24 their geographic locations are shown in Figure 1g. It has to be pointed out that about 75 stations from EPN also belong to IGS.
clude improper data and this leads to 237 stations in all. The geographic locations of those stations are shown in Figure 1e.
The EUREF Permanent Network (EPN) is now made up of more than 280 continuously operating GNSS reference stations distributed over Europe, the North Atlantic, and Greenland. Since 2001, the EPN analysis centers routinely archive and distribute the tropospheric zenith delay products in time interval of an hour [31]. These ZTD data are available from https://igs.bkg.bund.de/root_ftp/EUREF/products/, accessed on 10 December 2019. Applying the above-mentioned criteria leads to 256 stations in all and their geographic locations are shown in Figure 1g. It has to be pointed out that about 75 stations from EPN also belong to IGS.  (b,d,f,h) show normalized histograms of height differences between ERA surface data grids and observation stations of IGRA, IGS, CMONOC, and EPN, respectively. In (b,d,f,h), the black lines refer to results for the 2.5° × 2.5° grid while the red lines refer to those for the 0.75° × 0.75° grid. The percentages of height differences within ±300 m are also displayed, black for the resolution of 2.5° × 2.5° and red for the resolution of 0.75° × 0.75°.
An overview of these data sets is shown in Table 1. Figure 1 also displays histograms of the height differences between ERA surface data grids and the observation stations of four validation data. The x axis refers to grid height minus station height. In most cases, the height differences lie within ±300 m. When the 2.5° × 2.5° grid is applied, the percentages of height difference within ±300 m are found to be 80.8, 75.5, and 82.0% for IGRA, IGS, and EPN stations, respectively; using a 0.75° × 0.75° grid will increase the percentages  (b,d,f,h) show normalized histograms of height differences between ERA surface data grids and observation stations of IGRA, IGS, CMONOC, and EPN, respectively. In (b,d,f,h), the black lines refer to results for the 2.5 • × 2.5 • grid while the red lines refer to those for the 0.75 • × 0.75 • grid. The percentages of height differences within ±300 m are also displayed, black for the resolution of 2.5 • × 2.5 • and red for the resolution of 0.75 An overview of these data sets is shown in Table 1. Figure 1 also displays histograms of the height differences between ERA surface data grids and the observation stations of four validation data. The x axis refers to grid height minus station height. In most cases, the height differences lie within ±300 m. When the 2.5 • × 2.5 • grid is applied, the percentages of height difference within ±300 m are found to be 80.8, 75.5, and 82.0% for IGRA, IGS, and EPN stations, respectively; using a 0.75 • × 0.75 • grid will increase the percentages by about 3 to 4%. The percentages are relatively lower for the CMONOC stations, which is due to the fact that those stations are mainly built to monitor crustal deformation and movement, so they tend to be located on regions with high elevation fluctuations. This will also provide good opportunity to test our methods. The propagation delay of a radio signal caused by the neutral atmosphere depends on the refractivity along the travel path. The delay in the zenith direction is expressed as: where n is the refractive index, N is the refractivity and N = 10 6 (n − 1). H is the height of the station in m and z is the signal path in the zenith direction in m. The expression of N is as follows: where P d and e are the partial pressures of dry air and water vapor in hPa. T is the temperature in K. The constants k i (i = 1, 2, 3) were determined experimentally. Here, we adopt the values given by [32]. The refractivity N can be split into a hydrostatic term N h and a nonhydrostatic term N w , N = N h + N w , where ρ is the air density and ρ w is the density of water vapor; R d and R w are the specific gas constants of dry air and water vapor, R d = 287.053 J/kg/K, R w = 461.525 J/kg/K.
Then, the zenith tropospheric hydrostatic delay and wet delay are obtained by integrating N h and N w respectively.
When calculating ZTD, ZHD, and ZWD from atmospheric profiles such as radiosonde observations and NWM pressure level data, trapezoidal numerical integration is often employed. However, since the refractivity usually follows an exponential decline with increasing height, trapezoidal integration would cause slightly positive biases, especially when there is a large height gap in real observations. It is more accurate to use the logarithmic trapezoidal integration method [14], as shown below.
where the subscript i refers to the ith level. Since N = dN/(dN/N) = dN/dInN, dN/dInN can be used to obtain the mean N value between two successive levels as in the Remote Sens. 2021, 13, 3813 7 of 24 above formulas. For exponential variation parameters, it is a more precise representation than a linear average.

Calculation of Atmospheric Parameters from ERA Surface Data
This section introduces the method of deriving the values of P, T, e, and PWV at the target station. The procedure consists of three steps. First, vertical reduction is performed to acquire the properties at the same height with observation station. Then, horizontal interpolation is applied using bilinear interpolation between four surrounding grid points. The final step is a linear trend interpolation in the time domain. The last two steps are apparently easy to implement, so vertical reduction is the primary concern of this section. In the following analysis, subscript H refers to properties at the height of the station, while subscript s refers to those at the height of surface data, e.g., the grid of ERA surface data.
In the troposphere, the temperature decreases with height at a nearly constant rate. Thus, the vertical reduction of temperature is normally performed based on the temperature lapse rate β, as shown in Equation (7). In previous studies [7,33,34], the variation of water vapor pressure with height is assumed to be related to the variation of total pressure using the factor λ, as shown in Equation (8). In this study, both equations are applied with the values of β and λ derived from gridded file of the GPT2w model.
T s , H s , P s , and e s refer to corresponding parameters at the grid of ERA surface data. The reduction of pressure is realized through applying assumptions on vertical atmosphere structure. For example, the GPT2w model uses a vertical reduction function obtained under the assumption of isothermal atmosphere, i.e., where g is the gravity acceleration in m/s 2 . Tv is the virtual temperature in K, and it is the temperature of dry air that has the same density at the same pressure as moist air of temperature T. All relationships for dry air are valid for moist air if we replace temperature T by the virtual temperature Tv. As is known, the atmosphere in the lower troposphere should be treated as moist air, especially for middle and lower latitudes, or it will lead to obvious bias. In comparison to the isothermal atmosphere, the polytropic atmosphere, i.e., a model atmosphere assuming hydrostatic equilibrium with a constant nonzero lapse rate, reflects the realistic troposphere better. Under this assumption, the pressure in the moist atmosphere is calculated as: where Γ is the lapse rate of the virtual temperature in K/m. Equation (10) is adopted for pressure reduction in this study. Before using Equation (10), the expression of Γ is needed. From the definition of Tv, where q is the specific humidity: q = εe P = 0.622 e P . The formula of Γ can be deduced as: Remote Sens. 2021, 13, 3813 8 of 24 Using Equation (8) and the hydrostatic equation ∂P ∂z = −ρg while taking the derivative of q with respect to z, we obtained the following expression: However, in the polar regions, the structure of the temperature profile sometimes significantly deviates from the above constant lapse rate assumption, for the actual temperature in the lower troposphere shows an increase with height. For example, the mean temperature lapse rate values of GPT2w are mostly negative in the higher latitudes. Thus, in those regions isothermal atmosphere assumption is still adopted for pressure reduction, i.e., Equation (9).
The vertical reduction of PWV is achieved according to its definition: where ρ 1 is the density of liquid water, ρ 1 = 1000 kg/m 3 . The integration term of the right side of Equation (14) can be calculated using logarithmic trapezoidal integration method, as shown below.
In most GNSS applications, the user location has dramatic height difference with the base station, so vertical reduction of meteorological parameters is always necessary to account for errors caused by height difference. Currently, many analyses resort to simple models with global uniform coefficients, for example the Berg function [35] for pressure reduction (see Equation (10) in [36]). The GPT2w model is capable of providing accurate lapse rate data with high spatial resolution in the vicinity of earth's surface; thus, vertical reduction of parameters relates to atmospheric effect on GNSS data analysis can be conveniently accomplished by these lapse rate values. Obviously, the methods discussed in this section can also be applied for other surface atmospheric data.

Calculation of Tropospheric Delay from P, T, e, PWV
After the P, T, e, and PWV at the target location have been retrieved, ZHD, ZWD, and ZTD are calculated from those parameters using three methods. The first method utilizes the Saastamoinen model [37] with P, T, and e as input parameters. The expressions for ZHD and ZWD are as follows: ZTD is the sum of ZHD and ZWD. The second method is based on the formulas presented by Askne and Nordius [34]: wherein the values of parameter β and λ are extracted from the grid file of the GPT2w model. g m is the gravity acceleration at the atmospheric column centroid in m/s 2 and can be computed as: where ϕ is the latitude. Equation (19) is highly accurate because it is a theoretical solution of ZHD (the integration of N h ) under the assumption of hydrostatic equilibrium, which is normally satisfied in a realistic environment. The third method differs from the second method only in ZWD calculation, i.e., T m is obtained from the GPT2w model. From the above, it is seen that methods 1 and 2 are suitable for the condition that surface pressure, temperature, and humidity are acquired while method 3 is applicable when the PWV value is also available. Table 2 presents the global mean accuracies of P, T, and e retrieved from ERA surface data, using surface values of IGRA observations as references. It lists validation results for three ERA surface data sets: analysis fields on a grid of 2.5 • × 2.5 • and on a grid of 0.75 • × 0.75 • , and forecast fields on a grid of 2.5 • × 2.5 • . Overall, the retrieved P, T, and e exhibit good accuracy for all three data sets. The meteorological parameters determined from the 2.5 • × 2.5 • grid surface analysis data show a mean RMS of 1.16 hPa for pressure, a mean RMS of 1.95 K for temperature, and a mean RMS of 1.76 hPa for water vapor pressure. However, increasing the resolution from 2.5 • to 0.75 • only leads to marginal accuracy improvements, about 0.10 hPa, 0.14 K, and 0.07 hPa for P, T, and e, respectively. All parameters derived from the forecast fields show relatively higher RMS errors than those from the analysis fields. To illustrate more detailed results, Figure 2 displays station-wise bias and RMS values of pressure, temperature, and water vapor pressure derived from ERA surface analysis data (blue triangle: 2.5 • × 2.5 • ; red cross: 0.75 • × 0.75 • ) over IGRA stations. Here, characteristics of the errors from the forecast data are not shown in the figure, because they are very similar to those from the analysis data.

Validation of P, T, e
As can be seen from Figure 2a, the biases for pressures derived from ERA surface 2.5 • × 2.5 • analysis data vary from −6 to 6 hPa on a global scale, showing no obvious latitude-dependent variation. Pressure biases are within ±2 hPa over more than 96% of the stations, and RMS values are less than 2 hPa over more than 90% of the stations, which indicates comprehensive good performance around the globe. It is seen from Figure 2d that in general pressure demonstrates higher RMS values in the higher latitudes than in the lower latitudes, which is consistent with the latitudinal variation of the magnitude of surface pressure. Remarkable offsets in pressure are observed over a few stations located in the polar regions and middle latitudes of the Northern Hemisphere (NH), showing mostly positive values above 3 hPa. As found in Figure 2a As can be seen from Figure 2a, the biases for pressures derived from ERA surface 2.5° × 2.5° analysis data vary from −6 to 6 hPa on a global scale, showing no obvious latitudedependent variation. Pressure biases are within ±2 hPa over more than 96% of the stations, and RMS values are less than 2 hPa over more than 90% of the stations, which indicates comprehensive good performance around the globe. It is seen from Figure 2d that in general pressure demonstrates higher RMS values in the higher latitudes than in the lower latitudes, which is consistent with the latitudinal variation of the magnitude of surface pressure. Remarkable offsets in pressure are observed over a few stations located in the polar regions and middle latitudes of the Northern Hemisphere (NH), showing mostly positive values above 3 hPa. As found in Figure 2a,d, over a few polar stations, the absolute biases can be apparently reduced by about 1.5 to 1.9 hPa when using the 0.75° × 0.75° grid instead of the coarser 2.5° × 2.5° grid and the RMS values can be reduced by about 1.0-2.1 hPa. However, improving horizontal resolution does not seem to have an obvious impact on the performance of the lower and middle latitude regions.
As seen in Figure 2b, the biases of temperatures derived from ERA surface 2.5° × 2.5° and 0.75° × 0.75° analysis fields both vary from about −4 to 4 K, with values falling into the range of ±2 K over about 96% of the stations. More pronounced biases occur in the Antarctic, which is possibly attributed to its unique vertical structure of atmosphere (temperature) that is usually dominated by polar vortex. Those features cannot be realistically and accurately represented by the GPT2w model, for only seasonal variations are considered in GPT2w. The temperature RMS values obviously show an increase as latitude increases, with lower values (<2 K) in the tropics and more scattered distribution (from 1 to 5 K) outside the tropics. RMS values are within 3 K over more than 91% of the stations for temperatures derived from the 2.5° × 2.5° grid and over more than 93% of the stations for those derived from the 0.75° × 0.75° grid.
It can be seen from Figure 2c,f that the accuracy of derived water vapor pressure is highly latitude-dependent, showing more remarkable bias and RMS values in the lower latitudes and smaller values (less than 1 hPa) in the higher latitudes. However, unexpectedly, the most significant errors appear in the subtropical regions of both hemispheres. For example, a few coastal stations of the Arabian Peninsula, eastern Mediterranean, South Mexico, and several high elevation stations over the Tibetan Plateau, according to our further investigations. We also find that those coastal stations are basically consistent As seen in Figure 2b, the biases of temperatures derived from ERA surface 2.5 • × 2.5 • and 0.75 • × 0.75 • analysis fields both vary from about −4 to 4 K, with values falling into the range of ±2 K over about 96% of the stations. More pronounced biases occur in the Antarctic, which is possibly attributed to its unique vertical structure of atmosphere (temperature) that is usually dominated by polar vortex. Those features cannot be realistically and accurately represented by the GPT2w model, for only seasonal variations are considered in GPT2w. The temperature RMS values obviously show an increase as latitude increases, with lower values (<2 K) in the tropics and more scattered distribution (from 1 to 5 K) outside the tropics. RMS values are within 3 K over more than 91% of the stations for temperatures derived from the 2.5 • × 2.5 • grid and over more than 93% of the stations for those derived from the 0.75 • × 0.75 • grid.
It can be seen from Figure 2c,f that the accuracy of derived water vapor pressure is highly latitude-dependent, showing more remarkable bias and RMS values in the lower latitudes and smaller values (less than 1 hPa) in the higher latitudes. However, unexpectedly, the most significant errors appear in the subtropical regions of both hemispheres. For example, a few coastal stations of the Arabian Peninsula, eastern Mediterranean, South Mexico, and several high elevation stations over the Tibetan Plateau, according to our further investigations. We also find that those coastal stations are basically consistent with the regions where large standard deviations of λ from GPT2w are observed, as shown in Figure 1 of [7], which possibly explains the large errors over those regions. For both 2.5 • × 2.5 • and 0.75 • × 0.75 • grid analysis data, the biases of water vapor pressure are between −2 hPa and 2 hPa over about 90% of the stations, and the RMS values are below 3 hPa over about 90% of the stations. Moreover, using the 0.75 • × 0.75 • grid instead of the 2.5 • × 2.5 • grid only leads to marginal improvement, which probably implies that the resolution of 0.75 • × 0.75 • is still not sufficient for accommodating the spatial distribution of water vapor pressure, especially for regions with significant topographic variations or along the seaside.
Unlike pressure and temperature, the derived water vapor pressures exhibit positive offsets in terms of both overall performance (Table 2) and results for most stations (Figure 2c). This is likely attributed to the determination of parameter λ. In GPT2w and related research [7,34], it is a common approach to invert the zenith wet delay expression (Equation (19)) with the surface parameters (ZWD, T m , e) obtained beforehand to the λ.
Hence, the value of λ is actually consistent with the surface ZWD and describes an average effect of e decreasing with height through the whole troposphere. Then it is sometimes not capable of capturing accurate vertical variation within small height range. In addition, we have made further analysis and noticed that for cases with large height differences between grid and station (>300 m), the water vapor pressure biases are mostly negative when the stations are situated higher than the surrounding grid points, and vice versa.
In addition, we also find that the errors partly stem from the discrepancy between ERA analysis data and radiosonde observations. Stations showing large errors are chosen for further analysis and we find that for some stations, similar remarkable biases against radiosonde observations are still seen when using the ERA pressure level data instead of the surface data to derive the meteorological parameters.

Effects of Vertical Reduction of Pressure
For pressure reduction from grid height to the target station height, the power function proposed by Berg [35] is often recommended [38,39]. Actually, the Berg function is an approximation of the barometric formula of polytropic atmosphere (i.e., Equation (10)), which employs the standard (constant) atmospheric parameter values at mean sea level (i.e., 1013.25 hPa for pressure, 288.15 K for temperature, and 6.5 K/km for temperature lapse rate) and ignores regional and temporal climate variations. Some research has applied a first-order approximation, with parameters from the GPT model or standard atmosphere [40]. Unfortunately, those methods cannot yield realistic results when there is a significant height difference. Figure 3 depicts the station-specific errors of pressures obtained by using the Berg function for height correction. In comparison to Figure 2a,d, pressure reduction with the Berg function presents much poorer performance, especially for the NH midlatitudes and polar regions, showing RMS values exceeding 5 hPa for lots of stations; while in Figure 2d the corresponding results are mostly below this level. The error induced by the Berg function is mainly related to two factors: the height difference between the grid and target station and the discrepancy between local climate and the standard atmosphere. The resultant pressure error would be more significant when both the large height difference and bad representation of standard atmosphere are met. According to Equation (10), if the target station is situated lower than the height of surrounding grid points and meanwhile the standard temperature is warmer than the local value, it will cause negative offset for the derived pressure. This is reflected by those polar stations with remarkable biases less than −5 hPa and RMS values reaching 5-15 hPa, as shown in Figure 3. Comparison results about the performance of the Berg function and our method for pressure reduction are shown in Table 3. Statistics are conducted respectively for two data sets obtained in terms of absolute height differences between ERA surface data grid and IGRA stations: less than 300 m and exceeding 300 m. It is seen that when the height dif- Comparison results about the performance of the Berg function and our method for pressure reduction are shown in Table 3. Statistics are conducted respectively for two data sets obtained in terms of absolute height differences between ERA surface data grid and IGRA stations: less than 300 m and exceeding 300 m. It is seen that when the height differences are within 300 m, the differences in pressure RMS errors between two methods are relatively small; our method is obviously superior when the absolute height differences are above 300 m, showing a mean improvement of 2.08 (1.31) hPa and better performance over about 90% (84%) of the stations when using ERA 2.5 • × 2.5 • (0.75 • × 0.75 • ) grid data. From Figure 3 and Table 3, it is also noticed that the Berg method is more sensitive to grid resolution for increasing the resolution from 2.5 • to 0.75 • can apparently improve the pressure reduction results, which is probably caused by smaller height differences for denser grid. However, the method used in this study is generally more stable, as seen from Figure 2. Table 3. Performance comparison of the Berg function and our method for pressure reduction. The bias and RMS differences (in hPa) refer to the results obtained from the Berg function minus those from our method. The percentage of stations at which our method outperforms the Berg function is also presented.

ERA Surface Data
Height It was pointed out that due to the hydrostatic/wet mapping separation effects, a 10 hPa error in total pressure at the station will result in a 2.4 mm error in height solution when the cutoff elevation angle is set to 5 • [36]. Therefore, in order to achieve mm level accuracy, we recommend applying the more complicated vertical reduction method for pressure proposed in this study instead of the Berg function when there is a remarkable height difference.

Validation of Tropospheric Delays over the Globe
Tropospheric delays, including ZHDs, ZWDs, and ZTDs, are calculated via functions shown in Section 3.3, with the atmospheric parameters extracted from ERA surface data as inputs. In this procedure, three methods were utilized: the Saastamoinen model with P, T, and e, the Askne and Nordius model with P, T, and e, and tropospheric delays derived with P, T, e, and PWV. In the following analysis, the Askne and Nordius model is labeled as A&N, and the Saastamoinen model is labeled as Saas. Table 4 shows the global mean accuracies of ZHDs, ZWDs, and ZTDs derived from the ERA surface analysis data, through comparisons with reference values calculated from IGRA data. For a better illustration, corresponding results obtained from actual meteorological observations (i.e., surface values of radiosonde soundings) are also shown in Table 4, as well as results of two empirical tropospheric models, GPT2w and IGGtropSH, which was developed by us [41]. As can be seen from Table 4, the global mean RMS errors are about 3.26 and 3.35 cm for ZTDs computed with the A&N and Saas models using P, T, and e derived from 2.5 • × 2.5 • ERA surface analysis data, respectively, which are very close to the corresponding values (3.20 and 3.41 cm) computed using measured surface parameters. It is indicated that atmospheric parameters derived from ERA surface data exhibit basically the same level of accuracy with actual observations in terms of ZTD estimation. It is also noticed that the A&N model slightly outperforms the Saas model when either derived or measured atmospheric parameters are as inputs. Unfortunately, no improvement is observed when higher spatial resolution ERA surface data are used, although the derived P, T, and e all exhibit marginally higher accuracy, as seen in Table 2. This is not a surprise result considering the limited relation between ZWD and water vapor pressure at the surface.

Overall Performance
When the PWV parameter is added for ZTD estimation, the performance is significantly enhanced, showing a mean RMS value of 1.52 cm for the 2.5 • × 2.5 • grid analysis, which is due to the high correlation between PWV and ZWD. Besides, it is noted that using higher resolution ERA PWV data has a positive effect on ZTD estimation, since the mean RMS error is reduced by more than 2 mm.
For ZHDs, the global mean RMS error is about 0.27 cm for both the A&N and Saas models using pressures derived from 2.5 • × 2.5 • ERA surface data, and this value decreases to 0.24 cm when the 0.75 • grid data is utilized, as seen in Table 4. Obviously, those results are fully consistent with the errors of derived pressure shown in Table 2. Nevertheless, given measured station pressures, both models estimate the ZHD at submillimeter level, which means hydrostatic equilibrium is always satisfied in the troposphere. The accuracy of ZHD obtained from ERA surface data is greatly superior to that obtained from the empirical model GPT2w, which has no auxiliary input, and shows a mean RMS value of 1.47 cm. In the GNSS-PWV retrieval procedure, the uncertainty in PWV due to the uncertainty in ZHD (σ ZHD ) is calculated as Π×σ ZHD (Π is the proportionality of PWV to ZWD and has a typical value of 0.15 [42]). Thus, the accuracy of ZHD obtained from ERA surface data can satisfy mm level PWV retrieval and it is convenient to utilize this approach when there is no collocated meteorological sensor.
For ZWDs obtained from surface P, T, and e data, the global mean RMS error is above 3 cm for either model. According to the results obtained using the in situ P, T, and e measurements as input, the Saas model tends to underestimate the ZWD, showing a global mean bias of about −0.49 cm; while the A&N model tends to overestimate the ZWD, showing a mean bias of about 0.39 cm. However, such overestimation becomes more pronounced as the A&N model is applied with P, T, and e derived from ERA surface data, since an overall bias of up to about 1.0 cm is seen on a global scale. This is generated by combination of two error sources: positive bias of derived e values (see Table 2) and positive ZWD bias obtained with the A&N model. Similarly, the relatively lower mean bias for ZWDs obtained from ERA surface P, T, and e data with the Saas model can be explained by two factors having opposite effects: the negative bias of the Saas model itself and the positive bias of the input parameter e. In terms of ZWDs derived with ERA surface P, T, e, and PWV data, the globally averaged bias is below 2 mm, reflecting no obvious offset between our approach and radiosonde observations. Table 5 shows similar results with Table 4 but obtained with ERA surface forecast data. It is noted that ZTDs calculated with the ERA surface forecast P, T, and e data demonstrate slightly smaller absolute bias and RMS than the corresponding values from analysis data, which is due to the smaller bias of e derived from the surface forecast data, as shown in Table 2. In terms of ZTDs obtained from ERA surface P, T, e, and PWV data, the RMS error is about 1.80 cm on average when the forecast data is utilized, about 3 mm greater than those from analysis data.

Spatial Behaviors
In addition to overall performance, more detailed statistics on the spatial behavior of tropospheric delay estimation accuracy is presented in the following section. Figure 4 shows the station-specific bias and RMS values of ZTDs derived from 2.5 • × 2.5 • ERA surface analysis data over the IGRA stations: (a) and (c) refer to the results obtained from surface P, T, e, and PWV data; (b) and (d) refer to those obtained from surface P, T, and e data with the A&N model.
As seen from Figure 4b, over the globe, the biases of ZTDs derived from ERA surface P, T, and e data range from about −1. It is noted from Figure 4a,c that the approach using ERA surface P, T, e, and PWV data provide highly accurate ZTD estimations and approximately homogeneous performance over the entire globe. The ZTD biases range roughly between −2.1 and 3.4 cm, and according to our analysis, at about 91% of the stations the biases are within ±1 cm and at about 98% of the stations the biases are within ±2 cm. The most pronounced over-prediction of ZTD appears over three stations located in northwest China's Xinjiang Province and Aswan (southern Egypt). As seen in Figure 4c, the RMS values are generally below 2 cm for most middle and higher latitude regions while slightly larger values appear over Xinjiang and the southeastern part of China, as well as a few stations located on the tropical islands and coasts. The RMS values are less than 2 cm over more than 80% of the stations, whereas this percentage can increase to 90% as the higher-resolution data (0.75 • ) is utilized, according to our investigation. We also find that some of the large biases result from the discrepancy between the two data sources for biases still exist in statistical comparison between ZTDs calculated from ERA pressure level data and radiosonde observations. about 98% of the stations the biases are within ±2 cm. The most pronounced over-prediction of ZTD appears over three stations located in northwest China's Xinjiang Province and Aswan (southern Egypt). As seen in Figure 4c, the RMS values are generally below 2 cm for most middle and higher latitude regions while slightly larger values appear over Xinjiang and the southeastern part of China, as well as a few stations located on the tropical islands and coasts. The RMS values are less than 2 cm over more than 80% of the stations, whereas this percentage can increase to 90% as the higher-resolution data (0.75°) is utilized, according to our investigation. We also find that some of the large biases result from the discrepancy between the two data sources for biases still exist in statistical comparison between ZTDs calculated from ERA pressure level data and radiosonde observations.  The effect of spatial resolution of ERA surface data on ZTD estimation is explored through displaying the station-specific differences of RMS values between two different resolutions in Figure 5. The symbol refers to the results obtained from 2.5 • × 2.5 • grid data minus those from 0.75 • × 0.75 • grid data. As found in Figure 5a, when using ERA surface P, T, and e data, the differences caused by different resolutions are within ±1 cm over most stations; higher spatial resolution seems not to yield superior performance as expected. Figure 5b indicates that, when ERA surface PWV data are also utilized for ZTD retrieval, slightly better accuracy is obviously acquired with higher spatial resolution data for almost all the latitudes. Such enhancement is more significant at the lower and middle latitude stations where water vapor is abundant. In addition, the three stations in Xinjiang Province that show relatively larger errors in Figure 4c, now exhibit an improvement of 1-2 cm in RMS errors with 0.75 • grid data, as indicated by the red marks in Figure 5b. stations; higher spatial resolution seems not to yield superior performance as expected. Figure 5b indicates that, when ERA surface PWV data are also utilized for ZTD retrieval, slightly better accuracy is obviously acquired with higher spatial resolution data for almost all the latitudes. Such enhancement is more significant at the lower and middle latitude stations where water vapor is abundant. In addition, the three stations in Xinjiang Province that show relatively larger errors in Figure 4c, now exhibit an improvement of 1-2 cm in RMS errors with 0.75° grid data, as indicated by the red marks in Figure 5b. In Figure 6, the station-specific bias and RMS values of ZHDs as a function of latitude are depicted, in which the red line refers to the results derived from ERA surface data and the black line refers to those derived from local pressure measurements. As seen from the red line, the ZHD RMS errors vary from 0.1 to 1.2 cm, with more than 92% of the values being less than 0.5 cm. In addition, Figure 6 also shows that ZHD calculated from local pressure measurements is accurate at sub mm level over each station, although slightly better performance is observed in the lower latitudes. In Figure 6, the station-specific bias and RMS values of ZHDs as a function of latitude are depicted, in which the red line refers to the results derived from ERA surface data and the black line refers to those derived from local pressure measurements. As seen from the red line, the ZHD RMS errors vary from 0.1 to 1.2 cm, with more than 92% of the values being less than 0.5 cm. In addition, Figure 6 also shows that ZHD calculated from local pressure measurements is accurate at sub mm level over each station, although slightly better performance is observed in the lower latitudes.  Figure 7 displays the latitude-dependent distribution of bias and RMS values of ZWDs calculated from different approaches. For the sake of clarity, zonal mean bias and RMS values for each 5-degree width latitude band are depicted instead of station-specific results. The two types of blue symbols present results of ZWDs calculated from the A&N (triangle) and Saas (square) models with realistic P, T, and e measurements. Since they only reflect the effect of empirical ZWD models, we explore their features firstly. From Figure 6. Distribution of station-specific bias (a,c) and RMS (b,d) values of ZHDs derived from ERA surface analysis data as well as local meteorological measurements. The reference ZHDs are calculated from radiosonde profiles over 530 IGRA stations. The red line refers to ZHDs obtained from the A&N model with P from ERA surface analysis data as input, and the black line refers to those obtained with in situ P measurements. Figure 7 displays the latitude-dependent distribution of bias and RMS values of ZWDs calculated from different approaches. For the sake of clarity, zonal mean bias and RMS values for each 5-degree width latitude band are depicted instead of station-specific results. The two types of blue symbols present results of ZWDs calculated from the A&N (triangle) and Saas (square) models with realistic P, T, and e measurements. Since they only reflect the effect of empirical ZWD models, we explore their features firstly. From Figure 7a, it is found that the in the lower latitudes there exists obvious biases for ZWD estimations derived with both A&N and Saas models, although with opposite signs. The zonal mean bias of the A&N model reaches about 2 cm between 15 • latitude north and south of equator, and then it decreases to approximate zero for other latitudes. The Saas model shows remarkable negative zonal mean biases up to about −3.5 cm at the equator, and as the latitude increases the magnitude of bias decreases, being mostly negative in the NH and positive in the southern hemisphere (SH); biases close to zero only occur in the Antarctic and a few latitude zones. It has to be mentioned that, for station-specific ZWD results, more pronounced biases are observed for both models than those shown in Figure 7, according to our investigation. ZWDs are calculated from radiosonde profiles over 530 IGRA stations. The red triangle refers to results obtained using P, T, and e from ERA surface analysis data with the A&N model, and the red square refers to those with the Saas model. The black asterisk refers to results obtained using P, T, e, and PWV from ERA surface analysis data. The blue triangle refers to results obtained using P, T, and e at the surface level of radiosonde data with the A&N model, and the blue square refers to those with the Saas model.
The above analyses elaborately reveal that the empirical Saas model causes unrealistic ZWD estimations for most regions of the world, which is due to its globally uniform and season-invariant coefficients; by using gridded model coefficients from GPT2w, the A&N model has shown obviously improved performance, but systematic offsets still exist over the tropics and further refinement on the parameter λ is still required. In addition, our statistics also find that, because of their significant biases, despite using in situ meteorological measurements as inputs the Saas and A&N models perform worse than blind models such as GPT2w and IGGtropSH for certain tropical and subtropical stations. We have made a preliminary attempt on model refinement, which investigated the inaccurate coefficients of the Saas model over the whole China region and established a station-spe- ZWDs are calculated from radiosonde profiles over 530 IGRA stations. The red triangle refers to results obtained using P, T, and e from ERA surface analysis data with the A&N model, and the red square refers to those with the Saas model. The black asterisk refers to results obtained using P, T, e, and PWV from ERA surface analysis data. The blue triangle refers to results obtained using P, T, and e at the surface level of radiosonde data with the A&N model, and the blue square refers to those with the Saas model.
The above analyses elaborately reveal that the empirical Saas model causes unrealistic ZWD estimations for most regions of the world, which is due to its globally uniform and season-invariant coefficients; by using gridded model coefficients from GPT2w, the A&N model has shown obviously improved performance, but systematic offsets still exist over the tropics and further refinement on the parameter λ is still required. In addition, our statistics also find that, because of their significant biases, despite using in situ meteorological measurements as inputs the Saas and A&N models perform worse than blind models such as GPT2w and IGGtropSH for certain tropical and subtropical stations. We have made a preliminary attempt on model refinement, which investigated the inaccurate coefficients of the Saas model over the whole China region and established a station-specific tropospheric model [43]. The improved model shows ZTD prediction residuals closer to a normal distribution with very low systematic bias. The TropSite model, which was proposed for the Galileo System and consists of station-specific coefficients respectively for hundreds of globally distributed stations, demonstrated a mean accuracy of about 2.5 cm when using realistic surface meteorological properties [44].
The two types of red symbols in Figure 7 depict the accuracy of ZWDs obtained from the A&N and Saas models with meteorological parameters derived from ERA surface data as model inputs. Therefore, they represent a combination effect of model error and error of inputs. As seen in Figure 7a, the red symbols exhibit similar latitudinal variation with the corresponding blue symbols expect for higher biases of about 0.5 to 1.5 cm, which is probably associated with the positive biases of the derived e values; the difference becomes much smaller for higher latitudes. Although these approaches show significantly different features in bias, the patterns of ZWD RMS values are similar, showing larger values in the tropics and then decreasing continuously with the increasing latitude, as found in Figure 7b,d. In addition, when driven with local meteorological measurements, the A&N model outperforms the Saas model generally over all latitude bands between 30 • S-30 • N; it is also superior for most latitude bands when driven with parameters obtained from ERA surface data.
It can be seen from Figure 7 that the approach (black symbol) that based on PWV data perform significantly better in terms of both bias and RMS. Constantly low bias around zero is seen for all the altitudes, and zonal mean RMS values vary from about 0.5 cm in the polar regions to around 2.0 cm in the tropics. Figure 7 also reveals that slightly better accuracy is achieved with higher resolution ERA surface data in almost all the latitude bands.
As seen in Table 4, the accuracy of ZWDs is highly correlated with that of ZTDs with almost identical bias and RMS values, suggesting that most of the ZTD errors stem from the mismodeling of the wet delays. Hence, Figure 7 can also indicate the features of ZTD errors.

Comparison with IGS GNSS-ZTDs
Statistical results of ZTD validation against another global data set, i.e., GNSS-ZTDs at 302 IGS stations, are summarized in Table 6. The results are comparable to those obtained by comparisons with IGRA data (Tables 4 and 5), except for slightly lower bias and RMS values for ZTDs calculated from ERA surface P, T, and e data shown in Table 6. The differences of ZTD estimation accuracies obtained from two validation data sets should be mainly due to the different geographical distributions of radiosonde and GNSS stations. A previous study about comparison of GNSS-ZTDs at global IGS stations and ray-traced delays from pressure level data of ECMWF had found that the agreement is about 0.5 cm in bias and 1.4 cm in standard deviation [45]. Recently, a similar study has found the mean bias to be 0.47 cm and standard deviation to be 1.40 cm [46]. In this study, the overall agreement between ZTDs obtained from ERA surface P, T, e, and PWV data and those from GNSS observations are about 1.55 cm for the 2.5 • × 2.5 • grid and 1.41 cm for the 0.75 • × 0.75 • grid, which is comparable to the results of [45,46]. Ray tracing through atmospheric profiles is normally considered to be the most accurate method for ZTD calculation; however, our study indicates that similar accuracy can also be achieved using surface data of NWM. In comparison to ray tracing, our approach is apparently easy to implement, because of its much simpler calculation procedure and requirement of less meteorological data. Therefore, for geodetic applications in which atmospheric effect must be processed precisely, surface data of NMW together with suitable vertical reduction method can be used instead of ray tracing.  Figure 8 shows latitudinal variations of bias and RMS values of ZTDs derived from ERA surface analysis data over global IGS stations. As mentioned, the feature of ZTD errors can be totally represented by that of ZWD errors; therefore, we can directly compare Figure 8 with Figure 7 here. Although in Figure 7 zonal means are displayed instead of station-specific values, we can still discern that biases obtained respectively from validations against two data sets demonstrate roughly consistent behaviors: noticeable positive biases for the A&N model and negative biases for the Saas model in the equatorial regions; smaller absolute biases in the higher latitude regions; more remarkable biases for the Saas model relative to the A&N model at most stations. Small difference is observed for RMS features: validation with IGRA data shows the largest RMS errors in the tropical and subtropical regions while validation with IGS data demonstrate smaller RMS values near the equator than subtropical regions. This is likely attributed to the much sparser coverage of IGS stations in the tropics.  Figure 8 shows latitudinal variations of bias and RMS values of ZTDs derived from ERA surface analysis data over global IGS stations. As mentioned, the feature of ZTD errors can be totally represented by that of ZWD errors; therefore, we can directly compare Figure 8 with Figure 7 here. Although in Figure 7 zonal means are displayed instead of station-specific values, we can still discern that biases obtained respectively from validations against two data sets demonstrate roughly consistent behaviors: noticeable positive biases for the A&N model and negative biases for the Saas model in the equatorial regions; smaller absolute biases in the higher latitude regions; more remarkable biases for the Saas model relative to the A&N model at most stations. Small difference is observed for RMS features: validation with IGRA data shows the largest RMS errors in the tropical and subtropical regions while validation with IGS data demonstrate smaller RMS values near the equator than subtropical regions. This is likely attributed to the much sparser coverage of IGS stations in the tropics.

Validation of Tropospheric Delays with Regional Data
This section describes the validations of tropospheric delays against GNSS-ZTDs over two regional networks. Based on denser coverage over regional areas, they will lead to more detailed assessment results compared to global networks. Table 7 shows the sta-

Validation of Tropospheric Delays with Regional Data
This section describes the validations of tropospheric delays against GNSS-ZTDs over two regional networks. Based on denser coverage over regional areas, they will lead to more detailed assessment results compared to global networks. Table 7 shows the statistics of mean accuracies of ZTDs calculated from ERA surface data over the CMONOC stations and EPN stations, respectively.
As mentioned above, the accuracy of ZTD estimation is impacted by two factors: the error of surface parameters and the error of the chosen empirical model. Before analysis, we introduce a simple way to estimate their respective influences roughly. According to Equation (17), we can deduce that an error of 1 hPa in water vapor pressure will result in a ZTD error of about 1 cm (here, T = 288 K is used, but the choice of T value has little effect on the estimation result). We can use this rule to roughly separate the ZTD error induced by e error from that induced by inaccurate tropospheric model. Firstly, performance over China is explored. As seen in Table 7, for the approach using the A&N model and derived P, T, and e, in comparison to its global statistics results ( Table 4) the overestimation of ZTD is more pronounced in China, showing a mean bias of 1.35 cm. We find that the mean bias of derived e over China (deduced from IGRA data) reaches up to 1.0 hPa, about 0.3 hPa larger than the global mean value. Such remarkable positive biases in e would account for most of the ZTD overestimations in China, although the CMONOC stations are not completely collocated with the radiosonde stations. Besides, the ZTDs obtained from derived P, T, e, and PWV also show slightly larger mean bias over China than corresponding global results. In terms of RMS errors, the overall performance is almost identical for the A&N model and the Saas model.
The station-specific validation results over China are shown in Figure 9, in which the biases and RMS errors of ZTDs obtained from the 2.5 • × 2.5 • ERA surface analysis data using three approaches are displayed. As seen in Figure 9a, over a few subtropical stations, ZTDs derived from P, T, and e with the Saas model demonstrate great negative biases, reaching about −4.5 to −3 cm. In fact, considering the fact that the mean bias of derived e values varies from about 1.6 hPa for latitudes below 30 • N to about 0.8 hPa for higher latitudes over China (obtained from radiosonde data), the actual negative bias caused by the Saas model should be more severe than that demonstrated in the Figure 9a. The A&N model performs much better than the Saas model in terms of bias over the subtropics.
It is seen in Figure 9b that RMS values of ZTDs derived from P, T, e, and PWV are less than 2.5 cm over most stations, except for a few stations within two latitude bands, 25 • N-31 • N and 38 • N-43 • N. Investigations find that those stations are located in the eastern Hengduan Mountains (in southwest China) and eastern Xijiang Province, respectively. From ZTD validation with radiosonde observations (Figure 4), relatively poor performance over Xinjiang Province is also found, but detailed features over the Hengduan Mountains are not available due to sparse coverage of radiosonde stations. Therefore, the CMONOC data facilitate a more comprehensive assessment of ZTD estimations in western China. The mean RMS errors of ZTDs derived from P, T, e, and PWV are about 1.54 cm and 1.43 cm for the 2.5 • × 2.5 • and 0.75 • × 0.75 • grid data, respectively. However, they seem relatively lower than the results of a few similar studies that also carried out for China: the mean RMS difference between the ZTDs provided by GNSS observations and ECMWF pressure level data is about 2.4 cm [13]; the mean RMS difference between ZWDs provided by radiosonde observations and ECMWF pressure level data is about 2.1 cm [14]. This is attributed to the relatively sparser sampling of CMONOC in the eastern and southern parts of China (see Figure 1e), where the climate is moister and thus causes higher discrepancy between different data sources.  Figure 1e), where the climate is moister and thus causes higher discrepancy between different data sources. Similar validation results over GNSS stations in Europe are presented in Table 7 and Figure 10. It is found that for each approach, even for empirical models such as GPT2w and IGGtrop, the overall ZTD estimation performance in Europe is better than that over the globe as well as in China. This is attributed to the mild weather in Europe and relatively higher latitudes of the EPN stations. It is noted that the mean latitude of stations is 48.8°N for EPN, 34.6°N for CMONOC, and 31.9°N for IGRA.
From Figure 10a, it is observed that between 25°N and 40°N, the Saas model overestimates the ZTD, with maximum biases reaching about 3 to 5 cm, while the A&N model shows superior performance due to its smaller absolute biases ranging between −1.5 and 2.5 cm. For latitudes higher than 50°N, ZTDs derived from the two models both demonstrate good agreement with the observations. The mean difference of water vapor pressures derived from ERA surface data and radiosonde observations is around 0.7 hPa in Europe, corresponding to a bias of 0.7 cm in ZTD estimations. This indicates that ZTD biases shown in Table 7 are mostly caused by the inaccurate estimation of water vapor pressure. Similar validation results over GNSS stations in Europe are presented in Table 7 and Figure 10. It is found that for each approach, even for empirical models such as GPT2w and IGGtrop, the overall ZTD estimation performance in Europe is better than that over the globe as well as in China. This is attributed to the mild weather in Europe and relatively higher latitudes of the EPN stations. It is noted that the mean latitude of stations is 48.8 • N for EPN, 34.6 • N for CMONOC, and 31.9 • N for IGRA.
From Figure 10a, it is observed that between 25 • N and 40 • N, the Saas model overestimates the ZTD, with maximum biases reaching about 3 to 5 cm, while the A&N model shows superior performance due to its smaller absolute biases ranging between −1.5 and 2.5 cm. For latitudes higher than 50 • N, ZTDs derived from the two models both demonstrate good agreement with the observations. The mean difference of water vapor pressures derived from ERA surface data and radiosonde observations is around 0.7 hPa in Europe, corresponding to a bias of 0.7 cm in ZTD estimations. This indicates that ZTD biases shown in Table 7 are mostly caused by the inaccurate estimation of water vapor pressure.
shows superior performance due to its smaller absolute biases ranging between −1.5 and 2.5 cm. For latitudes higher than 50°N, ZTDs derived from the two models both demonstrate good agreement with the observations. The mean difference of water vapor pressures derived from ERA surface data and radiosonde observations is around 0.7 hPa in Europe, corresponding to a bias of 0.7 cm in ZTD estimations. This indicates that ZTD biases shown in Table 7 are mostly caused by the inaccurate estimation of water vapor pressure.

Conclusions
This study has investigated the methods of deriving tropospheric parameters, including P, T, e, PWV, ZTD, ZHD, and ZWD, from ERA-Interim surface data with the gridded lapse rate data of the GPT2w model. Complicated formulas have been deduced for the vertical reductions of pressure and PWV, based on realistic assumptions for atmosphere structure. The global mean RMS errors of derived P, T, and e from 2.5 • × 2.5 • grid surface analysis data are about 1.16 hPa, 1.95 K, and 1.76 hPa, respectively. Both P and T demonstrate lower accuracy in the higher latitudes and some places of middle latitudes of NH. The derived water vapor pressure exhibits positive biases over most stations. Vertical reduction approach of P proposed in this study is superior to the Berg function for over 80% of the IGRA stations, especially when there exists great height difference between the grid point and the observation station. The average improvement is about 0.27 hPa for absolute height differences less than 300 m while it increases to about 2.08 hPa for those above 300 m.
The ZTDs calculated from the ERA surface P, T, and e data (2.5 • × 2.5 • analysis) show a globally averaged accuracy of 3.26 cm and 3.34 cm respectively while using the A&N model and Saas model. They are comparable to the performance of those calculated from local meteorological observations. The RMS errors of ZTDs calculated from the ERA surface P, T, e, and PWV data range between 1.26 and 1.55 cm for the four reference data sets. The approach using the ERA PWV data together with suitable vertical reduction is of the same level of accuracy with respect to the results from NWM pressure level data.
Some common characteristics are noticed during validations using different reference data sets: for ZTDs calculated from P, T, and e, forecast data perform similarly to analysis data and higher spatial resolution does not result in better accuracy; for ZTDs calculated from P, T, e, and PWV, analysis data and higher spatial resolution both yield slightly but definitely better accuracy. Due to its globally and seasonally constant coefficients, the Saas model shows obvious biases in most zones, especially remarkable underestimations of ZWD over the tropics and the subtropics of NH. With the gridded lapse rates from GPT2w, the A&N model displays lower offset, while an overestimation of ZWD by about 2 cm is also seen over the tropics, implying requirement of further refinement.
Tropospheric parameters determined from NWM surface data is also capable of satisfying the accuracy requirement of many geodetic applications and can facilitate users without local meteorological instruments. Using NWM surface data instead of pressure level data cut down the budget for memory and transmission bandwidth and provides great convenience to remote users.