The authors present a unique and important new dataset of Earth thermoelastic vertical displacement at 15000 GPS stations globally. This is an important source of GPS signal, especially at annual and semi-annual periods, and yet there is no operational or even non-operational widely available product for users to apply to GPS bedrock time series in order to reduce noise and signal. The methods look solid and robust in the main. The work builds on previous work published with some methodological advances. The dataset and quality control is well described. I expect this will be a widely used dataset, although note there is no statement on if the dataset will be kept up to date - something that will terminate its usefulness as GPS time series get longer. Likewise, the decision to commence in 2000 may cut off the greatest value. During the sparsely observed (by GPS) 1990-2000 period the effects of signals like this will be even greater in some studies due to asymmetry of the global netowrk. it also leaves users wondering what to do with older data.

I have three major concerns and some other suggestions and points.

First, as with all previous studies of global thermoelastic deformation, the authors mishandle areas of the planet that are permanently and seasonally ice- or snow-covered. Take Antarctica for example. The permanent ice cover there insulates the sub-ice region from experiencing any annual thermal effects. This presumably affects the deformation of sites located next to the ice sheet or on rock outcrops through the ice sheet. This will include Greenland, and the many small glacier regions globally. As such, I suspect the modelling is unreliable in these regions and should be removed from the dataset until the model is revised.  To investigate this, I looked at two Antarctic sites CRDI and MAW1, which are adjacent to the Antarctic Ice Sheet, on small (~100x100m outcropping bedrock). The modelled series (see in the PDF) show significant annual signal and, in the case of CRDI, an apparent offset in 2010. My initial reaction is that these are likely not robust signals. However, maybe the thermal effect is of such small scale that even 100mx100m outcrops are enough to experience these thermal deformations (although presumably the very cold ice to the side of the sites also regulates the bedrock temperature considerably). I suspect the CRDI offset is not caught by the QC as there are few sites in Antarctica. But I also presume this is driven by an ERA5 temperature artefact. Also, there does not appear to be any data in the SoilsGrid paper under the permanent ice, so it is not clear how these regions are modelled anyhow.

Seasonal snow cover is more complex. This would be available in ERA5 I guess, in addition to temperature. Again, snow will insulate the soil/bedrock from air temperatures.

Second, I was not sure how ocean and lake areas are considered and what resolution any mask has. There are many coastal GPS sites and presumably this will matter.

Third, the data itself has no version control. I would expect there will be revisions and hence the dataset should be versioned. The data instructions are not clear that the trend component is not to be trusted.

Finally, it is a shame that the authors only released the vertical components. It would have been interesting to see the horizontals. Also, it would be useful in some cases to consider subdaily corrections for thermal expansion, and hence, the 6h series would have been useful, not just the daily ones. Maybe that will be a focus of future work, including validating the subdaily signals. Given that they are discussed here, it would be appropriate to include a map of subdaily amplitudes and phases, although I guess these are also seasonally modulated.

Minor comments

L15 I am not sure the current TED is reproducible. That would require code and workflows. I think you can say comprehensive.

L24 add how many stations had independent site characterisations

L32 please express what % of stations not just vague 'most'. Here and later, the authors merge the effect of TED into the combined effect with loading products. The authors should stick to just the TED effect here in the abstract. it is no more than +-10% at most sites. That is ok.

L49 obscure *other* subtle

L66 this is an important point. you could say a little more. there is no CM movement in this case, even if there is CF and CE movement.

L69 This tidal aliasing type effect could reference Penna and Stewart 2003 GRL as an e.g.

L103 this is where I would like to see versioning

L114 in this section it would be good to say what is not modelled. Buildings and monuments above ground, such as pillars. These may cause significant effects, and the authors could cite such works.

L149 gives the horizontal resolution of ERA5 forcing

L189 give the horizontal resolution of the FEM

L230-233 I didn't understand this sufficiently well to be able to reproduce it. Why not just take daily averages?

Results. consider adding a map of the trends. I think more is required in terms of discussing hte correlation with ERA5 trends or where else trends are coming from. A map of trends of the two (forcing and response) would be helpful to know if these may be real or a modelling artefact of kinds.

L367 this paragraph, please explicitly cite the panels of Fig 6 as relevant

Fig 6 caption. explain all panels explicitly

L388 'weak' downplays it without showing us the trends. Delete weak and let the reader decide that. change long-term to trend-like.

L437 I did not see Neff defined.

L438 I think a new subheader is appropriate here.

Section 3,3 As also noted for the abstract, the merging of the TED product with surface loading products, some with much larger importance, is not satisfactory. Please commence with the TED effect on WRMS. You could reasonabl remove first the other signals from the series.

Figure 14 labels the y-axis as %. Note the different colour scale for NTAL and HYDL

L646 new para at 'in practical'. Give TED only results not the merged results with NTAL etc.

L657 add also above-ground pillars.

L662 new paragraph

L668 new paragraph

Matt King, March 6 2026

This paper models the effect of the thermoelastic component of ground displacement driven by ERA5 temperatures down to 2 m of depth and soil properties also down to 2 m using a soil data base.  It presents the result of a thermoelastic model at > 15,000 sites and it provides time series for the model output for each station in a companion website.  It is evident that the sub millimeter annual amplitudes of the displacements for the layered model are different than those for a half-space model as the physics would predict and might warrant inclusion into GNSS post-processing analyses.  The paper represents a tremendous amount of research. 

90% of the manuscript is dedicated to comparing the half space model with the layered model.  It is clear that the physics of the layered model represents an improvement over the half-space model.  But this is not tested with real observations.  At the end of the paper, there is a comparison of the layered model with GNSS observations.  The improvement of the model to the WRMS of the GNSS data is very minor.  But this is where the half-space model and the layered model should be compared.  The WRMS reduction of the GNSS observations using both models should be compared.

The paper is generally well written and well referenced. It could be shortened substantially that would help to make it more focused.  There are just too many comparisons of the half space model with the layered model.  What the reader wants to know is if the correction reduces the RMS of the GNSS observations in any substantial way.

The question remains are submillimeter annual signals sufficiently significant compared to model limitations and other uncertainties (e.g. monument thermoelastic effects, model errors/simplifications, GNSS uncertainties) to justify including this model as standard GNSS post-processing practice?  In this respect, the authors could provide arguments of scale to justify the correction.

1. The biggest concern is a robust validation of the model. Does the RMS of the GNSS up time series decrease across the majority of stations when the model is applied?  The reader does not learn the answer to this question until page 22-23.  The abstract states ‘When TED corrections are applied together with non-tidal mass loading correction, the residual vertical dispersion decreases at most stations, with vertical scatter reduced by up to 70% at selected sites’.  This analysis does not test the thermoelastic model in isolation.  It tests the thermoelastic model + large amplitude loading effects.  We cannot determine whether the TED corrections or the loading models cause the change in the RMS.   This comparison is misleading.  Near the end of the manuscript, the authors finally report that ‘TED provides a more ubiquitous but generally modest benefit: dispersion reductions are typically a few percent (often ~5%)…’
2. The paper does not present any order of magnitude comparisons for the reader. How big are these Earth thermoelastic effects compared to monument thermoelastic effects?  How big are these thermoelastic effects with respect to loading (there is a figure late in the paper but a figure showing the global effects of a half-space model versus the loading early in the paper would be useful. And the figure should be a global figure not only a few stations.) How sensitive is the model to sediment thicknesses and thermal/mechanical properties? A single plot that shows the Earth thermoelastic effects over the continents at equal grid spacings would be extremely helpful for the reader to interpret the significance of the results and to determine the spatial (latitudinal variability).  In this vein, it would be useful to see a plot of sediment thicknesses over the continents.  When I compare the differences between the half space model and the Global 1-km Gridded Thickness of Soil, Regolith, and Sedimentary Deposit Layers from ORNL DAAC there is a lot of correlation.

Additional Comments

1. Soil model: The authors use sediment characteristics down to 2 m. For greater depths they assume a half-space model assumed to be bedrock. Do you stop at 2 m because the thermoelastic effect diminishes at this depth or because you do not have model data. 
   1. It would be useful if they presented the results for an order of magnitude calculation. To what depth are thermoelastic effects important? Of course this varies by material.  But the authors could provide a half-space model for two end points of high thermally conductive and low thermally conductive material to demonstrate how deep the effect can go.
   2. ALSO, what does the soil model look like for the edges of ice covered regions where the temperature range is usually damped by the low thermal conductivity of the oceans?
   3. Can you include a representative figure of your soil model?
   4. Can you provide a description of the half-space model? Sorry if I missed this. but I’ve gone through the paper twice trying to find this description.
   5. It appears soil moisture content of the soils are not considered. Could you provide an amplitude analysis and a discussion of how the results would change with changes in soil moisture content?  The authors mention vegetation as a controlling factor. It isn’t the vegetation; it is most likely the soil moisture.
2. TED model. Please instead of showing the model results at the locations of the GNSS sites, calculate the annual model at equally spaced grid points over the continents.  This would help the reader to visualize the spatial dependence of the temperature forcing over the continents.  It is difficult to glean this information by plotting only the model results at the GNSS locations.
3. P9; L 277: where are these results presented?
4. Sensitivity analysis: where are the results shown? What are volcano plots?
5. Figure 3:  where are these stations and what climatological regimes do they represent?  Are these daily model solutions? Please plot using a consistent y-axis in both panels.  Is the trend in most plots real global warming or is there a modeling issue in ERA5? Please estimate the amplitude of the monument thermoelastic signal at these sites for comparison.
6. Sensitivity analysis: where are the results shown? What are volcano plots?
7. Figure 5 peak to peak differences looks very similar to Figure 4 annual differences. Please comment.
8. Figure 6: Caption is inadequate.
   1. Figure 6a is peak to peak displacement. It is not clear what is plotted in figure 6b.  The mean of the peak to peak displacement or something different. 
   2. Would the latitude means (6b) look any different (smoother) in the southern hemisphere if you calculated the TED over an equally spaced grid over the continents? The number of sites in Asia, Africa, South America and India are sparse.
   3. You should show the zonal mean temperature for comparison in Figure 6b.
   4. Figure 6e is not explained. What are the implications of the results?
9. Linear trends in the TED model. A reader cannot interpret figure 7 because we cannot see the lengths of the vectors when there are many sites.  Again, it would useful to calculate TED over an equally spaced grid on the continents and use a color legend.
   1. ERA5 is expected to contain trends due to global warming.  The trends seem to range from +- 2 deg. C per decade (see Evaluating global and regional land warming trends in the past decades with both MODIS and ERA5-Land land surface temperature data, Wang et al, 2022).  A discussion of the trends in the soil temperatures in the ERA5 data should be discussed. 
   2. The discussion should include a discussion of the uncertainty in the soil temperature trends.
10. Figure 8 does not add anything to the paper. The authors are comparing the half-space model and the layered model in a slightly different way. However, at this point the reader gets the point that the models are different and produce different results.  What we want to know is the validity of any of the models.
11. Figure 9 and its discussion is about validating the SoilGrids in the FEM_FL It is not about evaluating the robustness of the FEM_FL model. In both cases you are using the FEM_FL model and comparing output using 1) SoilGrids as input and 2) site geology information as input. 
12. No comments on figures 10-11. I feel like these figures do not help the reader gain insight into the physics of the problem. 
13. Figure 12 is valuable to show the relative amplitudes of the loading versus the thermoelastic signals.
14. The discussion starting at L556 regarding the WRMS reduction is important and should be highlighted in the abstract. L570-571 in particular are very important and should be stated in the abstract.  This is what the readers want to know.
    1. And it would be good to know how the WRMS reduction from the half-space model compares with the layered model.
    2. Figure 14 only the upper left panel showing the WRMS reduction from TED is important in this figure.
    3. Figure 15 is misleading and should be removed.
    4. I feel that removing the loading models just muddles the conclusions. The GNSS data WRMS change should only consider TED.