General Comments

This manuscript presents a significant and timely contribution to the field of satellite gravimetry by introducing HUST-Grace2026s, a new static gravity field model that consistently integrates over 20 years of data from both the GRACE and GRACE-FO missions, thereby bridging the gap between these two satellite gravimetry eras. The authors apply advanced processing strategies—including stochastic modeling based on postfit residuals and rate-term estimation—to derive the model, and provide a thorough validation against independent data sources. The work is well-structured, technically sound, and offers valuable insights into the challenges and opportunities of extending static gravity field solutions into the GRACE-FO era. However, several aspects of the manuscript would benefit from clarification and revision to enhance its clarity, consistency, and impact. The following comments are intended to help the authors strengthen the presentation and address a few issues that may affect the interpretation of their results.

Comment on key outcome 1 (P26;)

The first conclusion states that "processing chain refinements can play a more critical role than the simple addition of GRACE-FO data." While this finding is supported by the analysis of the current dataset (GRACE-FO data up to March 2025 at ~500 km altitude), the statement as written lacks the temporal qualifiers that the authors themselves include in the main text (e.g., in Section 3.1.2, they note that GRACE-FO "currently operates at a higher orbital altitude of about 500 km"). To maintain consistency and avoid overgeneralization, it is recommended that the authors explicitly qualify this conclusion in the same manner, noting that it applies to the current orbital phase of GRACE-FO and that this assessment may change as the orbit decays in later mission phases.

Comment on Contribution Matrices (Figures 10 and 11)

The contribution matrices in Figures 10 and 11 are presented as revealing the "relative importance" of each data source. However, these matrices are derived from normal equations that are themselves influenced by the chosen parameterization (e.g., rate-term truncation at degree 60), and VCE weighting. They reflect the modeler's decisions as much as the data's inherent information content. The claim that GRACE-FO LRI data "provides the principal constraint on zonal coefficients" may be an artifact of the specific processing choices rather than a fundamental property of the data. A more cautious interpretation would strengthen the manuscript by acknowledging that these contributions are conditional on the processing framework.

 Minor Comments

 Minor Comment on Section 3.1.2

In the second paragraph of Section 3.1.2 (Page 16), the authors state that they attribute the limited improvement from GRACE-FO data to "two main factors," but then proceed to list three factors (numbered 1, 2, and 3). This inconsistency should be corrected. Please revise either the introductory phrasing to read "three main factors" or adjust the numbering to match the stated count.

 Minor Comment on Section Numbering

In Section 3.1, the subsection titled "Combination with KBR1B data" is numbered as 3.2.2. This appears to be a numbering error, as it should be 3.1.2 to be consistent with the hierarchical structure. Please correct this numbering.

 Comment on Figures with Multiple Curves (e.g., Figures 4, 5, 6, 8, 9, 12, 15)

The authors present multiple curves in these figures to compare different models or processing strategies. However, the current color scheme, particularly the use of similar shades of blue and green, makes it difficult to distinguish between curves. This issue is compounded by the fact that many of these curves exhibit only small differences, which are already challenging to visualize. To improve readability and allow readers to better appreciate both the similarities and subtle distinctions among the compared solutions, it is recommended to use a more distinguishable color palette (e.g., colorblind-friendly options from ColorBrewer) and/or incorporate distinct line styles (solid, dashed, dotted). These adjustments will not only enhance clarity but also help convey that certain models are nearly identical in performance.

The authors describe HUST-Grace2025s, a long-term gravity field model based on GRACE and GRACE-FO data, including GRACE-FO LRI. The model is available in a standard file format and hosted on ICGEM, the central hub for (static) gravity field models. 

The manuscript is well structured and starts with an overview on current comparable gravity field models, presents the methods used to derive the gravity field model, then describes the model characterics and follows with internal and external evluation and validation. Some figures can be made clearer so the information can be parsed easier. In my opinion, as this is a data description article, some parts of the methods and results section can be shortened while the model evaluation and validation can be expanded so that potential users get a clear picture of the model performance and characteristics.  There are a number of experiments on certain aspects on the processing chain in the paper. It should be clear to the user what processing choices the final, published model comprises - here a concise description would help. The authors also do not present a clear use case of their model (i.e. why should we use their model for a specific application compared to already established data products). There are also some textual inaccurcies that need to addressed to provide proper context.  

From a methodological point of view, I have two major issues which need to be discussed:

- The first one is the low maximum-degree of 60 for the co-estimated time variable coefficients. The authors even state that there is information in the observation data beyond degree 60. Truncating at such a low degree applies a hard constraint on coefficients above degree 60 which are subsequently fixed to the used background models. This introduces spatial aliasing and makes the model dependent on the background fields. I would argue that such an implicit regularization needs to be communicated more transparently.

- Second, since the empirical covariance model has a major impact on the gravity field solution in the authors' approach, a more precise treatment of this processing chain component is crucial. Empirically derived (post-fit) covariance models depend on the parametrization, which should then obviously be consistent throught the processing chain. Especially time-dependent parameters such as N-CPR drastically affect the shape of the (stationary in time) covariance function. I realize that in practice some approximations are necessary, however, it would greatly improve the confidence in the model if this is presented in a strict and traceable way. 

Next to these points, please find below some detailed comments on parts of the manuscript: 

P2, L28/L29: GOCO25 -> GOCO2025s

P4, Table1: EOT11a is quite an old tide model, have you evaluated more recent ones (e.g. FES2022, GOT 5.6, EOT20) and do you use admittance functions? I suspect that a modern tide model will have a large impact. Since you use 4 second LRI data: did you perform any low-pass filtering before down sampling? The maximum degree 60 for the parametrized temporal variations seems quite low. Have you conducted a post-fit residual analysis to see whether residual time-variable signal is present? Osulation -> Oscillation (also missing period/frequency); I guess the kinematic orbit reference should be https://doi.org/10.1016/j.asr.2022.03.014 as it is in the Data Availability section

P5, L22/L23: How are the sensitivities computed? Please check the math notation: https://www.earth-system-science-data.net/submission.html#math

P6, Figure1: y-label: Orbital altitude (singular); also please either align the y-limits or even better put GRACE/GRACE-FO into a single panel so the altitude can be directly compared.

P7: Please check the math notation: https://www.earth-system-science-data.net/submission.html#math

P8, L3 - L16 1) computational efficiency is very implementation specific:

I would not formulate such statements as general fact, but relate them to your approach/implementation. 

2) post-fit residuals depend on the parametrization. You use a time-varible gravity field solution up to d/o 120 to determine the post-fit residuals. Your trend/annual parametrization goes up to d/o 60 beyond where time variable gravity field signal is certainly present. This mean that your stochastic model will be too optimistic and certainly not an optimal fit.

3) Using a correlation length of one month: A correlation length of one month is equivalent to a frequency spacing of ~4e-7. Is there a reason why such a frequency spacing is necessary? Shorter correlation lengths (coarser frequency spacing) and averaging within a certain time period (e.g. one month, one year, ...) increases the redundancy of the autocorrelation estimates and is certainly preferable.

P10, Figure3: Please make the correlation matrices square and use a divergent colormap with symmetric axis limits [-1, 1]. It would also help to have the lag in a time unit and a horizontal line at 0.

P11: Table 2: I find it hard to parse the numbers in the table. If you want to present the relative improvement, I suggest adding a panel to Figure 4 where this can be plotted for each degree and not only selected ones.

P11, L10: The term "Rate-terms" is somewhat misleading as you not only parametrize a rate but also an oscillation. I suggest something along the lines of "time-variable".

P11, L11: Calling a gravity field representation with trend and oscillation "static" is a little bit confusing, since purely "static" (i.e. only a long-term mean field) models also exist. This sentence should be rephrased.

P12: please check the math notation

P12, L6: Most (all?) recent models contain time-variable coefficients - please rephrase

P12, L9/L10: see comment above. computational burden is subjective since the parameter space does not expand drastically

P12, L10-L12: in the ITSG and GOCO models, all parameters are solved in one single step - not multi-tiered

P12, L12-L14: this sounds as though this is the representation of the GOCO models, which is not correct

P12, L20 - L23: as stated above the empirical covariance matrix depends on the parametrization. If you introduce or change parameters the model no longer fits, especially with time-domain parameters such as N-CPR

This is critical since the impact of the covariance model is substational (cf. Fig 4)

P13, L1: Using centimeters would make this more readable.

P13 L10 - L13: Since you use an empirically derived stoachastic model, AOD errors are also accounted for in your autocorrelation function. The issue is that a stationary (in time) covariance model cannot fully capture the AOD error (which are more stationary in space, see Kvas et al. 2019). 

Note that the same is true for ocean tide errors (Abrykosov et al. 2023)

It would be good to discuss this either in the section on the autocorrelation determination or here.

P13, L15: since "rate term" is both used in KBR parametrization and time-variable gravity paramtrization, please introduce a distinct nomenclature (e.g., "time-variable" for gravity field parametrization as mentioned above)

P13 L22 - L23: 1) this is not correct (e.g. Kvas et al. 2019 and Abrykosov et al. 2023 use model uncertainty estimates) 2) the four presented methods work differently

 

P14, Figure 6: Please state the reference gravity field in the caption.

P15 L18: One issue I see here is that KBR and LRI (also kinematic orbits) contain common noise sources (accelerometer and background models). VCE in this conditions will not give you an optimal solution. This should be discussed.

P16 L13-L14: it may also lead to degradation due to spatial aliasing

P16, FIgure 8: Combine -> Combined

P18, Figure 10: since the "static-only" solution is not the final published model, this section can be condensed. I would rather see an analysis of the final published model characteristics (e.g., a contribution plot for the trend/annual coefficients)

P18, Figure 11: It would be interesting to see the contribution to the time-variable components (trend, annual). Also: a diverging colormap is not ideal for a sequential quantity, please adapt.

P20, Figure 12: Please add the reference model to the figure caption.

P21 L8: please specify extactly what external information is used

P21, L9: Here it is worth to note that truncating the time-variable coefficients at 60 applies a hard constraint to coefficients above 60 which are then fixed to the background model (which makes your model dependent on the background model quality). This gives used a clearer picture of what is happening here.

P21, L14 - L15: See comment above: truncating at an (in my oppinion) too low degree constrains all coefficients to the background model.

P22, FIgure 13: Please do not use the Jet colormap (https://www.earth-system-science-data.net/submission.html#figurestables)

P23, Figure 14: using a circular colormap is not ideal here. since these are differences with a well-defined mean, a diverging colormap is more suitable.

P26, L24: the static coefficients of ITSG-Grace2018s are not regularized, please be more specific

P27 Data Availability Statement: some links do not work, please check