Reviewer 1:

This manuscript describes a new capability for evaluating viscoelastic Love Numbers that has been built within the ISSM coupled ice sheet modeling framework. The new functionality allows to model the solid Earth response to surface loads up to very short wavelengths, and supports different rheological laws including the EBM transient model. The manuscript is very well written, it is technically sound and discusses in detail a comprehensive set of benchmarks with community reference solutions for incompressible LNs in order to validate the numerical approach. The authors also devoted considerable effort in improving the code efficiency, in view of its possible applications to ensemble modeling or Bayesian inversion frameworks, hence I believe that this new tool will be of great relevance for the GIA community. I have just one major point that in my opinion should be addressed before recommending publication, which is illustrated below.

Specific comments

My main (and only) concern about the manuscript is about the stability of the approach discussed in the paper in the case of compressible models. A compressible self-gravitating model is known to be affected by gravitational instabilities both in the elastic and viscoelastic regimes (see e.g. Hanyk et al., 1999 or Vermeersen and Mitrovica, 2000) which, from a mathematical point of view, manifest themselves as singularities in the Laplace-transformed solution on the real positive axis. Within the classic framework of viscoelastic normal modes those correspond to exponentially growing modes, which are usually dropped since their timescale is generally much larger than those of interest for GIA, and only stable modes associated with singularities on the real negative axis are kept to evaluate the time-domain solution. However, if the time-domain LNs are computed according to eq (35), there is effectively no way to filter out those modes (since the individual modes are not identified at all, which is one the main advantages of the proposed method), and therefore it is possible that the Laplace solution is sampled near a singular point, which would break the inversion scheme. I do not have any specific suggestion about how this problem can be circumvented, but I think that it is an important aspect that would merit some discussion. 

[Authors’ response]: We thank the reviewer for bringing up this topic with relevant references. We agree that this is an important point and have added a new paragraph (see section 3.5) to discuss this potential downside of the Post-Widder transform. Our discussion builds on Hanyk et al and Vermeersen and Mitrovica to suggest mitigation strategies for the user.

Minor points:

Page 9, equation 17: For a tidal forcing, the gravity potential is computed without considering the direct term (i.e. on the r.h.s. ‘k’ is used instead of ‘1+k’). Is it to be compliant with the GRACE convention for the definition of potential, in which the direct term is not included? If this is the case, I suggest making this point explicit.

[Authors' response]: This is not related to GRACE, but to the definition of the system. We assume that the geoid is based on the gravity potential of the Earth, reflecting both changes of mass in the interior and at the surface (hence why 1 is included in the surface loading case). Since tidal forcing originates from the Moon, Sun or other celestial bodies (except for the pole tide case), including the direct term would reflect the gravity field of the combined system instead of the geoid. For the rotational feedback (aka the pole tide) and for GRACE applications, we have a dedicated process that uses the PMTF and degree-2 tidal love numbers to compute the resulting deformation. In the end, the specifics depends on the intended application by the user, the scope of this paper is to provide them with the values of k. We have revised the paragraph following eq 17 to explain this choice.

Page 17, lines 335-340: the “pseudo-fluid Love numbers” and the “degree n_vmax where the viscous part of the LNs becomes negligible” are first mentioned here, but the meaning of those concepts is made explicit only in section 5 below. Maybe it could be beneficial to the reader if a few more words are given here about those two points, anticipating that they will be discussed in more detail in section 5.

[Authors' response]: We have added a few sentences in section 4 as suggested, and referenced the appropriate section for more information.

Technical corrections

Page 3, line 75: Missing period after the reference to Kierulf et al.

Page 4, line 92: maybe βtⁿ should read βt^1/n, or n ≅ 3 should be n ≅ 1/3 ?

Page 7, line 122: “an mantle” -> “a mantle” (or “a solid mantle”?)

Page 12, line 123: my understanding is that rⁿ should be r^-n 

Page 15, line 290: “This methods” -> “This method”

Page 21, line 357: “theses” -> “these”

Page 25, lines 439-440: “point to with” -> “point to”

Please note that some entries in the bibliography are missing the DOI.

[Authors' response]: All of these points have been corrected. We thank the reviewer for a thorough examination of the manuscript.

Reviewer 2:

Review: Love number computation within the Ice-sheet and Sea-level System Model (ISSM V4.24) – Caron et al.

Summary:

The manuscript of Caron et al. presents the theory for computing Love numbers for Maxwell, Burgers, and Extended Burgers rheologies. In this discussion they present a general overview of the parallelization strategy and demonstrate their optimization strategy. Furthermore, they compare their Love numbers to those of Spada et al. (2011), a benchmarking paper in the GIA community, and find good agreement for a Maxwell rheology. This body of work and the associated code is an important contribution to GIA modeling and coupled GIA/ice sheet modeling.

The code may also have utility for other deformation problems given suitable modifications. Overall, this manuscript is well written and in fairly good shape. Below I outline a few major items and several minor items that would help improve the quality of the manuscript. None of these should take very long to implement and thus, I suggest this manuscript be accepted following minor revisions.

[Authors' response]: We thank the reviewer for their thorough examination of both the content and presentation of the manuscript. Their comments have significantly improved the clarity and quality of the paper.

Major requests:

• Abstract: Please add the something about the Love number approach being for radial 1D (i.e., spherical symmetric earth models) so the limitations of this approach is more clear to the general audience. You could also add a sentence on the benefits of this approach (e.g., that it is computationally quick compared to codes with 3D viscosity structures) to balance out the negative.

[Authors' response]: We have clarified the abstract as requested.

• The dot used to show the divergence is at the bottom “.” as opposed to the middle “∙”. Please fix this in equation 2, 18. The dot product in equation 12 is at the bottom “.” and should be in the middle “∙”

[Authors' response]: Corrected.

• Please make sure all variables and super-/subscripts are defined within the text. Table 1 is a helpful guide, but it doesn’t contain everything and it makes for an unpleasant read to hunt for variable definitions that could easily be stated in the text. For example: eq 11-14: SH degree n and SH order m are not explained nor are the bounds for their summations; eq 19 – what is E; eq 29 – explain 𝜌# in text please; etc

[Authors' response]: We have modified the text throughout the manuscript to label variables as they are introduced.

• Line 166-171: There is a disagreement in the ordering of the text and yi in comparison to the equations 11-14. For a savvy reader this will not be a problem, but for others this could be confusing. I suggest first presenting equations 11, 12, 13, and 14 in order. So, introduce the full displacement vector field, the full traction vector field, etc. Then introduce “yi” and state the order of the vector elements following the convention of Alterman et al.

[Authors' response]: We moved this sentence after the y_i equations, as suggested, and clarified explicitly which one related to each physical variable therein.

• Equation 17 – Please state what N(\theta,\phi,t) is (i.e., how is it different from N^{m}_{n}, the

normalization factor or the number of layers N_{L}) and also consider a different variable to avoid

confusion

[Authors' response]: Noted, we changed the number of layers to N_L to \mathcal{L}, the spherical harmonics Normalization N^{m}_{n} to \mathcal{N}, and kept N for the notation of the geoid, as standard in the field.

• Is there a standard you followed for using ~, ‘, and f_{i} for the different factors in 38. It would be helpful if there was and if you could please explain it.

[Authors' response]: We have added a sentence (section 5.1.1) to make our approach to notations explicit: $f_{1..11}$ refer to the precomputed terms in order in which the appear in eq. (39), terms with a prime symbol $'$ are used for non-dimensionalized quantities (e.g. $\mu'$ is the non-dimensionalized shear modulus) while terms with a tilde $~$ are used for other scaling constants.

• Figure 7: Please add a bounding box to left figure to show the location of the right figure.

[Authors' response]: Done, this is also referenced in the caption now.

• Figure 8: Please add something to the caption describing where the h love numbers are plotting.

[Authors' response]: We have amended the caption to indicate that the curve for h tends to be identical to that of k (red line).

• Multi-panel figures: Please choose a consistent method for referencing multi panel figures. For example, Figure 8 uses (a) and (b), but none of the other figures use letters.

[Authors' response]: We have removed the letters from Figure 8 to make all figures consistent.

• Multi-panel Figures: Please ensure there is a description of each panel in the caption. For example, Figure 2 does not explain both plots. Second please modify the captions so that it is easier to find where the descriptions of each panel is. Bold and italic font might be helpful, particularly when using “Right”, “Top Row”, etc.

[Authors' response]: We revised most of the captions in this way, highlighting text pointing to rows and columns in bold. For figure 1 and 3, there was little to add to the information provided by the figure frame titles so we opted not to.

• Section 5.2: As described, it is unclear how experiment 1 & 2 are changing the mesh discretization. Both are described as “increasing the radial grid resolution ….”. The term layer is also used, but in this context, it is unclear if the authors are referring to the layers shown in Figure 1 or if they mean a spherical shell. A figure showing the change to the mesh for each experiment would be helpful.

[Authors' response]: We mean neither the layers themselves nor a spherical shell. For the sake of performing the propagation of the yi system within each layer (from bottom to top boundary), we perform the iterative procedure explained in section 3.3, updating the system after each radial step Δr. This subdivision creates a 1D grid from r_{inner} to the surface. The question asked in comparing experiments 1 and 2 is twofold: 1) for a given number of grid points (equivalent resources), is it better that Δr be uniform across layers, or that each layer have the same number of propagation steps? 2) how fine should this grid be, in order to solve for surface Love numbers accurately, before yielding diminishing returns?

We have added a short paragraph at the start of section 5.2 to better explain what get across these points.

Minor requests:

• Line 28: “… model disparities related to solid Earth structure and mantle Rheologies.” Please add “1D” to the solid Earth structure to be clearer.

[Authors' response]: Added “radial”.

• Line 75: missing a period at the end of sentence “ … Kierulf et al., 2022).”

[Authors' response]: Corrected.

• Line 79: Table 4 is introduced before Table 1, 2, and 3. Consider revising the text of table locations/labels so things appear in the order in which they are discussed.

[Authors' response]: Done.

• Line 83: Since you are discussing changes in Earth’s shape would \dot{J}_2 be more appropriate than just J_2?

[Authors' response]: Corrected.

• Line 103: “… Heaviside function forcing …” -> “… Heaviside forcing function …”, maybe easier to read?

[Authors' response]: Agreed, changed.

• Line 161: “… and material …” -> “… and then material …”

[Authors' response]: Changed to “and the material isotropic.”.

• Line 205: What is meant by “resolved”?

[Authors' response]: Changed to “solved”.

• Equations 24, 25, 26: The 0 on the right-hand side should be written as a vector and not a scalar.

[Authors' response]: Fixed.

• Table 3: Please add units to your times. Currently one must assume they are in seconds like the precision you report in the table caption.

[Authors' response]: Changed “Computation time (s)” to “Computation time (seconds)” at the beginning of the caption to make this unit harder to miss. We would prefer to not include it in each of the table cells in order to avoid clutter.

• Line 440: delete “with”

[Authors' response]: Done.