Anderson, E. and 10 others. 1999. LAPACK User’s Guide. Third Edition. Philadelphia, PA, Society for Industrial and Applied Mechanics.
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.
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.
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.
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.
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.
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.
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
The diffusivity D splits into a deformation component Ddef
and a sliding component Dslid
and their discretizations lead to
Finally, the discretized version of the ice-thickness Equation (2) can be written
where k = (ρg3/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 t1, t2 or t3. 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
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 to
+ Δ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 = to
+ Δt/2 and t
3 = to
For each of the lines (terms of constant index j), this leads to a system of Ny
equations with three unknowns (Ny
being the number of rows) written in a matrix form implying a tridiagonal matrix. After inverting all of the Nx
resulting matrices (Nx
being the number of lines), an intermediate surface at t
+ Δ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 to
+ Δ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
Appendix C. SI Scheme
In this scheme, all neighbouring points are considered implicitly, that is, directly at t0 + Δ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 NxNy x NxNy
matrix resulting from NxNy
equations with five unknowns (central and all neighbouring points evaluated at to
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 NxNy
gridpoints (i, j) of the domain.
Reshuffling Equation (A10) gives the entries of the A matrix of Equation (C1). A is a Nx Ny x Nx Ny matrix, indexed by ind = (j − 1) Nx
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
where c = 2(ρg)3/5. Spatial integration of ∂qx/∂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), (Nx, j) and (Nx, 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 ∂qx/∂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 Si,j terms (first line of Equation (A10)) have been considered to be implicit for both half time-steps, whereas Si,j
[1 + k(Di,j
+ Di—1j +D
) could be considered implicit and Si,j
+ Di–1j + Di–1
)] left explicit. Over the first half time-step, this leads to the same expression for ∂qx/∂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.