A New Grid Zenith Tropospheric Delay Model Considering Time-Varying Vertical Adjustment and Diurnal Variation over China

: Improving the accuracy of zenith tropospheric delay (ZTD) models is an important task. However, the existing ZTD models still have limitations, such as a lack of appropriate vertical adjustment function and being unsuitable for China, which has a complex climate and great undulating terrain. A new approach that considers the time-varying vertical adjustment and delicate diurnal variations of ZTD was introduced to develop a new grid ZTD model (NGZTD). The NGZTD model employed the Gaussian function and considered the seasonal variations of Gaussian coefficients to express the vertical variations of ZTD. The effectiveness of vertical interpolation for the vertical adjustment model (NGZTD-H) was validated. The root mean squared errors (RMSE) of the NGZTD-H model improved by 58% and 22% compared to the global pressure and temperature 3 (GPT3) model using ERA5 and radiosonde data, respectively. The NGZTD model’s effectiveness for directly estimating the ZTD was validated. The NGZTD model improved by 22% and 31% compared to the GPT3 model using GNSS-derived ZTD and layered ZTD at radiosonde stations, respectively. Seasonal variations in Gaussian coefficients need to be considered. Using constant Gaussian coefficients will generate large errors. The NGZTD model exhibited outstanding advantages in capturing diurnal variations and adapting to undulating terrain. We analyzed and discussed the main error sources of the NGZTD model using validation of spatial interpolation accuracy. This new ZTD model has potential applications in enhancing the reliability of navigation, positioning, and interferometric synthetic aperture radar (InSAR) measurements and is recommended to promote the development of space geodesy techniques.


Introduction
The troposphere, which ranges from the ground to the beginning of the stratosphere, is a layer of atmosphere near the Earth.Air in the troposphere accounts for 75% of the total atmospheric mass.The troposphere's height over China is approximately 10 km.The global navigation satellite system (GNSS) provides high-precision three-dimensional coordinates and speed.Tropospheric delay limits the accuracy of space geodesy [1].GNSS technology inevitably generates errors in its applications owing to the tropospheric delay, which seriously affects the accuracy of navigation, positioning, and interferometric synthetic aperture radar (InSAR) [2].The zenith tropospheric delay (ZTD) consists of the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD).Owing to the particularities of the troposphere, the ZHD with a stable variation law can be accurately modeled.Unlike ZHD, ZWD, which is caused by precipitable water vapor (PWV), is difficult to model accurately [3].In addition, ZTD can be a useful signal for PWV inversion [4].Therefore, research on ZTD modeling is of considerable significance for navigation, positioning, InSAR monitoring, and PWV inversion [5].
Existing ZTD models include meteorological and non-meteorological parameter models.The meteorological parameter models rely on measured meteorological parameters, whereas the non-meteorological parameter models only require the input of spatiotemporal information.The Hopfield [6], Saastamoinen [7], and Black models [8] are the main meteorological parameter models.Hopfield established a ZTD model called the Hopfield model using data from 18 radiosonde stations worldwide; this model requires meteorological parameters.Saastamoinen established the Saastamoinen model based on the U.S. Standard Atmospheric Model (SAM), which also requires input values such as temperature, pressure, and station location information.The Black model was developed as a new-generation model that follows the Hopfield model, which also requires the input of temperature and pressure.The accuracy of the abovementioned models can reach a centimeter level, provided that the measured meteorological parameters are available [9].However, the measured meteorological parameters can only be obtained from stations equipped with meteorological sensors.The distribution of these stations with low spatial resolution is uneven, and there is a time delay, which considerably limits the ability of meteorological parameter models to realize real-time applications.Therefore, real-time calculation of ZTD using non-meteorological parameter models has become a considerable challenge that must be solved, and it has also become a hot research topic for scholars [10,11].The TropGrid and GPT series models were established to realize a higher temporal resolution [12][13][14].Owing to the GPT2 model only calculating certain meteorological parameters, Böhm et al. [15] established the GPT2w model using monthly ERA-Interim data, which adds two output values: water vapor decline rate and atmospheric weighted mean temperature.The global pressure and temperature 3 (GPT3) model is a new-generation model that follows the GPT2w model and improves the empirical mapping function coefficients [16][17][18].The GPT3 model, which has the ability to output comprehensive meteorological parameters, can calculate ZTD by combining the Saastamoinen and Askne models.The GPT3 model provides two horizontal resolutions of 1 • × 1 • and 5 • × 5 • , which need to be further improved.To improve the applicability of the ZTD model globally, studies have been conducted on high-precision global ZTD modeling, which considers multiple factors.Li et al. [19] established an IGGtrop model that considers latitude, longitude, elevation, and annual period, which attained a further improvement in the global ZTD model's accuracy.In addition, Huang et al. [20] proposed a new global ZTD grid model based on a sliding window algorithm with different spatial resolutions, a novel model which shows excellent performance compared to the GPT3 and UNB3m models and significantly optimizes model parameters.However, the above models cannot capture the diurnal variation of ZTD, and the temporal resolution of these models needs to be further improved.
Recently, the calculation of the ZTD using the fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA5) and the Second Modern-Era Retrospective Analysis for Research and Applications (MERRA-2) has received much attention [21][22][23].These atmospheric reanalysis data exhibited high-precision results when they were validated by other reference data [24].Therefore, the ERA5 and MERRA-2 datasets are expected to be used widely in the future.Because these atmospheric reanalysis grid data exhibit a high spatiotemporal resolution, the grid point data around the target point can be interpolated to calculate the data at the user's position with high precision.However, ZTD varies much more vertically than horizontally.Direct interpolation produces large errors in regions with an undulating terrain.Therefore, to solve the aforementioned problems, it is necessary to develop a high-precision model to adjust the ZTD vertically.The selection of a fitting function for the variation in layered ZTD is an important research topic.Therefore, extensive studies have been conducted on ZTD vertical profile functions [25].The vertical variations of the ZTD are typically modeled using polynomials [26,27] and negative exponential functions [28][29][30].Zhu et al. [31] developed a segmented global ZTD vertical profile model (GZTD-P) with different spatial resolutions that considers the time-varying characteristics of the ZTD vertical variation factor to address the limitations of a single function in expressing the ZTD vertical profile, which shows better performance compared to the GPT3 model.However, the GZTD-P model can only vertically adjust ZTD from the starting height to the target height and cannot directly calculate ZTD at the target position, which limits its application.Sun et al. [32] proposed a global ZTD model, the GZTDS model, which considers delicate periodic variations by adopting a nonlinear function.The GSTDS model was developed as a new-generation model that follows the GZTDS model, and the performance of the new model deteriorated as the zenith angle increased.Hu and Yao [33] adopted the Gaussian function to fit the vertical ZTD and then established a ZTD vertical profile model that considers seasonal variation, which shows good performance in both time and space.The model achieved good results on a global scale.The model was developed using the monthly mean ZTD with horizontal resolutions of 5 • × 5 • provided by ERA-Interim.Its adaptability needs to be further improved in regions with greatly undulating terrain and complex climates.Zhao et al. [34] proposed a high-precision ZTD model that considers the height effect on ZTD after analyzing the relationship between the ZTD periodic residual term and the height of the GNSS station at different seasons.Although the aforementioned models have demonstrated their respective advantages, it is necessary to conduct further research on the more delicate vertical and temporal variations in the ZTD over China.
China is characterized by greatly undulating terrain, large latitudinal spans, various climate types, and large diurnal atmospheric differences [35][36][37].Under such conditions, the ZTD exhibits complex variations in both time and space.Existing ZTD models have difficulty meeting the requirements of applications in China with the aforementioned characteristics.Therefore, a detailed investigation of the spatiotemporal variations of the ZTD and the selection of better functions are of considerable significance for applications in GNSS and InSAR.Our aim was to develop a NGZTD model that considers the time-varying vertical adjustment and delicate diurnal variations of ZTD.To attain this objective, this paper is organized as follows.Section 2 introduces the data and methods for calculating ZTD and analyzes the spatiotemporal characteristics of Gaussian coefficients and surface ZTD.Section 3 develops a NGZTD-H model to vertically adjust ZTD, and the vertical interpolation accuracy of the model is validated using ERA5 and radiosonde data.In addition, the NGZTD model is developed to directly calculate ZTD, and its accuracy is validated using GNSS and radiosonde data.China is characterized by greatly undulating terrain, large latitudinal spans, various climate types, and large diurnal atmospheric differences [35][36][37].Under such conditions, the ZTD exhibits complex variations in both time and space.Existing ZTD models have difficulty meeting the requirements of applications in China with the aforementioned characteristics.Therefore, a detailed investigation of the spatiotemporal variations of the ZTD and the selection of better functions are of considerable significance for applications in GNSS and InSAR.Our aim was to develop a NGZTD model that considers the timevarying vertical adjustment and delicate diurnal variations of ZTD.To attain this objective, this paper is organized as follows.Section 2 introduces the data and methods for calculating ZTD and analyzes the spatiotemporal characteristics of Gaussian coefficients and surface ZTD.Section 3 develops a NGZTD-H model to vertically adjust ZTD, and the vertical interpolation accuracy of the model is validated using ERA5 and radiosonde data.In addition, the NGZTD model is developed to directly calculate ZTD, and its accuracy is validated using GNSS and radiosonde data.Section 4 discusses the main error sources of the NGZTD model.Section 5 contains the Conclusions.The research framework is shown in Figure 1.The proposed model compensates for the limitations of existing ZTD models in China.The new model is expected to be applied to precise point positioning (PPP) and InSAR atmospheric correction to improve their monitoring accuracy.

Data
ERA5 is the latest product provided by ECMWF, which is characterized by high temporal and spatial resolution.ERA5 data consist of 37 vertical pressure layers.The ERA5 datasets are generated via assimilation schemes that include observation data from GPS occultation observations, satellite altimeters, satellite radiation, and inversion equipment [38].Jiang et al. [39] systematically evaluated the effectiveness of ZTD data calculated using ERA5 data.The datasets are used as important resources for research in GNSS meteorology and space geodesy [40,41].In this paper, ERA5 data were used to analyze the spatiotemporal characteristics of Gaussian coefficients and surface ZTD.They also provided data sources for developing NGZTD-H and NGZTD models.In the vertical interpolation validation of the NGZTD-H model, ERA5 data were used as the reference values to validate the accuracy of the model.
Radiosonde stations comprise converters, sensors, and radio transmitters, which can provide measured meteorological parameters near the Earth.Radiosonde data are obtained by releasing balloons equipped with sensors and are widely used as reference data to test other data and tropospheric models [42].Radiosonde data were used as reference values to validate the accuracy of the NGZTD-H and the NGZTD model.
The ZTD from GNSS stations can reach an accuracy of 4 mm [43].It is recognized as a truth value and is widely used to evaluate the ZTD derived from other data and ZTD models [44].GNSS-derived ZTD was used as reference value to validate the accuracy of the NGZTD model.It can be used to validate the NGZTD model's ability to capture diurnal variations of ZTD, as the temporal resolution of GNSS data is 1 h.

The Method of Calculating ZTD
The layered ZTD were calculated using the integral method for atmospheric refractivity.The equation used to calculate atmospheric refractivity is expressed as follows:

NdH
(1) where e is the water vapor pressure (hPa); N is the total refractivity; sh is the specific humidity; T is the temperature (K); H is the height (km); h low and h top are the lowest and topmost heights; k 1 , k 2 , k 3 are the refractive constants.The Saastamoinen model was employed to estimate the ZHD at the top of the troposphere, which was used to correct the integral result of the ZTD [29].The equation is expressed as follows: where ZHD top is the ZHD at the top of troposphere (m), P top is the pressure (hPa) at the top of troposphere, and φ is the latitude ( • ).

Analysis of the Gaussian Coefficient
A detailed investigation of the spatiotemporal variations of the ZTD provides important references for developing ZTD models.Direct interpolation results in large errors owing to the inconsistent heights of the four grid points.It is necessary to use appropriate functions to express the variation in ZTD with height.Hu and Yao [33] found that the Gaussian function is superior to the exponential function.Therefore, we used the Gaussian function to fit the vertical variation in ZTD, as expressed by where h is the height (m), a is the ZTD (m) at the mean sea level, and b and c (m) are the height-scale factors (also called the Gaussian coefficients).
To verify the effectiveness of this function in expressing the vertical variation in ZTD, four representative grid points in China were selected to determine the ZTD-layered profile information from ERA5 and fitted vertical variations using the Gaussian function.
From Table 1, the mean bias of the four grid points was negative, indicating that a certain systematic error occurs in the Gaussian function to fit the ZTD profiles, but the value was small: approximately −1.00 mm.The root mean squared errors (RMSE) of the four grid points were 14.56, 10.24, 10.47, and 8.86 mm, respectively.From the above results, the Gaussian function attained high accuracy and stability to fit the ZTD profiles.Analyzing the temporal variation characteristics of the Gaussian coefficients b and c is critical for developing a ZTD vertical adjustment model.Therefore, the annual mean values and amplitudes of Gaussian coefficients b and c were calculated using the ERA5 ZTD to analyze the distribution characteristics of Gaussian coefficients b and c over China.The results are shown in Figure 2. To calculate these values, the ZTD-layered profiles calculated using ERA5 atmospheric reanalysis data from 2014 to 2017 were fitted using Equation (5) to obtain b and c for all grid points.Second, the Gaussian coefficients b and c for each grid point were fitted by considering the annual and semi-annual periods using Equation (6) [20].Finally, the annual mean values and amplitudes of Gaussian coefficients b and c were obtained using Equations (7) and (8).(6) where i represents the i-th grid point, co f i represents Gaussian coefficients b or c, α i 0 is the mean value of Gaussian coefficients, α i 1 to α i 4 are the period amplitude coefficients, and doy is the day of the year.
where amp a represents the annual period amplitude (m) of b and c and amp s represents the semi-annual period amplitude (m) of b and c.As shown in Figure 2, the characteristics of the Gaussian coefficient b varied across different regions.The annual mean values of the Gaussian coefficient b were larger in northwest China and smaller in southeast China, showing a decreasing trend from northwest to southeast China.The annual period amplitude gradually decreased from north to south, and there was a sudden increase in the Tibetan Plateau.This phenomenon may be caused by climatic differences between northern and southern China and differences between the ocean and land.The regional variation characteristics of the amplitude of the semi-annual period amplitude were subtle.In addition, the regional distribution characteristics of the annual mean values and amplitudes of Gaussian coefficient c were similar to those of Gaussian coefficient b.The annual mean values of Gaussian coefficient c showed an apparent regional variation law, as did the annual period amplitude, whereas the regional variation law of the semi-annual period amplitude was also subtle.Because the period characteristics of Gaussian coefficients b and c varied considerably across different regions, the delicate periods of Gaussian coefficients b and c must be considered for each grid point when developing a ZTD vertical adjustment model.
Remote Sens. 2024, 16, 2023 6 of 21 the regional variation law of the semi-annual period amplitude was also subtle.Because the period characteristics of Gaussian coefficients b and c varied considerably across different regions, the delicate periods of Gaussian coefficients b and c must be considered for each grid point when developing a ZTD vertical adjustment model.

Analysis of the Surface ZTD
Analyzing the delicate temporal variations of ZTD contributes to the provision of important references for optimizing the ZTD model.The surface ZTD has significant annual and semi-annual periods [34].Two representative grid points in China were selected to determine the diurnal and semi-diurnal periods by adopting the fast Fourier transform.
As shown in Figure 3, the surface ZTD exhibited significant diurnal and semi-diurnal period variations.The diurnal variation trends of the surface ZTD for the two grid points were different.For the first grid point, the surface ZTD increased and then decreased throughout the day, whereas the second grid point showed the opposite result.The diurnal period variation characteristics of the first grid point were significantly stronger than those of the semi-diurnal period, whereas the diurnal and semi-diurnal period variation characteristics of the second grid point were comparable.Therefore, annual, semi-annual, diurnal, and semi-diurnal periods must be considered when developing a ZTD grid model.
As shown in Figure 4, the annual mean surface ZTD was larger on the Tibetan Plateau and smaller elsewhere, which may be due to altitudinal differences.The annual period amplitude of the surface ZTD showed a decreasing trend from southeast to northwest, which may be because the northwest is far from the ocean and the southeast is close to the ocean, resulting in a large difference in the type of climate.The semi-annual period amplitude of the ZTD did not exhibit an apparent variation law.The distributions of the diurnal and semi-diurnal period amplitudes were similar.Therefore, these results further demonstrate that the delicate period of the ZTD must be considered when developing a ZTD grid model.As shown in Figure 3, the surface ZTD exhibited significant diurnal and semi-diurnal period variations.The diurnal variation trends of the surface ZTD for the two grid points were different.For the first grid point, the surface ZTD increased and then decreased throughout the day, whereas the second grid point showed the opposite result.The diurnal period variation characteristics of the first grid point were significantly stronger than those of the semi-diurnal period, whereas the diurnal and semi-diurnal period variation characteristics of the second grid point were comparable.Therefore, annual, semi-annual, diurnal, and semi-diurnal periods must be considered when developing a ZTD grid model.
As shown in Figure 4, the annual mean surface ZTD was larger on the Tibetan Plateau and smaller elsewhere, which may be due to altitudinal differences.The annual period amplitude of the surface ZTD showed a decreasing trend from southeast to northwest, which may be because the northwest is far from the ocean and the southeast is close to the ocean, resulting in a large difference in the type of climate.The semi-annual period amplitude of the ZTD did not exhibit an apparent variation law.The distributions of the diurnal and semi-diurnal period amplitudes were similar.Therefore, these results further demonstrate that the delicate period of the ZTD must be considered when developing a ZTD grid model.

Development of the NGZTD-H Model
According to the above analysis, both Gaussian coefficients b and c have annual and semi-annual period characteristics.Therefore, these coefficients of each ERA5 grid point were fitted by considering the annual and semi-annual periods.Five coefficients were calculated using Equation ( 6), which were stored in grid points with a horizontal resolution of 0.25° × 0.25°.After obtaining Gaussian coefficients b and c, the ZTD was adjusted vertically using the following formula [33]: where ℎ and ℎ are the starting and target heights (m), respectively, and ZTD and ZTD are the ZTD values (m) of the starting and target heights.

Development of the NGZTD-H Model
According to the above analysis, both Gaussian coefficients b and c have annual and semi-annual period characteristics.Therefore, these coefficients of each ERA5 grid point were fitted by considering the annual and semi-annual periods.Five coefficients were calculated using Equation ( 6), which were stored in grid points with a horizontal resolution of 0.25° × 0.25°.After obtaining Gaussian coefficients b and c, the ZTD was adjusted vertically using the following formula [33]: where ℎ and ℎ are the starting and target heights (m), respectively, and ZTD and ZTD are the ZTD values (m) of the starting and target heights.

Development of the NGZTD-H Model
According to the above analysis, both Gaussian coefficients b and c have annual and semi-annual period characteristics.Therefore, these coefficients of each ERA5 grid point were fitted by considering the annual and semi-annual periods.Five coefficients were calculated using Equation ( 6), which were stored in grid points with a horizontal resolution of 0.25 • × 0.25 • .After obtaining Gaussian coefficients b and c, the ZTD was adjusted vertically using the following formula [33]: where h r and h t are the starting and target heights (m), respectively, and ZTD r and ZTD t are the ZTD values (m) of the starting and target heights.
The ZTD vertical profile model (NGZTD-H model) consists of Equations ( 6) and ( 9).This model requires the following simple steps for application.First, after extracting the five model coefficients of the given grid point and inputting the doy, the Gaussian coefficients b and c were calculated using Equation (6).Second, after inputting the starting and target heights, the ZTD of the target height was calculated using Equation (9).The NGZTD-H model can achieve vertical adjustment of the ZTD from the starting height to the target height before horizontal interpolation, which reduces the interpolation errors.
To validate the NGZTD-H model, the effectiveness of vertical interpolation was evaluated using ERA5 and radiosonde data.

Accuracy Validation in ZTD Vertical Interpolation Using ERA5 Data
The proposed vertical profile model (NGZTD-H model) and the GPT3 model were evaluated using the ZTD layered profile information from ERA5 in 2018 over China.The UNB3 model can adjust the ZHD and ZWD from sea-level height to the target height.The GPT3 model provides the meteorological parameters required for the UNB3 model.Therefore, combining the GPT3 and UNB3 models can achieve a vertical adjustment for the ZTD.As the ERA5-layered ZTD recognized as a reference value does not start from sea-level height, it is necessary to derive the UNB3 model such that it can achieve vertical adjustment for ZHD and ZWD at any starting height.The derived formulas of the UNB3 model are as follows: where the subscript letters t and r represent the target and the starting heights, β is the lapse rate of temperature, T 0 is the thermodynamic temperature (K), R d is the dry gas constant with values of 287.0538J•kg −1 , g represents the ground gravity acceleration, and λ ′ is the water vapor lapse rate.
The NGZTD-H model was used to adjust the ZTD vertically from the surface to the ERA5 layer height.Because the GPT3 model cannot vertically adjust the ZTD directly, the surface ZHD and ZWD derived from the ERA5 data at the grid point need to be adjusted vertically and then the ZHD and ZWD must be added together to obtain the ZTD.Bias and RMSE were recognized as precision criteria.The bias and RMSE results of the NGZTD-H and GPT3 models are presented in Table 2.As shown in Table 2, the mean bias of the NGZTD-H model and the GPT3 model were 0.45 and −3.02 cm, which indicated that the GPT3 model exhibits a larger systematic error.Furthermore, the mean RMSE of NGZTD-H model was 1.70.Compared with the GPT3 model, the accuracy of NGZTD-H model improved by 10% over the ERA5 ZTD profile.Therefore, the NGZTD-H model exhibits better stability for vertically adjusting the ZTD.
The two models were used to calculate the hourly layered ZTD at the ERA5 grid points and then to calculate the annual mean bias and RMSE to analyze the geographic distribution of the accuracy.A is shown in Figure 5.The bias of the NGZTD-H model was approximately 0 cm in most regions, negative in high-elevation regions, and positive in Remote Sens. 2024, 16, 2023 9 of 21 marine regions.The bias of the GPT3 model was positive at high latitudes and negative at low latitudes, especially in the marine regions, where it showed a large negative bias of approximately −5 cm.In addition, the RMSE of the NGZTD-H model was approximately 2 cm and no more than 3 cm in most regions, attaining a better performance in Southwest China.The GPT3 model exhibited low accuracy both in land and marine regions.The two models were used to calculate the hourly layered ZTD at the ERA5 grid points and then to calculate the annual mean bias and RMSE to analyze the geographic distribution of the accuracy.A is shown in Figure 5.The bias of the NGZTD-H model was approximately 0 cm in most regions, negative in high-elevation regions, and positive in marine regions.The bias of the GPT3 model was positive at high latitudes and negative at low latitudes, especially in the marine regions, where it showed a large negative bias of approximately −5 cm.In addition, the RMSE of the NGZTD-H model was approximately 2 cm and no more than 3 cm in most regions, attaining a better performance in Southwest China.The GPT3 model exhibited low accuracy both in land and marine regions.To analyze the effectiveness of the NGZTD-H model in different atmospheric pressure layers and latitudes, the ERA5-layered ZTD was also recognized as a reference value in 2018 over China.Five pressure levels were selected to evaluate the model.In addition, the Chinese region was divided into four latitude zones at 10° intervals.Finally, the accuracy of the two models were validated for the selected pressure layers and latitude bands.
As shown in Figure 6, the NGZTD-H model exhibited a small bias value in all pressure layers, positive values in two pressure layers at 800 and 200 hPa, and a negative value in the rest of the pressure layers, which were all approximately 0 cm.The GPT3 model exhibited a large value in all pressure layers; the value reached −5 cm in the pressure layers with 600 and 400 hPa.The bias of the GPT3 model was negative in all pressure layers, indicating a negative systematic error.Moreover, the RMSE of the NGZTD-H model was approximately 2 cm for all the pressure layers.It attained a better performance in 200 and 10 hPa pressure layers, which indicated that the NGZTD-H model is more precise at larger vertical interpolation heights.However, the GPT3 model exhibited larger errors in most pressure layers, with an RMSE of approximately 5 cm.In addition, the bias of the NGZTD-H model was within 1 cm for all latitude bands.The bias of the GPT3 model reached −4 cm at lower latitudes.This finding may be due to the variations of meteorological parameters in the troposphere being more complicated at low latitudes, which causes difficulties in capturing their variations using the GPT3 model.The bias of the GPT3 model also To analyze the effectiveness of the NGZTD-H model in different atmospheric pressure layers and latitudes, the ERA5-layered ZTD was also recognized as a reference value in 2018 over China.Five pressure levels were selected to evaluate the model.In addition, the Chinese region was divided into four latitude zones at 10 • intervals.Finally, the accuracy of the two models were validated for the selected pressure layers and latitude bands.
As shown in Figure 6, the NGZTD-H model exhibited a small bias value in all pressure layers, positive values in two pressure layers at 800 and 200 hPa, and a negative value in the rest of the pressure layers, which were all approximately 0 cm.The GPT3 model exhibited a large value in all pressure layers; the value reached −5 cm in the pressure layers with 600 and 400 hPa.The bias of the GPT3 model was negative in all pressure layers, indicating a negative systematic error.Moreover, the RMSE of the NGZTD-H model was approximately 2 cm for all the pressure layers.It attained a better performance in 200 and 10 hPa pressure layers, which indicated that the NGZTD-H model is more precise at larger vertical interpolation heights.However, the GPT3 model exhibited larger errors in most pressure layers, with an RMSE of approximately 5 cm.In addition, the bias of the NGZTD-H model was within 1 cm for all latitude bands.The bias of the GPT3 model reached −4 cm at lower latitudes.This finding may be due to the variations of meteorological parameters in the troposphere being more complicated at low latitudes, which causes difficulties in capturing their variations using the GPT3 model.The bias of the GPT3 model also showed a decreasing trend with increasing latitude.Additionally, The NGZTD-H model exhibited good stability in all latitude bands, and its RMSE ranged from 1.5 to 2 cm.The RMSE of the GPT3 model reached 5 cm at lower latitudes.Therefore, compared with the GPT3 model, the NGZTD-H model is more suitable for the vertical adjustment of the ZTD over China.
Remote Sens. 2024, 16, 2023 10 of 21 showed a decreasing trend with increasing latitude.Additionally, The NGZTD-H model exhibited good stability in all latitude bands, and its RMSE ranged from 1.5 to 2 cm.The RMSE of the GPT3 model reached 5 cm at lower latitudes.Therefore, compared with the GPT3 model, the NGZTD-H model is more suitable for the vertical adjustment of the ZTD over China.

Accuracy Validation in ZTD Vertical Interpolation Using Radiosonde Data
The layered ZTD at radiosonde stations in 2018 over China were considered as reference values to further validate the effectiveness of the NGZTD-H model.The ZTD at the radiosonde station was the starting value for the NGZTD-H model.The surface ZHD and ZWD at the radiosonde station were the starting values for the GPT3 model.The ZHD and ZWD, adjusted vertically using the GPT3 model, were added together to obtain the layered ZTD, which was used for comparison with the NGZTD-H model.The Gaussian coefficients b and c required for the two models at the radiosonde stations were calculated via interpolation from four ERA5 grid points around the radiosonde stations.Finally, the two models were used to calculate the layered ZTD for radiosonde stations over China at UTC 00:00 and 12:00 in 2018.
Table 3 shows the mean bias of the NGZTD-H and GPT3 models were −0.32 and −2.66 cm, respectively, indicating that the ZTD from radiosonde stations was larger than the ZTD calculated by the GPT3 model.The mean RMSE of the NGZTD-H and GPT3 models were 3.10 and 3.97 cm.Compared with the GPT3 model, the accuracy of the NGZTD-H model enhanced by 0.87 cm (22%).The NGZTD-H model had a higher accuracy for the vertical adjustment of the ZTD.Compared with the validation results using the ERA5 profile ZTD, the NGZTD-H model exhibited lower accuracy, which may be due to the interpolation errors of Gaussian coefficients from four grid points to the radiosonde stations.

Accuracy Validation in ZTD Vertical Interpolation Using Radiosonde Data
The layered ZTD at radiosonde stations in 2018 over China were considered as reference values to further validate the effectiveness of the NGZTD-H model.The ZTD at the radiosonde station was the starting value for the NGZTD-H model.The surface ZHD and ZWD at the radiosonde station were the starting values for the GPT3 model.The ZHD and ZWD, adjusted vertically using the GPT3 model, were added together to obtain the layered ZTD, which was used for comparison with the NGZTD-H model.The Gaussian coefficients b and c required for the two models at the radiosonde stations were calculated via interpolation from four ERA5 grid points around the radiosonde stations.Finally, the two models were used to calculate the layered ZTD for radiosonde stations over China at UTC 00:00 and 12:00 in 2018.
Table 3 shows the mean bias of the NGZTD-H and GPT3 models were −0.32 and −2.66 cm, respectively, indicating that the ZTD from radiosonde stations was larger than the ZTD calculated by the GPT3 model.The mean RMSE of the NGZTD-H and GPT3 models were 3.10 and 3.97 cm.Compared with the GPT3 model, the accuracy of the NGZTD-H model enhanced by 0.87 cm (22%).The NGZTD-H model had a higher accuracy for the vertical adjustment of the ZTD.Compared with the validation results using the ERA5 profile ZTD, the NGZTD-H model exhibited lower accuracy, which may be due to the interpolation errors of Gaussian coefficients from four grid points to the radiosonde stations.As is shown in Figure 7, the bias of the NGZTD-H model was small at all radiosonde stations and approximately 0 cm in most regions.Both the NGZTD-H and GPT3 models exhibited a large positive bias at higher elevations in southwest China, which may be due to the undulating terrain.The bias of the GPT3 model reached −8 cm at certain stations, indicating that serious systematic errors occur in the GPT3 model.In addition, the RMSE of the NGZTD-H model was stable without abrupt variations.At higher latitudes, the difference in RMSE between the NGZTD-H and GPT3 models was small.The GPT3 model exhibited a larger RMSE at low latitudes for the reasons noted earlier.As is shown in Figure 7, the bias of the NGZTD-H model was small at all radiosonde stations and approximately 0 cm in most regions.Both the NGZTD-H and GPT3 models exhibited a large positive bias at higher elevations in southwest China, which may be due to the undulating terrain.The bias of the GPT3 model reached −8 cm at certain stations, indicating that serious systematic errors occur in the GPT3 model.In addition, the RMSE of the NGZTD-H model was stable without abrupt variations.At higher latitudes, the difference in RMSE between the NGZTD-H and GPT3 models was small.The GPT3 model exhibited a larger RMSE at low latitudes for the reasons noted earlier.To investigate the adaptability of the NGZTD-H model to seasonal variation, the diurnal mean bias and RMSE were determined at two representative radiosonde stations (Hailar and Hangzhou).As shown in Figure 8, at the Hailar radiosonde station, the GPT3 model, which had poor seasonal adaptation, was seriously affected by seasonal variations, and its bias fluctuated throughout the year.However, the NGZTD-H model could resist the influence of the season on the vertical adjustment of the ZTD, and its bias was stable throughout the year.The RMSE of the new model was 2 cm.At the Hangzhou radiosonde station, two models exhibited fluctuations.The bias of the NGZTD-H model fluctuated at approximately 0 cm, and its RMSE was approximately 2 cm.However, the GPT3 model exhibited a more evident negative bias.Therefore, the NGZTD-H model is more suitable for capturing seasonal variations than the GPT3 model.
Because the vertical adjustment of the ZTD is related to elevation, the radiosonde stations over China were divided into five intervals to validate the effectiveness of the NGZTD-H model at different elevations.Figure 9 shows the bias and RMSE statistics in each interval for the NGZTD-H and GPT3 models.The GPT3 model exhibited a larger positive bias in the interval of 3-4 km and a larger negative bias in the intervals of 0-3 and >4 km.The bias of the NGZTD-H model was smaller in each interval.A larger positive bias occurred in the interval of 3-4 km, indicating that the NGZTD-H model had a larger value than the reference value in this interval.Furthermore, the RMSE of the GPT3 model To investigate the adaptability of the NGZTD-H model to seasonal variation, the diurnal mean bias and RMSE were determined at two representative radiosonde stations (Hailar and Hangzhou).As shown in Figure 8, at the Hailar radiosonde station, the GPT3 model, which had poor seasonal adaptation, was seriously affected by seasonal variations, and its bias fluctuated throughout the year.However, the NGZTD-H model could resist the influence of the season on the vertical adjustment of the ZTD, and its bias was stable throughout the year.The RMSE of the new model was 2 cm.At the Hangzhou radiosonde station, two models exhibited fluctuations.The bias of the NGZTD-H model fluctuated at approximately 0 cm, and its RMSE was approximately 2 cm.However, the GPT3 model exhibited a more evident negative bias.Therefore, the NGZTD-H model is more suitable for capturing seasonal variations than the GPT3 model.
Because the vertical adjustment of the ZTD is related to elevation, the radiosonde stations over China were divided into five intervals to validate the effectiveness of the NGZTD-H model at different elevations.Figure 9 shows the bias and RMSE statistics in each interval for the NGZTD-H and GPT3 models.The GPT3 model exhibited a larger positive bias in the interval of 3-4 km and a larger negative bias in the intervals of 0-3 and >4 km.The bias of the NGZTD-H model was smaller in each interval.A larger positive bias occurred in the interval of 3-4 km, indicating that the NGZTD-H model had a larger value than the reference value in this interval.Furthermore, the RMSE of the GPT3 model was larger in the intervals of 0-1 and 2-3 km.The above analysis shows the NGZTD-H model exhibited more stable performances in vertical adjustment for ZTD over China.
was larger in the intervals of 0-1 and 2-3 km.The above analysis shows the NGZTD-H model exhibited more stable performances in vertical adjustment for ZTD over China.

Development of the NGZTD Model
From the above analysis of spatiotemporal characteristics, the surface ZTD exhibited regular seasonal and diurnal variations.Therefore, the development of a ZTD model must consider the delicate periods of the surface ZTD.This equation is expressed as follows: where ZTD is the surface ZTD (m),  is the annual mean ZTD,  to  are the period amplitude coefficients, and ℎ is the hour (UTC) of day.
The surface ZTD of all ERA5 grid points over China from 2014 to 2017 was used to calculate nine coefficients using the least squares method.The new ZTD grid model (NGZTD model), which directly calculates the ZTD, consists of Equation ( 12) and the NGZTD-H model proposed earlier.The method for using the NGZTD model is as follows: First, the user's coordinates are determined to search the ERA5 grid points around the user.After the model coefficients of these grid points have been obtained, as shown in

Development of the NGZTD Model
From the above analysis of spatiotemporal characteristics, the surface ZTD exhibited regular seasonal and diurnal variations.Therefore, the development of a ZTD model must consider the delicate periods of the surface ZTD.This equation is expressed as follows: where ZTD is the surface ZTD (m),  is the annual mean ZTD,  to  are the period amplitude coefficients, and ℎ is the hour (UTC) of day.
The surface ZTD of all ERA5 grid points over China from 2014 to 2017 was used to calculate nine coefficients using the least squares method.The new ZTD grid model (NGZTD model), which directly calculates the ZTD, consists of Equation ( 12) and the NGZTD-H model proposed earlier.The method for using the NGZTD model is as follows: First, the user's coordinates are determined to search the ERA5 grid points around the user.After the model coefficients of these grid points have been obtained, as shown in

Development of the NGZTD Model
From the above analysis of spatiotemporal characteristics, the surface ZTD exhibited regular seasonal and diurnal variations.Therefore, the development of a ZTD model must consider the delicate periods of the surface ZTD.This equation is expressed as follows: where ZTD i s is the surface ZTD (m), α i 0 is the annual mean ZTD, α i 1 to α i 8 are the period amplitude coefficients, and hod is the hour (UTC) of day.
The surface ZTD of all ERA5 grid points over China from 2014 to 2017 was used to calculate nine coefficients using the least squares method.The new ZTD grid model (NGZTD model), which directly calculates the ZTD, consists of Equation ( 12) and the NGZTD-H model proposed earlier.The method for using the NGZTD model is as follows: First, the user's coordinates are determined to search the ERA5 grid points around the user.After the model coefficients of these grid points have been obtained, as shown in Equation ( 12), the surface ZTD can be calculated by inputting doy and hod using Equation (12).Second, the NGZTD-H model is employed to calculate ZTD at the height as the user after obtaining the Gaussian coefficients b and c.Finally, the user's ZTD is calculated using inverse distance interpolation.The NGZTD model was developed based on high-resolution ERA5 data, employs the high-precision Gaussian function to adjust the ZTD vertically, and considers seasonal and delicate diurnal variations in the ZTD, which ensures the effectiveness of the model.Furthermore, the NGZTD model has the advantages of efficiency and simplicity because measured meteorological parameters are not required, and the model coefficients are fewer in number and are stored in regular grid points.

Accuracy Validation Using GNSS Data
The effectiveness of the NGZTD model for estimating the ZTD was validated using the GNSS-derived ZTD in 2018 over China.Because the NGZTD model estimated the ZTD based on the ERA5 grid points, the GNSS and ERA5 data were unified with the same elevation system.The ZTD with a temporal resolution of 1 h at the GNSS stations was determined using both the NGZTD and GPT3 models, and the results were compared with the GNSS-derived ZTD.The accuracy statistics for the two models are presented in Table 4 and Figure 10.Equation ( 12), the surface ZTD can be calculated by inputting  and ℎ using Equation (12).Second, the NGZTD-H model is employed to calculate ZTD at the height as the user after obtaining the Gaussian coefficients b and c.Finally, the user's ZTD is calculated using inverse distance interpolation.The NGZTD model was developed based on highresolution ERA5 data, employs the high-precision Gaussian function to adjust the ZTD vertically, and considers seasonal and delicate diurnal variations in the ZTD, which ensures the effectiveness of the model.Furthermore, the NGZTD model has the advantages of efficiency and simplicity because measured meteorological parameters are not required, and the model coefficients are fewer in number and are stored in regular grid points.

Accuracy Validation Using GNSS Data
The effectiveness of the NGZTD model for estimating the ZTD was validated using the GNSS-derived ZTD in 2018 over China.Because the NGZTD model estimated the ZTD based on the ERA5 grid points, the GNSS and ERA5 data were unified with the same elevation system.The ZTD with a temporal resolution of 1 h at the GNSS stations was determined using both the NGZTD and GPT3 models, and the results were compared with the GNSS-derived ZTD.The accuracy statistics for the two models are presented in Table 4 and Figure 10.As shown in Table 4, the mean bias of the GPT3 model was −0.38 cm and 0.16 cm.The ZTD calculated by the NGZTD model was generally smaller than the GNSS-derived ZTD but generally larger than the reference value for the GPT3 model.The bias of the NGZTD model was −0.38 cm and ranged from −3.98 to 1.06 cm, and the negative bias was larger, which indicated that some of the ZTD calculated by the NGZTD model were much smaller than the GNSS-derived ZTD.The NGZTD and GPT3 models had a similar RMSE range from 1 cm to 8 cm.The mean RMSE of the two models were 3.38 and 4.33 cm, respectively.As shown in Table 4, the mean bias of the GPT3 model was −0.38 cm and 0.16 cm.The ZTD calculated by the NGZTD model was generally smaller than the GNSS-derived ZTD but generally larger than the reference value for the GPT3 model.The bias of the NGZTD model was −0.38 cm and ranged from −3.98 to 1.06 cm, and the negative bias was larger, which indicated that some of the ZTD calculated by the NGZTD model were much smaller than the GNSS-derived ZTD.The NGZTD and GPT3 models had a similar RMSE range from 1 cm to 8 cm.The mean RMSE of the two models were 3.38 and 4.33 cm, respectively.Compared with the GPT3 model, the NGZTD model achieved an improvement of 0.95 cm (22%).
As shown in Figure 10, the bias of the GPT3 model was positive in the western and southwestern regions and negative in the eastern region, indicating that the performance of this model was unstable.The NGZTD model performed better than the other regions in northern China.Especially in western and southwestern China with undulating terrain, the NGZTD model exhibited better stability, indicating that the NGZTD model is more adaptable to the influence of complex geographical conditions.The NGZTD model exhibited a stable performance with an RMSE of approximately 2 cm.However, the GPT3 model exceeded 2 cm in most regions of China and reached 6 cm in the eastern region.The two models do not perform well in the eastern region, which may be affected by atmospheric differences between the sea and land.
To investigate the season adaptation of the NGZTD model, two representative GNSS stations (Kuqa and Delingha) were selected as validation stations.The diurnal mean bias and RMSE results are shown in Figure 11.The GPT3 model exhibited positive bias at both GNSS stations Kuqa and Delingha, whereas the NGZTD model exhibited smaller and more stable bias during the entire year.The NGZTD model performed a better result with an RMSE below 4 cm, indicating that the new ZTD model can cope with the influences of season.
Compared with the GPT3 model, the NGZTD model achieved an improvement of 0.95 cm (22%).As shown in Figure 10, the bias of the GPT3 model was positive in the western and southwestern regions and negative in the eastern region, indicating that the performance of this model was unstable.The NGZTD model performed better than the other regions in northern China.Especially in western and southwestern China with undulating terrain, the NGZTD model exhibited better stability, indicating that the NGZTD model is more adaptable to the influence of complex geographical conditions.The NGZTD model exhibited a stable performance with an RMSE of approximately 2 cm.However, the GPT3 model exceeded 2 cm in most regions of China and reached 6 cm in the eastern region.The two models do not perform well in the eastern region, which may be affected by atmospheric differences between the sea and land.
To investigate the season adaptation of the NGZTD model, two representative GNSS stations (Kuqa and Delingha) were selected as validation stations.The diurnal mean bias and RMSE results are shown in Figure 11.The GPT3 model exhibited positive bias at both GNSS stations Kuqa and Delingha, whereas the NGZTD model exhibited smaller and more stable bias during the entire year.The NGZTD model performed a better result with an RMSE below 4 cm, indicating that the new ZTD model can cope with the influences of season.ZTD is affected by the season as well as delicate diurnal variations.Two GNSS stations, Kuqa and Guilin, were selected to investigate the effectiveness of the NGZTD model to capture diurnal variations in the ZTD.The hourly ZTD over a five-day period in 2018 was estimated by both the NGZTD and GPT3 models, and their results were compared with the GNSS-derived ZTD and are shown in Figure 12.At the Kuqa GNSS station, the GNSS-derived ZTD ranged from 2.08 to 2.09 cm, exhibiting a diurnal period during those five days.The ZTD estimated by the NGZTD model also exhibited a diurnal period and ZTD is affected by the season as well as delicate diurnal variations.Two GNSS stations, Kuqa and Guilin, were selected to investigate the effectiveness of the NGZTD model to capture diurnal variations in the ZTD.The hourly ZTD over a five-day period in 2018 was estimated by both the NGZTD and GPT3 models, and their results were compared with the GNSS-derived ZTD and are shown in Figure 12.At the Kuqa GNSS station, the GNSS-derived ZTD ranged from 2.08 to 2.09 cm, exhibiting a diurnal period during those five days.The ZTD estimated by the NGZTD model also exhibited a diurnal period and agreed well with the reference value, indicating that the NGZTD model can capture the delicate diurnal variation of the ZTD.However, the ZTD determined using the GPT3 model was the same for one day and was larger than the reference value.At the Guilin GNSS station, although the ZTD calculated by the GPT3 model was close to the GNSS-derived agreed well with the reference value, indicating that the NGZTD model can capture the delicate diurnal variation of the ZTD.However, the ZTD determined using the GPT3 model was the same for one day and was larger than the reference value.At the Guilin GNSS station, although the ZTD calculated by the GPT3 model was close to the GNSSderived ZTD, this model was not able to capture the diurnal variation in the ZTD.Therefore, the proposed model can capture both seasonal and diurnal variations in the ZTD.

Accuracy Validation Using Radiosonde Data
Layered ZTD at radiosonde stations in 2018 over China were considered as reference data to validate the accuracy of the NGZTD model.To investigate the contribution of considering the seasonal variations of Gaussian coefficients b and c to the accuracy of estimating ZTD, the proposed model, which considered the seasonal variations of Gaussian coefficients and took the constants b and c as −7000 and 30,000 m, respectively, was used to calculate the layered ZTD at radiosonde stations at 00:00 and 12:00 and compared the results with those of the GPT3 model.
As shown in Table 5, the mean bias of the NGZTD model and the model with constant Gaussian coefficients were −1.53 and −4.63 cm, respectively, and the ZTD estimated by the NGZTD model with constant Gaussian coefficients showed an evident difference from the reference value, which indicated that seasonal variations in Gaussian coefficients need to be considered.The GPT3 model exhibited a more stable performance with a mean bias of 0.32 cm.The mean RMSE of the NGZTD model, the model with constant Gaussian coefficients, and the GPT3 model were 3.62, 5.91, and 5.25 cm, respectively.The NGZTD model achieved the highest accuracy, with an improvement of 2.29 (39%) and 1.63 cm (31%) compared with the other models.Therefore, these results further indicated that considering seasonal variations in the Gaussian coefficients contributed to the effectiveness of the NGZTD model.

Accuracy Validation Using Radiosonde Data
Layered ZTD at radiosonde stations in 2018 over China were considered as reference data to validate the accuracy of the NGZTD model.To investigate the contribution of considering the seasonal variations of Gaussian coefficients b and c to the accuracy of estimating ZTD, the proposed model, which considered the seasonal variations of Gaussian coefficients and took the constants b and c as −7000 and 30,000 m, respectively, was used to calculate the layered ZTD at radiosonde stations at 00:00 and 12:00 and compared the results with those of the GPT3 model.
As shown in Table 5, the mean bias of the NGZTD model and the model with constant Gaussian coefficients were −1.53 and −4.63 cm, respectively, and the ZTD estimated by the NGZTD model with constant Gaussian coefficients showed an evident difference from the reference value, which indicated that seasonal variations in Gaussian coefficients need to be considered.The GPT3 model exhibited a more stable performance with a mean bias of 0.32 cm.The mean RMSE of the NGZTD model, the model with constant Gaussian coefficients, and the GPT3 model were 3.62, 5.91, and 5.25 cm, respectively.The NGZTD model achieved the highest accuracy, with an improvement of 2.29 (39%) and 1.63 cm (31%) compared with the other models.Therefore, these results further indicated that considering seasonal variations in the Gaussian coefficients contributed to the effectiveness of the NGZTD model.Figure 13 shows that the NGZTD model exhibited a negative bias at most stations, whereas the GPT3 model exhibited a positive bias.The RMSE of the GPT3 model exceeded 5 cm in most regions.In particular, in regions with high altitudes and undulating terrain, the NGZTD model is more advantageous.In addition, Figure 14 shows the accuracy percentage for the NGZTD and GPT3 models.The NGZTD model had the largest proportion of RMSE in the interval of 2-3 cm, accounting for 30% of the RMSE, whereas the largest proportion of RMSE for the GPT3 model was in the interval of 5-6 cm, accounting for 69% of the RMSE.In the interval in which the RMSE was <2 cm, the NGZTD model accounted for 2%, whereas the GPT3 model accounted for 0%.According to the above results of the accuracy validation, the NGZTD model exhibited stable and high accuracy in terms of capturing seasonal and diurnal variations, adaptability to terrain, and adaptability to different elevations.The new ZTD model can directly provide reliable and high-precision ZTD for users in China.
Figure 13 shows that the NGZTD model exhibited a negative bias at most stations, whereas the GPT3 model exhibited a positive bias.The RMSE of the GPT3 model exceeded 5 cm in most regions.In particular, in regions with high altitudes and undulating terrain, the NGZTD model is more advantageous.In addition, Figure 14 shows the accuracy percentage for the NGZTD and GPT3 models.The NGZTD model had the largest proportion of RMSE in the interval of 2-3 cm, accounting for 30% of the RMSE, whereas the largest proportion of RMSE for the GPT3 model was in the interval of 5-6 cm, accounting for 69% of the RMSE.In the interval in which the RMSE was <2 cm, the NGZTD model accounted for 2%, whereas the GPT3 model accounted for 0%.According to the above results of the accuracy validation, the NGZTD model exhibited stable and high accuracy in terms of capturing seasonal and diurnal variations, adaptability to terrain, and adaptability to different elevations.The new ZTD model can directly provide reliable and high-precision ZTD for users in China.Figure 13 shows that the NGZTD model exhibited a negative bias at most stations, whereas the GPT3 model exhibited a positive bias.The RMSE of the GPT3 model exceeded 5 cm in most regions.In particular, in regions with high altitudes and undulating terrain, the NGZTD model is more advantageous.In addition, Figure 14 shows the accuracy percentage for the NGZTD and GPT3 models.The NGZTD model had the largest proportion of RMSE in the interval of 2-3 cm, accounting for 30% of the RMSE, whereas the largest proportion of RMSE for the GPT3 model was in the interval of 5-6 cm, accounting for 69% of the RMSE.In the interval in which the RMSE was <2 cm, the NGZTD model accounted for 2%, whereas the GPT3 model accounted for 0%.According to the above results of the accuracy validation, the NGZTD model exhibited stable and high accuracy in terms of capturing seasonal and diurnal variations, adaptability to terrain, and adaptability to different elevations.The new ZTD model can directly provide reliable and high-precision ZTD for users in China.

Discussion
As shown in Tables 2, 4 and 5, the mean RMSE of the vertical interpolation for the NGZTD-H model was 1.70 cm when the ERA5 data were used as the reference value, whereas the mean RMSE of the NGZTD model for the direct calculation of the ZTD exceeded 3 cm when both the GNSS data and radiosonde data were used as the reference value.Unlike the vertical interpolation validation for the NGZTD-H model in Section 3.1.1,the starting ZTD for the NGZTD model in Sections 3.2.1 and 3.2.2. was obtained by a model considering the seasonal and diurnal variations of ZTD, whereas the starting ZTD in Section 3.1.1.was obtained by the integral method.In addition, the NGZTD model required horizontal interpolation to calculate the ZTD at GNSS and radiosonde stations in Sections 3.2.1 and 3.2.2.Based on the above analysis, the errors of the NGZTD model may be attributed to the starting ZTD and horizontal interpolation.
To validate the above error sources, the GNSS-derived ZTD in 2018 over China was used to validate the spatial interpolation accuracy for the NGZTD-H model.First, the Gaussian coefficients b and c were determined.Second, the ZTD on the ERA5 grid points was vertically adjusted to the GNSS station height using the NGZTD-H model.It should be noted that the starting ZTD on the ERA5 grid points was obtained by the integral method.Finally, the adjusted ZTD was interpolated horizontally to the GNSS station employing an inverse distance interpolation method.A comparison of the interpolated results with the GNSS-derived ZTD validated the effectiveness of the NGZTD-H model.Because the elevations of GNSS stations are geodesic heights and the elevations of ERA5 grid points are geopotential heights, it is necessary to unify the elevations of these two products before vertically adjusting the ZTD.In this study, we employed the EGM2008 model [45] to convert elevations of the ERA5 data to geodesic elevations.
The spatial interpolation accuracy of the NGZTD-H model was compared with that of the GPT3 model.The mean bias of the GPT3 model was −1.14 cm (Table 6).This indicated the GPT3 model had certain systematic errors.The bias intervals of the two models were similar; both ranged from −3 to 1 cm.The mean RMSE of the NGZTD-H and GPT3 models were 1.48 cm and 1.79 cm.Compared with the GPT3 model, the accuracy of the NGZTD-H model improved by 0.31 cm (17%).Therefore, the NGZTD-H model has higher accuracy for the spatial interpolation of the ZTD. Figure 15 shows the stable results for NGZTD-H model at each GNSS station in China.Its bias was approximately 0 cm in most regions and reached −3 cm at individual stations in Southwest China, which indicated that the model struggles to express the vertical variation of ZTD in these regions.However, the result of the GPT3 model shows a large negative bias in most regions.The RMSE of the NGZTD-H model was smaller at high latitudes.However, the RMSE of the GPT3 model was greater than 2 cm at most GNSS stations.Therefore, compared with the GPT3 model, the NGZTD-H model has better stability in spatial interpolation.
The mean RMSE of spatial interpolation for the NGZTD-H model was 1.48 cm, which was apparently smaller than that of the NGZTD model in Sections 3.2.1 and 3.2.2.The validation of spatial interpolation accuracy for the NGZTD-H model used ZTD obtained by the integral method as the starting ZTD and required horizontal interpolation.This indicated that the errors of the NGZTD model mainly derive from the starting ZTD rather than from horizontal interpolation.polation in Section 3.1.1.This further indicated that the errors caused by horizontal interpolation were relatively small.According to the above analysis, the accuracy of ZTD vertical adjustment was higher in flat-terrain regions and lower in western China with undulating terrain.The GNSS stations are fewer in western China, whereas the distribution of ERA5 grid points is uniform, which resulted in a higher number of GNSS stations with smaller errors.Therefore, the accuracy of the spatial interpolation for the NGZTD-H model was higher than that of the vertical interpolation in Section 3.1.1.

Conclusions
To overcome the limitations of existing ZTD models that lack an appropriate vertical adjustment function and are unsuitable for use in China in light of its complex climate and greatly undulating terrain, this study proposed an NGZTD model considering time-varying vertical adjustment based on the Gaussian function and diurnal variation.First, a ZTD vertical adjustment model, NGZTD-H, was developed that can adjust the ZTD from the starting height to the target height.Then, an NGZTD model that considers the seasonal and delicate diurnal variation in the ZTD was developed to estimate the ZTD directly.The vertical interpolation accuracy of the NGZTD-H model was validated, and the model exhibited stability for different geographic regions, pressure layers, latitudes, seasons, and elevations.The NGZTD model was validated using GNSS-derived ZTD and layered ZTD at radiosonde stations, and the results show that the NGZTD model performed better than the GPT3 model in time and space dimensions.In particular, in terms of capturing diurnal variation, the NGZTD model can effectively detect diurnal variation in the ZTD.Because the NGZTD model adopted the Gaussian function to adjust the ZTD vertically, the new As shown in Tables 2 and 6, the spatial interpolation accuracy (the mean RMSE was 1.48 cm) of the NGZTD-H model is higher than the vertical interpolation accuracy (the mean RMSE was 1.70 cm).The validation of spatial interpolation accuracy for the NGZTD-H model required horizontal interpolation, whereas this is not required for vertical interpolation in Section 3.1.1.This further indicated that the errors caused by horizontal interpolation were relatively small.According to the above analysis, the accuracy of ZTD vertical adjustment was higher in flat-terrain regions and lower in western China with undulating terrain.The GNSS stations are fewer in western China, whereas the distribution of ERA5 grid points is uniform, which resulted in a higher number of GNSS stations with smaller errors.Therefore, the accuracy of the spatial interpolation for the NGZTD-H model was higher than that of the vertical interpolation in Section 3.1.1.

Conclusions
To overcome the limitations of existing ZTD models that lack an appropriate vertical adjustment function and are unsuitable for use in China in light of its complex climate and greatly undulating terrain, this study proposed an NGZTD model considering time-varying vertical adjustment based on the Gaussian function and diurnal variation.First, a ZTD vertical adjustment model, NGZTD-H, was developed that can adjust the ZTD from the starting height to the target height.Then, an NGZTD model that considers the seasonal and delicate diurnal variation in the ZTD was developed to estimate the ZTD directly.The vertical interpolation accuracy of the NGZTD-H model was validated, and the model exhibited stability for different geographic regions, pressure layers, latitudes, seasons, and elevations.The NGZTD model was validated using GNSS-derived ZTD and layered ZTD at radiosonde stations, and the results show that the NGZTD model performed better than the GPT3 model in time and space dimensions.In particular, in terms of capturing diurnal variation, the NGZTD model can effectively detect diurnal variation in the ZTD.Because the NGZTD model adopted the Gaussian function to adjust the ZTD vertically, the new model exhibited a more stable performance than the GPT3 model in western and southwestern China, where the terrain is undulating.Seasonal variations in Gaussian Section 4 discusses the main error sources of the NGZTD model.Section 5 contains the Conclusions.The research framework is shown in Figure 1.The proposed model compensates for the limitations of existing ZTD models in China.The new model is expected to be applied to precise point positioning (PPP) and InSAR atmospheric correction to improve their monitoring accuracy.Remote Sens. 2024, 16, 2023 3 of 21 a segmented global ZTD vertical profile model (GZTD-P) with different spatial resolutions that considers the time-varying characteristics of the ZTD vertical variation factor to address the limitations of a single function in expressing the ZTD vertical profile, which shows better performance compared to the GPT3 model.However, the GZTD-P model can only vertically adjust ZTD from the starting height to the target height and cannot directly calculate ZTD at the target position, which limits its application.Sun et al. [32] proposed a global ZTD model, the GZTDS model, which considers delicate periodic variations by adopting a nonlinear function.The GSTDS model was developed as a newgeneration model that follows the GZTDS model, and the performance of the new model deteriorated as the zenith angle increased.Hu and Yao [33] adopted the Gaussian function to fit the vertical ZTD and then established a ZTD vertical profile model that considers seasonal variation, which shows good performance in both time and space.The model achieved good results on a global scale.The model was developed using the monthly mean ZTD with horizontal resolutions of 5° × 5° provided by ERA-Interim.Its adaptability needs to be further improved in regions with greatly undulating terrain and complex climates.Zhao et al. [34] proposed a high-precision ZTD model that considers the height effect on ZTD after analyzing the relationship between the ZTD periodic residual term and the height of the GNSS station at different seasons.Although the aforementioned models have demonstrated their respective advantages, it is necessary to conduct further research on the more delicate vertical and temporal variations in the ZTD over China.

Figure 1 .
Figure 1.The research framework.Figure 1.The research framework.

Figure 1 .
Figure 1.The research framework.Figure 1.The research framework.

Figure 2 .
Figure 2. Distributions of the annual mean value and period amplitudes for Gaussian coefficients b and c.(a) The annual mean value of b.(b) The annual period amplitude of b.(c) The semi-annual period amplitude of b.(d) The annual mean value of c. (e) The annual period amplitude of c. (f) The semi-annual period amplitude of c.

Figure 2 .
Figure 2. Distributions of the annual mean value and period amplitudes for Gaussian coefficients b and c.(a) The annual mean value of b.(b) The annual period amplitude of b.(c) The semi-annual period amplitude of b.(d) The annual mean value of c. (e) The annual period amplitude of c. (f) The semi-annual period amplitude of c. 2.4.Analysis of the Surface ZTD Analyzing the delicate temporal variations of ZTD contributes to the provision of important references for optimizing the ZTD model.The surface ZTD has significant annual and semi-annual periods [34].Two representative grid points in China were selected to determine the diurnal and semi-diurnal periods by adopting the fast Fourier transform.As shown in Figure3, the surface ZTD exhibited significant diurnal and semi-diurnal period variations.The diurnal variation trends of the surface ZTD for the two grid points were different.For the first grid point, the surface ZTD increased and then decreased throughout the day, whereas the second grid point showed the opposite result.The diurnal period variation characteristics of the first grid point were significantly stronger than those of the semi-diurnal period, whereas the diurnal and semi-diurnal period variation characteristics of the second grid point were comparable.Therefore, annual, semi-annual, diurnal, and semi-diurnal periods must be considered when developing a ZTD grid model.As shown in Figure4, the annual mean surface ZTD was larger on the Tibetan Plateau and smaller elsewhere, which may be due to altitudinal differences.The annual period amplitude of the surface ZTD showed a decreasing trend from southeast to northwest, which may be because the northwest is far from the ocean and the southeast is close to the ocean, resulting in a large difference in the type of climate.The semi-annual period amplitude of the ZTD did not exhibit an apparent variation law.The distributions of the diurnal and semi-diurnal period amplitudes were similar.Therefore, these results further demonstrate that the delicate period of the ZTD must be considered when developing a ZTD grid model.

Figure 4 .
Figure 4. Distribution of annual mean value and period amplitudes for surface ZTD.(a) The annual mean value.(b) The annual period amplitude.(c) The semi-annual period amplitude.(d) The diurnal period amplitude.(e) The semi-diurnal period amplitude.

Figure 4 .
Figure 4. Distribution of annual mean value and period amplitudes for surface ZTD.(a) The annual mean value.(b) The annual period amplitude.(c) The semi-annual period amplitude.(d) The diurnal period amplitude.(e) The semi-diurnal period amplitude.

Figure 4 .
Figure 4. Distribution of annual mean value and period amplitudes for surface ZTD.(a) The annual mean value.(b) The annual period amplitude.(c) The semi-annual period amplitude.(d) The diurnal period amplitude.(e) The semi-diurnal period amplitude.

Figure 5 .
Figure 5. Distribution of vertical interpolation accuracy for NGZTD-H model and GPT3 model using ERA5 profile ZTD in 2018.(a) The bias of NGZTD-H.(b) The RMSE of NGZTD-H.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 5 .
Figure 5. Distribution of vertical interpolation accuracy for NGZTD-H model and GPT3 model using ERA5 profile ZTD in 2018.(a) The bias of NGZTD-H.(b) The RMSE of NGZTD-H.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 6 .
Figure 6.Distribution of vertical interpolation accuracy for NGZTD-H and GPT3 models in the selected pressure layers and latitude bands using ERA5 profile ZTD in 2018.(a) The bias of pressure layers.(b) The RMSE of pressure layers.(c) The bias of latitude bands.(d) The RMSE of latitude bands.

Figure 6 .
Figure 6.Distribution of vertical interpolation accuracy for NGZTD-H and GPT3 models in the selected pressure layers and latitude bands using ERA5 profile ZTD in 2018.(a) The bias of pressure layers.(b) The RMSE of pressure layers.(c) The bias of latitude bands.(d) The RMSE of latitude bands.

Figure 7 .
Figure 7. Distribution of vertical interpolation accuracy for the NGZTD-H and GPT3 models using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias of NGZTD-H.(b) The RMSE of NGZTD-H.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 7 .
Figure 7. Distribution of vertical interpolation accuracy for the NGZTD-H and GPT3 models using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias of NGZTD-H.(b) The RMSE of NGZTD-H.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 8 .
Figure 8. Distribution of vertical interpolation accuracy for NGZTD-H model and GPT3 model in different seasons using the ZTD-layered profiles at radiosonde stations in 2018.(a) Hailar.(b) Hangzhou.

Figure 9 .
Figure 9. Distribution of vertical interpolation accuracy for the NGZTD-H and GPT3 models at different heights using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias.(b) The RMSEs.

Figure 8 .
Figure 8. Distribution of vertical interpolation accuracy for NGZTD-H model and GPT3 model in different seasons using the ZTD-layered profiles at radiosonde stations in 2018.(a) Hailar.(b) Hangzhou.

Figure 8 .
Figure 8. Distribution of vertical interpolation accuracy for NGZTD-H model and GPT3 model in different seasons using the ZTD-layered profiles at radiosonde stations in 2018.(a) Hailar.(b) Hangzhou.

Figure 9 .
Figure 9. Distribution of vertical interpolation accuracy for the NGZTD-H and GPT3 models at different heights using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias.(b) The RMSEs.

Figure 9 .
Figure 9. Distribution of vertical interpolation accuracy for the NGZTD-H and GPT3 models at different heights using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias.(b) The RMSEs.

Figure 10 .
Figure 10.Distribution of accuracy for NGZTD and GPT3 models using the GNSS-derived ZTD at GNSS stations in 2018.(a) The bias of NGZTD.(b) The RMSE of NGZTD.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 10 .
Figure 10.Distribution of accuracy for NGZTD and GPT3 models using the GNSS-derived ZTD at GNSS stations in 2018.(a) The bias of NGZTD.(b) The RMSE of NGZTD.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 13 .
Figure 13.Distribution of accuracy for the NGZTD model, the model with constant Gaussian coefficients, and the GPT3 model using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias of NGZTD.(b) The RMSE of NGZTD.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 14 .
Figure 14.The percentage results of the RMSE for the NGZTD and GPT3 models using the ZTDlayered profiles at radiosonde stations in 2018.(a) NGZTD model.(b) GPT3 model.

Figure 13 .
Figure 13.Distribution of accuracy for the NGZTD model, the model with constant Gaussian coefficients, and the GPT3 model using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias of NGZTD.(b) The RMSE of NGZTD.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 13 .
Figure 13.Distribution of accuracy for the NGZTD model, the model with constant Gaussian coefficients, and the GPT3 model using the ZTD-layered profiles at radiosonde stations in 2018.(a) The bias of NGZTD.(b) The RMSE of NGZTD.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 14 .
Figure 14.The percentage results of the RMSE for the NGZTD and GPT3 models using the ZTDlayered profiles at radiosonde stations in 2018.(a) NGZTD model.(b) GPT3 model.

Figure 14 .
Figure 14.The percentage results of the RMSE for the NGZTD and GPT3 models using the ZTDlayered profiles at radiosonde stations in 2018.(a) NGZTD model.(b) GPT3 model.

Figure 15 .
Figure 15.Distribution of spatial interpolation accuracy for NGZTD-H and GPT3 models using the GNSS-derived ZTD at GNSS stations in 2018.(a) The bias of NGZTD-H.(b) The RMSE of NGZTD-H.(c) The bias of GPT3.(d) The RMSE of GPT3.

Figure 15 .
Figure 15.Distribution of spatial interpolation accuracy for NGZTD-H and GPT3 models using the GNSS-derived ZTD at GNSS stations in 2018.(a) The bias of NGZTD-H.(b) The RMSE of NGZTD-H.(c) The bias of GPT3.(d) The RMSE of GPT3.

Table 1 .
The bias and RMSE of ZTD fitting curves using the Gaussian function (unit: mm).

Table 2 .
Statistics of vertical interpolation accuracy for NGZTD-H model and GPT3 model using ERA5 profile ZTD in 2018 (unit: cm).

Table 3 .
Statistics of vertical interpolation accuracy for the NGZTD-H and GPT3 models using ZTDlayered profiles at radiosonde stations in 2018 (unit: cm).

Table 3 .
Statistics of vertical interpolation accuracy for the NGZTD-H and GPT3 models using ZTD-layered profiles at radiosonde stations in 2018 (unit: cm).

Table 4 .
Statistics of accuracy for NGZTD and GPT3 models using GNSS-derived ZTD in 2018 (unit: cm).

Table 4 .
Statistics of accuracy for NGZTD and GPT3 models using GNSS-derived ZTD in 2018 (unit: cm).

Table 5 .
Statistics of accuracy for the NGZTD model, the model with constant Gaussian coefficients, and the GPT3 model using the ZTD-layered profiles at radiosonde stations in 2018 (unit: cm).

Table 5 .
Statistics of accuracy for the NGZTD model, the model with constant Gaussian coefficients, and the GPT3 model using the ZTD-layered profiles at radiosonde stations in 2018 (unit: cm).

Table 6 .
Statistics of spatial interpolation accuracy for NGZTD-H and GPT3 models using GNSSderived ZTD in 2018 (unit: cm).