The preprint under consideration presents the free-surface stabilization scheme (FSSA), a stabilization scheme for the free-surface Stokes equation. The free surface Stokes equation is a time-dependent equation that is generally solved with explicit time-stepping schemes, due to the difficulty in solving the free-surface evolution and the Stokes equations simultaneoulsy. A consequence of using explicit schemes is that time steps have to be small. The FSSA has been shown to enable higher time steps without reducing the numerical accuracy signficantly. In this preprint, the authors study the impact of using the FSSA in the context of marine ice sheets, where the free surface detaches from the bedrock at the grounding line. 

The preprint is well written and contains two well-documented numerical tests that demonstrate the power of the FSSA in the context of marine ice sheets. My main criticism of this paper is that it is very similar to other existing preprints and publications, including 

- Josefin Ahlkrona, A. Clara J. Henry, and André Löfgren (2025), A Fully Implicit Second Order Method for Viscous Free Surface Stokes Flow – Application to Glacier Simulations (under consideration for publication in GMD),

- Tilda Westling Dolling, A. Clara J. Henry, and Josefin Ahlkrona (2025), A numerical stabilization scheme for the shallow shelf approximation (arXiv:2510.02943).

Since the need for stabilization schemes like this one is urgent in the glaciological community, I recommend this article for publication once the following comments are addressed. 

General comments

- As indicated above, this paper is closely related to a series of papers on the FSSA (and its depth-integrated variants) written by some of the authors. I think the introduction should contain a more exhaustive description of what exactly is novel about this publication (its application to grounding line dynamics) and in what way it is related to existing publications. 

- Section 4. For the MISMIP experiment, you only solve the unstabilized system ($\theta = 0$) for one timestep. Does the scheme work for higher timesteps if $\theta = 0$. If so, it would be important to include these results so that the reader understands the effects of not including the FSSA. If the scheme breaks down for higher timesteps when $\theta = 0$, please indicate this. 

- The numbering of figures throughout the text seems to be incorrect. See for example reference to "Figure 4.2" in Line 281 and "Fig. 5.1" in Line 300. 

Minor comments

- Line 3. "The restrictive time step size of ice-sheet simulations (∼ ∆t = 0.01−0.5 years) has led to the routine use of approximate models that compromise physical complexity compared to full Stokes models." In my opinion, approximations of the Stokes equations are used in glaciology because the Stokes equations, as opposed to e.g. depth-integrated models, are very expensive to solve. Please rewrite along the lines of ""The restrictive time step size of ice-sheet simulations (∼ ∆t = 0.01−0.5 years) is one of the reasons why... "

- Figure 1. This figure is not referred to in the text. Please refer to to it when introducing the notation for the domain $\Omega$ and its boundary. 

- Line 81. Clarify that, despite the Stokes equations appearing to be stationary in time, time-dependence enters the problem because the domain Ω changes in time.   

- Lines 105 and 110. Writing $v \in \mathbb{R}^3$ and $q \in \mathbb{R}$ is very confusing. Please rewrite as e.g. $v: \Omega \to \mathbb{R}^3$ and $q:\Omega \to \RR$.

- Figure 2. Please refer to this figure in the text. 

- Line 137. "However, due to the elliptical nature of the Stokes equations, the coupled system of the Stokes and the free-surface equations cannot be fully implicit in the velocity." I do not understand this sentence. Are you not solving this coupled system implicity in the paper "A Fully Implicit Second Order Method for Viscous Free Surface Stokes Flow - Application to Glacier Simulations"? Please clarify and explain relationship to this paper. 

- Line 185. What does the velocity with a tilde, $\tilde{u}$, refer to? It seems to be like this variable has not been defined. 

- Line 221. Write "Table 1" as opposed to "Table (1)". Same for "Table (2)" in the next line.  

 - Line 245. Why are Picard iterations used? Is a Newton solver not implemented in Elmer/Ice? The p-Stokes equations are the optimality conditions for differentiable convex optimization problem and therefore work very well, in general, when Newton's method is applied.

- Line 308. "The isothermal ice fluidity, A, and the friction coefficient, C, were tuned..." Please indicate how they were tuned. Did you solve an inverse problem to find the corresponding scalar values?

- Equation (28). Is this choice of boundary condition for the ice divide common? Please reference other works where this is used. What value is chosen for $C_n$? I can imagine that, as $C_n$ becomes very large, an impenetrability boundary condition is effectively enforced (this is equivalent to penalizing u.n at the boundary). Is this true? If so, please indicate this in the text to facilitate understanding for the reader. 

- Line 343. "Figs 5.1 and 5.1". Please correct.

- Line 363. "a time-step size of ∆ = 50 years". Please correct.

- Line 382. "We hypothesis...". Please write verb correctly.

The manuscript “Grounding-line dynamics in a Stokes ice-flow model (Elmer/Ice v9.0): Improved numerical stability allows larger time steps” by Henry et al. extends a previous method developed by Löfgren et al. (2022) for stabilizing ice-sheet free surfaces (FSSA) to the case of a marine-terminating ice sheet. The FSSA method uses the Reynolds transport theorem to approximate the gravitational body force on the domain at the next time step, introducing new terms to the weak form of the Stokes equations. This stabilization complements the previous “sea-spring” method of Durand et al. (2009) that applies similar implicit approximations to the water pressure boundary condition. The method is clearly described and motivated, and will make a substantial contribution to computational glaciology by providing a roadmap for implementation in other ice sheet models beyond Elmer/Ice. However, the presentation requires improvement prior to publication, and I have several specific comments and suggestions that should be considered.

1. Line 16: Specify that the model domain is the Ekström Ice Shelf.

2. Line 17: In the abstract, I think you need some caveats about slow retreat/advance, because a time step of ten years will not be universally applicable to more rapidly evolving systems driven by ocean melting (e.g., Thwaites). 

3. Line 28: Regarding that statement “grounding line can be subject to MISI”, consider citing some of the recent papers by Sergienko et al.:

a. Sergienko, O., & Haseloff, M. (2023). ‘Stable’and ‘unstable’ are not useful descriptions of marine ice sheets in the Earth's climate system. Journal of Glaciology, 69(277), 1483-1499.

b. Sergienko, O. V. (2022). No general stability conditions for marine ice-sheet grounding lines in the presence of feedbacks. Nature Communications, 13(1), 2265.

4. Line 39: “Only a small number of models support full Stokes” Please specify which models, for example I know that ISSM supports full Stokes. 

5. Line 49: I would rephrase this is as “extend FSSA to the ice-ocean interface” rather than “introduce FSSA at the ice-ocean interface”.

6. Line 56: Suggest motivating the choice for Ekström Ice Shelf here, as it is a very particular choice.

7. Line 114: It is not clear whether the iteration j is referring to the time step or nonlinear solve step. In either case, how is the coefficient initialized?

8. Line 130: A description of the algorithm for grounding line migration is needed, especially because the in-text description in Durand et al. (2009) is unclear. I am suggesting using standard algorithm notation to say exactly how the contact conditions are resolved during the nonlinear iterations. This could help illuminate how the stabilization scheme(s) influence grounding line migration.

9. Line 142 and elsewhere: I think it is worth emphasizing that the sea-spring and FSSA really complement each other by introducing implicit approximations of the gravity-driven forcings. That is, it makes sense to use both methods so that the forcings are consistently approximated.

10. Line 168: Suggest providing some intuition behind why Reynolds transport is only applied to the gravitational terms, following some of the descriptions from Löfgren et al. (2022) , e.g., “This makes the FS equations ‘aware’ of the spatially evolving domain by estimating the impact of the force of gravity at the next time step.”

11. Line 176: Are the mass balance terms a_s and a_b evaluated at timestep k or k+1?

12. Line 179: remove “on” and reference previous weak form equation (13)

13. Line 182: \Gamma_{c,z<0} is confusing notation. Why not just modify the forcing term to be zero outside the submerged portion?

14. Line 182: You say the calving front does not move, which is a simplification. This should be mentioned earlier like in section 2.2. Somewhere you need to describe if it is possible to extend the method to the general case where a calving law is implemented, because modelers will want that capability in the future.  

15. Equation 21: Why does the u have a tilde and k+1 here?

16. Line 187: It seems like these details about the mass matrix are not used anywhere, so I suggest omitting them.

17. Line 195: In this paragraph, it seems like you are describing that the FSSA and sea-spring balance out in a way, in line with my previous comment about how they complement each other.

18. Line 205: “size of the largest stable time step size” redundant size

19. Line 235: The mesh is updated so that the refined area tracks the grounding line position. Specific details about how this is implemented would be helpful.

20. Line 250: This paragraph described a lot of numerical details about how you are solving the free surface equations. Consider writing down the discretization details to show exactly what you are solving, as you did with the Stokes problem.

21. Figure 4: There are five colors and three linestyles in this figure. It is impossible to parse the results with all the closely spaced and overlapping colors and textures. Please split these line plots into multiple panels for the different timestep sizes, with the reference solution plotted in each panel.

22. It’s never explained why the referenced grounding line solution in Figure 4 has a mesh resolution of 200m. I found this confusing because Figure 4b and the results section seems to use the reference solution to measure the “accuracy” (Line 264) of the stabilization methods, but it’s not clear to me why this should be the case. Wouldn’t a higher resolution simulation be better to measure accuracy?

23. Figure 5: Same comment about too many colors and linestyles. I cannot easily interpret this figure. Also, in the caption, specify that theta=0 is the black line.

24. Line 286: Maybe this is mentioned in the Discussion (Line 383), but how does FSSA consistently having a further oceanward grounding line position relate to the previous discussion of changes in the predicted ice load or water pressure (or thickness) in advancing or retreating cases?

25. Figure 7 is missing units in the axis labels.

26. Line 296: This sentence is confusing because it’s unclear if “respectively” is referring to the timestep sizes or spatial resolution.

27. Line 308: How exactly were the friction coefficient and ice fluidity tuned?

28. Line 323: Is this Neumann condition like a weak penalty condition? Do these numerical details belong in Section 2 or Section 3?

29. Line 340: Convergence tolerances have already been stated previously in Section 4, so I suggest removing unless they are really necessary here.

30. Section 5: The evolution of the grounding line position is not shown for the 3D results. I know this might be more tedious for the 3D simulation, but it seems necessary as the paper is about grounding line dynamics. 

31. Line 347: How does buttressing relate to the thickness statistics in the previous sentences? I don’t see the connection 

32. Line 349: clarify what “good convergence” means here. Converged in a small number of iterations?

33. Figure 9: Include horizontal ticks on the upper panels

34. Line 363: Delta is missing a t

35. General suggestion to discuss applicability to shorter time scales, e.g., fast evolving or vulnerable systems like Thwaites, and 100-yr projections. The discussion on Lines 393-396 needs to appear or be mentioned much earlier like in the introduction.

36. Line 370: What does “correction of datasets” mean? Are you referring to bed topography?

37. Line 382: Change “hypothesis” to “hypothesize”

38. The attached code repository needs a readme that describes the different files, dependencies, how to run the Elmer code, etc.....