2.1. A family of similarity solutions with non-zero accumulation

We now consider certain exact, time-dependent and radial solutions to Equations (1) and (2) with H ≥ 0, called similarity solutions. (See Reference BarenblattBarenblatt (1996) for examples of similarity solutions in other fields.) It turns out that such a solution exists if accumulation is proportional to thickness and is inversely proportional to time. In particular, as will be shown below, there is a similarity solution H = Hχ to the equation

where χ is an accumulation parameter explained below. Such solutions are useful for verifying a numerical code’s ability to respond to a time-dependent accumulation (source) term.

We seek a solution of the similarity form (ansatz)

where r is the radial distance r ² = x ² + y ². (Regarding derivatives in radial directions, recall that Δ _g = (dg/dr)ř and ∇ • (gř) = r^{— 1}d(rg)/dr if g = g(r) and ř is the unit radial vector.) The constants α, 3 and the function φ are to be determined.

Following Reference HalfarHalfar (1983), one might try to start the derivation by supposing H(r, t) = aF²(t)G(brF(t)) for unknown constants α, b and unknown functions F, G. This leads to ordinary differential equations for F and G by separation but presupposes α = 23 which excludes a nonzero-accumulation solution. Similarly, the elegant derivation of Halfar’s solution in Nye (2000) starts from the assumption that the volume is constant;this is not true of the λ = 0 solutions here.

Instead, substitute H = H _λ and h =(1 – f)H_λ into Equation (3). The course of the derivation that follows is best seen by noting that the conditions for eventual success are, most importantly,

Note that these equations completely determine α and 3 as functions of n and λ.

The role of the first of conditions (5) is to allow the elimination of t from the equation found by substituting the similarity form. The result is an ordinary differential equation

where (•)⁰ denotes d/ds. Equation (6) can be integrated once using the second of conditions (5) to write (α + λ)sφ + 3s² φ’ = 3(s²φ)’; the constant of integration is seen to be zero because φ(0) and φ’(0) are finite. If φ’ < 0, i.e. if flow is in the positive radial direction, then

This equation is separable. Note φ (s) < 0 for some s > 0 implies 3 > 0.

For the next integration, suppose 3 > 0 and suppose that at time t ₀ > 0 the dome thickness at r = 0 is H₀. Integrating Equation (7) yields

Let R ₀ be the margin radius at time t ₀, that is, such that We have the relation

between t ₀, H ₀ and R ₀ (cf. equation (19) in Nye (2000)). The time t₀ gives the characteristic scale for substantial changes in the similarity solution.

In summary,

is a solution to Equation (3) where

and t ₀, H ₀ and R ₀ satisfy Equation (9). The dome height is H _d(t) = H _λ(0, t) = H ₀(t/t ₀)^—α as a function of time while the margin radius is R_m (t) = R ₀(t/t ₀)³. The accumulation function corresponding to H_λ is M _λ(r, t) = λt ^–1 H _λ(r, t), with dome value M_d (t) = λH ₀ 1 1^–1 (t/t₀)^—α.

Formula (9) for the absolute time-scale, t₀, is related to the ‘t = 0’ value of our solution. By choosing H₀ and R₀ we have determined the t! 0+ volume of the ice, for λ in the appropriate range (see Table 1) and also the time t₀ for the dimensions to evolve to H₀ and R₀. For the purposes of numerical verification, it is convenient to parameterize by H₀, R₀ and accept t₀ as the time-scale. (Compare comments by Nye (2000) and Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson (2004) on the modelling significance of t₀ for the Halfar solutions.)

Table 1.

Initial and ultimate conditions for solutions H_λ (r,t); see Equations (10) and (11). Equation (13) defines H_d and M_d

To understand the relation between the value of λ and the initial (t →0+) and ultimate (t →+∞) conditions of these similarity solutions, we non-dimensionalize t ~ t₀, r ~ R₀, H ~ H ₀ and write Equation (10) as

in the new coordinates. Also

are dome height, margin radius and dome accumulation rate, respectively, in the non-dimensional coordinates.

From Equations (11) and (13) we have the various cases in Table 1.

Seeking only realistic solutions, we recall Equation (7), and note that if φ’(s) < 0 for some s > 0 then 3 > 0 and λ > —(2n + 1)^–1. It follows that R_m (t) → 0+ as t → 0+ and R_m(t) →+∞ as t ! +1 A ‘delta function initial condition’ description is only appropriate for the solutions H _A with — (2n + 1)^–1 <A< 2(n + 1)^–1; this range includes the Halfar A = 0 solution, of course.

Similarity solutions of Equation (3) for —∞ <A < – (2n + 1)^–1 do exist. They satisfy φ’> 0 and β < 0. Equations (7), (8), (9), (10) and (12) must be changed appropriately. These solutions have infinite ice-sheet extent (thus R ₀ is not defined) and infinite ice volume at all times. Flow is towards r = 0 into an interior ‘cup’ (except when β = 0, then thickness is constant in r). In any case, note that when A < 0 there is increasing ablation with increasing height. These curiosities are, perhaps, not useful for modelling or verification.

Figure 1 shows the profiles of H _A for n = 3 at t = 0.11₀, t₀ and 10t₀ and for A = –1.7, 0 and 1.5. The A = 0 solutions have time-dependent volume;they are all new.

Fig. 1.

Three cases of the similarity solution (Equations (10) and (11)). (a) The Reference HalfarHalfar (1983) solution (A = 0) with constant volume and zero accumulation. (b) The A = –1/7 solution with height-dependent ablation and decaying volume; the margin does not advance. (c) The A = 1.5 solution with height-dependent accumulation and increasing volume; the margin advances rapidly.

An interesting similarity solution is obtained by setting A = 5. In this case the dome thickness H_d(t) grows linearly in time, from zero at t = 0, and the dome accumulation M_d(t) for some sliding coefficient μ > 0. Various physical motivations for such laws exist (Reference MacAyealMacAyeal, 1989; Reference PatersonPaterson, 1994). Simple modification of the derivation below allows construction of exact solutions with basal sliding given by any reasonable function μ_b = F(τ_b). is constant as a function of time. Test C in section 3 is based upon this case.

The concatenated solution H = H_A>₅ for 0 < t < t ₀ and H = H_A<₀ for t > t ₀ forms a simple analytical ‘complete lifecycle’ model for an isothermal shallow ice sheet. It starts with zero ice and experiences a period of linear-in-surface-elevation accumulation. The accumulation switches off at t = t ₀ (when A = 0) or switches to ablation (when A < 0) and the sheets decay to zero thickness as t →+∞. There is no discontinuity of H at t = t ₀ because H _A(r, t₀) is independent of A.

2.2. Solution with basal sliding in an ice-stream-like region

In this subsection we show that compensatory accumulation functions can be defined to create essentially any profile we choose. The resulting exact solutions can be used, for example, to test the ability of numerical codes to deal with a spatially varying accumulation source term and/or basal sliding (for a broad range of sliding laws).

The Reference BodvarssonBodvarsson (1955)-Reference VialovVialov (1958) radially symmetric ice sheet with Vialov-Nye margin H = 0 at r = L and constant accumulation M ₀ (Equation (17) below) is a well-known exact steady solution to the isothermal shallow-ice equation. In this subsection we construct an exact steady solution with basal sliding in an ice-stream-like sector, based upon the Bodvarsson-Vialov solution.

In particular, consider the equation

where q _f is given by Equation (2) as before and ub is the basal sliding velocity, assumed proportional to a power p of the basal shear stress (τ_b = (τ_xz, τ _yz) =—pgHVh in the shallow-ice approximation):

Our exact solution will have non-constant sliding parameter μ = μ(r, θ), in polar coordinates, varying continuously from zero in the regions of frozen base to a maximum μ_max within the ice-stream-like region. In particular, we suppose the region is a sector q < r < r₂, θ₁ < θ < θ₂, and we suppose such a sector appears in each quadrant of the ice sheet by reflection across x = 0 and y = 0. (Such symmetry produces a convenient verification test (see section 4) but is inessential to the derivation here.) Supposing r = L is the position of the margin, we require q > 0 and r ₂ < L for non-singularity of certain quantities in what follows. In test E (section 3) we choose

within the region, and μ = 0 otherwise, for (relative) simplicity. Note that μ((q + r₂)/2, (θ ₁ + θ₂)/2) = μ_max.

Figure 2 shows the region of non-zero sliding and contours of the sliding parameter. Note θļ = 10, θ₂ = 40, r1 = 200 km and r2 = 700 km, where L = 750 km in this map plane view; these are the parameters used in test E. The angles θ1 and θ2 have been chosen to produce a relatively generic sliding sector, not specially aligned to the finite-difference grid in section 4.

Fig. 2.

Contours of sliding parameter, μ, for an exact solution with basal sliding.

For the figures and for verification we choose p = 1 and μ_max = 2.5 × 10^–11 Pa^–1 ms^–1. These values produce, for example, a velocity of 80 m a^–1 at a basal shear stress of 100 kPa (1 bar).

To construct an exact solution to Equation (14) we find a compensatory accumulation M_b as follows. First, H(ŕ) is given by a predetermined function, in this case the Bodvarsson–Vialov profile

with M₀ a positive constant. Now let

Then H_V solves Equation (14) with M = M₀ + M_b and u _b given by Equation (15) with p = 1. This can be seen by substitution. Note that Mb is bounded and piecewise-differentiable but discontinuous; it is well known that discontinuous source terms pose no difficulties for diffusion equations (Reference EvansEvans, 1998).

The main point, however, is simply that HV is the exact solution to Equation (14) because M_b has been specifically chosen to make it so. Figure 3 shows H and M = M₀ + M_b along the center line of the sector. The accumulation is greater than M0 upstream in the region, and is negative (ablation) downstream, in order to compensate for the basal sliding. That is, faster flow from basal sliding exports ice from the upstream region and piles it downstream; accumulation must compensate. Details of the calculation are addressed in the Appendix; test E in section 3 is based on H_V and M_b.

Fig. 3.

Profile of ice thickness, H, (solid curve) and accumulation, M, (dashed curve) along the center line of the sector.

2.3. Time-dependent solution with compensatory accumulation

The technique used in section 2.2 to generate an exact solution with basal sliding has more general applicability. Many exact solutions suitable for testing, including the time-dependent solution we now describe, can be found this way. This solution, which has profile H_p (Equation (24)), is a perturbation of an apparently new steady-state solution H_s (Equation (21)). We construct a new steady solution because no existing smooth exact solution has its margin in an ablation zone. We seek a solution subject to a natural (ablation zone) boundary condition because, in the following verification test, margin drift is an indication of numerical error.

(Our derivation will produce a compensatory accumulation from the derivative of the flux, which therefore should be smooth. For the 2-D steady case there is an exact solution with piecewise constant accumulation/ablation (Reference PatersonPaterson, 1994, p. 245) which has been used in verification (Reference Calvo, Diaz, Durany, Schiavi and VazquezCalvo and others, 2002b), but that solution does not generalize to the radial case with a convenient expression.)

Thus, for the moment, we seek a steady radial solution to the isothermal shallow ice-sheet equation with margin at r = L situated in an ablation zone. With the scalar flux is Assuming Integrating from an arbitrary radius r to a margin at r = L, where H(L) = 0, we have

The purpose of this recapitulation is to note that an analytical expression for H(r) exists if the nth root of the flux function q(r) has an analytical integral. Thus we are invoking a technique for producing myriad steady analytical solutions. Note, however, that if r = 0 is in an accumulation zone and if the margin r = L is in an ablation zone then q(r) will grow linearly at r = 0 from a value of zero and will decay linearly at r = L to zero. One possible flux function q(r) with all of these properties is

for some C, where s = r/L. A dome thickness H₀ at r = 0 determines The corresponding thickness function is

where s = r/L. From Equation (20) we find an accumulation function M_s(r) = r^–1 (d/dr)(rq(r)) by differentiation, that is,

Note M_s (0⁺) =2C/L and M_s (L^—) = —C/L by limit arguments.

The particular functions H_s and M_s together form a new steady solution with dome in accumulation region and margin in ablation region (see Fig. 4). Both H_s and M_s are smooth on 0 < r < L, and they satisfy 0 = M_s + exactly. They will be incorporated into test D in section 3.

Fig. 4.

The profile (solid curve) of a steady exact solution with margin in an ablation zone. Accumulation function shown (dashed) with location of equilibrium line noted. Dimensions and constants as in test D; see section 3.

We now return to the construction of an exact time-dependent solution. Suppose that we add to the steady profile H _s a sinusoidal-in-time perturbation in the vicinity of the equilibrium line. Suppose, for example, a perturbation

if 0.3L < r < 0.9L, and P(r, t) = 0 otherwise; period of the perturbation and Cp its magnitude. Tp is the Let

Now, H = H _p is an exact solution to the time-dependent equation

where h_p = (1 – f)H_p. That is, one makes the constructed function H_p an actual solution to Equation (25) by

with compensating accumulation defined by

where h_p = (1 – f)H_p. That is, one makes the constructed function H_p an actual solution to Equation (25) by with compensating accumulation defined by compensating with precisely that additional accumulation which makes it a solution.

Figure 5 shows thickness profiles and accumulation over one period T_p . Note that M is bounded but discontinuous because P has a continuous first, but not second, spatial derivative. During verification, the accumulation function is extended by –0.1 ma^–1 in the region beyond the exact margin position so as to induce mild negative feedback for margin drift (see section 4).

Fig. 5.

(a) The envelope of H_p (Equation (24)), a perturbed steady solution. (b) The corresponding accumulations M = M_s + M_c (Equations (22) and (26)). Profiles are shown at t = 0, Tp/8,... 7Tp/8, that is, at eight equally spaced times in one period of the perturbation.

In a map plane view of the ice sheet the perturbation and corresponding non-zero values of M_c are confined to an annulus away from r = 0 and r = L. Because of the large variation of M over time, and because of the discontinuous nature of M, the resulting numerical verification test is quite challenging. Details of the calculation of M_c, as needed for test D in section 3, are given in the Appendix.

4.1. Numerical method summary

The most general form of the equation of continuity considered in this paper is

where q _f is given by Equation (2) and u _b is given by Equation (15) with p = 1. For simplicity we use f = 0 so H = h.

We suppose a rectangular grid (xj, yk, t/) with fixed spacing ∆x, ∆y, ∆f. Let Hjy approximate H(xj, yk, t_l) where H(x, y, t) is the exact solution to a particular initial/boundary-value problem. Our explicit finite-difference scheme for Equation (27) is

where, defining and recalling that μ is the spatially varying sliding coefficient,

The quantity in square brackets in Equation (29) is the

numerical ‘diffusivity’. The staggered grid values of H and α

respectively, with appropriate modifications at the other staggered gridpoints. The ‘staggered grid’ flux Qy+₁//₂ is computed analogously. The deformational part of the Q^x , Q^y formulas is exactly as in Reference MahaffyMahaffy (1976) and is described as ‘type I’ in EISMINT I (see also Reference Hindmarsh and PayneHindmarsh and Payne, 1996).

To implement the natural boundary condition H > 0, tentative thicknesses Hj, k, _l+₁ are found from Equation (28)

using values at t = t_l . If Hj, _k, _l+₁ < 0 then Hj, _k, _l+₁ is set to zero, otherwise the value is kept (Reference Hindmarsh, Greve, Ehrentraut and WangHindmarsh, 2001). In the case of a Vialov-Nye boundary condition, those gridpoints predetermined as outside the ice region are set as H = 0 at every time-step.

Overall dimensions -L_x < x < L _x, -L_y < y < L _y and run time T are supposed fixed. In our verifications the spatial domain is square so we suppose L_x = L _y, N = N_x = N _y, ∆x = 2L_x/N, ∆y = 2L_y//N and ∆t = T/K; there are (N + 1) gridpoints in each spatial direction.

Explicit methods, while easy to implement, are well known to be subject to stability restrictions (Reference Hindmarsh and PayneHindmarsh and Payne, 1996). To analyze numerical stability, one must temporarily linearize the problem. In particular, one must fix the diffusivity D =(1 - f)ГHⁿ+²1Vh^n-1 at some time-step so that it no longer evolves in time. Steps of our explicit scheme for ∂H/∂t = M +V ■ (DVH) are then subject to the restriction (maxD)∆t/∆x² < C for some constant Cof order one (Reference Morton and MayersMorton and Mayers, 1994). In a fixed total time run, a doubling of N therefore requires an increase in the number of time-steps Kby a factor of 4. On the other hand, the work of a run is proportional to N²K for this, and any other, explicit method. Therefore the work increases by a factor of 16 when N is doubled. Note that a fastest possible ‘refinement path’ (Reference Morton and MayersMorton and Mayers, 1994) in ∆x, ∆t space will therefore have ∆t ≈ C∆x ² for some C’ (see Fig. 6).

Fig. 6.

Step sizes ∆t and ∆x used for each test; ‘refinement paths’. Values of ∆x correspond to N = 30, 60, 120 and 240 grids. Values of ∆t chosen as large as allowed by stability.

This restriction (max D)∆t/∆x² < C does not directly apply to the actual (non-linear) evolution, but it has practical value. For instance, in test A we empirically determine C ≈ 0.15 from several runs while the linear stability analysis of Reference Hindmarsh, Greve, Ehrentraut and WangHindmarsh (2001) predicts that C = (2n)^–1 = 1/6 from the slab flow. The restriction is reliable enough to be made the basis of an adaptive time-stepping scheme, though in this paper we used a fixed time-step in each run. Tests B and C would show distinct benefits of adaptation: in the first case, ∆t would lengthen as the sheet became shallower;∆t would shorten in the second as the sheet thickened.

To compute ice-sheet volume we have chosen the numerical method

This method is a reasonable choice because it is exactly the two-variable trapezoid rule. However, since H has unbounded gradient at a margin, all volume methods based upon polynomial approximation (trapezoid, Simpson’s, Gauss’...) will suffer from volume errors at the margin. The trapezoid rule at least makes minimal presumptions on differentiability.

It turns out that volume method (31) is compatible with the flux conservation of our explicit method and with our implementation of the H > 0 natural boundary condition. By ‘compatible’ we mean that this numerical scheme conserves V_num to within the rounding error (less than one part in 10¹⁴) in accumulation-free cases. See the verification with test B below.

4.2. Verification and grid refinement study

Because the tests in this paper have radial symmetry, we have saved time by performing all computations in the quadrant 0 < x < L _x, 0 < ỵ < L _ỵ. (In test E there is symmetry of accumulation across both x = 0 and ỵ = 0.) We impose reflection boundary conditions along the lines x = 0, ỵ = 0; the origin (x, ỵ) = (0,0) is the center (‘dome’) of each ice sheet.

In each time-dependent test, the horizontal extent Lx = Lỵ was chosen large enough so that exact and numerical ice sheets never advance to within one gridpoint of the outside boundaries x = L _x, ỵ = L_ỵ of the computational domain. In particular, the natural boundary condition H > 0 in tests B-D was always imposed at gridpoints away from the x = L _x, ỵ = L_ỵ boundaries.

The runs were done with fixed ∆t, which was adjusted by hand to be close to the stability limit for the particular problem. From the grid refinement study we are convinced that, for time-dependent problems, choosing ∆t small enough to satisfy stability requirements makes the O(∆t) error of our first-order-in-time explicit method negligible relative to errors in spatial derivative approximations.

Runs using N = 30, 60, 120 and 240 were performed for each test;N = 30 is the EISMINT I choice and N = 60 is the EISMINT II choice. Figure 6 shows the refinement paths, in ∆t, ∆x space, used for each test.

Figures 7 and 8 summarize the results. We show maximum errors, i.e. max|H_num – H_exact| over all gridpoints (x, ỵ) = (xj, ỵk) at the final time t = t_f, and dome errors |H_num(0, 0, t_f) – H_exact(0,0, tf)| as functions of N, respectively. Best-fit lines, i.e. linear regression in (logN)× log(error) plane, and their slopes (‘O(N ^—p)’), are added to suggest the rate of convergence, discussed below.

Fig. 7.

Maximum error for each test; behavior of the errors under grid refinement. Errors for tests A and E are indistinguishable at this scale.

Fig. 8.

Dome error for each test; behavior under grid refinement.

Dome errors are surrogates for the ill-defined concept of ‘away-from-the-margin errors’. For example, Figure 9 shows contours of the errors for tests A and B with N = 60. The largest errors are confined to the vicinity of the margin, while error is relatively uniform and small in the interior. Dome errors are substantially smaller than maximum errors and they converge to zero faster. An explanation appears in section 5.1.

4.3. Further remarks on test A

Test A is a steady-state solution for which a run of 25 000 model years is sufficient for convergence to numerical equilibrium from the exact profile as an initial value (cf. EISMINT I). The error H_num – H_exact therefore compares numerical to exact equilibrium.

We believe our numerical method is completely typical of EISMINT I (see section 5.3) but the maximum errors for test A in Figure 7 are larger than those suggested by figure 2 in EISMINT I. In EISMINT’s figure an error apparently on the order of 300m is shown at the last interior gridpoint (near the margin);only a 26m dome error is actually reported. (We necessarily refer to the EISMINT I 2-D results;three-dimensional (3-D) errors cannot be known without an analytical solution to the EISMINT I 3-D experiments.) As pointed out earlier, our radial tests are much more challenging than the EISMINT I steady-state experiment, in that circular margins are not aligned to the grid.

Figure 9a shows the spatial distribution of the error in test A. From Figure 7 we have a suggested O(N^{–0 33}) rate of maximum error convergence. By contrast, the centered difference approximation of the V • qf term has O(N^–2) local truncation error (found from Taylor’s theorem assuming four continuous derivatives (Reference Morton and MayersMorton and Mayers, 1994)), and thus the rate O(N^–2) is the best possible global (discretization) error achievable by a numerical scheme using the centered formula. Note that in easier problems than ice flow, local errors tend to accumulate as one moves in from the boundary to give larger global errors. For ice flow with a grounded margin we see the opposite situation: errors are large near where H is exactly imposed (the boundary condition) and are small at the dome. In this case, unbounded derivatives obviously degrade the best behavior (see sections 5.1 and 5.2).

Fig. 9.

Spatial dependence of the errors in tests A and B when N = 60. Error is defined as H_num – H_exact, which is positive everywhere using our type I scheme. Only one-quarter of the sheet is shown owing to symmetry. (a) Test A: maximum error of 650 m near margin; interior of sheet has errors in range 30–70m. (b) Test B: maximum error of 170 m near margin; interior of sheet has errors in range 1–5 m.

4.4. Further remarks on test B

Test B is the zero-accumulation similarity solution of Reference HalfarHalfar (1983), describing the spreading of a delta mass of ice driven solely by gravity. As shown in Figures 7 and 8, the errors are substantially smaller than in test A and converge to zero faster. They are smaller, in part, because the dome thickness at the final time is smaller (2283.4 m as opposed to 3600m), but primarily because the margin gradient singularity is milder (see section 5.1). Figure 9b shows the spatial dependence of the error at the final time when N = 60.

In Figure 10 we examine the numerical mislocation of the margin. We show gridpoints where the exact and numerical solutions disagree about whether ice is present at the final time. It turns out that, regardless of the value of N in test B, the width of the arc of mislocated-ice gridpoints is always three. That is, the arc can always be crossed by ‘stepping on’ at most three of them. We do not know the explanation for this particular number, but we may conclude that the method is doing an acceptable job of locating the margin. The location error is linear in ∆x, that is, O(N^–1), while the O(N ^_0•³⁹) thickness error near the margin decays much more slowly.

Fig. 10.

(a) Gridpoints on an N = 30 grid where the exact and numerical solutions disagree about whether ice is present at the final time in test B. (b) The corresponding result for N = 240. The arc is three gridpoints wide in both cases. The arc is also three gridpoints wide when N = 60 and N = 120 (not shown).

The exact solution in test B has constant volume over time. The finite-difference method used here also yields constancy of the numerical volume computed by Equation (31). Indeed, when the initial and final numerical volumes were compared, the absolute errors were 2, 0, 0.5 and 7m³ for N = 30, 60, 120 and 240, respectively. The numerically computed volume in question is approximately 4 × 10⁶ km³ and thus these errors are all smaller than 10^–14 in a relative sense. Note that the numerical volume (Equation (31)) is not a perfect estimate of the exact volume. For instance, formula (31) approximated the exact (analytical) volume in test B only to within 0.07% for N = 60. Rather, the numerical choices made here (finite-difference method, numerical boundary condition and numerical volume method) are consistent in terms of volume conservation. Local volume conservation (flux conservation) for the finite-difference method is reasonably obvious;the property of our numerical scheme illuminated by this experiment is the (numerical) conservation of volume at the margin with a natural H > 0 boundary condition.

4.5. Further remarks on test C

Test C is a new non-zero-accumulation similarity solution with zero initial condition and linearly increasing dome height over time. The accumulation function has a constant (in time) maximum value of about 1.2 m a^–1. Because M is proportional to H, the region of non-zero accumulation grows as the ice sheet advances.

Because diffusivity D depends on the n + 2 power of H, as thickness increases stability can ‘suddenly’ become an issue at some time-step. This test is therefore well suited to testing adaptive time-stepping (although we did no such test here).

The accumulation function M is irregular in the sense that it goes to zero with unbounded derivative at the margin. We use the final volume to measure the degree to which the numerical method accurately accumulates ice over time. The result is that the final numerical volume V_num converges to the exact value at a relatively rapid O(N^–2–41) rate.

In this paper, we have not taken full advantage of the verification opportunities offered by test C, namely, to verify a shallow numerical model of the ‘mass-balance altitude feedback effect’ (Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson, 2004) in the case where accumulation is linear in surface height (thickness). In our tests we regarded the exact M(r, t) as a predetermined function of space and time. However, since M = \t ^–1H, one can instead compute the accumulation M_num = 5t^H_num at each gridpoint based on the approximate value H_num. Such a calculation requires using the exact solution H(r, t) at positive t for an initial condition;an ice sheet of zero thickness does not grow if accumulation is linear in thickness.

4.6. Further remarks on test D

Test D is a time-dependent, sinusoidal perturbation of a new steady-state solution. After each period of T_p = 5000 model years the exact solution returns to the same thickness profile and thus to the same volume. It has margin in an ablation zone and the margin boundary condition is simply H > 0. The numerical margin is allowed to drift. We impose a 0.1 ma^¹ ablation rate in the region outside the exact margin for mild negative feedback.

Now, because of oscillatory and irregular accumulation (Fig. 5), this test is ideal for measuring the degree to which the numerical method does not conserve volume. Figure 11a shows that as the grid is refined the numerical volume at the final time of 25 000 model years seems to converge to the exact value. Figure 11b shows the numerical and exact volumes over time; the N = 30 and N = 240 curves are omitted for clarity. Note, in particular, that no apparent phase error has occurred, and we believe this is because our time-steps are relatively small (for stability reasons).

Fig. 11.

Test D. (a) Relative volume error at 25 000 model years. (b) Exact (solid curve) and numerical volumes over time for N = 60 (dotted curve) and N = 120 (dashed curve).

4.7. Further remarks on test E

Test E is a steady-state solution which adds linear basal sliding and compensatory accumulation, in a sector, to test A (see Fig. 2). Again, a run of 25 000 model years was regarded as sufficient to produce numerical equilibrium, given the exact profile as an initial value.

The errors are maximum near the margin as in test A, and indeed dominate the error made in approximating basal sliding. We conclude that our numerical scheme for basal sliding is completely compatible, in terms of numerical difficulties, with difficulties already present in the frozen base case.

A contour map of the test E thickness error appears in Figure 12. It shows that the two processes of basal sliding and compensatory accumulation do, as designed, essentially cancel each other. By comparison with Figure 9a, we see that a slight asymmetry of the error contours, especially the 30 m contour, is the only visible indication of numerical non-cancellation.

Fig. 12.

Error H_num – H_exact in test E with N = 60. Near-margin and interior errors are nearly identical to those in test A (cf. Fig. 9a). Compensatory accumulation and basal sliding essentially cancel as intended.

5.1 On the approximation of margins

We now redo the regression on the maximum error graph (Fig. 7) by discarding the N = 30 case as an ‘outlier’. The result (Fig. 13) is that the three margin shapes - steady accumulation zone/Vialov-Nye (A and E), moving similarity/natural (B and C), and steady ablation-zone/natural (D) – determine the rates of convergence of the maximum error everywhere on the sheet. Unfortunately, we also see by the poor linear fit that for the Vialov–Nye margin cases (tests A and E) the rate of convergence may be slower than these linear fits suggest. It is conceivable, at least, that the error does not approach zero as O(N^—p) for any p > 0, even in theory. In practice, it is clear that no achievable grid exists for our finite-difference method that would give a maximum error of, say, 50 m in test A. By contrast, a maximum error of 50 m is achievable in test B using N = 960, i.e. two more doublings of N = 240, if the convergence rate continues.

Fig. 13.

The convergence of maximum error is influenced by the type of boundary condition imposed. Here the linear regression is to the N = 60, 120 and 240 values only. Maximum errors in tests A and E coincide at the resolution of this figure.

The source of large maximum errors is the difficulty of numerically approximating ice-sheet margin shape. It turns out, in other words, that large maximum errors do not come from some elaborate non-linear mechanism or from the misapplication of numerical boundary conditions, but rather from the ‘simple’ issue of singular margin gradient.

Certainly the shallow approximation is invalid as a precise continuum description of ice flow as one approaches the margins, where neglected stress components become important. However, the shallow-ice approximation is both a good approximation for ice-sheet-scale behavior and is more practical at ice-sheet scales than the full system equations (Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson, 2004), so the extent of disagreement between continuum solutions and their numerical approximations is important. In particular, we will explain why large numerical errors at margins do not ‘corrupt’ accuracy over the interior of the ice sheet.

The known (asymptotic) margin shapes correspond to three specific exponents. The mildest are ablation-zone steady margins with shape H(x) ~ x^1/2 where x is a generic coordinate describing distance into the ice sheet from the margin (Reference FowlerFowler, 1997). (These are ‘mild’ in the sense that x^p is of increasing differentiability at x = 0 as p increases.) The next mildest profiles appear in the similarity solutions H = Hχ of section 2.1. The shape is independent of the accumulation parameter χ: H_χ (x) ~ x^n/(2n+Г). Finally, accumulation-zone steady margins (Vialov–Nye margins) have the steepest shape H(x) ~ x^n/(2n+²) (Reference FowlerFowler, 1997).

We assume other time-dependent margin shapes are possible. We also suppose, however, that the x^n/(2nþ2) power of a Vialov–Nye margin represents the most singular shape of interest to modellers.

The importance of these powers can be seen by evaluating the flux, given by Equation (2), as x →0^þ, in these three cases. The result is that q(x) ~ x¹, x^n/(2n+1), x⁰, respectively. In the last case, q is not even continuous at x = 0 because, presumably, q(x) = 0 for x outside the ice.

All finite-difference methods which depend on the standard O(∆x²) centered-difference formulas for derivatives can be regarded as being implicitly based upon global piecewise linear approximation (although they are derived by Taylor series arguments). This can be seen clearly in certain examples, for instance, the perfect correspondence between the standard finite-difference and finite-element treatments of the simplest Poisson problem on a rectangle (section 1.2 of Reference JohnsonJohnson (1987)). Finite-element methods are, of course, explicitly based upon piecewise linear (or, at least, piecewise polynomial) approximations.

Consider the piecewise linear approximation of an ice-sheet margin with asymptotic shape H(x) = x^p , 0 < p < 1, as in Figure 14. Suppose that a gridpoint coincides with the exact position of a margin. The next gridpoint is ∆x up-glacier and the error between the exact H and its linear approximation is E(x) =x^p — (line) =x^p—(∆x)^p–1x. The maximum E_max = max E(x), for 0 < x < ∆x, is at x = ∆x • p^{1/(1— p)} and has magnitude

for 7 = p^p/(1—p) – p ^1/(1—p). A more complete version of this analysis shows that even as the position of the last on-ice gridpoint is varied the approximation error is still O(N^—p ).

Fig. 14.

Linear approximation of a margin H(x) = x^p gives error E (x) = C (∆x )^p .

Our hypothesis is that Equation (32) predicts the numerical error in the vicinity of the margin, and thus the maximum error on the whole ice sheet (which we hypothesize always occurs at a point near the margin). Table 4 compares the prediction of this hypothesis to observed global maximum errors in the tests of section 4. The correspondence is close for the (steady) ablation-zone and similarity margins, but for the Vialov–Nye (steady) accumulation-zone margin the numerical errors are even larger than predicted by Equation (32). For natural margins there is a high-quality approximation of a continuum quantity evolving ‘in the background’, specifically of the transformed thickness η = H⁽²ⁿ⁺²’^)/n , considered in section 5.2, and Equation (32) really determines the error. For the Vialov–Nye margin something is going wrong even for the η-evolution, namely that η is not continuously differentiable at the margin.

Table 4.

Margin shape predicts global error decay rate; Glen exponent n = 3

In any case, we believe this is strong evidence that margin shape, and not more optimistic local truncation error predictions from Taylor’s formula (e.g. Reference WaddingtonWaddington, 1981), controls the global discretization error in a shallow-ice-sheet simulation. As noted, such Taylor-series arguments predict O(N^–2) local error which accumulates as one moves in from the margin.

5.2. Regularization of the margin by a transformation

Following Reference Calvo, Diaz, Durany, Schiavi and VazquezCalvo and others (2002b), we let η = H^(2n+2)/n for Glen exponent n. Equations (1) and (2), with u _b = 0 and f = 0 for simplicity, are equivalent to for Г = [n/(2n + 2)]ⁿГ. Equation (33), with its natural boundary condition η > 0, describes an evolution equivalent to the original ice equation.

The η-margin is much less steep than the H-margin. To see this, let n = 3 for simplicity. Margin shapes H(x) ~ x^1/2, x³¹⁷ and x^3/ŭ are transformed to η(x) ~ x^4/3, x⁸⁷⁷ and x¹,

The n-margin is much less steep than the H-margin. To see this, let n = 3 for simplicity. Margin shapes H(x) ~ x ^1/2, x ^3/7 and x ^3/8 are transformed to Z(x) ~ x ^4/3, x ^8/7 and x ¹ respectively. That is, the ablation-zone steady and similarity margins are transformed to functions η with continuous derivative. The accumulation-zone (Vialov-Nye) margin is transformed to a η-margin with bounded but discontinuous derivative. For illustration, the η-analog of the Reference HalfarHalfar (1983) solution, that is, the 8/3 power of H_λ (r, t) in Equation (10) with λ = 0, is shown in Figure 15.

Fig. 15.

Profiles of η for the λ = 0 similarity solution at times t₀ (dotted curve) and 3t₀ (solid curve). Note there is no discontinuity of ∂η/∂r at the margin; shape is η = x^8/7 if x is distance from margin. (Cf. Fig. 1a.)

We rerun test B doing the calculation in terms of η. The numerical method used for η (Equation (33)) is the direct analog of the explicit finite-difference method described in section 4.1. Figure 16 shows the rate at which the (relative) maximum and (relative) dome errors converge as N increases. (‘Relative maximum’ error is ‘relative to the value of η at the dome’.) For the evolution of η the maximum error and dome error now decay at comparable rates. The maximum error goes to zero at rate O(N^–1•⁰⁷), instead of O(N^–0Л4) for H, because the η-margin is well behaved. The significance of this good performance is, we believe, that errors in H are much greater near the margin than in the interior only because H is the 3/8 power of η and not because of more complicated non-linear effects. The convergence rate of the H dome thickness closely matches that of η. ‘Corruption’ from margin errors does not extend into the interior for either the H or η computation, except to reduce the order of interior convergence from the best conceivable rate O(N^–2) to about O(N^–15).

Fig. 16.

Convergence of relative maximum (solid lines) and dome (dashed lines) errors under grid refinement, for both H-evolution and η-evolution. (Cf. Figs 7, 8 and 13.)

5.3. Methods for improving the EISMINT benchmarks

Having verified our numerical method relative to exact solutions, we consider the continuum problems used by EISMINT as benchmarks. An important criticism of EISMINT I is that the values reported as benchmark results are not numerical results of ‘benchmark quality’. For several examples of benchmark-quality results for problems in fluid dynamics without an analytical/exact solution, see Reference RoacheRoache (1998). We will show how such benchmark-quality results can be extracted for these particular problems. Nonetheless, for verification the use of exact solutions is easier and the sources of error are more understandable than if benchmark-quality numerical solutions are used.

We have reproduced three of the six EISMINT I experiments, namely the ‘fixed-margin experiment in steady state’, the ‘fixed-margin experiment forced by sinusoidal boundary conditions with a period of 20 kyr,’ and the ‘moving-margin experiment in steady state’, to show how the benchmarks can be improved upon. We compute on the EISMINT I grid with N = 30 and then on refined grids with N = 60, N = 120 and N = 240.

Table 5 reports the EISMINT I results and our results for dome/divide thickness H(0), ‘midpoint’ flux (which is at 400 km away from the dome;the margin in the fixed case is at position 750 km) and margin position itself (in the moving-margin experiment). The table gives the EISMINT I 3-D value and range, that is, the mean and standard deviation of the 11 results from 3-D codes, and also the value and range of results from those five codes which, like ours, use ‘type I’ diffusivity/flux finite-difference approximations. We report our results on the N = 30 and N = 240 grids (‘X₃₀’ and ‘X₂₄₀’);though N = 60 and N = 120 results were calculated, they are omitted from the table for clarity. In each case our N = 30 values are in the EISMINT type I range.

Table 5.

A sampling of EISMINT I benchmarks values (columns 3 and 4), our values on the N = 30 rough grid (column 5) and versions improved by grid refinement and Richardson extrapolation (columns 6 and 7). See section 5.3 for last column

Consider the dome/divide thickness H(0) in the fixed-margin and steady case. There exists a correct value for the continuum problem; unfortunately the EISMINT I square geometry means that it is not known analytically. However, we can use grid refinement and Richardson extrapolation (Reference Burden and FairesBurden and Faires, 2001) to get a much more accurate estimate than from the N = 30 grid, or indeed from any fixed grid. Figure 17 shows H(0) values from the N = 30, 60, 120 and 240 grids. These are associated with decreasing values of ∆x, of course. Following the technique of Richardson extrapolation, we suppose that the numerical approximation to H(0) is a polynomial function of ∆x; the convergence results of section 4 suggest limited validity to this assumption as dome errors converge as O(∆x^q) for 1 < q < 2 (Fig. 8). We nonetheless choose the degree of that polynomial to be two (and ignore the N = 30 case;see Fig. 13).

Fig. 17.

Richardson extrapolation gives a better estimate of the correct dome/divide thickness for the EISMINT I fixed-margin, steady-state experiment.

Using data for three values of ∆x (N = 60, 120 and 240), we find the coefficients a, b and c in X = a + b∆x + c∆x² where X is the computed quantity. Evaluation of this polynomial X(∆x) at the ultimately refined grid spacing ∆x = 0 gives the Richardson extrapolation value reported in the table as ‘X_RE’ and shown in Figure 17. The column ‘|diff|’ in Table 5 gives values of |X₂₄₀ - X_RE|. Richardson extrapolation does not come with an error bound, but this difference is a (probably conservative) estimate of the precision of X_RE.

Two details in Table 5 require comment. First, the position of the margin in the moving-margin experiment is a discrete quantity and extrapolation is not appropriate. The exact (analytical) value is available, 579.81 km, and our N = 240 value (581.25 km) is, as expected, better than the EISMINT I value of 600 km. Second, the midpoint flux value in the moving-margin experiment did not behave with sufficient regularity under grid refinement to produce a good extrapolation. This value is sensitive to discrete jumps of margin position. These cases show some of the limitations of Richardson extrapolation of numerical results for producing benchmark values, but codes can nonetheless be tested against such extrapolated results if relevant exact solutions are not available.