In this paper the authors investigate the stability and performance of well known depth-integrated ice sheet models (SSA, SSA+SIA, L1L2, DIVA). The stability is studied theoretically on 1d problems with the assumption of uniform viscosity, and for explicit time integration schemes. The theoretical results are confirmed, for the most part, by numerical experiments. The authors also compare the stability and performance of these different models using the ice sheet codes CISM and Yelmo for modeling the Greenland ice sheet. Among these models, DIVA has better stability properties, and consequently performs better.



Overall the paper is well written and I think its contribution is very important. The theoretical and numerical stability results provided therein will certainly guide modelers and developers in choosing what depth-integrated models to use/implement.



I'm concerned about the discrepancy of their stability analysis and the results obtained using the L1L2 model implemented in CISM and Yelmo. I think it might have to do with the peculiar discretization chosen -see detailed review below- and I wonder whether they could get a better agreement if they changed the discretization used in the analysis.

I also suggest that the authors develop the theoretical analysis in an unbounded domain, which would make it simpler and cleaner.





Detailed review:

1. Title: I think this paper is about stability more than performance, so I would suggest changing it into "A comparison of the stability and performance of ..."
2. Abstract: please mention that the stability analysis is performed for an explicit time discretization scheme.
3. page 1, line 20. Saying that the solution of Stokes problem is still infeasible is too strong of a statement. It would require a lot of resources, large clusters and parallel scalable implementations, but I would not say it's infeasible nowadays.
4. Eq. (6) maybe it's worth pointing out that these equations implicitly define the velocity, as the viscosity depends on the velocity. Explicit formulas can be obtained.  Also, it is possible to account for a sliding term in the SIA model by adjusting the SIA velocity with the term -ρ g H ∇ s / β. I guess this is the formulation you use in the Greenland runs.
5. Section 2.2: SSA solver relies on the hypothesis that the velocity is uniform in the vertical direction. So the velocity, at any z-coordinate is the same as the depth-averaged velocity. I think it would be clearer in equations (8), (9) and (10) to substitute the velocity and the basal velocity with the averaged velocity.
6. Eq. (14). I would point out that the generalized integrals depend on the velocity (except from the n=1 case).
7. Eq. (19). One can go directly to (17) to (19) byt taking the limit for \beta that goes to infinity, which gives in fact the no-slip condition.
   
   After eq. (19). The two-steps recipe to compute the velocity only holds for the case n=1, otherwise one needs to compute the vertical velocity in order to compute the viscosity needed for computing F_2 in eq. (19). Please explain how the problem is solved in the general case  n >= 1.
8. page 10, line 1: For the stability analysis you can keep infinite domains. Just define x = n Δ x, with n being any integer (not just natural). I think this simplifies the analysis since you do not have to worry about boundaries, and in particular you can properly invert the circulant matrix arising from the discretization. It's understood that when solving the problem in practice you'll work on finite domains and have boundary effects.
9. Eq. (43): Only the solution with the negative sign is acceptable.
10. Section 3.2 I see a potential issue in the way the stability analysis is conducted. I would think that in most codes, when solving the thickness equation (59), the SSA and SIA velocities are discretized in the same way, whereas in the proposed analysis the SIA velocity term is expressed as a function of the thickness and treated differently. I think it would be best to consider the full velocity (u_sia + u_ssa) in (59) as a single object, take the derivative of H (u_sia + u_ssa), to get the advective and divergence terms as done in the DIVA case, and discretize them as done in the Diva case. In this way, the effective hybrid velocity should be the same as the depth-averaged hybrid velocity in (61).
11. Section 3.3, Similarly, I would not treat separately the terms of the L1L2 depth-averaged velocity in eq. (71), and I would avoid reformulations in (72) but simply discretize the thickness equation as for the Diva method. This might be the source of discrepancy seen when solving the problem with the CISM and Yelmo discretizations.
12. Eq. (72) It took me a long time to figure out how the equations have been derived. I would add some more steps to make it clearer.
13. Section 3.4 I think the instability in the L1L2-SIA model has been observed by BISICLES developers as well. To alleviate this, I believe BISICLES adopts a modification of the L1L2 method, consisting of avoiding the step in eq. (28), and taking the velocity to be uniform in the vertical direction and equal to the basal velocity.  It might be worth mentioning this. (In your analysis, this modified L1L2 method would be indistinguishable from the SSA method, because you consider a uniform viscosity and n=1.)
14. Section 4. Can you comment on the fact that the fitted exponent p goes from 1.5 to 2.8 instead of from 1 to 2 as in the theory. Does it depend on the Glen's law exponent n?
15. Page 28, Line 22. There are two possible values for r, say r₊ and r_-. We have that r₊ > 1, which implies, see (91), that the effect of a point x on the solution at y, grows exponentially with the distance between x and y. This is clearly not physical, so r_- is the only physical solution. The two solutions exist because we are considering an infinite domain. Anyway, as I mentioned before I think that there is no need to do the analysis for a finite domain. That complicates the analysis in the Appendix without really adding any value.