An improved global zenith tropospheric delay model 1 GZTD 2 considering diurnal variations 2

Abstract. The zenith tropospheric delay (ZTD) is an important atmospheric parameter in the wide application of global navigation satellite systems (GNSS) technology in geoscience. Given that the temporal resolution of the current global zenith tropospheric delay model (GZTD) is only 24 h, an improved model, GZTD2, has been developed by taking the diurnal variations into consideration and modifying the model expansion function. The data set used to establish this model is the global ZTD grid data provided by Global Geodetic Observing System (GGOS) Atmosphere spanning from 2002 to 2009. We validated the proposed model with respect to ZTD grid data from GGOS Atmosphere, which was not involved in modeling, as well as International GNSS Service (IGS) tropospheric product. The obtained results of ZTD grid data show that the global average bias and root mean square (rms) for the GZTD2 model are 0.2 and 3.8 cm, respectively. The global average bias is comparable to that of the GZTD model, but the global average rms is improved by 3 mm. The bias and rms are far better than the EGNOS model and the UNB series models. The testing results from global IGS tropospheric product show the bias and rms (−0.3 and 3.9 cm) of the GZTD2 model are superior to that of GZTD (−0.3 and 4.2 cm), suggesting higher accuracy and reliability compared to the EGNOS model, as well as the UNB series models.


Introduction
Radio space-based geodesy techniques suffer from atmosphere propagation delays, of which the ionospheric delay can be largely eliminated by iono-free, carrier-phase combination techniques (Spilker, 1980), and then the tropospheric delay becomes the main error source.In general, we project the slant delay to zenith direction with mapping function in GNSS navigation and positioning, so modeling the ZTD is a common method to reduce the tropospheric influence on signal traveling.In order to improve the accuracy and efficiency of the application in earth science based on space geodesy techniques, a reliable tropospheric delay model is required.
Some tropospheric delay models are developed to mitigate the tropospheric delay.The traditional models like the Hopfield model (Hopfield, 1969), Saastamoinen model (Saastamoinen, 1973), and Black model (Black, 1978) require realtime meteorological data to reach a correction accuracy better than 10 cm.Given the location and time information, the UNB series models (Collins andLangley, 1997, 1998;Leandro et al., 2006Leandro et al., , 2008) ) and the EGNOS model (Dodson et al., 1999;Penna et al., 2001;Ueno et al., 2001) use the empirical meteorological parameters in the form of the latitude band table to estimate the ZTD with an accuracy of about 5 cm, while the IGGTrop model (Li et al., 2012) is based on the empirical three-dimensional parameters in the form of grids to calculate the ZTD with an accuracy of about 4 cm.However, the IGGTrop model needs a large number of parameters.Then, Li et al. (2015) developed the new versions of IGGtrop named IGGtrop_ri (i = 1, 2, 3) by simplifying the algorithm and lowering the resolution, which substantially reduce the required numbers with a similar accuracy.Krueger (2004Krueger ( , 2005) ) and Schüler (2014) obtained the annual and diurnal coefficients for underlying parameters by fitting every grid point's meteorological parameters' time series of the National Centers for Environmental Prediction (NCEP) atmospheric data, and established two global tropospheric delay models -TropGrid and TropGrid2.The correction accuracy of TropGrid2 is 3.8 cm.Böhm et al. (2015) proposed the Global Pressure and Temperature 2 wet model (GPT2w) as an extension to GPT2 (Lagler et al., 2013) with an improved capability to determine zenith wet delays in blind model.The GPT2w model accounts for the annual and semiannual variations of meteorological parameters, and the validation with IGS data and an extended validation with raytraced delays (Möller et al., 2014) show a high accuracy of about 3.6 cm for GPT2w.However, GPT2w has numerous parameters for storage like the above grid models such as IGGTrop series models and TropGrid series models.Yao et al. (2013) established a global non-meteorological parameters tropospheric delay model GZTD (Global Zenith Tropospheric Delay) based on spherical harmonics using the global zenith tropospheric delay grid data provided by Global Geodetic Observing System (GGOS) Atmosphere.The harmonic function including three terms (mean, annual, and semiannual) is used to fit the ZTD time series from 2002 to 2009 for each grid, then the fitted coefficients of all the grids are expanded with 10th-order and 10-degree spherical harmonics.Its modeling approach was very simple, and the overall accuracy of 4.2 cm was similar to the IGGtrop on a global scale, but the required parameters were reduced greatly to about 600.The GZTD model is constructed by global daily average ZTD grid data and the model param- eters were expanded with a low order spherical harmonics, whose temporal resolution is only 1 day in theory and spatial resolution is low.
In this paper, using the ZTD grid data provided by the GGOS Atmosphere, the diurnal variations in ZTD were analyzed to prove the practical necessity for temporal resolution improvement of the GZTD model.Then, on the basis of the GZTD model, and by taking the diurnal variations into consideration and modifying the expansion function, we developed an improved global non-meteorological parameters ZTD model -GZTD2.The data set used to establish this model is the global ZTD grid data provided by the GGOS Atmosphere from 2002 to 2009.Using ZTD grid data obtained from GGOS Atmosphere and tropospheric product (Byun et al., 2009) provided by IGS for model validation, the accuracy of the GZTD2 model is superior to that of the GZTD model, and this model performs much better than other commonly used models such as the EGNOS model and UNB series models.

The new tropospheric delay model
The GGOS Atmosphere is a project that aims to establish atmospheric models, which has been carried out at Vienna University of Technology and has been funded by the Austrian Science Fund (Böhm and Schuh, 2013).It provides grid data of global zenith delays (including zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD)) with a temporal resolution of 6 h (00:00, 06:00, 12:00, 18:00 UTC) and spatial resolution of 2.5 • × 2 • (long × lat), which are derived from the reanalysis data (Uppala et al., 2005) provided by the ECMWF.The ZTD grid data can be obtained by simply  We can see clearly from Fig. 1 that the ZTD estimates of the GZTD model can almost be fitted with a straight line parallel to the time axis which only varies about 1 mm in a single day.The real variations of GGOS grid ZTD data are mostly up to centimeter level, which is 1 order larger than the variations of the GZTD model estimates.Furthermore, we calculated the mean diurnal ZTD values of these six GGOS grid points over the entire year 2010 (Fig. 2), and the significant signal of diurnal variation can be seen at all these six grid points.We can draw a conclusion that the GZTD model could not reflect the characteristic of diurnal variations in ZTD, so the model estimations have nearly no difference when calculating real value or corresponding integer value of the input DOY.Therefore, it is necessary to improve the temporal resolution of the GZTD model to reflect diurnal variations.It should be noted that Jin et al. (2009) has investigated the diurnal and semidiurnal variations in ZTD which were obtained from a decade of global GPS observations, and thought that the atmospheric tides were the major driver of these variations after finding the general similarities of diurnal variations between ZTD and pressure.However, the semidiurnal variations could hardly be described because of the low temporal resolution (6 h) of GGOS ZTD data, so we did not consider the semidiurnal components of ZTD in modeling in the following section.

Establishment of the GZTD2 model
According to the previous research conducted by Jin et al. (2007) and Yao et al. (2013), ZTD decreases exponentially with increasing height, is featured by 1-year periodicity and half-year periodicity, and has a strong correlation with lati- Here, In Eq. ( 1), DOY is the day of the year; hod is the UTC time; h is the height (altitude); a 0 is the annual mean of ZTD on the mean sea level (MSL); a 1 is the annual variation amplitude of ZTD; a 2 is the initial phase of annual variation; a 3 is the semiannual variation amplitude of ZTD; a 4 is the initial phase of semiannual variation; a 5 is the diurnal periodic variation amplitude of ZTD; a 6 is the initial phase of diurnal variation; β = −0.00013137 is the constant to reduce the ZTD in height compared to the MSL, which was determined by Yao et al. (2013) by fitting the global GGOS grid ZTD via exponential function with respect to height; P nm are the Legendre polynomials; ϕ is the latitude of the grid point; λ is the longitude of the grid point; A i nm and B i nm are the coefficients of spherical harmonics determined by least-square optimization.
For each grid-point-specific ZTD time series derived from GGOS Atmosphere, we used Eq.(1) to fit them to temporal coefficients at MSL.Our previous GZTD model only accounts for the annual and semiannual variations of ZTD, whose first equation is similar to Eq. ( 1) but without the fourth term (diurnal term) on the right of Eq. (1).However, there are seven coefficients for each grid, which need large storage space on global scale.Then, referring to the idea of   2) to express the temporal coefficients (mean, annual terms, etc.) of all grids as a function of location (latitude, longitude, and height), thus reducing the parameters.In contrast with the GZTD model established using daily average global ZTD data, we utilized the ZTD time series data of four moments per day (00:00, 06:00, 12:00, 18:00 UTC) from 2002 to 2009, provided by GGOS Atmosphere, to fit ZTD values to obtain temporal variation parameters via Eq.( 1), then expanded these parameters with an 18th-order and 18-degree spherical harmonic function (Eq.2), respectively.The expansion equation of the GZTD model is a 10th-order and 10-degree spherical harmonic function which is 8 orders and degrees less than Eq. ( 2).We used this spherical harmonic function instead of the 10th-order and 10-degree function adopted in the GZTD model because it is not sufficient to apply the previous 10th-order function for the expansion of the temporal variation parameters with relatively high resolution.The number of orders and degrees of spherical harmonics determine the horizontal resolution of model.However, higher order and degree bring more parameters for model.The 10 spherical harmonics adopted by GZTD result in a resolution of about 18 degrees, which is too low for the GZTD2 model to reflect the small diurnal variations.To keep a balance between the resolution and number of parameters, we used 18 spherical harmonics for GZTD2 whose resolution is about 10 degrees.
Figure 3 shows the global distributions of the annual mean of ZTD on MSL and amplitude parameters after fitting by Eq. ( 1).As can be seen from Fig. 3a, the coefficient a 0 in low latitudes, especially near the equator, is significantly larger than in high latitudes, and the distribution in the Southern Hemisphere is more uniform than in the Northern Hemisphere; these results are mostly in agreement with the results of Li et al. (2012) and Yao et al. (2013).For the sawtooth shape in the 40 • N-40 • S region, Yao et al. (2013) found that this shape appears in coastal areas and is consistent with the directions of equatorial trade winds, so they assumed that the distributions of ZTD are effected by some physical impacts such as terrains and heat circulation.Compared with the previous discovery, the sawtooth shape in Fig. 3a is more evident, indicating that the GZTD2 model incorporates these physical impacts.Figure 3b and c show the global distributions of annual amplitude and semiannual amplitude, respectively, both of which are more uniform in the Southern Hemisphere than in the Northern Hemisphere, which is probably due to the fact that most parts of the Southern Hemisphere are covered by oceans, while the Northern Hemisphere has many seacoast regions which lead to relatively complex spatial variation.
Figure 3d shows the global distribution of diurnal variation amplitudes.It can be seen that diurnal variation amplitudes are less than 3 mm in most parts of the world, but up to 1 cm in some low-latitude equatorial areas such as Central America, South America, central Africa, and tropical Asia, indicating notable diurnal variations in these areas.The distribution characteristics of diurnal variation amplitudes is similar to the results of Jin et al. (2009).So taking these diurnal variations into consideration in the GZTD2 model is quite reasonable and necessary.
The GZTD2 model only needs DOY, UTC time, latitude, longitude, and height as input parameters in practical application.GZTD2 uses Eq. (2) to derive temporal parameters a 0 , a 1 , a 2 a 3 , a 4 , a 5 , and a 6 , which are then entered into Eq.( 1) together with the DOY to get the ZTD at MSL.The realization of the GZTD2 model is simple with a few parameters, and the calculation is convenient without inputting any real-time meteorological parameters.Table 1 summarizes the main improvements and features of the newly suggested model compared to the GZTD model.3

Validation and analysis of the GZTD2 model
To analyze the effectiveness and reliability of the new model and verify its accuracy and stability on global scale, as well as to compare it with the GZTD model, this section will exploit some data sources to conduct model validation.Two kinds of data sources are used here: the first is ZTD grid data from GGOS Atmosphere which is not used in modeling, and the other is tropospheric product data provided by IGS.The accuracy is characterized with the average deviation (bias) and root mean square (rms) which are usually used for model validation (Yao et al., 2013;Li et al., 2015;Böhm et al., 2015).The expressions of bias and rms are bias where ZTD M i is the value estimated by the model and ZTD 0 i is the reference value.
As can be seen from Table 2, for the total 13 104 points involved in the global validation, the GZTD2 model's mean bias is 0.2 cm with a maximum of 6.2 cm, and the average of rms is 3.8 cm with a maximum of 8.3 cm, significantly better than the EGNOS and UNB series models, and the rms is reduced by 3 mm compared with that of the GZTD model.The UNB3m model's accuracy is about 1 cm better than the UNB3 and EGNOS models, so we only chose the UNB3m as the representative of commonly used model in our following comparison analysis.Figure 4 shows the global distributions of bias and rms of the three models.
As can be seen from Fig. 4, compared with the other two models, the new model has better accuracy on the worldwide scale, and the accuracy of the areas where larger errors appear improves significantly.Compared with the GZTD model, the GZTD2 model improves the accuracy in the equator area.Obviously, all these three models have suffered large errors in the Pacific Ocean near the equator and Indian Ocean.These areas are near the equator where the deep moist convection effects related to the change of ZTD are more intense (Trenberth et al., 2005;Pramualsakdikul et al., 2007), so the weather changes in these areas are more complex compared with other areas, resulting in difficulty for modeling tropospheric delay.In addition, the GZTD2 and GZTD models are comparable in the Northern and Southern hemispheres, but the UNB3m model's accuracy is obviously lower in the Southern Hemisphere; this is because the UNB3m model is based on the assumptions that tropospheric delay is symmetrical with the equator (Leandro et al., 2006).In fact, this assumption is not reasonable enough and the modeling data sources are derived from North America, so the accuracy of the model is higher in the Northern Hemisphere, especially in North America.

Validation with IGS tropospheric delay data
IGS has provided final troposphere products with a temporal resolution of 5 min since 1998.In 2010, some IGS sites have the severe problem of ZTD data missing.For a convinced validation, only the IGS sites with at least 120 days (approximately a third of the year) of tropospheric delays are selected.Consequently, there are 362 IGS sites selected in 2010 to verify the accuracy of the GZTD2 model, and the distribution of IGS sites is shown in Fig. 5.The uncertainties of the ZTD products are very small (see Fig. 5) with a mean value of 1.5 mm, indicating high quality of the ZTD products.
Considering the ZTD products of IGS sites as true value, we tested and analyzed the ZTD estimates of the GZTD2 model, the GZTD model, the EGNOS model, and the UNB series models.The bias and rms statistical results are shown in Table 3.
As can be seen from Table 3, in terms of the results of accuracy and stability testing for all IGS sites throughout the year, the GZTD2 model performs with the best average rms, and then the GZTD model follows.Global correction accuracy of the new model reaches centimeter level: bias average value is −0.3 cm, average rms is 3.9 cm.Compared with the GZTD model, the range of bias of the GZTD2 model reduces by 2.4 cm and the maximum rms of the GZTD2 model decreases by 0.2 cm, indicating that the new model has a higher  stability.Bias and rms of the EGNOS model are very close to those of the UNB3 model and both are worse than the UNB3m, which is similar to the results of Li et al. (2012).To display the correction effects of different models in a more intuitive way, we computed the distributions of bias and rms of all IGS stations.Figure 7 shows the histograms of bias and rms for the three models.
As can be seen from Fig. 7, the bias of the GZTD2 model concentrates in the range of [−3 cm, 3 cm], while the main distribution range of the bias of the GZTD model is 1 cm larger, and the bias for the UNB3m is distributed with the range of more than 8 cm.It indicates that the GZTD2 model and the GZTD model have small systematic deviations compared with IGS data on a global scale, with the former performing better than the latter, but problematic systematic deviations exist in the UNB3m model within some special areas.Figure 7 also shows that the rms of the GZTD2 model is mostly around 4 cm, and its distribution is more concentrated compared to the GZTD model, indicating the GZTD2 model has higher stability than GZTD.The rms of the UNB3m model is mainly around 5 cm and exceeds 9 cm at many sites, which further suggests the existence of systematic deviations in certain areas in the UNB3m model.
To further analyze the accuracy of the different models varying with location, Fig. 8 shows the global distributions of bias and rms calculated from different models for IGS sites.As can be seen from Fig. 8, the GZTD2 and GZTD models largely eliminate the effects caused by latitude and longitude variations, and the former is more stable than the latter in terms of global distribution of bias and rms in spite of a few sites with large relative error, of which most sites are located in the ocean and coastal areas.A more clear comparison, in terms of rms, between GZTD and GZTD2 is shown in Fig. 9.The reduction for rms can be found at most sites (the number is 273) when moving from GZTD to GZTD2, which accounts for 75.4 % of all sites.The significant improvements of rms are found at the sites in low-latitude areas such as the Pacific Ocean, the coast of South America, and the coast of west Africa where the diurnal variations are notable (see Fig. 3d).This result proves the reasonability of adding diurnal variations in GZTD2.For the UNB3m model, as it is presented in Fig. 8, biases are negative in most parts of the Northern Hemisphere and positive in most parts of the Southern Hemisphere with significantly larger deviations, and rms are smaller for areas in the latitudes higher than 30 • , again suggesting that the correction effect of the UNB3m model is regional.
Figure 10 shows the global distribution of bias and rms with respect to height for the GZTD2 model, the GZTD model, and the UNB3m model.As can be seen, the bias and rms are larger with a height less than 500 m for all three models.Between 500 and 2000 m height, the bias and rms of the GZTD model and the GZTD2 model perform better than those of the UNB3m model, and the overall correction effects of the GZTD and GZTD2 models are also better than that of the UNB3m model.Due to the same exponential function and reducing constant for height, the distribution patterns of the bias and rms of the GZTD and GZTD2 models with respect to height are roughly similar, but the latter is obviously superior to the former.
For a more comprehensive analysis of the relationship between model stability and height, Fig. 11 presents the global distribution of relative rms for three models with respect to height.The relative rms is the ratio of the rms to the annual mean ZTD at the site.Basically, a relative accuracy between 1 and 2.5 % can usually be stated for the majority of the sites from the GZTD2 model, and the relative accuracy is less than 3 % for the GZTD model, showing that both perform better than the UNB3m model.
Figure 12 illustrates the comparisons between IGS ZTD data and ZTDs determined by the UNB3m, GZTD, and GZTD2 models over the year 2010 at site KOUR and TWTF.During the entire year 2010, the ZTD values estimated by the GZTD2 model show the best agreement with the IGS data, which is better than that of the GZTD model without diurnal terms.The ZTDs determined by the UNB3m model vary slightly throughout the year 2010, thus resulting in poor performance.The results in Fig. 12 indicates that the GZTD2 model has a temporal stability for correction accuracy.
From the above analysis, we can conclude that the overall accuracy of the GZTD2 model is up to centimeter level.The GZTD2 model is substantially superior to other commonly used models in terms of bias and rms, and the accuracy improves significantly compared with the GZTD model, thus performing a higher reliability and stability.

Conclusions
In this paper, using the time series data of global tropospheric zenith delays provided by GGOS Atmosphere, we analyzed the diurnal variation in the ZTD which is neglected in the previous GZTD model, then we modified the model function to develop an improved model named GZTD2.We conducted external validation testing with GGOS ZTD grid data which were not involved in modeling, and with IGS tropospheric product.The testing results of GGOS ZTD grid data show that the global average bias and rms for the GZTD2 model are 0.2 and 3.8 cm, respectively.The global average bias is comparable to that of the GZTD model, but the global average rms has been reduced by 0.3 cm.Both the bias and rms are far better than the EGNOS model and the UNB series models.The testing results of global IGS tropospheric product show that the bias and rms for the GZTD2 model are −0.3 and 3.9 cm, superior to those of GZTD (−0.3 and 4.2 cm), indicating higher accuracy and reliability compared to the EGNOS model and the UNB series models.
Overall, compared to the GZTD model, the GZTD2 model improves the temporal resolution and spatial resolution by considering diurnal periodic variations and modifying the expansion function, further completing and optimizing the theory of model establishment.The reliability and stability for the GZTD2 model are much better than other commonly used models.However, like other empirical models such as the UNB3m, the GZTD2 model would be inaccurate in extreme weather events.The Saastamoinen model is recommended if the real-time meteorological observations are available under extreme weather events.Moreover, the GZTD2 model does not consider the semidiurnal variations due to the temporal resolution of GGOS data.In order to build a global tropospheric model with high accuracy, ZTD data with high quality and resolution are required, and the diurnal and semidiurnal variations as well as the subtle secular variation trend of ZTD need more detailed and further study.

Figure 1 .
Figure 1.The GZTD model estimates (blue circles) and corresponding GGOS grid values (green triangles) at the first DOY of 2010.

Figure 2 .
Figure 2. Mean diurnal ZTD values of GGOS grid points with error bars denoting the standard deviations from the average over the year 2010.

Figure 3 .
Figure 3.The global distribution of the annual mean ZTD on MSL (a), the annual variation amplitude (b), the semiannual variation amplitude (c), and the diurnal variation amplitude (d).

Figure 5 .
Figure 5. Distribution of global IGS sites involved in validation.

Figure 6 .
Figure 6.Histogram of uncertainty of ZTD at selected IGS sites.

Figure 7 .
Figure 7. Histograms of bias and rms for three models.

Figure 10 .
Figure 10.Global distributions of bias and rms for different models with respect to height.

Figure 11 .
Figure 11.Relative rms for different models with respect to height.

Table 1 .
Improvements of GZTD2 with respect to GZTD.

Table 2 .
Modeling errors of different models validated by GGOS data.

Table 3 .
Error of different considered models versus IGS data.