Dear Daniel Moreno-Parada,

thank you very much for you review! I appreciate your comments very much in times when most of the reviewers comment only on the way how the introduction is written. On the other hand, it looks as if you want a paper different from what I want, which makes it a bit complicated. Since you made several comments, I uploaded my preliminary repsonses as an attachment.

Best,

Stefan 

Dear Thomas Zwinger,

thank you very much for your thorough review! 

For me, the central point is the "applicability" of the SIA to fine grids since it defines whether or not the entire exercise makes sense. Therefore, I want to discuss this point here, while I uploaded  the detailed responses as a separate document. 

Your statements about the SIA at fine resolutions confuse me. One could even interpret them in such a way that the SIA became increasingly wrong by itself if we use finer grids. I guess that this  is not what you want to say.

The mathematical analyses, typically using a vertical/horizontal scaling factor epsilon refer to the situation that we expand or compress the domain (and the topography) horizontally. If we compress it and make the valleys narrower, horizontal velocity gradients increase and the SIA becomes worse. So far clear, but we cannot identify the spatial resolution of the grid with the horizontal length scale in the denominator of epsilon. If we cover the same topography with a finer grid (without introducing additional roughness), nothing happens.  

For a given topography, the SIA is good if the ice thickness is small (what the SIA is made for). It becomes worse with increasing thickness and finally becomes better again when the ice overflows the small valleys. If the valleys are not wide enough in relation to their depth, horizontal stresses become relevant in nature. Then the SIA overestimates the velocity in the middle of the valley and underestimates it at the sides. This effect is basically independent of the spatial resolution for fine grids, but the accuracy does not automatically become better if we decrease the resolution so that we cannot resolve the valleys any more. Therefore, I disagree to your statement that an SIA appropriate mesh resolution would be > 1000 m. From my point of  view, the typical test examples include regions (valleys) where the SIA is sufficient and regions where it is inaccurate.

As a central point, however, I think that we must distinguish properly between the physical deficiencies of the SIA (narrow valleys, strong sliding, ...) and numerical issues. The SIA leads to some kind of nonlinear diffusion equation (although the linearized diffusivity is 3D instead of D, see detailed responses). Diffusion equations dampen spatial undulations exponentially. If L is the wavelength of the undulation, the respective decay (damping) constant is D/(2*pi*L)^2. So short-wavelength undulations are damped very rapidly by the SIA itself. Tentatively, I would say that the assumption of hydrostatic pressure even makes damping too fast (since it affects the entire vertical column), but this is practically not relevant.

Numerically, however, it is opposite. Explicit schemes cannot capture the fast damping, but overshoot and generate growing oscillations from small undulations. The hydrostatic approximation makes these oscillations even worse. But in any case, this is a deficiency of the numerical schemes and not a deficiency of of the SIA itself,  as I would read from your comments.

Summarizing the first point, I strongly disagree to your points that the SIA becomes worse with increasing spatial resolution (still thinking that this is not what you wanted to say) and that the numerical issues with the SIA are a deficiency of the SIA itself and not of the numerical scheme. Therefore, I still believe that the new numerical scheme brings the solution closer to the exact solution of the SIA (with all its problems in narrow valleys and if sliding is strong).

There will probably be an extension by horizontal stresses in the future, but I do not know so far what this extension will look like.

As you also raise an issue about the correction for steep slopes in your major points, let me briefly comment on this here. The equations for the hydrostatic pressure and for the horizontal shear stress are the same for the original shallow water equations and for the SIA. To my knowledge, this is not coincidence, but arises because it makes no difference up to this  point whether we start from Navier--Stokes or Stokes and what the rheology looks like. As far as I know, Hutter (1983, doi 0.1007/978-94-015-1167-4) already proposed a version for an inclined bed.

The central idea (beyond zero pressure at the free surface) behind this extension is that the velocity is parallel to the bed (not necessarily horizontal), which reduces the hydrostatic pressure compared to the original shallow water equations. It is good to know that this version is consistent with the slab solution (I only have the 2010 lecture notes of Ralf Greve and  not the textbook, but it should be similar). I wrote a short proof of this, as requested, in the detailled responses.

In essence, I still think that the version my correction is better for steep slopes (ice surface almost parallel to the bed), while the original form may be better for glacier fronts with an almost horizontal bed (although it does not comply with the condition defined by Hutter, 1983). So I found it the best option to leave it to the users which version to take.

Best regards,

Stefan