The paper by Yuan et al. presents the next step of the GIA modelling setup using CitcomS. While the method has been published previously (SLE solver in Zhong et al. and compressibility in A et al.), here they present a benchmark of the two combinations. The paper summarizes the applied methods very well. However, I felt sometimes a bit lost in the paper as you introduce a method/set-up and then something else is mentioned afterwards. You come back to the method/set-up at a later point. I would suggest splitting the paper in two parts: 1) introducing and benchmarking your model to simple harmonic loads (where no SLE is involved), and 2) the benchmark with the SLE and the ICE-6G_D ice model, which means that you introduce the SLE part after you have presented the first benchmark.

A few other things could be done to convince the reader that your GIA modelling setup works:

• Show a comparison of the geoid change in the benchmark with the SLE and ICE-6G_D.
• Compare your radial and horizontal displacement rates as well as geoid and RSL rates to the results by Peltier et al. (which you can find on his website).

Minor comments:

L.73 - Huang et al. is not a coupled spectral-finite element method as it is purely based on the finite-element method. In addition, the method by Wu is also based on the finite-element method and only the potential calculations are done in the spherical harmonic domain. Not sure, if this should be then called coupled spectral-FE method. If I understand you correctly, you are doing the same by "coupling" the spherical harmonic domain with the FE method.

L.77 - "center of mass".

L.144 - You mention that the Lamé parameters and viscosity can vary laterally and radially. Is this also true for the density? I couldn't find such a statement.

L.183 - "the applied potential, which is only relevant"

L.185 - "on the unknown incremental"

L.195-202 - Please move this paragraph somewhere else. It doesn't fit here.

L.212 - What "given threshold error tolerance" is implemented in CitcomSVE 2.1 & 3.0?

L.301 - Interesting approach. Please write it out in an equation how O(t) is calculated.

L.324 - Why do use the density and Lamé parameters from the mantle for the crust?

Section 3.2.1 - You refer only to figures in the supplement. Either move the text to the supplement and summarize the results in the main article or move figures from the supplement to the main manuscript. I would prefer the former.

Section 3.2.1 - Just out of curisosity: Could you please plot K2 in Fig. S1 as well? So, the results of the 2nd iteration based on the approach by Kendall et al. No need to do this for the other figures.

L.419 - In the abstract you say ICE-6G_D, but here it is ICE-6G now. Please use the same and correct abbreviations throughout the text.

L.546 - Please use more representative names for your model types, e.g. GIA_R1 could be GIA_135 to indicate that this model has a resolution of 135 km and the reader doesn't have to go back and check what R1 was referring to.

Table 3 - I find it difficult to compare the CPU hours when you use different amounts of cores. Please unify this to see the direct difference.

Fig.S2/S3 - Use the same abbreviations for both figures.

Figures – Copernicus recommends the authors to use colorblind-friendly colors/colormaps, which is stated on the GMD website (https://www.geoscientific-model-development.net/submission.html#figurestables). This includes avoiding the parallel usage of red and green as well as the avoiding jet and rainbow colormaps. The link above provides examples where you can find suitable colormaps and find colors for your line plots that are suitable.