Restoring mass conservation to shallow ice flow models over complex terrain

Numerical simulation of glacier dynamics in mountainous regions using zero-order, shallow ice models is desirable for computational efficiency so as to allow broad coverage. However, these models present several difficulties when applied to complex terrain. One such problem arises where steep terrain can spuriously lead to large ice fluxes that remove more mass from a grid cell than it originally contains, leading to mass conservation being violated. This paper describes a vertically integrated, shallow ice model using a second-order flux-limiting spatial discretization scheme that enforces mass conservation. An exact solution to ice flow over a bedrock step is derived for a given mass balance forcing as a benchmark to evaluate the model performance in such a difficult setting. This benchmark should serve as a useful test for modellers interested in simulating glaciers over complex terrain.


Introduction
Numerical simulation of glaciers and ice sheets is essential for understanding the cryospheric response to a changing climate and is increasingly an integral part of modern climate change projections.Although the vast majority of fresh water capable of causing sea level rise over the long-term lies in the Antarctic and Greenland ice sheets, arguably it is the mountain glaciers which are most susceptible to climate change in the near future.Most of them lie at moderate to high latitudes in mountainous terrain.It has been shown that these glaciers are the largest contributor to contemporary sea level rise and that they will contribute to sea level rise in the coming century (e.g.Radić and Hock, 2011;Marzeion et al., 2012).This importance of alpine glaciers creates a need to understand their behaviour in coming decades.One approach is to explicitly simulate glaciers at a sub-kilometer resolution over large, ice-covered regions of the globe.Such an approach demands models of ice dynamics capable of simulating mountain glaciers in computational domains containing many (e.g. 10 7 ) grid nodes over century-long model periods.Higher-order ice dynamical models are capable of simulating individual glaciers or large ice sheets, but presently their high computational demands restrict their use over domains required to simulate regional mountain glacier evolution.By reducing the complexity of the stresses that are simulated in a dynamical model, greater computational effort can be put into addressing large-scale problems at some cost to model accuracy.One such model is the vertically integrated, shallow ice formulation discretized using finite differences (e.g.Mahaffy, 1976).This approach to the shallow ice problem has been used for example to simulate mountain glacier complexes in the Sierra Nevada, USA during the last deglaciation (Plummer and Phillips, 2003) and glacier advances on the summit of Mauna Kea, Hawaii during the last deglaciation (Anslow et al., 2010).
One major problem with standard numerical solvers for shallow ice models is a tendency not to conserve mass in regions where thin ice is draped over steep bed topography.Many of the spatial discretization schemes developed for ice sheets (e.g.Hindmarsh and Payne, 1996;Huybrechts and Payne, 1996), where steep bed topography can be an issue as well (e.g.East Greenland), will spuriously create mass at such troublesome locations.Where surface gradients are large, this mass creation can lead to very large errors Published by Copernicus Publications on behalf of the European Geosciences Union.
A. H. Jarosch et al.: Restoring mass conservation to SIA models in modelled steady-states.This paper describes the application of a second-order flux-limiting spatial scheme to the finite difference solution of a shallow ice model that ensures mass conservation.Furthermore, we describe a benchmark test case along with an exact solution upon which models can be tested for mass conservation when such situations arise.Confirming that a shallow ice model can meet the benchmark described here along with the benchmarks for the transient simulation of a growing ice sheet described by Bueler et al. (2005) is strongly recommended prior to conducting simulations of glaciers over rough topography.

Standard shallow ice models and numerical methods
In a Cartesian coordinate system with the xy-plane oriented horizontally, the continuity equation, together with Glen's flow law (Glen, 1958), are the equations solved for the isothermal shallow ice model (Fowler and Larson, 1978;Morland and Johnson, 1980) where s(x, y, t) is ice surface elevation, h(x, y, t) = s(x, y, t) − b(x, y) is ice thickness, b(x, y) is bed elevation and A and n are the rate factor and power law constants in Glen's flow law, while ρ and g are ice density and acceleration due to gravity, respectively, and ṁ is surface mass balance.∇ is the 2-D gradient operator.Importantly, Eq. (1) holds only where there is ice (h > 0).Ice geometry evolution models are intrinsically free boundary models in which parts of the domain may be ice-free.A complete formulation of the ice flow problem must therefore incorporate a means of evolving ice-covered and ice-free parts of the domain geometrically.In ice-free parts of the domain, h = 0 and s = b.Ice will grow if ṁ > 0, but not otherwise.Taken together, this implies that, when h = 0, Negative ice thicknesses are never realized.In addition, at the ice margin (the free boundary between regions where h = 0 and h > 0), mass must be conserved and in addition we expect the surface s to be at least continuous.This implies at this free boundary, with n normal to the free boundary in the xy-plane.The formulation of Eqs. ( 1)-( 3) is known mathematically as an obstacle problem (Evans, 1998).Using the inequality constraints one can re-write the problem in its so-called weak form as a variational inequality, which allows various theoretical advances to be made, mostly in demonstrating the well-posedness of the shallow ice problem and analyzing the convergence of finite element discretizations (e.g.Calvo et al., 2002;Jouvet et al., 2011;Jouvet and Bueler, 2012).
Our aim here is more practical.We address shortcomings in widely used numerical methods for solving shallow ice problems.A frequently used approach is to treat Eq. (1) as a parabolic (i.e.diffusion) problem, writing it in the form where is a diffusion coefficient.This underlies the numerical methods first developed in Mahaffy (1976) and described in more detail in Hindmarsh and Payne (1996) or Huybrechts and Payne (1996).An often used time-stepping method updates ice surface elevation s i (x, y) = s(x, y, t i ) by using a lagged diffusivity D i = D(h i , |∇s i |) and solving for an unconstrained updated ice surface elevation si+1 through where t = t i+1 − t i .This corresponds to the semi-implicit time-stepping scheme described by Hindmarsh and Payne (1996).The actual ice thickness is then updated by truncating this solution anywhere the unconstrained ice surface elevation corresponds to negative ice thickness using a "naive" projection step such that A slightly more self-consistent approach to the inequality constraint h > 0 which governs Eq. ( 1) would be to apply Eq. ( 6) only where si+1 > b, and to demand instead that where si+1 = 0 (this essentially being the discretized version of inequality (2)), while not allowing negative si+1 at all.This is mathematically equivalent of finding the updated ice thickness by minimizing the functional subject to s ≥ b, where is the entire domain, i.e. the union of ice-covered and ice-free regions (Kinderlehrer and Stampacchia, 1980;Evans, 1998).In other words, s i+1 = arg min s≥b J (s), and this can be solved numerically using projected successive over-relaxation (PSOR) methods (Glowinski, 1984) that are similar to solving Eq. ( 6) with the projection step Eq. ( 7).Importantly, however, the continuum formulation of Eqs. ( 6)-( 7) as well as of Eq. ( 9) is misleading: D i = 0 anywhere the ice thickness h i is zero, suggesting that ice flow alone should not be able to expand the ice covered area, when clearly this should be possible.In the methods described above, a spatial discretization must be applied first, and the nature of this spatial discretization is crucial.
In particular, spatial discretization schemes designed for diffusion equations of the form of Eq. ( 4) with bounded diffusivities D may spuriously generate negative ice thicknesses.
In fact, such methods may not be appropriate at all in settings where bed topography is steep.The easiest way to understand this is to re-write Eq. (1) as a conservation law for ice thickness h = s − b: where h > 0, with an analogous inequality to Eq. ( 2) holding where h = 0.In steep terrain, the gradient term ∇(b+h) may now be dominated by bed slope ∇b, leading approximately to the hyperbolic problem (see also Fowler and Larson, 1978) In the absence of a surface mass balance term (i.e. when ṁ = 0), this hyperbolic equation in its continuum form preserves positivity, i.e. given non-negative initial conditions on h, negative h will never be generated.Spatial discretizations appropriate for hyperbolic equations will maintain this property.However, discretizations designed for parabolic problems, including the symmetric centered difference schemes described in, e.g.Huybrechts and Payne (1996), may not preserve positivity for h, and can therefore spuriously generate negative ice thicknesses.The projection step Eq. ( 7) of course will then set ice thickness back to zero where this occurs.However, in the process, this causes the numerical scheme to create mass, which can severely affect its results.Two of the most widely used discretizations in ice sheet models are those referred to as "method 2" (abbr.M2) and "method 3" (abbr.M3) by Hindmarsh and Payne (1996) or as "type I" and "type II"1 by Huybrechts and Payne (1996).All these schemes are appropriate for the primarily diffusive case of small bed slopes.Essentially, we can view these as finite volume discretizations on a regular mesh, with ice surface elevation piecewise constant on each cell.The location of cell centers are (x k , y l ) on a grid with uniform spacing such that x = x k+1 − x k and y = y l+1 − y l .We label the cells by indices (k, l) and denote the normal component of flux on cell boundaries such that the y-component of flux on the cell edge between cells (k, l) and (k, l + 1) is q y k,l+ 1 2 and the xcomponent of flux on the cell edge between cells (k, l) and (Fig. 1).
The M2 and M3 schemes both relate these fluxes to differences in surface elevation through where are the diffusivities evaluated on the cell boundaries.The fully discretized version of Eq. ( 1a) is then and the projection step Eq. ( 7) is applied cell-wise.Hindmarsh and Payne's M2 and M3 schemes only differ in how they handle the diffusivities . We define the norm of the surface gradient at the cell boundary (k + 1 2 )2 in accordance with Mahaffy (1976) as A. H. Jarosch et al.: Restoring mass conservation to SIA models and can now write M2, which uses an averaged ice thickness at the cell, as M3 uses an average over the factor h n+2 that appears in the definition of D, With these discretizations, the projection step scheme (cf.Eq. 7) has performed well in many ice sheet models, and in particular, has reproduced a number of known exact solutions outlined in Bueler et al. (2005Bueler et al. ( , 2007)).However, these exact benchmarks all refer to the case of a flat bed, for which we have h = s.Our aim here is to explore a number of complications that arise precisely when this is not the case.That is, we wish to study complications that are typically associated with bed undulations, and which become particularly relevant for modelling mountain glaciation.

Mass conservation problems in projection step schemes
One simple yet problematic case is the one of a mountain glacier sitting in a u-shaped valley.The projection step with a M2 or M3 diffusivity can generate a spurious mass flux out of the bare rock sidewall of a u-shaped valley into a glacier in the bottom of that valley.Here, we have a cell in which h i k,l = s i k,l − b k,l = 0 adjacent to a cell in which h i k+1,l > 0 and yet we also have s i k,l > s i k+1,l , as displayed in Fig. 2a 3 .That is, the ice-free cell has a higher surface elevation than the ice-covered cell.Consequently we expect that q x k+ 1 2 ,l > 0 and either an M2 or M3 scheme above predicts ice flowing from the ice-free into the ice-covered cell.If ice does flow out of the ice-free cell, then the time-stepping scheme (Eq.13) will predict a negative ice thickness for the respective cell after a single time step.The projection step (Eq.7) then sets the actual surface elevation s i+1 k,l back to the bed elevation b k,l .In terms of mass conservation, we have just extracted mass from the ice-free cell (k, l) and transferred it to the icecovered cell (k + 1, l).In the projection step (Eq.7), we have added that mass back into the cell (k, l) in a bid to avoid unphysical negative ice thickness.Formulated this way, the projection step scheme therefore creates mass.
This mass conservation issue was previously recognized by Plummer and Phillips (2003), who proposed a slightly modified scheme that prevents such a mass violation.In particular, Plummer and Phillips (2003) to either the 3 Subsequently we focus on mass conservation problems along the x-axis, but they can equally be generated in three dimensions.
M2 or M3 values suggested above except at cell boundaries that correspond to a glacier-rock wall boundary.These can be recognized as boundaries with indices (k + 1 2 , l) for which we have The first of these statements says that ice thickness is greater in the cell that is at a lower elevation, while the second statement says that one of the cells has zero ice thickness (which must therefore be the one with the greater surface elevation).For these cell boundaries, Plummer and Phillips (2003) set = 0 and they also apply an analogous scheme for D i k,l+ 1 2 at cell boundaries parallel to the y-axis.A slightly different scheme with essentially the same properties was developed independently by Haseloff (2009).
There is however another possible complication that is not captured by this adjustment of diffusivities.This can occur when a relatively thin glacier flows over a steep bedrock step, as in an icefall.Figure 2b shows the situation we have in mind.Here we can generate a significant ice flux out of the upstream grid cell (k, l) at the top of the ice fall, simply because of the large surface slope between the upstream cell (k, l) and the downstream cell (k + 1, l).This large flux can then lead to more ice flowing out of the grid cell (k, l) in a single time step from t i to t i+1 than was present at time t i .The updated ice thickness value si+1 k,l − b k,l becomes negative after a single time step and the projection step (Eq.7) sets it s i+1 k,l back to zero.Ice mass is created in the process.At time t i+2 , the upstream cell is likely to acquire a non-zero ice thickness s i+2 k,l again either from a positive mass balance, ṁ > 0, or through inflow from the cell (k − 1, l).Such a situation is possible in Plummer and Philips' scheme.After time t i+2 , we can therefore return to the same situation as time t i , with a thin ice cover in cell (k, l) and a steep surface slope into cell (k + 1, l).Mass can therefore be created on alternating time steps, causing the resulting error to grow over time.
The main reason why the M2 and M3 schemes above are able to create mass in this way is that they do not limit the flux across a cell boundary as the ice thickness in the upstream cell goes to zero.Consider a very small ice thickness h k,l in cell (k, l) whose bed elevation is greater than the surface elevation in the next cell downstream, b k,l > s k+1,l .Using h k,l h k+1,l and s k,l = b k,l + h k,l ≈ b k,l , an M2 scheme approximately gives the flux q x k+ 1 2 ,l across the (k + 1 2 , l) cell boundary as As h k,l → 0 the expression on the right does not go to zero.Thus a finite amount of mass can be extracted from a cell with vanishingly small mass content.This occurs because the diffusivity on the cell boundary is dominated by the non-zero Bedrock in grey and ice in light blue.
ice thickness in the downstream cell (k + 1, l).An analogous observation applies to M3 schemes.Below, we will illustrate this shortcoming of M2 and M3 discretizations further by showing that they fail to reproduce certain exact steady-state solutions to the shallow ice equations.Before we do so, we propose an alternative scheme for computing diffusivities that restores the mass conservation.

A mass-conserving scheme
The difficulties of conserving mass with both M2 and M3 schemes are all rooted in the computation of the diffusivities and D k,l+ 1 2 .These numerical artifacts stem from the evaluation of the ice thickness term h n+2 in the definition of diffusivity (Eq.5).In both schemes, h on the (k + 1 2 , l) cell boundary is evaluated numerically as an average over the ice thicknesses in the adjacent cells.Consequently, the diffusivity on the cell boundary does not go to zero when the ice thickness in just one of these cells goes to zero.
When there is an advancing ice margin, it is important that the diffusivity should not go to zero at a cell boundary adjoining the ice-free cell.Otherwise ice could never flow from an ice-covered cell into an ice-free cell, and the ice margin could never advance due to flow.However, we need to avoid the reverse situation in which too much ice flows from a barely ice-covered cell into another cell with lower surface elevation.
To do this, a flux-modification scheme is required.We adapt one of the flux-limiting schemes from the conservation law literature, namely a second-order Monotone Upstreamcentered Schemes for Conservation Laws (MUSCL, e.g.van Leer, 1979;Gottlieb and Shu, 1998) for the ice flux discretization.To do so we propose and use a new factorization of ice flux such that with as used in Eq. ( 4).This allows us to think of the shallow ice, mass continuity equation, Eq. ( 1), vaguely as a Burger's equation in the form of Note that Eq. ( 22) has the appearance of a hyperbolic conservation law, which is not actually the case as the gradient term ∇s that appears in the definition of ω depends on the gradient of h, and the flux is at least in part diffusive.However, for steep terrain, ∇s may be dominated by the bed slope ∇b, and in this limiting case Eq. ( 22) does become hyperbolic as previously explored by Fowler and Larson (1978).Our main interest in applying a flux-limiting scheme to the shallow ice problem defined by Eq. ( 22) lies in the positivity-preserving property of such schemes.This prevents the spurious extraction of excessive volumes of ice from grid cells in steep terrain that lies at the heart of the mass conservation problem identified in the previous section.
A distinct feature of MUSCL schemes is the separation of flux at the cell boundary (k + 1 2 , l) into two components, the (k + 1 2 + , l) and (k + 1 2 − , l) term, which we define below for our application.Ice thickness h at the cell boundary is once again the dominant term in the flux discretization, so we can www.the-cryosphere.net/7/229/2013/The Cryosphere, 7, 229-240, 2013 with the ratio of downstream to upstream ice thickness change and φ(r k,l ) being the flux-limiting function.We investigate the usability of two flux limiters in our study, the minmod limiter φ mm (r) and superbee limiter φ sb (r) (Roe, 1986): Using the ice thickness estimates from Eqs. ( 23) and ( 24), we can define two flux terms at the cell boundary and D i k+ 1 2 − ,l by using Eq. ( 23) instead of Eq. ( 24).To limit the flux at the cell boundary, one defines a minimum and maximum diffusivity such that and constructs a diffusivity for the (k + 1 2 , l) cell boundary as The diffusivities , and D i k,l− 1 2 can be constructed in a similar manner.Note that the local surface slopes are used to identify the upstream direction, which is needed in a MUSCL scheme to assign the correct limited flux terms.We recall that our initial equation was which can be discretized in time, as a simpler alternative4 to Eq. ( 6), explicitly using a forward Euler scheme All that is left to do is to define the gradient of the flux in its fully discretized form: The value for the time step t used is crucial in this forward scheme to provide numerically stable solutions.A stability condition can be used to automatically calculate a suitable value as Hindmarsh ( 2001) analysed time-stepping stability criteria and reports for explicit time-stepping schemes c stab < 1 2n for one-dimensional and c stab < 1 2(n+1) for two-dimensional configurations.In case of n = 3 this leads to c stab < 0.166 6 and c stab < 0.125, respectively.

One-dimensional steady-states
A good way to test a shallow ice code is to compare results with exact solutions (Bueler et al., 2005(Bueler et al., , 2007)).Below we construct such a steady-state solution which includes bed topography and a prescribed accumulation rate which is a function of position only.
In one dimension with the assumption of steady-state, the shallow ice model (Eq. 1) can be written in the form where the subscript "x" denotes an ordinary derivative, and To simplify matters, we assume that accumulation rate ṁ depends only on position x and is such that there is a continuous ice region for the interval 0 < x < x m .Here x m is the margin position, which must be determined as part of the solution.At x = x m , and we have In addition, we assume that there is no inflow of ice at the fixed upstream boundary x = 0.In that case, ice flux q can be found explicitly as a function of position for any 0 < x < x m : To simplify our notation, we write Q(x) = x 0 ṁ(x )dx .Given ṁ(x), Q(x) is then a known function of position.The unknown margin position is then determined implicitly by the second condition in Eq. ( 40), Given x m , ice thickness h must then be found as a function of position through solving the differential equation q = Q(x), or subject to the first condition in Eq. ( 40), h(x m ) = 0.There are no general methods for solving Eq. ( 44) analytically.To get around this, we restrict our choice of bed topography to generate a tractable problem.Our objective is to develop a test for numerical shallow ice codes that incorporate bed topography.Consequently, we do not wish to put b ≡ 0. On the other hand, Eq. ( 44) is easiest to deal with for a flat bed, in which case s x = h x .To make use of this, we consider a bed for which b is a step function, where b 0 and x s are constants, and we assume that 0 < x s < x m .
In the interval 0 < x < x s and x s < x < x m , this allows us to write Eq. ( 44) as which we can re-write as where we have expanded according to Eq. (1c) for clarity.Integrating using the boundary condition h(x m ) = 0, and subsequently solving for h, we get At the bedrock step at x = x s , we can therefore define an ice thickness just downstream of the step as In order to extend the solution to the interval 0 < x < x s , we can then integrate Eq. ( 47) backwards from x = x s : where h s− = lim x→x − s h(x).To close this solution, it remains to determine the ice thickness h s− at the top of the bedrock step.
In general, we expect the surface elevation s to be continuous.But s = h + b, so this implies This must be substituted in Eq. ( 50) with h s+ given by Eq. ( 49).It is, however, possible that h s− computed in this way is negative.Specifically, this occurs when h s+ computed in Eq. ( 49) is less than b 0 .In that case, Eq. ( 51) cannot hold, as h will be negative just upstream of the bedrock step, and will therefore violate the condition given in Eq. ( 41).A more acceptable solution can instead be obtained in the case that h s+ < b 0 if we put h s− = 0 in Eq. ( 50).
Allowing for this possibility, the required exact steadystate solution is given by Eqs. ( 50) and ( 48), with h s− determined by This solution, with a discontinuity in surface elevation, may seem an unnatural test for a shallow ice model.However, it can be shown that the solution we have given is in fact the correct limit of a solution with a continuous but steep bedrock step as the width of that bedrock step goes to zero.
In numerical simulations with finite grid size, steep steps in bed topography may not be well resolved, and it is desirable to have a numerical scheme that remains robust when this is the case.Crucial to the mechanism for mass creation described in Eq. ( 18) was a setting in which bedrock elevation b k,l in the upstream cell (k, l) is greater than ice surface elevation s k+1,l in the downstream cell, and this is precisely the situation realized when there is a sufficiently tall bedrock step h s− = 0 as advocated above.In fact, we show explicitly in Appendix A that there are conditions under which M2 and M3 schemes cannot reproduce such steady-states.The solutions in Eqs. ( 48) and ( 50) are still given in terms of the general flux Q(x) = x 0 ṁ(x )dx .Next, we give a choice Q(x) that allows us to compute h explicitly, and which is such that the corresponding accumulation rate function ṁ = Q x (x) is sensible (in particular, which satisfies the obvious requirement that Q(0) = 0 and which is such that ṁ is negative for x > x m , so that there is indeed a single ice body in steady-state).This is given by with a corresponding accumulation rate function Q satisfies Q(x m ) = 0, and x m can indeed be identified with the steady-state margin position.In addition ṁ < 0 if x > x m , so there is no ice outside the margin x m .With this choice of ṁ and Q, we have and hence the ice surface profile in Eqs. ( 48) and ( 50) can be computed as (56) for x s < x < x m , and (57) for 0 < x < x s .Here h s+ and h s− are given through the calculations 7 Numerical benchmark experiments

Cliff benchmark
To demonstrate the performance of our newly introduced scheme and to showcase our exact solution as a benchmark for numerical ice flow schemes in mountainous regions, we numerically implement Eqs. ( 12) and ( 13)5 .We test three different schemes.First, we use the diffusivity from Eq. ( 15) and refer to the corresponding solution as "M2" results.Secondly, we compute results using the diffusivity from Eq. ( 30) along with the superbee flux limiter, Eq. ( 29).We label this solution "MUSCL superbee".Using the minmod flux limiter, Eq. ( 28), gives slightly different results, which we will discuss below.Thirdly, we use a simple upstream scheme based on writing the evolution Eq. (1a) in the form where u = − h n+1 |∇s| n−1 ∇s.Applying a simple upwinding to h in Eq. ( 60) leads to the following scheme (e.g.Aðalgeirsdóttir, 2003): define an upstream ice thickness through and so the diffusivity becomes We refer to this case as "Upstream" in Fig. 3.For temporal evolution, we solve Eq. ( 35) with a sufficient stability condition as mentioned earlier.
We start our numerical solutions with an initial condition of zero ice.We assume that the continuum solution should evolve toward the single steady-state solution which we have found exactly.The results of our numerical computations for the M2 scheme, the MUSCL superbee scheme, and the upstream scheme are displayed in Fig. 3 in comparison with a result computed with Eqs. ( 56) and (57).We plot numerical results in 1000 yr intervals for a 50 000 yr evolution of the models.The MUSCL superbee scheme (blue lines in Fig. 3) converges towards the steady-state solution (orange line in Fig. 3. Comparison of the "MUSCL superbee" scheme (blue lines) with a classical "M2" scheme (red lines), an upstream scheme (green lines), and a solution computed with Eqs. ( 56) and (57) (orange line).For all numerical schemes, solutions are plotted at 1000 yr intervals for a 50 000 yr evolution.Fig. 3), whereas the classical M2 scheme and the upstream scheme fail to do so and create a large amount of spurious mass.
We compare volume estimates between the model outputs and the explicit solution.Integrating the steady-state solution results in a target 2-D volume of 4.539371 × 10 6 m 2 .After 50 000 yr of evolution, the "MUSCL superbee" scheme ends with a volume of 4.399017 × 10 6 m 2 , or a relative error of −3.092 %.The classical M2 scheme leads to a volume of 10.93219 × 10 6 m 2 , or a relative error of 140.830 % and the upstream scheme to a volume of 7.586381 × 10 6 m 2 or a relative error of 67.124 %.Convergence of the MUSCL scheme with both flux limiters as well as the M2 scheme and the upstream scheme for different x towards the explicit solution is demonstrated in Table 1.Note that the relative error of the M2 scheme is increasing with decreasing x.The upstream scheme displays a quite different behaviour.As long as the horizontal resolution, x, is sufficiently larger than half the vertical bedrock step height, b 0 , the upstream scheme performs reasonably well in the benchmark, even though not as well as either of the MUSCL schemes.As soon as x ≤ b 0 /2 the relative volume error of the upstream scheme increases dramatically and the scheme fails the benchmark (cf.Table 1).Therefore we limit the remaining comparisons between numerical schemes in the manuscript to the M2 and MUSCL schemes.
The mass conservation problem in the projection step scheme, described in Sect.3, has been tested with our new scheme as well.We create a setup similar to the one displayed in Fig. 2a, with as spatial resolution of x = 200 m and let it evolve for 50 000 yr with ṁ = 0.The result is dis-Table 1. Relative volume errors, RE vol = (V numerical − V exact )/V exact • 100, for our schemes, the M2 scheme, and the upstream scheme for different spatial resolutions.V exact = 4.546878 × 10 6 m 2 in the 2-D case described in Sect.7 with results plotted in Fig. 3 Fig. 4. Comparison of the "MUSCL superbee" (blue lines), with a classical "M2" scheme (red lines) for the mass conservation problem described in Sect.3. The initial surface is displayed as a magenta line.For both numerical schemes, solutions are plotted at 1000 yr intervals for a 50 000 yr evolution.
played in Fig. 4. We monitor the changes in ice volume, which should be zero as ṁ = 0.After 50 000 yr, the solutions with the MUSCL superbee scheme (blue lines in Fig. 4) as well as the MUSCL minmod scheme (not shown) conserve mass whereas the M2 scheme has a relative volume error of −9.5 % in comparison with the initial volume.The earlier described modification, Eq. ( 17), to the M2 scheme has not been applied in this comparison, which demonstrated that both of our schemes have no mass conservation difficulties in this test as well.Thus a correction step as described in Eq. ( 17) is not required when using our schemes.

Bueler C benchmark
The Bueler C benchmark (Bueler et al., 2005)   grow an ice dome over 15 208 yr (as required by the benchmark) after which the numerical solution is compared to the exact one.We take error in central dome ice thickness, e dome =|h exact (0, 0) − h num (0, 0)|, and maximum ice thickness difference in the whole domain, e max = max(|h exact − h num |), as our performance measures.Figure 5 displays the decrease in e dome and e max with decreasing values of x for the same benchmark setup as displayed in Figs.7 and 8 in Bueler et al. (2005).
In both cases we demonstrate that the MUSCL scheme performs better than the M2 scheme for smaller grid sizes ( x ≤ 6250 m) if the right flux limiter is chosen, i.e. the superbee limiter (cf.Eq. 29).This is an anticipated result as the MUSCL scheme is second-order and thus more accurate than the M2 scheme, but it is surprising that an unfortunate choice of flux limiter, i.e. the minmod limiter (cf.Eq. 28), makes the MUSCL scheme perform worse in comparison to the M2 scheme.

Conclusions
After revisiting a well-known mass conservation problem of finite difference models for glacier flow in mountainous regions, we have identified another complication which arises with very steep topography.In that case, several widely-used numerical schemes will extract excess mass from cells with thin ice cover, and subsequently add mass to these cells again to avoid negative ice thicknesses, thereby violating mass conservation.
To overcome both problems, we propose using a secondorder flux-limiting spatial discretization for the diffusion term in the standard shallow ice equation.In this contribution we have investigated the applicability of a MUSCL scheme with two different flux limiters, the minmod and the superbee.
As a benchmark to evaluate the performance of the MUSCL scheme in comparison to M2 and upstream schemes in such steep topographies, we have derived an exact solution to ice flow over a bedrock step for a given mass balance forcing.Using this newly derived exact solution in combination with the well-established exact solutions of Bueler et al. (2005), we find the MUSCL scheme in combination with the superbee flux limiter a very suitable spatial discretization for mountain glacier flow models, which has no difficulties with the abovementioned mass conservation issues.
Our newly developed exact solution for ice flow over a bedrock step adds another exact solution-based benchmark to the existing ones (Bueler et al., 2005(Bueler et al., , 2007)), against which numerical ice flow models should be evaluated.If shallow ice flow models are to be applied in mountainous regions with complex topography, we anticipate that our proposed scheme and benchmark will help significantly to improve and evaluate such models.
Let q x 1 2 = 0, so the cell boundary to the left of the cell k = 1 is a domain boundary with no inflow, corresponding to x = 0 in the continuum solutions above.Assuming that cells 1, 2, . . ., k are ice covered, Eq. (A1) then shows that which is analogous to the statement that q(x) = x 0 ṁ(x )dx in the continuum problem above.
Suppose that there is a single ice mass between x = 0 and the margin x = x m , and that there is no ice for x > x m .Let discrete margin position be the cell boundary k m + 1 2 , so that h k > 0 for k ≤ k m but h k = 0 for k > k m , and similarly q k+ 1 2 > 0 for k ≤ k m but q k+ 1 2 = 0 for k > k m .Equation (A3) of course holds only for k ≤ k m .It can be shown from the projection step (cf.Eq. 7) that the margin location k m in steady-state is then given by the value of k m that satisfies both of the following inequalities: which are equivalent to the statement that x m 0 ṁ(x )dx = 0 in the continuum formulation above; it is easy to show that the margin location defined by these inequalities converges to the continuum solution of x m 0 ṁ(x )dx .Given a discrete margin location k m , ice thicknesses h k for k ≤ k m can then be computed recursively, starting with ice thickness just upstream of the margin at k = k m .For each k ≤ k m , we have flux q k+ 1 2 explicitly through Eq. (A3).To take an example, consider a M2 discretization for diffusivity, though the argument below can also be applied in slightly modified form to a M3 discretization.With s k = h k + b k , we then have This nonlinear equation must then be solved for h k , given ice thickness h k+1 at the next grid cell downstream, as well as the bed elevations b k and b k+1 .This procedure is started with k = k m , for which we have h m+1 = 0. Problems arise in this procedure if at some value of k we have b k > h k+1 + b k+1 .This occurs when surface elevation in the (k + 1)th cell is lower than bed elevation in the k-th cell.If we also demand that h k ≥ 0, then one can show that the expression on the left-hand side of Eq. ( A5) is bounded below by a quantity q min that depends only on bed geometry and on ice thickness downstream from the current cell, Hence no non-negative solution for h k can be computed from Eq. (A5) if q k+ 1 2 < q min (h k+1 , b k , b k+1 ).
In this situation, the assumption we have made in arriving at Eq. (A5) must break down.In particular, the assumption of a single connected ice mass in which h k > 0 for k ≤ k m must fail for the discrete solution even if it holds for the continuum problem, and the discrete solution will not approximate the continuous solution.Again, this occurs because the flux q k+ 1 2 does not go to zero even as the ice thickness h k in the upstream cell does.

Fig. 2 .
Fig. 2. The valley glacier case in (a) and the icefall case in (b).Bedrock in grey and ice in light blue.

Fig. 5 .
Fig. 5. Results of the Bueler C benchmark (cf.Sect.7.2) for increasing spatial resolution on a log-log scale.The maximum error in the whole domain, e max , is displayed in (a), and (b) shows the central dome height error e dome .

Jarosch et al.: Restoring mass conservation to SIA models
is an ideal test case to compare finite difference discretization schemes with a time evolving exact solution.In this benchmark, a time evolving mass balance is given to the flow code to www.the-cryosphere.net/7/229/2013/The Cryosphere, 7, 229-240, 2013 A. H.