Please see the attention for the detail response.

My response:

Thanks for the reviewer’s suggestions and we have made changes for all the suggestions in the new manuscript accordingly. The detailed modifications are as follows.

1) The motivation of the study includes:

1. a) Improving the blind models used in real-time retrieval of GNSS-PWV under the condition that meteorological sensor are not equipped at the GNSS stations.
2. b) Providing an alternative method in consideration of the inconvenience in acquiring forecast pressure from numerical weather models. The forecast data need to be downloaded in advance, which increases the complexity of data processing, not to mention the fact that the forecast data may not be obtained due to various reasons, e.g. problems of some servers or agencies.
3. c) Coming up with a new method for the determination of the ZHD in the retrieval of GNSS-PWV.

2) Line 54 has been corrected.

3) There are no topography and interpolation used in Figure 8. When ZHD-ERA5 was used as the reference, the ZHD-ERA5 at a grid point was calculated by the integral from the lowest level (1000 hPa) to the highest level (1 hPa), and the ZHD-GPT3 at the same grid point was calculated using the Saastamonien model with the input pressure derived from GPT3 at the lowest level. Note that the geopotential height in ERA5 was converted to ellipsoidal height using the method described in the following paper:

Wang, X., Zhang, K., Wu, S., Fan, S., and Cheng, Y.: Water vapor-weighted mean temperature and its impact on the determination of precipitable water vapor and its linear trend: Water Vapor-Weighted Mean Temperature, J. Geophys. Res. Atmos., 121, 833–852, https://doi.org/10.1002/2015JD024181, 2016.

The reason for the large biases in Figure 8 may be that the GPT3 model performed poorly in the regions where the height is significantly different from the earth surface. The height of the grid points from ERA5 was above the earth surface in most regions but below the earth surface in mountainous regions, since the pressure at the lowest level of the grid points from ERA5 is almost constant (1000 hPa) but the heights at the lowest level in different grid points are different.

4) Figure 3 caption has been corrected, as advised.

5) Lines 110, 277 and 305 have been corrected, as advised.

6) In addition to the above suggestions, we also found a new problem: when the integral method was used to calculate the ZHD, geoidal height or ellipsoidal height should be used, so the geopotential height contained in ERA5 and sounding data at each pressure level needs to be converted to geoidal or ellipsoidal height. However, in the previous manuscript, the height was not converted. Although this fault might cause large biases in the ZHD, it is not necessarily to largely affect the RMSE of the ZHD. This is due to that the effect of this fault changes biases from positive to negative, but the values are about equal. In the new manuscript, all the problems have been corrected.

Please see teh attention for the detail response.

My response:

Thanks for the reviewer’s suggestions and the modifications made accordingly in the new manuscript are as follows.

1) The GZHD model was evaluated using out-of-sample data: the radiosonde and ERA5 data during the one-year period of 2020 (in which no data were used to train the model). Results showed that the performance of GZHD was almost the same as that in the training period (please see Section 3.3 in the new manuscript).

2) The effect of the GZHD model was tested on real-time GNSS-PWV and results showed that GZHD outperformed GPT3 (please see Section 3.4 in the new manuscript). Note that the real-time ZTD used to retrieve GNSS-PWV was for 154 days in 2020, and it was processed by a modified BNC software package, more details can be found in the following paper:

Sun, P., Zhang, K., Wu, S., Wang, R., and Wan, M.: An investigation of real-time GPS/GLONASS single-frequency precise point positioning and its atmospheric mitigation strategies, Meas. Sci. Technol., https://doi.org/10.1088/1361-6501/ac0a0e, 2021b.

3) A reference of GPT has been added in line 52.

4) The color in Figure 1(d) agrees very well with the terrain, e.g. the Tibet Plateau. This is because the penetration depth of a profile has to be above the earth surface. In addition, there are over 3,000,000 profiles in this figure, thus some of the profiles are covered by others.

5) In addition to the above suggestions, we also found a new problem: when the integral method was used to calculate the ZHD, geoidal height or ellipsoidal height should be used, so the geopotential height contained in what data?? at each pressure?? level needs to be converted to geoidal or ellipsoidal height. However, in the previous manuscript, the height was not converted. Although this fault might cause large biases in the ZHD, it is not necessarily to largely affect the RMSE of the ZHD. This is due to that the effect of this fault changes biases from positive to negative, but the values are about equal. In the new manuscript, all the problems have been corrected.