1. Introduction

Since the 1970s, knowledge of the dynamical behaviour of ice sheets and their interaction with climate has greatly expanded through the development of ice-flow models, which simulate the evolution of ice under changing climatic conditions. All ice-sheet models are based on the principle of mass conservation; the ice flows under its own weight (e.g. Reference Oerlemans and van der VeenOerlemans and Van der Veen, 1984; Reference PatersonPaterson, 1994). However, a hierarchy of these models has evolved, with models differing in many aspects, such as applied rheology, treatment of the boundaries, coupling to the atmosphere and numerical details, such as the use of finite-element or finite-difference schemes. More fundamental differences come from how many elements of the full system of equations are actually solved. Besides solving the complete system, an approach called the shallow-ice approximation (SIA) is widely used (Reference HutterHutter, 1983). This approximation is based on the assumption that large ice bodies are mainly deformed by basal shear stress, and therefore neglects longitudinal stresses. Whereas solving the full system of equations obviously results in more realistic results (especially near the margin and ice divide), the SIA simulates the essential processes of ice flow quite well (Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson, 2004). The most important advantage of models based on the SIA is their computational speed. It is not yet possible to calculate the evolution of large ice sheets throughout complete ice-age cycles with the full system of equations, which is one of the reasons why the SIA is still widely used today (e.g. Reference MahaffyMahaffy, 1976; Reference OerlemansOerlemans, 1981a, Reference Oerlemansb; Reference Fastook and ChapmanFastook and Chapman, 1989; Reference HuybrechtsHuybrechts, 1990; Reference Calov and HutterCalov and Hutter, 1996; Reference GreveGreve, 1997; Reference Pollard and ThompsonPollard and Thompson, 1997; Reference Pfeffer, Sassolas, Bahr and MeierPfeffer and others, 1998; Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999; Reference Van de WalVan de Wal, 1999; Reference Van de Wal, Wild and de Wolde.Van de Wal and others, 2001; Reference Bintanja, van de Wal and OerlemansBintanja and others, 2002; Reference Näslund, Rodhe, Fastook and HolmlundNäslund and others, 2003; Reference Plummer and PhilipsPlummer and Philips, 2003). We are interested in the modelling of large ice sheets throughout ice ages. Therefore, we limit ourselves in this paper to a model using the SIA.

An attempt to compare and assess many different SIA ice models on continental scales was made as part of the European Ice-Sheet Modelling Initiative (EISMINT, funded by the European Science Foundation). Several models were evaluated for their physical properties and accuracy, using tests with simplified geometries (Reference HuybrechtsHuybrechts and others, 1996; Reference PaynePayne and others, 2000). These tests, however, did not focus on the numerical properties of the models. Reference Greve and CalovGreve and Calov (2002) investigated numerical behaviour. They focused on the use of different numerical schemes, whereas the present paper addresses the behaviour of a single numerical scheme as a function of varying spatial discretization. Reference Greve and CalovGreve and Calov (2002) only used two different gridpoint spacings. More importantly, neither EISMINT nor Reference Greve and CalovGreve and Calov (2002) combined any experiments for numerical properties with mass-balance parameterizations, which included the non-linear feedback with height. Instead, mass balance was parameterized as a function of distance only. This is an essential modelling element (Reference Oerlemans and van der VeenOerlemans and Van der Veen, 1984), as numerical errors can grow non-linearly over time due to this feedback mechanism. The mass-balance feedback with height is especially important for research which deals with ice-sheet inception, where the integration periods of model runs are very large and the climatic conditions remain relatively uncertain.

In this study we have assessed a two-dimensional, vertically integrated, finite-difference SIA ice model specifically for dependency of modelling results on varying gridpoint spacings (1–20 km) in cases where the mass balance was a function of altitude. When we say varying gridpoint spacings, we mean variations between model runs, not variations within one grid. All our experiments used spatially uniform grids. We stress that we did not attempt to assess the modelling errors by solving the SIA instead of the full system of equations. We focused on the numerical errors arising from inaccuracies in the numerical treatment of the SIA. As Reference Bueler, Lingle, Kallen-Brown, Covey and BowmanBueler and others (2005) discuss, these two approaches are fundamentally different. Although we focused on the effects of varying spatial discretizations, this did not include the resolving power of small-scale structures. We investigated the ability of each discretization to solve for the global behaviour of a large, smooth ice sheet as well as the influence of initial conditions and the sensitivity of the model results to the individual model parameters.

2. Theory

The ice model is based on the vertically integrated continuity equation (e.g. Reference HuybrechtsHuybrechts, 1992; Reference Van der Veen and BalkemaVan der Veen, 1999; Reference Oerlemans and BalkemaOerlemans, 2001)

(1)

where H is the ice thickness, U is the vertical mean horizontal velocity and B is the mass balance. Equation (1) can be rewritten as

(2)

where F is the horizontal ice flux. It is defined as

(3)

where h is the ice surface height and D is the diffusion. We used the SIA (e.g. Reference HutterHutter, 1983; Reference HuybrechtsHuybrechts, 1992; Reference Van der Veen and BalkemaVan der Veen, 1999), which assumes that the horizontal scale of the ice extent is much larger than the vertical scale. For a two-dimensional model, the diffusion D is then given as

(4)

with ρ the ice density, g the gravitational acceleration, and A and n the flow parameters. The value for n is set to 3. The parameter A is known to be temperature-dependent, but we set it to a constant value, for simplicity. The vertical mean horizontal velocity is given by

(5)

where U _d is the vertical mean horizontal deformation velocity (called ‘deformation velocity’, hereafter) and U _s the sliding velocity. We tried several parameterizations for the sliding velocity, but as long as the sliding velocity was smaller than the deformation velocity, the numerical properties of the model did not change. Therefore, we set Us to zero everywhere. For the two-dimensional model, the expression for the vertical mean horizontal velocity then reduces to the deformation velocity, which is given by (e.g. Reference Oerlemans and van der VeenOerlemans and Van der Veen, 1984; Reference HuybrechtsHuybrechts, 1992)

(6)

We also used a one-dimensional model, which is based on the two-dimensional model and is equivalent to a two-dimensional plane strain approximation. It is discussed in detail by Reference Oerlemans and BalkemaOerlemans (2001). The final expression for the vertical mean horizontal velocity is parameterized as

(7)

The flow parameter f _d (1.1 × 10^–24 Pa^–3 s^–1) represents the situation on a large ice sheet like Greenland (Reference Van de WalVan de Wal, 1999). The driving stress S _d is caused by the gravitational force of the ice. In the one-dimensional model, the diffusion D is given by

(8)

We prescribed the mass balance as a linear profile for which the only variable is height above a reference level

(9)

where β is the mass-balance gradient in a^–1 and E the equilibrium-line altitude (ELA) in metres. The mass balance is zero at E and increases with elevation until a positive upper limit B _max (in ma^–1 ice equivalent). Above that elevation, the mass balance remains constant with height.

3. Numerical Methods

3.1. Discretizations

The one-dimensional model is solved on a staggered grid (Fig. 1). At each time-step, between each two gridpoints (marked by circles in Fig. 1), F is calculated either as

Fig. 1. Staggered grid in one dimension.

(10)

with i the gridpoint and t the time-step, or as

(11)

Equations (10) and (11) correspond, respectively, to type I and II ice models of the EISMINT experiments (Reference HuybrechtsHuybrechts and others, 1996), which are widely used by ice-sheet modellers. The difference between the two types is whether the diffusion D is calculated at the staggered gridpoints (I) (marked by the circles in Fig. 1) or the regular gridpoints (II) with central differencing on the regular gridpoints (marked by crosses in Fig. 1). For type I the diffusion is calculated as

(12)

For type II the diffusion is given by

(13)

Type I models are numerically more precise, but require much smaller time-steps than type II models (Reference Hindmarsh and PayneHindmarsh and Payne, 1996; Reference HuybrechtsHuybrechts and others, 1996).

Finally, the time integration is calculated with forward differencing as

(14)

The numerical approach in the two-dimensional model is similar to that used in the one-dimensional model. All spatial derivatives are calculated with central differencing on a staggered grid (Fig. 2). However, the time integration is calculated differently, since explicit integration requires very small time-steps which results in a long computation time. Instead, an ‘alternating direction implicit’ method is used as described by several authors (Reference MahaffyMahaffy, 1976; Reference HuybrechtsHuybrechts, 1992). In this method, each time-step is divided into two steps. In the first step, the time integration for all the x-direction components of the equations is performed implicitly, whereas the y-direction components remain explicit. In the second step, this is reversed. All calculations were performed in double precision.

Fig. 2. Staggered grid in two dimensions.

3.2. Boundary conditions

At the lefthand side of the domain, we impose a boundary condition that corresponds to an ice divide, i.e. the surface gradient is zero:

(15)

At the righthand side, we impose that the ice thickness H is equal to zero. In the two-dimensional model, the ice thickness at all domain boundaries is set to zero.

At each time-step the new ice thickness is calculated from the old ice thickness and the mass balance throughout the grid, where all gridpoints are treated the same, whether the ice thickness is positive, zero or negative. Then, all negative ice thicknesses are set to zero. The initial ice thickness, i.e. the ice-thickness distribution at time-step t = 0, is free as long as the grid boundary conditions are fulfilled.

4. Results

We performed experiments with both the one-and twodimensional (type I and II) models. In each experiment, the ice margin evolved freely as a function of ice flux and mass balance. We calculated steady-state solutions by integrating over a period of 250 000 years for parameter sets listed in Table 1. For the two-dimensional model, a gridcell size of 2 × 2 km was the finest grid feasible. The ice sheet was too flat for the surface gradients to be calculated within the limits posed by the use of double precision for smaller gridcell sizes. The one-dimensional model did not show this problem.

Table 1. Parameter values for the experiments, with λ the bedrock slope, β the mass-balance gradient, B _max the cut-off value of the mass balance, E the ELA, and initial condition that of the ice sheet at the start of the experiment. The ice sheets used as initial conditions in experiments 2 and 3 are calculated using the original solution. The starting condition for experiment 2 is the result of experiment 1 for ∆x = 1 km. For experiment 3, instead of the result from experiment 1, the steady-state ice sheet for E = 150 m was used

We did not include any bedrock response, but prescribed the bedrock profile as

(16)

with b the bedrock height in metres, and b_max the maximum bedrock height, which we kept at 400 m for all experiments. The bedrock slope is given by λ, and r is the distance from the centre in kilometres, which in one dimension is given by x and in two dimensions by (x ² + y ²)^1/2. The combination of a sloping bedrock and an altitude-dependent mass-balance parameterization is similar to an ice sheet developing on a flat bed with a sloping equilibrium line. According to Reference WeertmanWeertman (1961), this kind of situation has at least one stable solution, as long as the slope is finite and non-zero. This means that under changing climate conditions the ice sheet is able to approach a new equilibrium.

4.1. Initial conditions

4.1.1. Experiment 1: growth to steady state

As a start, we evaluated ice growth towards a steady state with the parameters listed in Table 1 (experiment 1). We started the time integration without ice, and at time t = 0 we set the ELA E to 250 m and let the ice sheet evolve to steady state. Figure 3 shows the results. We scaled all calculated ice volumes with the volumes from the type II model associated with a gridcell size ∆x of 2 × 2 km in the two-dimensional case, such that the volume curve for ∆x = 2 km was always equal to 1. In the onedimensional case, all curves were scaled with the areas from the type II model associated with a gridcell size of 1 km.

Fig. 3. Top: Ice volume/area as a function of time for the parameter values listed in Table 1, experiment 1, for several gridpoint distances ∆x for (a) the two-dimensional type II model (volumes) and (b) the one-dimensional type I and II models (areas). The dashed lines represent the reference type II solutions, the solid lines the type II solutions calculated with a modified surface gradient (section 4.1.4), and the dash– dotted lines represent solutions calculated with a type I model. All volume curves are scaled such that the reference type II solution with ∆x = 1 km is equal to 1 for the one-dimensional case and ∆x = 2 km is equal to 1 for the two-dimensional case. The lower panel (c) shows steady-state profiles for the same discretizations with height in metres on the y axis and distance from the centre in kilometres, calculated with the one-dimensional type II model. The thick, black, dashed line is the ELA, in this case constant at 250 m.

The dashed lines in the upper panels of Figure 3 show the ice volumes as a function of time for both the twodimensional (a) and the one-dimensional (b) type II model for several discretizations. (The solid lines are discussed in section 4.1.4.) In the one-dimensional case, the curves were calculated as the integral over all ice thicknesses, and represent the area of a cross-section in an infinite plane. The dashed lines are based on the standard type II model as described in sections 2 and 3.

The dash–dotted lines in Figure 3b show the ice area for the one-dimensional type I model. We do not show the corresponding ice volumes for the two-dimensional type I model, since (see section 3.3) the required time-steps were so small that the calculation time was no longer feasible. Tests at coarser grids confirmed that the numerical behaviour in two dimensions was similar to the behaviour in one dimension. The type I models generally perform better than the type II models.

The lower panel (c) of Figure 3 shows the steady-state solutions for the different discretizations for the onedimensional type II model. The results in Figure 3 are striking: the results from the individual discretizations differ by tens of per cent. This is observed in both the onedimensional and the two-dimensional model, calculated with different numerical methods.

Since both models gave similar results, and since the onedimensional model required less computation time, all other experiments were performed with the one-dimensional model only. Even though type I models generally performed better than type II models, the fundamental behaviour was the same. For this reason, we performed the sensitivity tests in this paper only for type II models.

4.1.2. Experiment 2: maintaining steady state

In experiment 2 (see Table 1), we studied the dependence on the initial conditions by starting the test with the steady-state ice sheet for ∆x = 1km from experiment 1, instead of starting without ice. We kept the equilibrium line E at 250 m and continued the calculations for the different discretizations for the one-dimensional model. Figure 4 shows the results. In principle, the ice sheet should not change. As coarse grids introduce larger inaccuracies in the calculation of the surface gradient than fine grids, the modelled ice sheets adjust to slightly lower volumes for the different discretizations. The jumps in the curves are caused by the retreating ice margin. Each jump corresponds to the ice margin retreating over a distance of one gridpoint.

Fig. 4. Relative ice area as a function of time for several discretizations for experiment 2 (see Table 1). All curves are scaled relative to the reference solution calculated with ∆x = 1 km. The solid lines represent the solutions calculated with a modified surface gradient, and the dashed lines represent the type II reference solutions.

4.1.3. Experiment 3: retreat to steady state

In experiment 3 we studied the retreat of a steady-state ice sheet towards a smaller steady-state ice sheet corresponding to a warmer climate. The initial conditions were given by the steady-state ice sheet for E = 150 m, calculated with ∆x = 1 km (Table 1, experiment 3). At time t = 0, the value of E was set to 250 m. Figure 5 depicts the results. The jumps in the curves represent the moving ice margin. For the finest grids, the curves are practically smooth. The effects shown in Figure 5 are remarkably smaller in this retreat scenario than in the growth scenario and are of the same order of magnitude as the minor adjustments shown in Figure 4. This suggests that, for coarse grids, results are strongly dependent on initial conditions.

Fig. 5. Fig. 5. Same as Figure 4, but for experiment 3.

To further examine this, we calculated steady-state ice sheets for a wide range of initial conditions from no ice to very large ice sheets. Figure 6 summarizes the results; it shows the resulting area as a function of the ELA E for ∆x = 1 km (Fig. 6a) and for ∆x = 20 km (Fig. 6b). We only show stable solutions. By stable, we mean that when the steady states are locally perturbed while keeping the ELA constant, the model returns to the steady-state solution belonging to that ELA. Figure 6a shows three different regimes. The first is the regime for ELAs <400 m. In this regime there is only one stable solution possible for each ELA. The second regime is the regime for ELAs >800 m. In this regime there is also only one stable solution possible, namely no ice. In between these two regimes is a third regime where for each ELA there are two stable solutions. Depending on the initial conditions, the numerical integration of Equation (1) will result in either of the two stable steady states.

Fig. 6. Steady-state ice areas for two discretizations (a) ∆x = 1km and (b) ∆x = 20 km both for the type I (open circles) and the type II model (closed circles) as a function of ELA. The solid lines in (a) show branches of possible stable steady-state solutions. The grey areas in (b) show regions for possible stable steady-state solutions; lightest grey is the region for the type II model, the slightly darker grey that for the type I model and the darkest where the regions for both model types overlie each other.

Whereas for ∆x = 1 km type I and II models produced the same results, this is no longer true for ∆x = 20 km (Fig. 6b). The open and closed circles do not coincide. In addition, instead of clear lines with stable solutions, as shown in Figure 6a, there are now regions of possible steady-state solutions. Depending on initial conditions the numerical integration of Equation (1) will result in a stable steady state somewhere in the grey area. There are three different grey regions in the plot: the two lightest shades represent the regions with possible stable steady states for type I (darker of the two) or type II (lighter of the two) models; the third and darkest region shows where the regions for type I and II models overlie each other. This different behaviour for ∆x = 1 km and ∆x = 20 km is the result of numerical inaccuracies.

An interesting aspect is that the maximum possible steady states (indicated by the upper row of open and closed circles) for ∆x = 20 km are slightly too big for type I and too small for type II, compared to the solutions for ∆x = 1 km. We come back to this in section 5 where we discuss the reasons for the numerical problems.

4.1.4. A slope correction formula

One reason for the inconsistencies between the results from the different spatial gridcell sizes is the discontinuity in the surface gradient at the ice margin, as illustrated in Figure 7. The discretized surface gradient at the last point with ice is too small, as illustrated in Figure 7a by the difference between lines 1 and 2. Hence, the ice flux through the margin is underestimated. This hampers the growth of the ice sheet in the model. For one time-step, this is a small effect. However, for many time-steps, this error is significant since the feedback of height with mass balance is non-linear. Since the mass balance is limited to a maximum value, the ice sheet can only compensate up to a certain degree by thickening. This suggests that this value for the maximum mass balance is an important parameter. We examine this in section 5.

Fig. 7. Schematic drawing of the ice margin. (a) Line 1 represents the numerically calculated surface gradient at the last point before the ice margin. The thick dash–dotted line (2) represents the true surface gradient. (b) The lines represent the surface gradient at the last few points before the ice margin, except the last point, which is given by line 2 in (a).

Figure 7b shows the numerically estimated surface gradients at the last few gridpoints before the ice margin, which are not influenced by the discontinuity of the surface. Assuming a smooth ice surface, the difference between these surface gradients, which is the curvature of the surface, contains information about the surface gradient at the last gridpoint before the ice margin. So, as an alternative to the central differencing operator, we defined a new, forward operator for the type II model by using a Taylor expansion. For the one-dimensional case, this operator is given by

(17)

In this way, we avoid having to use any points beyond the ice margin, thus avoiding the discontinuity. We cannot use this for the type I model, because in that scheme the surface gradients are calculated on the staggered grid

(18)

It can be seen that in this case the curvature term is zero. Hence, it is not possible to use the curvature information, which was the basis for the new operator. The solid lines in Figures 3–5 show the results for this improved surface gradient method. While this approach clearly improves the solution, it does not remove the fundamental problems. The results with this new operator are not better than the type I model, because the flux through the margin is also calculated with the diffusion of the first point outside the ice margin, which is zero. The method implicitly assumes that the diffusion goes to zero in a linear way, which is of course not the case. This problem, however, is not easily solved.

The bifurcation curves in Figure 6 show a different effect. It seems that for a retreating ice sheet two numerical effects are competing. The first is the already described effect of the underestimated slope. The other effect is also the result of grid coarseness. When the ELA is moved up by 100 m, an ice sheet evolving on a fine grid will register that the ablation area is larger and that a part of the accumulation area receives less snow. On a coarse grid of ∆x = 20 km, however, the differences in surface elevation between the gridpoints are so large that the ablation area actually remains the same. Also, the first point in the accumulation area is already so high that it continues to receive the maximum value of the mass balance B_max. Hence, the ice sheet evolving on a coarse grid will hardly register any difference when the ELA is moved up by 100 m, and so will remain larger than an ice sheet evolving on a fine grid. Since type I models suffer less from the effect of the underestimated slope, this latter effect is detectable. For the type II model the latter effect is obscured by the underestimated slope effect.

4.2. Sensitivity of modelled ice sheets to individual parameters

Here we discuss the results of sensitivity tests to individual model parameters. In each experiment, we used both the original model, which we refer to as the reference model, and the improved gradient model, as defined by Equation (17). Again, all tests reported here were performed for the one-dimensional type II model only. Tests with the type I model showed similar behaviour. Table 1 (experiments 4–6) shows the specific parameters of the sensitivity tests that were carried out. In addition, we varied the ELA E. The results were not sensitive to the actual value of E as long as enough points covered the ice, which was only a problem at ∆x = 20 km. This, however, is a normal numerical effect and is not related to the effects studied in this paper. Again, all curves were scaled with the reference solution for ∆x = 1 km, such that the areas for ∆x = 1 km were always equal to 1. Each experiment started without ice. At time t = 0, the parameters were set to the values given in Table 1, after which the ice sheet freely evolved to steady state.

4.2.1. Experiment 4: bedrock slope, λ

In experiment 4, we varied X, the slope of the underlying bedrock, as shown in Table 1. The changes in mass balance between individual gridpoints are larger for a steeper slope than for a smaller slope. Figure 8a demonstrates that the larger the bedrock slope, the worse the results. As in the case of the varying ELA, the smaller the ice sheet, the fewer gridpoints actually cover the ice. However, this only affects the results for ∆x = 20 km.

Fig. 8. (a) Steady-state areas as a function of bedrock slope λ for several discretizations (Table 1, experiment 4). (b) Steady-state areas as a function of mass-balance gradient β (Table 1, experiment 5). (c) Steady-state areas as a function of maximum mass balance B_max (Table 1, experiment 6). All curves are scaled relative to the reference solution calculated with ∆x = 1 km. The solid lines represent the solutions calculated with a modified surface gradient, and the dashed lines represent the type II reference solutions. Colours of lines: see Figure 5.

For large slopes between the bedrock and the equilibrium line, the results are worse than for small slopes. However, with a slope equal to zero, the mass balance would be equal in the entire domain and no stable ice sheet would develop.

4.2.2. Experiment 5: mass-balance gradient, β

Here, the mass-balance gradient β was varied. A large mass-balance gradient increases the differences between the values of the mass balance at the individual gridpoints. Figure 8b shows that the larger the mass-balance gradient β, the larger the differences between the solutions for the different gridpoint distances. So, the solution turns out to be sensitive for the mass-balance difference between grid-points.

4.2.3. Experiment 6: maximum mass balance, B _max

The maximum mass balance B_max was varied. Figure 8c demonstrates that the problem with discretization only occurs when the maximum mass balance is lower than approximately 1 m a^–1 ice equivalent for this particular set of parameters, hence only in dry areas. The quantity B_max is strongly related to the total mass budget of an ice sheet, to which the solution is very sensitive. We clarify this further in section 5.

4.3. Ice inception

Experiments 1–6 used a linearly sloping bedrock, and the resulting solutions turned out to be sensitive to the specific slope. In practice, the bedrock slope is seldom linear. Therefore, the next and final experiment had a non-linear bedrock slope. Figure 9 shows the results. We superimposed a sinusoidal perturbation on the linear bedrock of experiment 1 (see Table 1), but now with b_max = 800 m instead of the previously used value of 400 m, thus creating a bedrock with smoothly varying slopes (Fig. 9b). Furthermore, we used an ELA of 130 m. For the first 5000 model years, the ice lies on relatively small local slopes, resulting in minor differences between the area curves (Fig. 9a). Once the ice reaches larger slopes, the differences between the area curves grow to the same order of magnitude as in the case of a purely linear bedrock. Finally, the ice sheets for the two finest discretizations meet and merge into a larger ice sheet. The coarser-discretization runs suffered from severely underestimated ice volumes, as a result of which the ice did not merge. Figure 9b shows the steady-state solutions for all discretizations; the different discretizations resulted in completely different ice histories. Hence, the reduced-volume effect, caused by the discretization, could have serious implications for the modelling of ice-sheet inception.

Fig. 9. (a) Area as a function of time for a sinusoidal bedrock and E = 130 m for several discretizations. The dashed curves represent the reference solutions, the solid lines the solutions calculated with a modified surface gradient. All curves are scaled relative to the reference solution calculated with ∆x = 1 km. (b) Steady-state solutions for the same discretizations with height in metres on the y axis and distance from the centre in kilometres.

5. Discussion

We performed a series of sensitivity experiments on the model parameters in a one-dimensional and a two-dimensional ice-flow model. From these experiments it emerges that the model’s numerical behaviour is very sensitive to the precise initial conditions, local geometry and climate. This may result in very unpredictable numerical behaviour and thus unreliable results for complex ice sheets. The cause of this behaviour is the sensitivity to the discretization of the ice margin. For a fine grid, the ice flux at the ice margin is large enough to overcome the negative mass balance at the next gridpoint in the model. The ice sheet will grow and continue to the next gridpoint and so forth. However, for a coarse grid, the same ice flux may not be large enough to overcome the more negative mass balance of the next gridpoint much farther away, hence ice growth in the model will stop. Because of the height–mass-balance feedback this seemingly small error can grow rapidly over time. The sensitivity tests in experiments 5 and 6 both affected the difference in mass balance between points. The most important parameter however, was the maximum value for the mass balance, B_max. This parameter affects the ability of the ice sheet to thicken and compensate.

We can explain the model’s numerical behaviour more fundamentally through the use of dimension analysis. Let us define new dimensionless variables according to Reference SaltzmanSaltzman (2002):

(19)

where L_v is a characteristic vertical length scale of the ice, a a characteristic value for the mass balance, U_c a characteristic velocity scale for the horizontal ice motion and LҺ a typical horizontal length scale. Whereas Reference SaltzmanSaltzman (2002) used these dimensionless variables to study the thermodynamic equations, the variables are equally suited to study Equation (1). This leads to a new dimensionless equation

(20)

The mass balance B is prescribed analytically. The ice flux Г ■ (HU’) therefore contains all numerical errors. The factor in front of the ice flux, henceforth called ‘flux number’, determines the relative effect of these numerical errors in the total ice-thickness change. On the one hand, if the flux number is small, the ice-thickness change is dominated by the mass balance. This implies that the results do not suffer much from the numerical problems as shown in Figures 3–9. On the other hand, if the flux number is high, the ice flow, rather than mass balance, dominates the ice-thickness change, which implies that the numerical errors may cause serious problems in the model results. For wet areas, for example, a is large and the mass balance is dominant. For dry areas, however, the numerical effects in the calculation of the ice flux dominate.

Figures 3–9 show that we need extremely small gridcells even when the modelled ice sheets are large and smoother than small glaciers. For small glaciers, the flux number L_vU_c/L_ha will be relatively small (compared with the value on large, growing ice sheets). For example, a typical ice velocity for a mountain glacier is approximately 100 ma^–1. Typical mass balances have an order of magnitude of 1–10 m a^–1, and a typical size for a mountain glacier is 10 km with a typical thickness of 100 m. This results in a typical flux number of 0.1.

We also need to distinguish between the centre and the margin of ice sheets. In the centre, the ice velocities are small, <1 ma^–1. The mass balance in the centre is small too, around 0.1 ma^–1. Together with a typical ice thickness of 1000m and horizontal dimensions between 100 and 1000 km, this results in values for the flux number between 0.01 and 0.1. This is comparable to or smaller than that for mountain glaciers. Thus, in the centre of the ice sheet the numerical effects in the ice-flux term are not important. However, near the ice margin, the typical velocity is of the order of 100–1000 m a^–1 and the mass balance around 1 ma^–1, resulting in values for the flux number between 0.1 and 10. This is larger than for a typical mountain glacier. Hence, compared with mountain glaciers, numerical effects in the ice-flux calculation are much more important in the ice margin of large ice sheets. These marginal zones are also the sensitive locations that determine ice growth. Note that for a typical size LҺ of 100 km, the flux number is higher than for a typical size of 1000 km. This confirms the test results that the numerical problems mostly occur during the growth phase of an ice-sheet model. Results of numerical tests for glaciers cannot be extrapolated to ice sheets; we must expect different numerical behaviour for ice sheets and glaciers.

6. Conclusions

We assessed the effects of varying spatial discretizations through a series of experiments with simplified geometries in a one-dimensional and a two-dimensional vertically integrated ice-flow model based on the SIA. This approximation is widely used by climate researchers, making it essential to understand its numerical behaviour. The focus was not on resolving fine structures, but on the global behaviour of large, smooth ice sheets, where we included the feedback of the mass balance with height. Note that we cannot extend our results directly to situations where the mass balance does not agree with such a parameterization. Note also that these results are only valid for stationary grids. We expected each of the discretizations that we used to perform well, given the size and geometry of the ice sheets. However, after testing for the sensitivities to initial conditions and the individual model parameters, we found that not all discretizations performed well: