Analysis of Systematic Biases in Tropospheric Hydrostatic Delay Models and Construction of Correction Model

.


Introduction
Propagation delay caused by troposphere is inevitable in space observing technologies which employ electromagnetic waves as the primary means of detection, such as Global Navigation Satellite System (GNSS), 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 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 GNSS 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 furtherly exert 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.
However, hydrostatic equilibrium will be broken if vertical wind acceleration occurs, which shall probably be 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 locations or 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 high-altitude angle, the estimates of technologies like GNSS can hardly be affected by this kind of bias.That's 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 on to the vertical direction by the rule of thumb (Boehm et al., 2006), which furtherly affects the accuracy of geodetic solutions.
In addition, biases of ZHD will be transferred to ZWD, which is often pursued by ZTD minus ZHD.Since ZWD is closely related to precipitable water vapor (PWV), roughly 0.15-0.25 times of ZWD (Yao et al., 2016), ZHD biases will cause PWV errors indirectly, which is undoubtedly unfavourable 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 analysed the characteristics of model ZHD biases from different aspects, and constructed a set of grid correction model to provide reliable reference for improving various solutions to space geodetic surveying data and obtaining precise meteorological parameters such as PWV.
Deleted: some 65 Deleted: to determining the total delay Commented [ft1]: Referee's comment: Line 37: "To some extend..." -I suggest moving this sentence to the end/middle of the paragraph, starting at line 56.It is unclear here why the ZHD is the key to the total delay determination, and you explained it there.

Commented [ft2R1]:
Thanks for your suggestion.We put an explanation before this sentence to make it more understandable.

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, which 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 st andard atmospheric model COESA 1976 (Minzner, 1977).
Given the low vertical resolution of ERA5 and COESA 1976, Equation ( 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).

T T T T h h h h
The subscript  in Equation ( 1) -(3) refers to the meteorological elements (temperature , total air pressure  and the water vapour pressure   );   ,   and  , represent the air pressure, temperature and water vapour pressure respectively at the interpolation height ℎ  , which locates between the  th and ( + 1) th layer;  , means the virtual temperature and its calculation follows Equation (4) (Hofmeister, 2016)., , 0.378

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 Equation ( 5).
Where   (ℎ), ℎ 0 ,  and  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.   and    are the hydrostatic component refractivity and refractive index respectively, and ℎ  denotes the thickness between subdivide layers.Therein,    can be calculated through Equation ( 6) (Davis et al., 1985). 1 Where  1 is a constant (77.6890 [K/hPa]) selected from the "best average" parameters of Rüeger (2002), while   and   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 ECMWF possesses high-quality meteorological parameter values of profile grids at both global and yearly scales, ZHD calculated using these data were treated as the reference value in subsequent statistical analysis and modelling.

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 Equation ( 7) Where  0 , ℎ 0 and  denote the air pressure (mbar), height (m) and geodetic latitude (°) of the study site, while the constant  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 within 1 mm.Two values of  were 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.

Location-specific analysis
Based on the data prepared in Section 2, we obtained the ZHD model biases by subtracting model ZHDs from reference ZHDs.Twelve sites, divided into 4 groups, were chosen to learn the time-varying characteristics, whose locations and climate types are depicted in Figure 1 and Table 1.Three-year (2019Three-year ( -2021) ) time series of biases of study sites were drawn Commented [ft3]: Referee's comment: Line 88: Put a space between "k" and "represent".
Commented [ft4R3]: Sorry for my negligence, and we revised it.

Deleted: hydrostatic component
Formatted: Font: 10 pt, Not Bold Deleted: Figure 1 Formatted: Font: 10 pt, Not Bold Deleted: Table 1 according to Equation ( 5) and ( 7), and two types of constant  were used (Figure 2).Here we note Equation ( 7) as model SAASD for  = 0.0022768 and model SAASZ for  = 0.0022794.Additionally, the ZHD mean absolute biases (MAB) and standard deviation (STD) of biases for the two traditional models in each study site are displayed in Table 2.The MAB is calculated as , where Bias  means the bias of the i th grid point, and  means the total number of grid points.Deleted: Climate type of each site classification.For more details please refer to Peel et al. (2007) and Beck et al. (2018).ST07 to ST12 are located in the 140 middle of ocean, which have no type attribution in Köppen-Geiger's classification, so a "-" was left in the table above.It can be seen from Figure 2 that if no correction applied, SAASD model would generate about 2.5-millimetre biases; by contrast, those from SAASZ seem to be restricted within 0.5 mm mostly (10 sites in 12 in Table 2).Yet, in fact, there exists unknown high-frequency residual information in both two ZHD models, which might be speculated that it should be related 150 to the vertical pressure gradient in each region.Also, it turns out that SAASZ only increase by ~ 2 mm overall, with the periodic trend still reserved, and the biases' STD of SAASZ keeps almost the same with that of SAASD (Table 2).Besides, from Figure 2, the STDs in marine areas are generally smaller than those on land, and the value of STD seems to be related to the complexity of climatic condition; for example, the sites in Group 1 show little variety in STDs, while those in Group 2 experience a quite different thing.155 Since the difference between the two models is not statistically significant, and SAASD has been used for much longer, we took SAASD for an example in the following part to analyse the characteristics in further.The correlation coefficient matrix of biases of 12 sites was shown in Figure 3, as well as their spectrum analysis diagram (Figure 4) and proportion of annual and semi-annual periodic terms in total energy spectrum (Figure 5).
It can be inferred from Figure 3 that the correlation coefficient gets larger if the two sites are closer.For instance, the 160 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 varies with the location, and the correlation between 165 them seems stronger than those between the sites with the same distance spacing on land.Combined with the spectrum diagram (Figure 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 site's location: Regions with high altitude appear more sensitive to seasons, such as ST04 and ST05, while situations for those in the ocean in south hemisphere are opposite, such as ST10-ST12 (Figure 5).170

Global analysis
To further analyse the global distribution characteristics of ZHD systematic biases, we calculated the 1 arc-degree grid values of ZHD biases during January-December, 2020, and depicted their monthly averaged distributions (Figure 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 averagely negative, while those in the southern hemisphere, especially in highaltitude areas, are positive.Then, with season shifting, the values in the northern hemisphere gradually become positive and those in the southern hemisphere winter go to the opposite way.This annual regression phenomenon can also be found by combining the MABs and STDs of the global monthly averaged biases in Table 3.
2) Combining with Figure 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 America), yet more negative ones in marine areas (such as the North Atlantic and the East Pacific) or low-altitude land (such as northwest of Eurasian continent and northeast of North America).However, this pattern is opposite in Antarctica and areas with high altitude on 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 SAASZ, 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) varies 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, which is latitudinally dependent, this probably leads to the latitudinal distribution.As a counterexample, that pattern is weak in the northern hemisphere.4) Overall, the MAB of SAASD over the whole year is ~0.77mm, not that large, but it changes significantly with season and geographical location, and the difference between the maximum and minimum values exceed 30 mm, especially in the northern hemisphere in summer (Table 3).Hence, such a kind biases should be fully considered in the field of high-precision measurement.
5) Looking back at the SAASZ model, it is only modified in terms of the coefficient of , which is equivalent to an overall enlargement of ~0.11% on the basis of the SAASD.Therefore, the SAASZ model will also experience biases similar to the SAASD, yet their variation with geographical location may differ from that of the latter model.For this reason, a correction is necessary for SAASZ as well.Commented [ft12R11]: Thanks for your question.We chose SAASD just for an example since the difference between the two models is not that large.And we added this explanation at Line 146.240

Grid correction model
Based on the periodic characteristics that the biases perform in the above experiments, we approximated them using a trigonometric function as Equation ( 8), which contains annual and semi-annual periodic terms.
Where  ZHD is the approximated ZHD bias,  0 - 4 represent annual average, annual and semi-annual periodic term coefficients, respectively, and  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 distribution of  0 , √  1 2 +  2 2 and √  3 2 +  4 2 were depicted by Figure 7, where √  1 2 +  2 2 and √  3 2 +  4 2 denote the amplifications of annual and semi-annual periodic terms, respectively.Because of the location-related feature (based on Section 3.1), a grid correction model plus space interpolation shall be applicable to recover the bias at any position at any time during year 202 0.
Five model parameters ( 0 - 4 ) will be needed for each grid point.Once the  ZHD is calculated, it will just serve as the correction of ZHD model.We call it ZHD_crct.
From Figure 7 we have conclusions: 1) In mid & low-latitude areas, annual average of ZHD biases clearly varies with terrain (or air pressure, since pressure generally decrease with height increasing), 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 & low-latitude marine areas are significantly smaller than those in high-latitude land areas.Therefore, if we don't ask for high accuracy, the annual variation of these areas can be directly neglected.However, on the land and in high-latitude areas, especially in Antarctica, Greenland, southern and western Asia, northeast and northwest Africa and other regions, this variation should be well considered, for the annual impact reaches ± 4mm or more.
3) The amplitudes of semi-annual periodic term are insignificant, by and large, yet in areas with 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 semiannual signal is outstanding from surroundings, which remains further study.

Internal coincidence examination
We collected the residual biases of ZHD during 2020 after ZHD_crct being added to check the model's internal coincidence.285 As the global STD in summer in the northern hemisphere reaches the maximum, here the monthly averaged biases of July were painted out for instance by Figure 8a, as well as the distribution chart of their frequency number before and after correction by Figure 9a.Also, the annual averaged MABs and their statistical feature can be referred to Figure 8b and Figure 9b.Comparing Figure 8 and Figure 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 290 and southern hemispheres; and from Figure 9, the MABs and STDs are sharply reduced by ~50% both over July and the whole year.However, biases larger than ± 2mm still exist in Antarctica, Andes Mountains and parts of Asia-Africa continent,  As can be seen from Figure 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 STDs decrease by about 50% in all periods of the year.Figure 10b shows that the improvement in highlatitude regions and those in southern hemisphere with mid & low latitudes is superior; by contrast, the reduction of STDs of 330 biases in the Antarctic region doesn't perform well enough compared with that in others.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 the low temperature (below 240K) areas (also mainly near the south pole), the improvement doesn't outstand others; although the MABs are 335 reduced by half, the STDs are not significantly improved compared with that 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 of ZHD by SAASD can be as large as 10 mm or worse; in the northern hemisphere in summer, in the regions with pressure lower than

External coincidence examination
In order to learn the external coincidence of this grid model, we collected the global biases over the year 2021 be fore and after correction using ZHD_crct, which was fed with coefficients solved according to 2020 bias series.Results of four 350 months' MAB were painted out for comparison (Figure 11), also the statistics of biases (Table 4) and the MAB/STD of biases in different dimensions (Figure 12) were exhibited.365  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 two years seems alike, which means that it is feasible to use historical meteorological data to predict the future bias corrections.

Conclusions
We made a detailed analysis of the biases from the traditional ZHD model.If no correction applied, SAASD, one of the most widely used models, would generate biases at millimetre level.Although the MAB over the whole year is ~0.77mm, biases differ significantly with season and geographical location, which reach ±10 mm or even larger, and the difference between the maximum and minimum values exceed 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 aeras.
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 well removed, and the model still works when used to rectify the biases in the next year, which can be inferred that this model may be possible to forecast corrections in quite a few years.Yet, unfortunately, our model only shows expert in detrending the seasonal impacts, with unknown high-frequency residuals still remained to be studied further, which is likely to be recovered by 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.

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

Figure 2 :
Figure 2: Bias time series of model SAASD and SAASZ.(a)-(l) represent the situations of the 12 sites respectively.Light red dot line and light blue dot line show the mean value of bias series from SAASD and SAASZ respectively.

Figure 5 :
Figure 5: Percentage of annual and semi-annual items in total energy for ZHD biases of model SAASD. 190

Formatted:
Font: (Asian) 宋体, 10 pt, Not Bold, (Asian) Chinese (China) Deleted: Figure 1 Commented [ft11]: Referee comment: Why did you consider only the SAASD for validation?What is the impact of your proposed correction model in SAASZ?

Figure 6 :
Figure 6: Global distribution maps of monthly averaged ZHD systematic biases between SAASD ZHD and integral ZHD.(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 colour bars in pictures.

Figure 7 :
Figure 7: Global distribution of coefficients in periodic function.Changes in the constant term, annual and semi-annual amplitudes are described in (a)-(c), respectively.
Figureand the improvement near equator appears not that good compared with other places at mid & low latitude, keeping similar to the situation of semi-annual amplitudes (Figure8b).

Figure 8 :Figure 9 :FormattedFigure 10 :
Figure 8: MAB over July (a) and the whole year (b) of 2020 after correction by grid model.

Formatted
or with the temperature lower than 260 K, or most land areas, if no correction is applied, the biases of ZHD model will be larger than 1 mm generally.
(a) MABs in January 2021 before correction (b) MABs in January 2021 after correction 355 (c) MABs in April 2021 before correction (d) MABs in April 2021 after correction Formatted: Font: pt, Not Bold, (Asian) Chinese (China) Deleted: Figure Formatted: Font: pt, Not Bold, (Asian) Chinese (China) MABs in July 2021 before correction (f) MABs in July 2021 after correction (g) MABs in October 2021 before correction (h) MABs in October 2021 after correction

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

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

Table 1
Geodetic Coordinates of Study Areas and Corresponding Climate TypesGroup Site abbreviation Latitude (Unit: °) Longitude (Unit: °) Altitude (Unit: m) Formatted: Font: 10 pt, Not Bold Deleted: Figure 2 Formatted: Font: 10 pt, Not Bold Deleted:Table 2 Commented [ft5]: Referee comment: 1.Table 1: Although you cite the climate type in the Note below the table, I don't understand why you put them in the table.If it is important to know the climate type of the Sites, please define them in a short paragraph, probably as a Note under the table.2.Table 1: There is no info for the Climate Type and Climate name for ST07 to ST12.Probably since they are in the middle of the ocean, this information is not available.If yes, please mention it in the text; otherwise, fill in the blanks in the table with the correct info.Commented [ft6R5]: 1. Thanks for your advice, and we made an explanation why the climate type matters in the study.2. You're right about it, and oceans indeed don't own any climate type in Köppen-Geiger's model.We explained that under the table.

Table 3
Statistics of Systematic Biases in Different Months and the Whole Year (Unit: mm)

Table 4
Statistics Bold numbers represent minimum absolute values in row comparisons.
of Systematic Biases Before/After Correction in Different Months and the Whole Year of 2021(Unit: mm)