General Comments 

The work presented in the manuscript address the question about long term stability in estimated atmospheric propagation delays using ground-based GNSS stations.

The part I find most interesting, and that may be worth to be published, is the assessment of estimated trends and how these depend on the used mapping functions and the elevation cutoff angle. I think this part is an important contribution to the community but it needs to be more critical. In the present version of the manuscript I think the results are overrated.

We cannot speak about an optimum elevation cutoff angle in general because it is station dependent, i.e. the time dependence of systematic errors in e.g. mapping functions and the multipath environment. Therefore, the presented results are not necessarily in contradiction with those presented by Ning and Elgered (2012) and Baldysz et al. (2018). There is no conflict between these results because the stations analysed in this manuscript have almost no overlap with those in the other studies. Ideally, without systematic errors,  the estimated trends shall be identical regardless of the elevation cutoff angle. This is different to the individual ZTD estimates where the geometry obtained for low elevation angles reduce the errors in the estimates (also for the estimated coordinates and especially the vertical). When trends are estimated individual errors are averaged out, if no systematic errors are present.

There ought to be a critical discussion about the uncertainties of the estimated trends as a base for a statement regarding which differences that are significant. For example, which differences seen between the estimated trends using the different elevation cutoff angles in Figure 8 are significant.  Noting the consequences from introducing changepoints as described, I think this shall be analysed in more detail.

Specific comments

Line (L)103: Radiosonde (RS) data were processed by Dai (2011) and are used as a reference. In the study data up to 2014 are used. This requires an explanation. How did you handle RS data acquired in the years thereafter?

L126: I assume that when mapping functions are compared in Table 3, all these solutions are carried out using an elevation cutoff angle of 7°. Can you mention this explicitly? Please also comment on to what extent you find the differences in Table 3 significant.

L 140: As I have understood the RMS is defined as the root-sum-squared of the standard deviation and the bias. But this is not the case in Tables 5 and 6. Please explain.

L155: The ABS method suffers from the fact that if an unusual cold and dry month is followed by an unusual warm and humid month a false detection is likely. This ought to be discussed and the different criteria used to identify a changepoint shall be stated.

L167: I agree with your conclusion that the REL method shall not be used when the goal is to compare "before" and "after" with the ERA5 (because an improvement is expected when the reference data set is used to add changepoints in the GPS time series, the agreement between the trends is of course improved.

The robustness of the trend results after adding changepoints can be assessed by studying subsets of the data and the stations.

Some suggestions related to Table 7:

(i) Apply changepoints only for the events that can be supported by the station log. 

(ii) Apply only those changepoints when Offset 1 and Offset 2 differ by less that a certain value. The fact that some of them are very different, as well as having opposite signs, I think is warning to be very careful.

(iii) A combination of (i) and (ii).

L182: Figure 5: The changepoints seen in the figure are not the ones in Table 7. Are not both of these carried out using an elevation cutoff angle of 7°? Furthermore, the ones in Table 7 are not supported by station logs. I think that if you present such results as in Figure 5 you should discuss them in more detail and arrive at some understanding why the two mapping functions result in such different trends. Can anyone of them be trusted?

Figure 8: Assuming that the work by Dai (2011) implied a significant improvement in the RS data, the results for the Raw comparison may be ignored. Adding that the introduction of changepoints seems to be a rather inaccurate method, the Dai and the ERA5 comparisons before homogenization are the most interesting. It is also worth noting that these two also give the best agreement for elevation cutoff angles of 20° and below. Using GPS satellites only (and not multi-GNSS) means that there are much less observtions for cutoff angles above 20°.

Technical Corrections 

Line(L) 6: "Homogenized atmospheric water vapor" sounds strange. To me it sound like something done in a chemistry lab.?

L6+: You use the American spelling of vapour, although ACP is a European journal?

L11: the word "latest" may not be true if and when the manuscript is accepted for publication.

L14: 0.3 mm/yr  -->  0.3 mm/year  (and a few more places in the manuscript. Note that there is no symbol for "year" in SI, although some use "a", for annual)

L23: 7%  -->  7 %  (see also line 131)

L25: There are more recent IPCC reports. Although it does not change the statement it would be more relevant with a more recent one.

L80: 300s  --> 300 s

Table 1: Perhaps it will be more clear if you note that the E5 solution is used both in the mapping-function comparison and in the elevation cutoff-angle comparison?

Table 2: The unit for the random walk shall not be in italic font

L110: Equation (1) would be informative to explain a bit more so that an overall understanding can be obtained without reading the reference. For example, are the i and j terms all possible combinations (where t_j > t_i) or adjacent values only? Please also define "hat x" in Equation (2).

Figure 2 (and Figure 5): Remove the text above the graphs and add it with an explanation in the figure captions?

Figure 4:  Should not the green bars to the right in the graphs be blue (rather than green). The way I interpret the text is that there shall be one green and one blue bar for each mapping function?

L189: Y-axis label is missing

Figure 6: top and bottom shall read left and right.

L216 (and other places):   homogenezation  -->  homogenization

L219: 30-yr  -->  30 years

L285: A doi link is missing, also for some other references and the established standard acronyms for journals are not used in all cases. Furthermore, sometimes they are given as "https://..." addresses and sometimes just as "doi:..."

General comments:

The paper „Impact of processing strategies on Long-term GPS ZTD” concerns investigating the influence of the selected GNSS observation processing strategies on the reliability of position and ZTD.  In general, the Authors compared the GPS processing approach, which differs in mapping function (GMF, GPT2, GPT3, VMF1, VMF2) and elevation cut-off angle (3°, 7°, 10°,15°, 20°, 25°, 30°). They have mainly focused on the ZTD time series but also provided some basic results regarding position repeatability. Although GNSS meteorology is a well-known concept, it still requires improving existing algorithms and validating possible/new solutions, including new mapping functions. In light of this, the general idea of the paper is justified. However, the complexity of the impact of individual observation processing elements on the reliability of the final solution is very high. Hence it requires a very detailed analysis, which in my opinion, has not been done by the Authors. 

Specific comments:

Firstly, there is no information about trend estimation uncertainties, which are significant when assessing various solutions. Some differences between different observation processing strategies are expected, but assessing their significance is the most important.

The Authors have analysed 46 IGS stations, while only 19 have presented results in Figures 4 and 7. There is no appendix to see what is happening with the rest of the stations. Figure 8 presents results for all stations?, but it is unclear. Additionally – we do not have any information about data quality. The data completeness probably varies for different stations and may affect final solution.

I am also wondering why the Authors have used 1995-2014? Before 2000, quite a poor quality of orbits and SA negatively affect GPS solutions.  The station selection is also questionable – 100 km is a lot and may result in different troposphere conditions. Here a table with exact differences between GPS and RS sites is necessary. Moreover, Dai et al. (2011) presented a homogenised dataset until 2011 (or at least that's what the text says). But if the GPS data were processed until 2014, what was the reference for the last three years?

Why did the Authors decide to verify different mapping functions using 7° cut-off angle?

The Authors wrote, “The method for calculating ZTD from ERA5 can be referred to Haase et al. (2003)” – please be more specific about whether used by the Authors method is the same as in Haase, or not. Additionally, what is a temporal and (even more important) spatial (vertical) interpolation between ERA5 and GPS site – there is no information about this.

The ERA5 homogenised dataset (according to Dai et al. 2011) should cover a time span until 2011. Please be more specific about the exact source of radiosonde data (the link given in Dai et al. 2011 does not exist at this moment). This also concerns ‘raw Radiosonde’. That would be helpful for the readers.  Additionally, please add the info on whether the Authors used exactly the Haase (2003) method for calculating ZTD from RS.  

I’m not sure why the Authors have focused on analysing position accuracy since there are no conclusions (just a description of the results) and, more importantly, the results from this part of the manuscript were not considered in any other part. The small variability of position is rather excepted and obtained differences are very small (hundredths of a millimetre).

There is also no specific conclusion from analysing bias, STD and RMS from differences between GPS ZTD and ER5 ZTD

There should be more discussion regarding the impact of homogenisation on long-term trends. It is clear that adopting various homogenisation approach influence the final solution the most (since homogenisation may ‘fix’ even distinct inhomogeneities resulting from adopting various processing strategies). Several papers concern different methods of GNSS time-series homogenisation. It is unclear from what the Authors wrote whether the changepoints they found are correct, better/worse than changepoints that may be found with other approaches.

Figures 4 and 5 make me worry about the reliability of the homogenisation process. After taking a closer look at e.g. JOZE station, we can see that the trends are similar before homogenisation, while after homogenisation there is a distinct difference between VMF1 and VMF3. These mapping functions rely on the numerical weather model and are very similar regarding the a and b coefficients. Therefore such differences are unexpected. BRMU station also looks interesting – before homogenisation all trends are similar, after homogenisation, there is a distinct division between climatological and discrete mapping functions. At this point, I would not worry about the comparison to the RS since it may even be 100 km away (there is no info about that).

Figure 5, in turn, makes me worry about the reliability of the GPS observation processing. Presented by the Authors monthly ZTD anomalies present a distinct shift in the case of GPT3 mapping function, while using VMF3 there is no such situation. The main problem is that GPT3 is a climatological mapping function and is therefore continuous. Therefore presented by the Authors shifts in this particular solution are not a problem of GPT3, but of the processing itself.  

I am also not sure why the Authors focus on ‘Raw radiosonde’ as a reference since they stated in the introduction that RS homogenisation is important. I am also not sure why the Authors focus on the un-homogenized GPS ZTD time-series and, based on them, assess various cut-off angles. “However, for other situations, i.e., taking Dai- or ERA5-derived ZTD trends as references for un-homogenized GPS ZTD evaluation….”. If we already know that GPS time series may be affected by various factors (e.g. antenna/receiver changes), why should we focus on un-homogenized ZTD, while comparing it to the reference set?

To all figures and tables – please change their description to make it possible to correctly understand the presented in them results, without looking for information in the manuscript’s main body.  Figures are often not well readable.

Overall it seems that the presented paper covers too many issues that are too briefly analysed. A proper analysis of each of its elements (i.e. the impact of processing strategy on ZTD, the impact of processing strategy on position and homogenisation on long-term ZTD reliability) is a big task. Therefore it is hard to find reliable outcomes from the conducted analysis.

More detailed comments:

Figure 1  - there is no ‘BOGO’ station in IGS

Page 6, Table 4 – please add info that all cut-off angles were tested using VMF3 (I know it was pointed out, but the table should be read correctly, without looking for further information in the manuscript body

Page 7 Tables 5 and 6 – same as above, but regarding cut-off angle, and mapping function

Page 7, Tables 5 and 6 – please add info that this is a difference

Page 10, Figure 5 – add y-axis description to the figure

Page 9, figure 2 – the colours are way too similar. Instead of the legend, I suggest you add the solution name to the axis

Page 11, line 208 – shouldn’t it be Baldysz et al.2016?

Page 12, Figure 8 – The figure description should be corrected (left/right instead of top/bottom)           

Page 13, lines 233-235 – this is rather expected. Since we estimate differences between GPS ZTD and ERA5 ZTD and then use these differences to correct GPS ZTD time series, the final GPS ZTD solution will be similar to the ERA5