Dear authors,

both reviewers have now seen the revised version of your manuscript and concluded that only some minor revisions are still needed. Please consider the following comments when improving your manuscript.

Kind regards

Heiko Goelzer

Referee #1
I would like to thanks the author for the revised version of their paper that account for most of my remarks. I have few more comments that I think could still improve the manuscript.

First, I think the paper could be improved by adding a systematic comparison (quality of the solution, cpu) of the 3 presented stabilisation methods for all examples. It is mentioned that all are done with the Crank–Nicholson method (page 11 line 1), but it seems it is not the case (discussion in 6.2.2).

Second, the oscillation issue is not really shown in the main text (only in Fig 11 and then 12 and 13). Would it be possible to show on the other examples this issue with dedicated plots? Also, clearer conclusions should be drawn on which method(s) should be preferred?

Then I have few other comments (the page number and line number are those of the last version with highlighted changes in the author response):

Eq. (2): again I don't see why this equation is needed for what follow. At least, it should be given in term of stress to avoid an incomplete definition of the effective viscosity

page 4, line 11: Eq. (4) instead of (3)?

In Eq. (5): put v always before (or after) the initial terms

In Eq. (8): not sure ":=" is needed ("=" is enough). Define what is F in (8). Is it a scalar or a vector?

page 4, line 19: n \times n (not n x n)

In Eq. (9): is operator F the same as in (8) (so why in bold?)?

In Eq. (11): "with" should not be in maths

In Eq. (12): there is a missing integral sign?

After Eq. (13) should be mentioned how F_stab is added to the functional?

page 5, line 15: in terms "of" stabilisation

In Eq. (16) the notation should be given

page 8, line 27: S^1 = S^0 + \Delta t is not homogeneous. There is a velocity missing in front of \Delta t to have a displacement

7.1.1: there is no 7.1.2


Referee #2
-The authors have responded to my full satisfaction to all my earlier comments.
I still don’t understand how they solve the inequality problem, but form eq 10 is seems they might be using what I know as the ‘active method’. There is now discussion however on how the iteration over the active set conduced. Anyhow, I think since the paper now has another title this is not a bit issue anymore. And there are many different ways of solving such problems.

-The paper fits nicely to the scope of the journal and it’s great to be able to read and learn about various tactical modelling details that one would typically not get from papers on other journals.

-$n x n $ -> $n \times n$ on page 4.
-‘harder test’ on page 12, maybe replace with ‘more stringent’ test

-Just a comment: With respect to applying elevation-dependent mass balance (page 12 and there about). I’ve always found that it helps if this can be done implicitly. That is the da/dz term is taken over to the left-hand side and one solves for velocities and thickness at the same time. (I understand that in a full Stokes situation this is going to be more difficult than in the situations that I usually deal with which as the SSA equations.)

-I’m not sure this is really needed at this stage. But because you have so many tests to describe it might have been a good idea to list them in a able, of give them some more descriptive names.

-Line 15 page 19. Rewrite first sentence. Maybe: Occasionally, for very steep advancing glacier fronts, unphysical oscillations developed.

-Page 19: Are you sure the ‘wave’ is a kinematic wave? It’s almost impossible for kinematic waves to develop on glaciers unless mean slopes very steep. But maybe these are finite-amplitude effects. See for example, Vieli, G. J. M. C. L., & Gudmundsson, G. H. (2010). A numerical study of glacier advance over deforming till. Cryosphere, 4(3), 359–372. https://doi.org/10.5194/tc-4-359-2010.

-I was just wondering if your SUPG stability parameter tau=h/2u should possible have been selected as 1/tau=1/taut + 1/taus ; where taut =dt/2 and taus=h/2v. This gives the limit tau->0 with dt->0 as is should be. However, I suspect that this might not actually help and that you would only get rid of the oscillations using some sort of forward/positive-backward/negative diffusion method.

Dear Authors,

thank you for your careful and detailed work in response to the reviewers' comments. I have only one additional comment that you may want to address, if you agree. I find the font choice in Figure 1 and Figure A1 a bit unusual and not very readable. The font types in your other figures would be much better suited, in my opinion.

Your paper can now advance to the typesetting, proofing and publication process, congratulations. Thank you for your patience and apologies for the lengthy review process, which fell unfortunately into the summer holiday season and a period of an ongoing global health crisis.

Best greetings

Heiko Goelzer