The article presents an ad-hoc-stabilized scheme for implicit time-evolution of the viscous free-surface Stokes equations. The idea relies on a relaxed explicit Euler discretization of the transport theorem that ties the discrete solution at two consecutive time steps. The main technical concern pertains to the introduction of a relaxation parameter controlling the size of the time-step appearing in the regularization. Its role is not discussed in full detail.

Since this parameter effectively reduces the time step size in the stabilization (but only in the stabilization and not in the general time-marching scheme)  it introduces a mismatch in the timescales (and time grids) involved in the discretization.  Thus, if the value of this parameter is not one, then it may introduce an inconsistency and lead to loss of convergence. The inconsistency, however, may only become dominant for much smaller values of the time step than the ones  feasible in the article's simulations. The information provided in the text is not nearly enough to determine if this indeed would be the case, or if some additional measures are taken to avoid the issue. As a side comment, this "reduced time step" could potentially be interpreted as an intermediate stage (as in a Runge-Kutta method). This technique is relatively common for stiff problems, but the numerical solutions at these intermediate steps should not be given any physical meaning (just like the stages of an R-K method are not). In that sense the value of the parameter theta could be taken to be of the form 1/n where n is the number of intermediate stages to be discarded, interpreting only every "n-th" step as a physical solution. Unfortunately, the manuscript does not discuss this parameter with enough details to determine if this interpretation or idea is correct.

Nevertheless, the authors provide numerical evidence and sufficient comparisons that seem to indicate stability and reliability of the technique in some realistic situations of simplified (2D) glaciar flows. 

The article is in general well written and, modulo the aforementioned discussion of the relaxation parameter, and the correction of several details scattered through the manuscript (see files attached), the referee can recommend the article for publication.

Annotations on the manuscript (see the attached file)

1.- limit

2.- Illustrate explicitly in the figure the domain \omega

3.- either "the boundary terms vanish" or "the boundary term vanishes"

4.- Illustrate the region \omega explicitly in figure 1

5.- "At the [free or top] surface... or even better, "on \Gamma_s" 

6.- Indicate explicitly "on \Gamma_s" in equation (5)

7.- "on"

8.- Indicate explicitly "on \Gamma_b" in equation (8)

9.- Are all the curly brackets really necessary?

10.- This expression and several related ones later on seem to be missing the terms corresponding to the y coordinate. In the numerical experiments at the end the problem is considered in only 2D, but the rest of the discussion pertains to the full 3D problem, and so the "y-terms" must be included or an explicit comment should be made at this point regarding why only a 2D slice is considered. Compare for example with eq (14), where the "y-terms" do appear and then eq. (15) when they disappear back.

11.- See previous remark

12.- See remark 10

13.- W^{1,2}

14.- The line numbers in the algorithm are missing

15.- This statement is purely speculative. Either provide a proof/stringer argument or supress the sentence.

16.- do

17.- "... on weak form, recalling that $\partial_t(\rho\mathbf g \cdot v) = 0$, which leads to"

18.- The first equality seems unnecessary

19.- The stabilization seems to consist of connecting the terms at times k and k+1 through an explicit Euler discretization. In particular, the substitution of eq (17) into the right hand side of (1). The Euler discretization comes with the addition of a factor $\theta$ controlling the size of the time step. This is referred to as "user defined parameter". However, unless theta is identically equal to one, this would result on an inconsistent discretization and therefore would lead to failure of convergence. The role of the parameter needs to be discussed in more detail, taking into account the loss of consistency that it introduces. 

20.- Just as in comment 15, this statement is purely speculative and it must be supported strongly or dismissed. In what specific sense do these solutions "mimic" the ones arising from an implicit discretization? Rather than this speculative statement it could be said that  "numerical experiments seem to indicate that this results in more stable discretizations" simply pointing to the experimental evidence without speculating.

21.- This line belongs in the previous paragraph

22.- Line numbers missing

23.- This expression is missing the "y-terms" (see 10)

24.- indeed

25.- missing "y-terms" (see 10)

26.- This is not a distance, but it is a difference

27.- better factoring out $\Delta t$

28.- "...side of (15) and the expression above..." 

29.- " , leading to: "

30.- this sentence should be right after the equation

31.- "is substituted into the Stokes equations"

32.- "A key observation is that, if the iteration converges,"

33.- this is not a distance, in this case "difference" can be used

34.- mark

35.- this is not a distance, but in this case "displacement" can be used

36.- two

37.- from

38.- "y-terms" missing, see (10)

39.- Is this the domain in Figure 3? If so, it should be mentioned

40.- This is not an error, unless you have an error estimator, it is simply the convergence tolerance. 

41.- This figure is not referenced explicitly in the text. It should.

42.- This is not measuring an error, so the term "accuracy" is incorrect. Do you mean that this is used to determine convergence of the iteration?

43.-  are

44.- this is not an error, it is the iteration tolerance

45.- Once again, this is not an error, unless there is an error estimator that is not mentioned anywhere

46.- "stabilization to be stable" is there such a thing as an unstable stabilization? Then it is no stabilization at all

47.- Do you mean "convergence of the second order method"?

48.- types

49.- remove comma

50.- Except for the first three points, all other errors are large to the point of making the numerical solution completely meaningless, so making a distinction between the intermediate and large step sizes seems a moot point

51.- Forward Euler

52.- This is not an error, it is simply the relative change between iterations

53.- Iteration tolerance

54.- convergence or iteration tolerance

55.- stabilization

56.- This point is quite important and is somewhat glossed over very lightly. It should be discussed in more detail, or at least references to a deeper discussion should be provided

57.- The notation is not clear. Is this a gaussian random variable and the first argument of octave( , ) is the mean and the second the standard deviation? Is it a different type of random variable? If so, what distribution does it follow and what do the arguments control?

58.- There is a double comma

59.- ms^{-2}

60.- convergence tolerance

61.- This is not an error, is simply relative change. An iteration can converge to the wrong solution, thus even id the relative change is 0 the error might be very large. Since there is no error estimator invovled, the term error should be avoided.

62.- scheme

63.- keep

64. - If I recall correctly there was an issue with the use of github. Should this be updated?

This is a very interesting and well-written manuscript that improves my earlier free surface stabilisation (FSSA) algorithm to allow for significantly larger timesteps in glacier and geodynamics codes that have a free surface. 

The authors demonstrate that a purely implicit timestepping algorithm doesn’t give much advantage compared to an explicit scheme with FSSA if a first-order timestepping scheme is employed.  This surprised me somewhat (I had thought otherwise), but is rigorously shown here. In addition, they show that using a second-order timestepping scheme can be made second-order by using a modified FSSA algorithm. 

This is initially described for a simple linear viscous setup, which is expanded to a more realistic ice flow simulation in section 7.

Overall, this is a very nice contribution that will be helpful for both the glaciology and geodynamics communities. I think it can be accepted with minor revisions, and I just have a few minor remarks:   

•  You seem to use Picard fixed-point iterations. If you have a good initial guess, a faster way to converge nonlinear iterations is Newton iterations, which should work particularly well for moderately nonlinear problems such as powerlaw viscous materials with an exponent of n=3. This would potentially also help the convergence if the implicit timesteps, even though it requires deriving the appropriate Jacobian. Is this not commonly done in the glaciology community, or are there other reasons you stick with Picard?
• eq. 29: Can you give a function or algorithm on how “octave_?” is computed? The way this is written is a bit puzzling to me, but this might be because I am not from glaciology...

Minor typos:

l. 299: kg m^{-3}

Fig. 8 caption: “stabilisation”