CC1

Please consider using different colour combinations, as red-green combinations are indistinguishable for many of your readers with colour-vision deficiencies, which excludes them from understanding most figures in the current pre-print.

As you would probably agree, figures should also be readable by everyone with colour vision deficiency (CVD; this affects ~5% of the global population). To achieve this basic standard of science figures, CVD-unfit colour combinations like red and green colour palettes have to be avoided and be replaced by accessible colour palettes. Amongst others, the Scientific colour maps (www.fabiocrameri.ch/colourmaps) are a freely available, citable, and validated resource of various colour palettes of all different types ensuring accurate and accessible data representation. Because the Scientific colour maps address known issues, they have been made compatible with most software frequently used to create graphics. Details and step-by-step instructions are given in the included user guide. The background of this tool and others, as well as the importance of addressing the broader issue, are described in Crameri et al. (2020), https://doi.org/10.1038/s41467-020-19160-7.



Finally, to check your figures for accessibility issues, there is open tools available, such as Coblis (https://www.color-blindness.com/coblis-color-blindness-simulator/).



Thanks for your work!

AC:

Dear Dr. Fabio Crameri,

Thank you for your suggestions. We have updated all our plots adopting the colormap provided by Colorbrew2 (https://colorbrewer2.org/), and the new plots are all colorblind safe. The following changes are made in the acknowledge section: “We appreciate the suggestions by Dr. Fabio Crameri and adopted the colorblind safe colormap provided by Colorbrew2 (https://colorbrewer2.org).”

RC1

The study was carried out properly and was well presented in terms of organization and use of language. I have some comments and suggestions in the upload pdf-file that should be clarified in revisions. The manuscript should be accepted for publication after minor revisions.

AC:

Dear reviewer,

thank you for your encouraging comments and helpful suggestions. We have adopted all your comments and revised our manuscript accordingly. In the following, we copied your comments in the PDF-file and provided our point-to-point response. These responses are also given in the PDF-file with comments.

R1C1 (line 20): what is "uncertainty"? what is "accuracy"?  what is uncertainty is quantified? In the empirical model, "formal error" is used for uncertainty later on. Better to use ZTD and its accuracy of the empirical model.

AC1: ‘Uncertainty’ here refers to the model’s formal error, while ‘accuracy’ refers to the agreement of modeled ZTD with external references, such as NWM and GNSS. We have changed to “formal error” in the revised manuscript.

R1C2 (line 22): Better express it as: Using GPS observations from xxx IGS stations processed in kinematic PPP mode.

AC2: Done.

R1C3 (line 23): is this import? providing uncertainty information is the most important issue instead of the form of the information. Or make it clear both the ZTD and the corresponding uncertainty models can be calculated without any meteo. data but epoch time and location.

AC3: As our modeling data is based on the similar Numerical Weather Model input that is used by the GPT series, our formal error information can also be adopted by the GPT series. We believe that for any model, formal error information is necessary. For the current empirical models, e.g., the GPT series, formal error is not provided. Therefore, we emphasize that our model has more broad applications. For empirical models such as GPT series and our model, it is acknowledged that no meteorological information is required and only time and location are necessary, thus we do not repeat this statement.

R1C4 (line 41): whether  Saas. and Hopfield models are empirical model or ?

AC4: Saastamoinen and Hopfield models require the meteorological data as input, such as pressure and temperature. If the meteorological data is provided by an empirical model, e.g., GPT, then Saastamoinen/Hopfield models are used empirically. On the other hand, if the meteorological data come from in situ observations or Numerical Weather Model (NWM), then it is not an empirical model.

R1C5 (line 65): if so, it should not be the accuracy but agreement with the reference used.

AC5: we agree with this comment and have changed it to “agreement” in the revised manuscript.

R1C6 (line 70): better than 1 cm?

AC6: modified.

R1C7 (line 41): zenith or slant path delay?

AC7: in this case we refer to slant delay, which is dependent on both the zenith delay and the mapping function.

R1C8 (line 91): the empirical model is derived based on data in the period of 2009 to 2018. Why data in 2017 and 2018 are used.

AC8: The modeling data covers from 2009 to 2018, so in the evaluation using NWM as reference we select the data of 2019-2021, which has no overlap with modeling period. When compared with GNSS ZTD, we would like to cover both the modeling and the prediction periods, so we select 2017-2020, where 2017-2018 covers the modeling period and 2019-2020 covers the prediction period.

R1C9 (line 110): fit all the ZTDs using the following function

AC9: corrected.

R1C10 (line 119): this is also because of higher altitude?

AC10: yes.

R1C11 (line 143): make it clear that this is the uncertainty information of the empirical model.

AC11: thank you and we’ve modified it.

R1C12 (line 143): this means using Res(t)^2 as accuracy. But It is difficult  to see the relationship with "minimizes the differences between the squared formal error and the squared ZTD residuals."  My suggestion is:

1. model the uncertainty to reflect (the distribution of ) the fitting residuals, for example to minimize the differences

\Sigma _i=1,ndays  \Delta_i ^2 =min

1. assumes \Sigma(t)^2 has the similar expression as ZTD, i.e. Eq(2).
2. the the coefficients are estimated according to the minimizing criteria by LSQ.

Is this reasonable?

AC12: thank you for the insightful suggestions regarding the fitting process. Given the residuals of ZTD model, we have a zero-mean time series randomly distributed. We can either fit the absolute or the squared values of the residuals. We found that when fitting the absolute residuals, the fitting RMS is not comparable to the RMS of the ZTD modeling residuals, and a scaling factor has to be applied. On the other hand, if we fit the squared values of the ZTD modeling residuals, the formal error modeling RMS can better represent the ZTD modeling accuracy.

R1C13 (line 156): it has a special meaning in GNSS, or you mean pseudo-obseravtion?

AC13: yes this is a typo. We have corrected it.

R1C14 (line 169): “Note that the presented amplitudes do not give the peak-to-peak values of the empirically modeled formal error, and only illustrate the strength of the annual and semi-annual signals”: not clear to me

AC14: when fitting both annual and semi-annual periodical terms, we cannot consider any of them as the peak-to-peak values of the time series, as the peak-to-peak value is a combined effect of both.

R1C15 (line 177): is there an explanation of the excellent agreement based on the fitting method?

AC15: as we fit the squared residuals using the function of a constant (R0) with annual signals, by definition the RMS of residual should be consistent with the constant term.

R1C16 (line 195): why? The improvement is too small compared to the accuracy of the empirical model? The estimation of Beta with Eq. (3) is anyway to consider the impact of the height differences among grids.

AC16: yes, the improvement of adopting the higher-order exponential function is only a few mm. For the empirical modeling of ZTD and formal error, the ZTD modeling error is around 3-4 cm. The higher-order exponential function is more significant when the epoch-wise NWM-based delays are used.

R1C17 (line 202): Is this the empirical model to be evaluated? If yes, I would emphasize the data interval (2009-2018?)used again.

AC17: Yes and thank you for your suggestions. We modified the manuscript accordingly.

R1C18 (line 258): “However, the Pole regions (latitude higher than 60°), especially the north one, show larger temporal noises than other regions.”: what are the possible reasons?

AC18: we believe this is due to the fact that in the higher latitudes, the annual and semi-annual variations of the ZTD modeling error is less significant than other regions, as shown in Figure 9.

R1C19 (line 273): why 2017 and 2018 are involved? From Fig.10-12, the evaluation is indeed carried out with data from 2017 to 2020. That means half for data involved in the model establishment and half for data not involved. What is the intension?

AC19: yes and we intentionally select the period of 2017-2020 when compared with GNSS ZTD. As the model is determined using NWM from 2009 to 2018, the evaluation using NWM is performed in the period of 2019-2021, i.e., external period. In the evaluation using GNSS, however, we would like to cover both the modeling and prediction periods, so we selected two years for both.

R1C20 (line 295): the axis for ZTD residuals, RMS, and \sigma is missing?

AC20: The axes of these variables are the same as that for ZTD values. However, we do apply a shift for the residual/RMS/sigma values for better visibility.

R1C21 (line 335): of which year? Each day there are 6 sessions and each year 120 days are selected, so each year 720 sessions for one station.

AC21: In the year of 2022. We agree that each station has 720 sessions per year. We have revised the manuscript.

R1C22 (line 344): This is not the same as JPL used (line 270)

AC22: No it is not. Based on our data processing experience, we believe that 5mm/sqrt(h) is a proper value. The temporal constraints also depend on other setup of each software. The value used by JPL (5x10^-8 mm/sqrt(s)) is 3 mm/sqrt(h), which is comparable to our value (5 mm/sqrt(h)), but we believe that their value (3mm/sqrt(h)) is too tight.

R1C23: question on the “an ideal and normal case”

AC23: we assume an ideal case with almost perfect observation geometry without any blocks, i.e., the 5 degree cut-off elevation angle. However, it is quite often that the low-elevation angle observations have poor quality and can be unavailable, e.g., blocked. Therefore, we also assume a "normal" case with the cut-off elevation angle of 15 degree.

R1C24: How about the result with 3  and 4 \sigma?  From the result with 5 degree cutoff-elevation, with 3 \sigma may still shorten the convergence time. Furthermore, what does it mean that the result with enlarged STD as constraint performs better? Are the uncertainty information  of the empirical model reliable or not?

AC24: we tested the results of 3 and 4 times of the sigma values, and the convergence time tends to perform comparable to that of solution “No”, i.e., no constraints (or extremely loose constraints). This is also expected. We only show the results of the 1 and 2 times of sigma as they are sufficient to illustrate our conclusion: the modeling of formal error can effectively improve the convergence time if it is used properly, i.e., with the scaling factor.

The uncertainly information of the empirical model, i.e., formal error, is reliable as it agrees well with the RMS of residuals when comparing to NWM and GNSS ZTD. However, when using the constraints in the GNSS PPP solution, we have to consider other effects, such as the temporal constraints. When both the temporal constraints and the absolute constraints of parameter (i.e., using the formal error of our model) are applied, the actual constraints applied to the parameter is tighter than the absolute constraint. We considered to perform a test with only absolute constraints on the parameter w/o any temporal constraints, which can better illustrate the effectiveness of our model. However, in that case we cannot conduct a fair comparison as the reference solution (baseline solution “No”) without the absolute constraint of the parameters still need the temporal constraints.

 RC2

Dear authors, I want to thank you for your interesting work and mostly well prepared manuscript. Still, I have some comments which should be addressed before the publication.

Response: thank you so much for your encouraging comments and the helpful suggestions. We have modified the manuscript accordingly. Please find the point-to-point response as follows.

Major comments

1, In the manuscript, you provide a broad set of figures showing various detailed results. In the text, you mention only a few statistical parameters. I would like to see more of these statistics in your manuscript. Not necessarily in every (sub-)section, but e.g. in:

a, section 2.1: provide not only mean RMS and range of RMS, but also other information. E.g.: is there any bias of the fitting residuals? What is their mean standard deviation?

Response: thank you for your advice. As we fit ZTD time series via the least-squares adjustment, the mean value of the fitting residuals at each grid point (one time series per point) is always zero, and the standard deviation is the same as the RMS value. Therefore, we present only the RMS values in section 2.

b, section 3: only some information about average RMS is given. Please, include more statistics for ZTD residuals (at least mean = bias, standard deviation). How much the statistics differ between individual years?

Response: in section 3 we evaluated the our model w.r.t NWM-based ZTD values in the period of 2009-2021. As the period of 2009-2018 is used for modeling, the average bias is zero and the RMS values are already in Section 2.

For the prediction period of 2019-2021, we added yearly statistics and the discussion, as shown in the following paragraph and table. The differences between STD and RMS statistics are mostly within 1 mm, thus we only give the Mean and RMS values.

Table 1 summarizes the yearly statistics of the Mean and RMS values in the prediction period (2019-2021). The year-to-year changes are rather stable, as the average bias varies within 1.5 mm and the average RMS varies within 1 mm. On the other hand, the maximum values can vary up to 10 mm, for example, the maximum bias of all grid points in 2019 is 41.2 mm and that in 2021 is 30.4 mm. This is expected as the minimum and maximum values are more relevant to local severe weathers and can change significantly. The average values indicate that the overall performance of the model is stable in different years.

Table 1 Statistics of the ZTD differences between modeled values and NWM-derived values in the prediction period (2019-2021), including the yearly values and the value over the whole period. For each grid point we calculated the Mean and RMS values, and then present the average, 95% of absolute, minimum, and maximum values of all grid points. Units are mm.

As the mean bias tends to be zero, the STD is similar to the RMS, and thus we do not repeat it. In section 4 where the mean bias is more prominent, we presented more statistics, including the mean, STD, and RMS values.

c, section 4: if possible, add more summarizing information about the differences between modelled ZTD and GNSS ZTD (I see only average bias and RMS values). E.g. results for individual years, range of bias/rms values at individual stations, etc.

Response: thank you for your suggestions and we have added the following paragraph and table.

The agreement of modeled ZTD w.r.t GNSS estimates is summarized in 2. The model shows no systematic biases as the averaged bias over all stations is -0.1 mm, despite that the maximum biases can reach up to 2 cm. The average absolute bias is 4.1 mm, which can be attributed to (1) the mismodeling effects at specific stations due to the deficiency of our model, which can be further improved by adopting a higher temporal resolution, and (2) the systematic biases of GNSS ZTD estimates due to the instrument effects (Ding et al. 2023). The RMS value varies from 19.2 to 63.8 mm, with an average value of 37.7 mm. Both the bias and RMS values agree well with previous investigations of the GPT3 empirical ZTD models and GNSS ZTD, for example, an RMS of 44.1 mm by Ding & Chen (2020), and an RMS of 38.56 mm reported by Yao et al. (2024).

Table 2 Agreement of the modeled ZTD w.r.t. GNSS estimates. For each station, we calculated the Mean, Median absolute error (MAE), and RMS values of the ZTD differences over the four years (2017.00-2021.00), then the average, median, mean absolute value, maximum, minimum, and 95% values of all stations are given. Units are mm.

d, section 5: you provide figure 14 and corresponding text which just generaly summarizes which solution is improved/not improved. Please, add also specific numeric information about the size of improvement. E.g. a table showing how much the solutions with 1/2 sigma uncertainties improved compared to the solution with no added uncertainties (a percentual improvement/worsening).

Response: thank you for your advice. However, constraining tropospheric delay with proper weighting mainly contributes to speeding up the convergence time, instead of improving the positioning accuracy. Therefore, we try not to focus on the “accuracy”, but instead, concentrate on “convergence”. We thus present the following table in this response letter, but not in the revised manuscript.

Table 3: Relative improvement (in %) of kinematic PPP accuracy w.r.t. the reference solution, i.e., w/o tropospheric delay constraints (solution “No”) in different arcs. A positive value indicates that the solution is improved w.r.t. the reference solution and a negative value (in red)  indicates that the solution is worse.

2, Although your Conclusion and Outlook section contains some elements of a scientific discussion, you should elaborate this topic more. Please, discuss more the quality and sufficiency of your methods and experiments. I also wonder if it would be possible to include the diurnal cycle into your model as it can be rather strong for ZWD in summer periods.

Response: thank you for your suggestions. We agree that the diurnal cycle can be strong for ZWD in summer periods at specific regions. However, the overall magnitude of diurnal cycle on a global scale is rather small, especially for the empirical delay models. The focus of our study is the uncertainty information, and ignoring the diurnal signals can hardly affect our model performance, especially considering the fact that our modeling data is the 6-hourly sampled time series. The diurnal signals (or sub-daily atmospheric tides) can be better described if we have higher temporal resolution, e.g., the 3-hourly ERA5 profiles. That is, however, not the focus of this work. On the other hand, we believe that to further improve the accuracy of empirical models, more sophisticated models, such as those AI-based and/or data-driven methods, should be adopted, as we discussed in the conclusion section.

Minor comments

L31: correction: replace ".. of water vapor in spatial and temporal" with "... of water vapor in space and time". Or add a noun after the "spatial and temporal" words.

Response: we appreciate your suggestions and have corrected it to ‘in space and time”

L35: in the first sentence you are writing about various techniques for tropospheric zenith delay estimation. Some of them are techniques or instruments, other models. However, you start your second sentence with words "Despite the relatively high accuracy of these models, ...". What do you mean with the term "models" in this sentence? Please, rephrase the text to clarify. Later, you write that "these methods require accurate meteorological information". E.g. radiosonde does not require any meteorological information as it is an instrument to measure them. Please, revise the paragraph (e.g. strictly separate measurement techniques/instruments and models). The current text can be confusing for reader.

Response: thank you and we agree with your comments. We have modified the manuscript as follows.

Tropospheric zenith delay can be derived using several approaches, including in situ meteorological observations combined with models such as the Saastamoinen or Hopfield function (Askne and Nordius, 1987; Hopfield, 1969; Saastamoinen, 1972); radiosonde measurements, which provide vertical profiles of meteorological data (Chen and Liu, 2016; Liou et al., 2001; Wu et al., 2019); and water vapor radiometers, which directly measure atmospheric water vapor content (Braun et al., 2003; Niell et al., 2001). Additionally, numerical weather models (NWMs) offer tropospheric delays on a global scale by integrating meteorological data vertically (Böhm et al., 2007; Kouba, 2007; Landskron and Böhm, 2017). Although these approaches can yield high accuracy, typically within 5 mm to 2 cm, they all require access to accurate meteorological information, which may not be available to real-time GNSS users.

L54: typo NECP -> NCEP

Response: done.

L105: I would suggest to provide two figures instead of one. In the first one, you would show the height of VMF3 grid points. In the second one, you would show position of GNSS stations. At the moment, the figure is a rather chaotic and e.g. in the central Europe, it is not possible to get any information about grid point heights/position of GNSS stations as everything is shown in white.

Response: thank you and we have replotted it.

Figure 1: Altitude (ellipsoidal height) of the VMF3 grid points (left) and distribution of GNSS stations (right). The GNSS stations used for ZTD comparison and PPP validation are given in blue and red dots, respectively.

L224: you write that in periods with peak of formal errors, the ZTD residuals could be extremely large. You write that reason of this situation is in "abundant water vapor and more likely happening extreme weather conditions". Please, where do you have any proof or support for this statement? Please, explain why do you think so.

Response: thank you for your question. As the empirical model is built using the long-term analysis of the ZTD, it represents the averaged (or smoothed) pattern of the local weather conditions. For the extreme weather events, the water vapor can be abnormally large or small and thus deviate from the usually cases. As a consequence, the ZTD residuals can be extremely large. In other words, the empirical models using only annual and semi-annual signals cannot describe the high-frequency local weather events, resulting in the large ZTD residuals.

L230: figure 7: please, increase size of the figures as they are hardly readable

Response: we have replotted the figures with larger font size.

L239: you are writing about 2017 to 2021 period, but do you mean just 2019 to 2021 period which you used for prediction evaluation? Similarly, in Figure 8 caption, you write about 2018-2021 period. Please, check and correct.

Response: sorry for the misleading information. These are two typos. We have corrected it.

L330: although the figure caption decribes yellow and green dots, I can see only orange and green dots. Please, correct.

Response: done.

L335: from which year you used GNSS data? I do not see this information.

Response: sorry that we forget to mention it. We use the data of 2020. Corrected.

L333: here you write about GNSS positioning, on L335 you mention usage of GPS observations. Please, clarify which GNSS systems were used in your data processing. If you used GPS-only solution, please explain this setup and discuss whether usage of multi-gnss would lead to different results (as nowadays a multi-gnss or at least gps+glonass processing is a rather standard, at least in scientific studies).

Response: we use “GNSS positioning” as a more general term as GNSS can include different constellations. In our experiment, we used GPS-only observations. For the multi-GNSS scenarios, the benefits of using external tropospheric delay model with weighting information is less significant, as the many satellites already provide a rather good observation geometry.