Analysis of systematic biases in tropospheric hydrostatic delay models and construction of a correction model

In the field of space geodetic techniques, such as global navigation satellite systems (GNSSs), tropospheric zenith hydrostatic delay (ZHD) is chosen as the a priori value of tropospheric total delay. Therefore, the inaccuracy of ZHD will definitely affect parameters like the wet delay and the horizontal gradient of tropospheric delay, accompanied by an indirect influence on the accuracy of geodetic parameters, if not dealt with well at low elevation angles. In fact, however, the most widely used ZHD model currently seems to contain millimeter-level biases from the precise integral method. We explored the bias of traditional ZHD models and analyzed the characteristics in different aspects on a global annual scale. It was found that biases differ significantly with season and geographical location, and the difference between the maximum and minimum values exceeds 30 mm, which should be fully considered in the field of high-precision measurement. Then, we constructed a global grid correction model, which is named ZHD_crct, based on the meteorological data of the year 2020 from the ECMWF (European Centre for Medium-Range Weather Forecasts), and it turned out that the bias of traditional models in the current year could be reduced by ∼ 50 % when the ZHD_crct was added. When we verified the effect of ZHD_crct on the biases in the next year, it worked almost the same as the former year. The mean absolute biases (MABs) of ZHD will be narrowed within ∼ 0.5 mm for most regions, and the SDs (standard deviations) will be within ∼ 0.7 mm. This improvement will be helpful for research on meteorological phenomena as well.

1 Introduction

Propagation delay caused by the troposphere is inevitable in space-observing technologies which employ electromagnetic waves as the primary means of detection, such as global navigation satellite systems (GNSSs), very-long-baseline interferometry (VLBI), satellite laser ranging (SLR), interferometric synthetic aperture radar (InSAR), and other technologies (Bevis et al., 1994; Boehm et al., 2006; Tang et al., 2017; Drożdżewski and Sośnica, 2021; Xiao et al., 2021). Affected by weather conditions, such as air temperature, pressure, and water vapor content, this kind of delay varies strongly over time, geographical location, and transmission path and cannot be eliminated by multi-frequency combined observing, which has become one of the main bottlenecks restricting the accuracy of space geodetic surveying (Mendes, 1999; Alizadeh et al., 2013; Younes and Afify, 2014; Yao et al., 2018; Ma et al., 2022).

To be facilitated, the tropospheric delay is usually divided into hydrostatic and non-hydrostatic (wet) components (Davis et al., 1985). As the main influence, the hydrostatic component accounts for ∼ 90 % of the total delay, which reaches 20 m or even larger at a low range of elevation angle and still remains ∼ 2 m in the zenith direction (Niell, 1996; Chen and Herring, 1997; Penna et al., 2001). In data processing of technologies such as GNSSs and VLBI, zenith hydrostatic delay (ZHD) is usually taken as the a priori value of the total zenith delay and the estimated term as the zenith wet delay (ZWD); finally, both of them are respectively multiplied by the corresponding mapping function and then summed up to obtain a slant-path one (Mendes et al., 2002). Therefore, if the ZHD contains bias, this error will probably be transmitted to the ZWD, which further exerts an influence on the total delay and the final solutions. To this extent, ZHD seems to be the key at the beginning of the processing.

Fortunately, ZHD models, such as those from Hopfield (1971) and Saastamoinen (1972), are capable of approximating the actual situation well. Moreover, with the help of subsequent researchers, uncertainties in model ZHDs have been limited to the sub-millimeter level, especially using closed formulae induced from the precise integral method based on hydrostatic equilibrium condition (Davis et al., 1985; Zhang et al., 2016). Thus, in fields including GNSSs and VLBI, calculated figures from the ZHD models are widely used even as true values (Wang et al., 2005; Tuka and El-Mowafy, 2013; Liu et al., 2017; Feng et al., 2020), which facilitate quite a few advanced models, such as Global Pressure and Temperature (GPT) model series, Vienna Mapping Function (VMF) model series (Boehm et al., 2006, 2007; Böhm et al., 2015; Landskron and Böhm, 2018), and extrapolation models (Li et al., 2018; Hu and Yao, 2019; Li et al., 2020; Wang et al., 2022).

However, hydrostatic equilibrium will be broken if vertical wind acceleration occurs, which is probably influential on the accuracy of traditional ZHD models depicted by closed formulae. In fact, it has been noticed that these models show some certain systematic biases when compared with the precise integral method, which vary with location and time yet are easily neglected (Liu et al., 2000; Chen et al., 2009; Yan et al., 2011; Dai and Zhao, 2013; Zhang et al., 2016; Feng et al., 2020). In the direction of zenith or at a high altitude angle, the estimates of technologies like GNSSs can hardly be affected by this kind of bias. That is because the hydrostatic mapping function is nearly equal to the wet one, and the ZHD and ZWD biases rightly offset each other. Nevertheless, when it comes to low elevation angles, the difference between mapping functions of the two components will lead to slant-path delay errors up to ∼ 10 mm (Fan et al., 2019) and coordinate errors of ∼ 2 mm when mapped onto the vertical direction according to the empirical criterion by Boehm et al. (2006), which further affects the accuracy of geodetic solutions.

In addition, biases of ZHD will be transferred to ZWD, which is often pursued by zenith tropospheric delay (ZTD) minus ZHD. Since ZWD is closely related to precipitable water vapor (PWV), roughly 0.15–0.25 times ZWD (Yao et al., 2016), ZHD biases will cause PWV errors indirectly, which is undoubtedly unfavorable for meteorological-phenomenon studies like the atmospheric water cycle, and forecasting disastrous weather, including rainstorms and typhoons.

From the point of view of situations above, we analyzed the characteristics of model ZHD biases from different aspects and constructed a grid correction model to provide a reliable reference for improving various solutions to space geodetic surveying data and obtain precise meteorological parameters such as PWV.

2 Data preparation

2.1 Profile meteorological data

The profile meteorological data were needed to calculate ZHD reference values based on the integral method, i.e., ray tracing. Data observed by radiosondes or meteorological products provided by international authorities are all suitable to achieve reference values. Our ZHD reference values were obtained mainly based on the multilayer grid meteorological data provided by the fifth-generation ECMWF reanalysis (ERA5) dataset; have been validated well in tropospheric delay calculations; and cover the whole world, even the middle of oceans (Abdelfatah et al., 2015; Bekaert et al., 2015; Graham et al., 2019; Sun et al., 2019; Dogan and Erdogan, 2022).

The ERA5 profile data are mostly valid at about 30–40 km high. As there still exists a small amount of air above, when the valid height of data was exceeded during experiments, they were complemented until 86 km using the American standard atmospheric model COESA 1976 (Minzner, 1977).

Given the low vertical resolution of ERA5 and COESA 1976, Eqs. (1)–(3) were used to interpolate air pressure, temperature, and humidity between adjacent layers (Nafisi et al., 2012), where the interpolating step follows the criterion described by Rocken et al. (2001).

$$\begin{array}{}\text{(1)}& {\displaystyle }{p}_{\mathrm{int}}=p_i\exp \left[-{\displaystyle \frac{(h_{\mathrm{int}}-h_i)g_i}{R_d{T}_{\mathrm{v},i}}}\right]\text{(2)}& {\displaystyle }{T}_{\mathrm{int}}=T_{i+\mathrm{1}}+{\displaystyle \frac{T_i-T_{i+\mathrm{1}}}{h_i-h_{i+\mathrm{1}}}}(h_{\mathrm{int}}-h_{i+\mathrm{1}})\text{(3)}& {\displaystyle }{p}_{\mathrm{w},\mathrm{int}}=\left\{\begin{array}{l}{p}_{\mathrm{w},i}\exp \left[\frac{h_{\mathrm{int}}-h_i}{h_{i+\mathrm{1}}-h_i}\log (p_{\mathrm{w},i+\mathrm{1}}/p_{\mathrm{w},i})\right],\\ \mathrm{if}{p}_{\mathrm{w},i}>\mathrm{0}\\ \frac{p_{\mathrm{w},i+\mathrm{1}}}{h_{i+\mathrm{1}}-h_i}(h_{\mathrm{int}}-h_i),\mathrm{if}{p}_{\mathrm{w},i}=\mathrm{0}\end{array}\right.\end{array}$$

The subscript i in Eqs. (1)–(3) refers to the meteorological elements (temperature T, total air pressure p, and water vapor pressure p_w); p_int, T_int, and p_w,int represent the air pressure, temperature, and water vapor pressure, respectively, at the interpolation height h_int, which is located between the ith and (i+1)th layer; T_v,i means the virtual temperature, and its calculation follows Eq. (4) (Hofmeister, 2016).

$$\begin{array}{}\text{(4)}& T_{\mathrm{v},i}={\displaystyle \frac{T_i{p}_i}{p_i-\mathrm{0.378}{p}_{\mathrm{w},i}}}\end{array}$$

2.2 Calculation of ZHD reference data

According to the principle of ray tracing (Thayer, 1967; Nafisi et al., 2012; Eriksson et al., 2014), ZHD can be accurately achieved by taking the integrals of refractive parameters in the sky above the study site, specifically as seen in Eq. (5).

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{ZHD}}_{\mathrm{reference}}& ={\int }_{h\mathrm{0}}^H\left[n^H\right(h)-\mathrm{1}]\mathrm{d}h\approx \sum _{i=\mathrm{1}}^{k-\mathrm{1}}(n_i^H-\mathrm{1})\mathrm{\Delta }{h}_i\\ & ={\mathrm{10}}^{-\mathrm{6}}\sum _{i=\mathrm{1}}^{k-\mathrm{1}}{N}_i^H\mathrm{\Delta }{h}_i,\end{array}\end{array}$$

where n^H(h), h₀, H, and k represent hydrostatic delay, the refractive index of hydrostatic delay at height h, geodetic height of the site, maximum tropospheric thickness, and the total number of divided tropospheric layers, respectively; $n_i^H$ and $N_i^H$ are the hydrostatic refractivity and refractive index, respectively; and Δh_i denotes the thickness between subdivided layers. Therein, $N_i^H$ can be calculated through Eq. (6) (Davis et al., 1985).

$$\begin{array}{}\text{(6)}& N_i^H=k_{\mathrm{1}}{\displaystyle \frac{p_i}{T_i}},\end{array}$$

where k₁ is a constant (77.6890 K hPa⁻¹) selected from the “best average” parameters of Rüeger (2002), while p_i and T_i represent the air pressure and air temperature between the corresponding height layers, respectively. Given accurate meteorological parameters, the integral results shall naturally be of high reliability, thus serving as reference or actual values in atmospheric-delay-effect analysis (Hobiger et al., 2008; Hofmeister and Böhm, 2017; Osah et al., 2021). Since the ECMWF possesses high-quality meteorological parameter values of profile grids at both global and yearly scales, ZHD values calculated using these data were treated as the reference value in subsequent statistical analysis and modeling.

2.3 Traditional ZHD model

According to findings of pioneers (Saastamoinen, 1972; Davis et al., 1985; Penna et al., 2001; Tregoning and Herring, 2006; Leandro et al., 2008), ZHD values can be determined by closed formula as Eq. (7)

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathrm{ZHD}}_{\mathrm{Empirical}}& (p_{\mathrm{0}},h_{\mathrm{0}},\mathit{\phi })\\ & =C\times {\displaystyle \frac{p_{\mathrm{0}}}{\mathrm{1}-\mathrm{0.0026}\cos \left(\mathrm{2}\mathit{\phi }\right)-\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{7}}{h}_{\mathrm{0}}}}\end{array},\end{array}$$

where p₀, h₀, and φ denote the air pressure (mbar), height (m), and geodetic latitude (^∘) of the study site, while the constant C in most cases is set as 0.0022768, refined by Davis et al. (1985), or it could be set as 0.0022794, modified by Zhang et al. (2016), which is claimed to be able to cut bias down to within 1 mm. Two values of C are both studied in the following part.

Obviously, traditional models are simpler and more practical than the integral method since only the pressure and position parameters of the site are required. In the experiments, we exploited the ECMWF surface data for meteorological inputs of traditional models.

3 Biases' characteristics and the correction model

3.1 Location-specific analysis

Based on the data prepared in Sect. 2, we obtained the ZHD model biases by subtracting model ZHDs from reference ZHDs. Twelve sites, divided into four groups, were chosen to learn the time-varying characteristics, whose locations and climate types are depicted in Fig. 1 and Table 1. Three-year (2019–2021) time series of biases of study sites were drawn according to Eqs. (5) and (7), and two types of constant C were used (Fig. 2). Here we note Eq. (7) as the model SAAS_D for C=0.0022768 and the model SAAS_Z for C=0.0022794. Additionally, the ZHD mean absolute biases (MABs) and standard deviations (SDs) of biases for the two traditional models at each study site are displayed in Table 2. The MAB is calculated as $\frac{\mathrm{1}}{N}{\sum }_{i=\mathrm{1}}^N\left|{\mathrm{Bias}}_i\right|$, where Bias_i means the bias of the ith grid point, and N means the total number of grid points.

Figure 1 The 10 arcmin global land topography. White circles are the locations of 12 sites, and the color bar represents global surface fluctuations. Topography data were abstracted from the ETOPO1 Ice Surface grid model (Amante and Eakins, 2009).

Figure 2 Bias time series of the model SAAS_D and SAAS_Z. Panels (a)–(l) represent the situations of the 12 sites, respectively. The light-red dotted line and light-blue dotted line show the mean value of bias series from SAAS_D and SAAS_Z, respectively.

Table 1 Geodetic coordinates of study areas and corresponding climate types.

Note: climate type determines the general weather conditions in this area; thus, it probably keeps close relation with ZHD. Because of this, the climate type of each site is displayed here for further analysis, which complies with Köppen–Geiger climate classification. For more details please refer to Peel et al. (2007) and Beck et al. (2018). ST07 to ST12 are located in the middle of the ocean and have no type attribution in the Köppen–Geiger classification, so a “–” was left in the table above.

Table 2 Mean absolute bias (MAB)/standard deviation (SD) of SAAS_D and SAAS_Z biases (unit: mm).

Note: bold numbers represent minimum absolute values in row comparisons.

It can be seen from Fig. 2 that if no correction is applied, the SAAS_D model would generate about 2.5 mm biases; by contrast, those from SAAS_Z seem to be restricted to within 0.5 mm mostly (10 sites of 12 in Table 2). Yet, in fact, there exists unknown high-frequency residual information in both ZHD models, and it might be speculated that it should be related to the vertical pressure gradient in each region. Also, it turns out that SAAS_Z only increases by ∼ 2 mm overall, with the periodic trend still reserved, and the biases' SD of SAAS_Z keeps almost the same as that of SAAS_D (Table 2). Besides, from Fig. 2, the SDs in marine areas are generally smaller than those on land, and the value of SD seems to be related to the complexity of climatic condition; for example, the sites in Group 1 show little variety in SDs, while those in Group 2 experience a quite different thing.

Since the difference between the two models is not statistically significant, and SAAS_D has been used for much longer, we took SAAS_D for an example in the following part to analyze the characteristics further. The correlation coefficient matrix of biases of 12 sites was shown in Fig. 3, as well as their spectrum analysis diagram (Fig. 4) and proportion of annual and semi-annual periodic terms in the total energy spectrum (Fig. 5).

Figure 3 Correlation coefficients of bias time series of the 12 sites.

Figure 4 Temporal frequency spectrograms for ZHD biases of model SAAS_D.

Figure 5 Percentage of annual and semi-annual items in total energy for ZHD biases of the model SAAS_D.

It can be inferred from Fig. 3 that the correlation coefficient becomes larger if the two sites are closer. For instance, the coefficients of the sites within the group are apparently larger than those between groups; also, as ST02 is away from the other two in Group 1, the correlation appears weaker. The same phenomenon can be observed from the relationship between ST06 and the other two in Group 2. When we compare the correlation between Group 1 and Group 2, it can be seen that the correlation is also generally larger in areas with the same climate category or at similar elevation. When it comes to the sites in Group 3 and Group 4, it turns out that the biases in marine areas also vary with the location, and the correlation between them seems stronger than those between the sites with the same distance spacing on land. Combined with the spectrum diagram (Fig. 4), it can be seen that the biases for most sites are generally influenced by annual and semi-annual signals, and the energy intensity varies significantly with the site's location: regions with high altitude appear more sensitive to season, such as ST04 and ST05, while situations for those in the ocean in the Southern Hemisphere are opposite, such as ST10–ST12 (Fig. 5).

3.2 Global analysis

To further analyze the global distribution characteristics of ZHD systematic biases, we calculated the 1 arcdeg grid values of ZHD biases during January–December 2020 and depicted their monthly averaged distributions (Fig. 6) with statistics listed in Table 3. Herein, we have some conclusions.

1. It confirms the seasonal characteristics, which is reflected in varying degrees throughout the world. For instance, the bias values in Northern Hemisphere winter are on average negative, while those in the Southern Hemisphere, especially in high-altitude areas, are positive. Then, with season shifting, the values in the Northern Hemisphere gradually become positive, and those in the Southern Hemisphere winter go the opposite way. This annual regression phenomenon can also be found by combining the MABs and SDs of the global monthly averaged biases in Table 3.

2. Combining with Fig. 1, we find that the biases have a different sign in regions with different altitude. In the Northern Hemisphere in winter, there are more positive values in high-altitude areas (such as the Qinghai Tibet Plateau in Asia and the Andes Mountains in South America), yet more negative ones in marine areas (such as the North Atlantic and the eastern Pacific) or low-altitude land (such as the northwest of the Eurasian continent and northeast of North America). However, this pattern is opposite in Antarctica and areas with high altitude in Greenland. The preliminary judgment may be subjected to the fact that those regions are in the polar high-pressure zone all year round, and the vertical pressure gradient force is usually less than the atmospheric gravity; at this time, the hydrostatic equilibrium is broken, and the closed formula deduced from this condition deviates. Looking back at SAAS_Z, if this model is used, biases in areas with high altitude will probably become worse.

3. By and large, biases in the Southern Hemisphere (except Antarctica) vary with latitude even in different seasons. Considering that the ocean area in the Southern Hemisphere is much larger than the land, and the meteorological condition is dominated by oceans, whose physical characteristics are affected by solar energy and closely related to latitude, this probably leads to the latitudinal distribution. As a counterexample, that pattern is weak in the Northern Hemisphere.

4. Overall, the MAB of SAAS_D over the whole year is ∼ 0.77 mm, not that large, but it changes significantly with season and geographical location, and the difference between the maximum and minimum values exceeds 30 mm, especially in the Northern Hemisphere in summer (Table 3). Hence, such kinds of biases should be fully considered in the field of high-precision measurement.

5. Looking back at the SAAS_Z model, it is only modified in terms of the coefficient of C, which is equivalent to an overall enlargement of ∼ 0.11 % on the basis of the SAAS_D. Therefore, the SAAS_Z model will also experience biases similar to the SAAS_D, yet their variation with geographical location may differ from that of the latter model. For this reason, a correction is necessary for SAAS_Z as well.

Figure 6 Global distribution maps of monthly averaged ZHD systematic biases between SAAS_D ZHD and integral ZHD. Panels (a)–(l) represent the monthly averaged biases from January to December 2020. In order to distinguish the variations between different months, we use the same scale color bars in the 12 panels above.

Figure 7 Global distribution of coefficients as a periodic function. Changes in the constant term and annual and semi-annual amplitudes are described in (a)–(c), respectively.

Figure 8 MABs over July (a) and the whole year (b) of 2020 after correction by the grid model.

Figure 9 Frequency distribution of MABs over July (a) and the whole year (b) of 2020 before and after correction.

Table 3 Statistics of systematic biases in different months and the whole year (unit: mm).

Note: bold numbers represent minimum absolute values in column comparisons.

3.3 Grid correction model

Based on the periodic characteristics that the biases perform in the above experiments, we approximated them using a trigonometric function as Eq. (8), which contains annual and semi-annual periodic terms.

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{ZHD}}& =a_{\mathrm{0}}+a_{\mathrm{1}}\cos (\mathrm{2}\mathit{\pi }t/\mathrm{365.25})+a_{\mathrm{2}}\sin (\mathrm{2}\mathit{\pi }t/\mathrm{365.25})\\ & +a_{\mathrm{3}}\cos (\mathrm{4}\mathit{\pi }t/\mathrm{365.25})+a_{\mathrm{4}}\sin (\mathrm{4}\mathit{\pi }t/\mathrm{365.25}),\end{array}\end{array}$$

where δ_ZHD is the approximated ZHD bias; a₀–a₄ represent annual average and annual and semi-annual periodic term coefficients, respectively; and t stands for the day of year (DOY) (Böhm et al., 2015; Landskron and Böhm, 2018). The parameters were estimated using the biases of each grid point on a global scale in 2020, and the distributions of a₀, $\sqrt{a_{\mathrm{1}}^{\mathrm{2}}+a_{\mathrm{2}}^{\mathrm{2}}}$, and $\sqrt{a_{\mathrm{3}}^{\mathrm{2}}+a_{\mathrm{4}}^{\mathrm{2}}}$ are depicted in Fig. 7, where $\sqrt{a_{\mathrm{1}}^{\mathrm{2}}+a_{\mathrm{2}}^{\mathrm{2}}}$ and $\sqrt{a_{\mathrm{3}}^{\mathrm{2}}+a_{\mathrm{4}}^{\mathrm{2}}}$ denote the amplifications of annual and semi-annual periodic terms, respectively. Because of the location-related feature (based on Sect. 3.1), a grid correction model plus space interpolation are applicable to recover the bias at any position at any time during the year 2020. Five model parameters (a₀–a₄) will be needed for each grid point. Once the δ_ZHD is calculated, it will just serve as the correction of the ZHD model. We call it ZHD_crct.

From Fig. 7 we have the following conclusions:

1. In mid- and low-latitude areas, the annual average of ZHD biases clearly varies with terrain (or air pressure, since pressure generally decreases with increasing height), which is maintained between −0.2–0.4 mm in ocean and plain areas and ranges from 1.0 to 10.0 mm in high-altitude areas like plateaus and mountains. Results seem to be reversed in cold high-latitude areas, such as Greenland (Arctic) and some of Antarctica, for annual average decreases to −4 mm or even lower.

2. The annual amplitudes in the mid- and low-latitude marine areas are significantly smaller than those in high-latitude land areas. Therefore, if we do not ask for high accuracy, the annual variation in these areas can be directly neglected. However, on the land and in high-latitude areas, especially in Antarctica, Greenland, southern and western Asia, northeastern and northwestern Africa, and other regions, this variation should be well considered since the annual impact reaches ±4 mm or more.

3. The amplitudes of the semi-annual periodic term are insignificant, by and large, yet in areas with a strong annual signal, the semi-annual effect is still obvious enough, even reaching ∼ 1.6 mm. In addition, in the ocean areas near the Equator, the semi-annual signal sticks out from its surroundings, which requires further study.

Figure 10 The MAB and SD of biases in different dimensions before and after correction over 2020. Panel (a) represents the global MABs and SDs on different DOYs, (b) represents the MABs and SDs at different latitudes, (c) represents the MABs and SDs in regions with different air pressures, (d) represents the MABs and SDs in regions with different temperatures, and (e) represents the MABs and SDs in regions with different heights.

Figure 11 Global MABs before (left ones) and after (right ones) correction by the grid model over January (a–b), April (c–d), July (e–f), and October (g–h) in 2021. Please note that the color bar of each month is different.

4 Validation

4.1 Internal coincidence examination

We collected the residual biases of ZHD during 2020 after ZHD_crct being added to check the model's internal coincidence. As the global SD in summer in the Northern Hemisphere reaches the maximum, here the monthly averaged biases of July are depicted for instance in Fig. 8a, as well as the distribution chart of their frequency number before and after correction in Fig. 9a. Also, the annually averaged MABs and their statistical feature can be referred to in Fig. 8b and Fig. 9b. Comparing Fig. 8 and Fig. 6g, it can be seen that the corrected biases are less than ±0.5 mm in most regions of the world, even in mountainous and plateau regions; positive and negative values are nearly evenly distributed in the Northern Hemisphere and Southern Hemisphere, and from Fig. 9, the MABs and SDs are sharply reduced by ∼ 50 % over both July and the whole year. However, biases larger than ±2 mm still exist in Antarctica, the Andes Mountains, and parts of the continents of Asia and Africa, and the improvement near the Equator appears to be not that good compared with other mid- and low-latitude places, similar to the situation of semi-annual amplitudes (Fig. 8b).

Figure 10 shows the MAB and SD of biases in different dimensions before and after correction, in which (a) represents the global MABs and SDs on different DOYs, (b) represents the MABs and SDs at different latitudes, (c) represents the MABs and SDs in regions with different air pressures, (d) represents the MABs and SDs in regions with different temperatures, and (e) represents the MABs and SDs in regions with different heights.

As can be seen from Fig. 10, the corrected biases have been significantly cut down in different dimensions. Figure 10a shows that the annual and semi-annual trends of biases are basically removed. The biases are weakened the most in winter and summer, and the SDs decrease by about 50 % in all periods of the year. Figure 10b shows that the improvement in high-latitude regions and in the Southern Hemisphere at mid-latitudes and low latitudes is superior; by contrast, the reduction in SDs of biases in the Antarctic region was not significant enough compared with that in other regions. Figure 10c shows that the reduction in the low-pressure (600–750 mbar) areas is weaker than that in other areas; these areas are mainly near the south pole, which is basically consistent with the previous conclusion. At the same time, the improvement in the high-pressure (above 1000 mbar) areas is also insignificant, which needs to be further studied. Figure 10d shows that in low-temperature (below 240 K) areas (also mainly near the south pole), the improvement does not stick out from the others; although the MABs are reduced by half, the SDs are not significantly improved compared with those before correction. Figure 10e shows that the improvement is more obvious in land areas with an altitude below 3000 m or over 5000 m.

In terms of the views above, it can be summarized that if no correction is made near the south pole, the error in ZHD by SAAS_D can be as large as 10 mm or worse; in the Northern Hemisphere in summer, in the regions with pressure lower than 750 mbar, with the temperature lower than 260 K, or most land areas, if no correction is applied, the biases of the ZHD model will be larger than 1 mm generally.

4.2 External coincidence examination

In order to learn the external coincidence of this grid model, we collected the global biases over the year 2021 before and after correction using ZHD_crct, which was fed with coefficients solved according to 2020 bias series. Results of 4 months' MABs are depicted for comparison (Fig. 11), and the statistics of biases (Table 4) and the MAB and SD of biases in different dimensions (Fig. 12) are also exhibited.

Figure 12 The MABs and SDs of biases in different dimensions before and after correction over 2021. Panel (a) represents the global MABs and SDs on different DOYs, (b) represents the MABs and SDs at different latitudes, (c) represents the MABs and SDs in regions with different air pressures, (d) represents the MABs and SDs in regions with different temperatures, and (e) represents the MABs and SDs in regions with different heights.

Table 4 Statistics of systematic biases before/after correction in different months and the whole year of 2021 (unit: mm).

Note: bold numbers represent minimum absolute values in row comparisons.

From the charts above, it is clear that the correction made according to the biases of 2020 still works in 2021, and the improvement over the 2 years seems alike, which means that it is feasible to use historical meteorological data to predict the future bias corrections.

5 Conclusions

We made a detailed analysis of the biases from the traditional ZHD model. If no correction is applied, SAAS_D, one of the most widely used models, would generate biases at the millimeter level. Although the MAB over the whole year is ∼ 0.77 mm, biases differ significantly with season and geographical location and reach ±10 mm or even larger, and the difference between the maximum and minimum values exceeds 30 mm. Therefore, it should be fully considered in the field of high-precision measurement, especially for the regions near the south pole and those with pressure less than 750 mbar, temperature less than 260 K, and most land areas.

In order to cut down the biases, a grid model, ZHD_crct, was constructed exploiting a periodic function based on the meteorological data of the current year from ECMWF. Generally, ZHD biases were reduced by ∼ 50 % after correction, with annual and semi-annual terms thoroughly removed, and the model still works when used to rectify the biases in the next year. It can be inferred that this model may be able to forecast corrections in quite a few years. Yet, unfortunately, our model only shows skill in detrending the seasonal impacts, with unknown high-frequency residuals still requiring further study, and this problem is likely to be conquered when we start our study of the vertical velocity of air.

In addition, nowadays profile meteorological parameters can be provided from various organizations, such as MERRA2 data products from the NASA Global Modeling and Assimilation Office, NCEP reanalysis data products, the Integrated Global Radiosonde Archive (IGRA) from the National Climatic Data Center of America (NCDC), and the University of Wyoming Atmospheric Science Radiosonde Archive. Hence, ZHD bias correction models could be further developed based on different datasets in follow-up studies.