## References

Anderson, E. and 10 others. 1999. LAPACK User’s Guide. Third Edition. Philadelphia, PA, Society for Industrial and Applied Mechanics.

Glen, J.W.
1955.The creep of polycrystalline ice. Proc. R. Soc. London, Ser. A, 228(1175), 519–538.

Greve, R. and Calov, R.. 2002.Comparison of numerical schemes for the solution of the ice-thickness equation in a dynamic/thermodynamic ice-sheet model. J. Comput. Phys., 179(2), 649–664.

Hindmarsh, R.C.A.
2001.Notes on basic glaciological computational methods and algorithms. In Straughan, B., R. Greve, H. Ehrentraut and Y. Wang, *eds. Continuum mechanics and applications in geophysics and the environment* Berlin, etc., Springer-Verlag, 222–249.

Hindmarsh, R.C.A. and Payne, A.J.. 1996. Time-step limits for stable solutions of the ice-sheet equation. *Ann. Glaciol*., 23, 74–85.

Hutter, K.
1983. Theoretical glaciology; material science of ice and the mechanics of glaciers and ice sheets. Dordrecht, D. Reidel Publishing Co./Tokyo, Terra Scientific Publishing Co.

Huybrechts, P.
1992.The Antarctic ice sheet and environmental change: a three-dimensional modelling study. Ber. Polar-forsch. 99.

Le Meur, E. and Vincent, C.. 2003. A two-dimensional shallow ice-flow model of Glacier de Saint Sorlin, France. *J. Glaciol*., 49(167), 527–538.

Le Meur, E.
Gagliardini, O.
Zwinger, T. and Ruokolainen, J.. 2004. Glacier flow modelling: a comparison of the Shallow Ice Approximation and the full-Stokes equation. *C. R. Phys*., 5(7), 709–722.

Van den Berg, J.
van de Wal, R.S.W. and Oerlemans, J.. 2006. Effects of spatial discretization in ice-sheet modelling using the shallow-ice approximation. *J. Glaciol*., 52(176), 89–98.

Vincent, C.
Vallon, M.
Reynaud, L. and E. Le Meur. 2000.Dynamic behaviour analysis of glacier de Saint Sorlin, France, from 40 years of observations, 1957–97. J. Glaciol., 46(154), 499–506.

Weertman, J.
1964.The theory of glacier sliding. J. Glaciol., 5(39), 287–303.

## Appendix A. Discretization of the Ice-Thickness Equation

Discretization of the ice-thickness Equation (2) requires the calculation of several derivatives according to

Therefore,

and

The diffusivity *D* splits into a deformation component *D*_{def}
and a sliding component *D*_{slid}
according to

and their discretizations lead to

Finally, the discretized version of the ice-thickness Equation (2) can be written

where *k* = (*ρ*g^{3}/5) x (*Δt* /Δ^{2}).

In the last equation, the lower indices for the surface *S* refer to the gridpoints. Depending on the position indices, *S* on the lefthand side is evaluated either at t_{1}, t_{2} or t_{3}. Respective values for these times depend upon the selected scheme and reflect its degree of implicitness, as detailed in Appendices B and C. The terms on the righthand side and the diffusitivity *D* are always evaluated explicitly, that is, at time t*
*_{o}
.

## Appendix B. ADI Scheme

Forward integration of Equation (2) is carried out in two half-steps during which an implicit direction is chosen. For instance, during the first half-step, the system is solved according to the lines x, *y* where the x, *y* directions are denoted by the indices *i, j*. This implies that the system variables can be reduced to the central point (i, *j*) and the two neighbouring points alongthe implicitdirection with indexes (*i −* 1, j)and (i +1, j). Therefore, upper surface points
are considered implicit (evaluated at *t*_{o}
+ Δ*t*/2), whereas remaining neighbouring surface points along the *y* direction remain explicit (i.e. evaluated at *t*
_{0}).That means *t*
_{1} = *t*
_{2} = *t*_{o}
+ Δt/2 and *t*
_{3} = *t*_{o}
.

For each of the lines (terms of constant index *j*), this leads to a system of *N*_{y}
equations with three unknowns (*N*_{y}
being the number of rows) written in a matrix form implying a tridiagonal matrix. After inverting all of the *N*_{x}
resulting matrices (*N*_{x}
being the number of lines), an intermediate surface at *t*
_{0}
*+ Δt/2* is obtained, and from this surface, the second step of the forward integration is now performed according to the rows (*y* direction) in order to provide the final surface at *t*_{o}
+ Δt. This procedure is exactly the same as above but with inversion of indices *i* and *j*. It then leads to *t*
_{1} = *t*
_{3} = *t*
_{0} +Δ*t* , *t*
_{2} = *t*
_{0} +Δ*t*/2.

## Appendix C. SI Scheme

In this scheme, all neighbouring points are considered implicitly, that is, directly at t_{0} + Δt, leading to *t*
_{1} = *t*
_{2} = *t*
_{3} = *t*
_{0} + Δ*t*. The area needs to be swept only once, leading to a system that can now be expressed with a pentadiagonal *N*_{x}N_{y} x N_{x}N_{y}
matrix resulting from *N*_{x}N_{y}
equations with five unknowns (central and all neighbouring points evaluated at *t*_{o}
+ Δ*t*),

where is the vector of all implicit surface elevations
*A* is a pentadiagonal matrix and is a vector of all explicit terms, including the previous step surface and mass balance at all *N*_{x}N_{y}
gridpoints (*i, j*) of the domain.

Reshuffling Equation (A10) gives the entries of the *A* matrix of Equation (C1). *A* is a *N*_{x} N_{y} x N_{x} Ny matrix, indexed by ind = (*j −* 1) *N*_{x}
+ *i*.

With *c*
_{1} = 0.5 x Δ*t*/Δ^{2}, the entries of the five diagonals of *A* can be written

The remaining terms are equal to zero. The values of the vector are given by

## Appendix D. Mass Conservation

### D1. ADI case

Once integrated over the first half time-step with the ADI scheme, the x and *y* derivatives of the flux along the x and *y* directions can be expressed as

and

where *c =* 2(*ρ*g)^{3}/5. Spatial integration of *∂q*_{x}/∂x gives the following equations for all values of *j*:

and finally goes to zero since the domain is chosen such that its outer fringe remains ice-free, leading to zero diffusion coefficients at points (0, j), (0, *j* – 1), (*N*_{x}, j) and (*N*_{x}, j – 1).

Spatial integration at the same half time-step of the term *∂qy/∂y* in the *y* direction gives

for all *i*. Because of remaining terms of the form

spatial integration of *∂qy/∂y* (and therefore that of *∂q*_{x}/∂x+ *∂qy/∂y*) over the first half-step is not equal to zero.

### D2. An alternative form of the ADI scheme

In order to guarantee as well as Δ *·*
at each half time-step, another form of the ADI scheme can be proposed. So far, all *S*_{i},j terms (first line of Equation (A10)) have been considered to be implicit for both half time-steps, whereas *S*_{i,j}
[1 + *k*(*D*_{i,j}
+ *D*_{i,j-1}
+ *D*_{i—1}j +D*
*_{i–1/}
·*
*_{j–1}
) could be considered implicit and *S*_{i,j}
[k(D *
*_{i,j}
+ *D*_{i,j}__{1}
+ *D*_{i–1}j + *D*_{i–1}
,*
*_{j–1}
)] left explicit. Over the first half time-step, this leads to the same expression for *∂q*_{x}/∂x as before, whereas the expression for *∂qy∂y* then reads

It can easily be shown that this leads to (similar to the approach used with Equation (D2)). An analogous calculation leads to the same result for the second half time-step.

Although it is inherently mass-conserving, under certain geometrical conditions this scheme is numerically unstable unless abnormally small time-steps are used. In a sensitivity test, we observed that some reasonable geometries (bedrock and ice surface slopes) require very small time-steps of the order of 0.001 years. This means that in order to avoid any numerical instability for the various geometries possibly encountered with a real glacier, the time-step has to be extremely small. It was therefore impossible for us to use an ADI scheme that was both inherently mass-conserving and numerically stable.