Boundary conditions

PDEs (1) and (2) require boundary and initial conditions to complete a (presumed) well-posed mathematical problem.

Because the margin is allowed to move, the boundary condition for Equation (1) used here (and in Reference Bueler, Lingle, Covey and BowmanBueler and others, 2005) is

which is really a constraint on the allowed thickness functions. In fact, it is essential to note that Equation (1) is solved (in general) as a free boundary problem. The map-plane region covered by ice is not predetermined. This free boundary problem is known to be well posed in the isothermal case (Reference Calvo, Dſaz, Durany, Schiavi and VasquezCalvo and others, 2002b). In particular, it is shown by Reference Calvo, Dſaz, Durany, Schiavi and VasquezCalvo and others (2002b) that at a free margin where H = 0 the flux Q is also zero in the correctly (weakly) formulated mathematical problem including inequality (14). Finding the location of the margin is part of the problem. To make up for the added unknown boundary location, at the boundary both the thickness and the flux are equal to zero. (These would be overdetermined boundary conditions if they were regarded as boundary values for a PDE in the classical sense.)

A numerical approximation of this free boundary problem is available in an appropriately formulated finite-element method (Reference Calvo, Dſaz, Durany, Schiavi and VasquezCalvo and others, 2002b). A principled approximation is much less clear for a finite-difference method such as the one used here. That is, it is harder to know in advance that a given finite-difference approximation treats the moving margin boundary condition correctly. We note, however, that straightforward enforcement of inequality (14) at each time-step and at each point of the finite-difference grid gives verifiably reasonable results in the current paper and in the isothermal case (Reference Bueler, Lingle, Covey and BowmanBueler and others, 2005).

The three-dimensional (3-D) region on which Equation (2) applies varies in time, but the boundary conditions are of the classic type. On the upper surface of the ice,

where r ² = x ² + y ² and (x, y) = (0, 0) is taken to be the center of the angularly symmetric ice sheet. In this paper, as in EISMINT II, the surface temperature is a specified function of the map-plane position. For construction of exact solutions (below) we take T _min = 223.15 K, the value used in experiment F of EISMINT II. Note that experiment F has the strongest basal temperature spokes of the EISMINT II nonsliding base experiments. (See Table 1 for values of constants not discussed in the text.) A constant geothermal flux is applied at the base of the ice, z = 0:

There is no bedrock thermal model.

Chosen thickness and temperature functions

We calculate two angularly symmetric exact solutions, denoted tests F and G, which are proposed additions to the suite of exact solutions in Reference Bueler, Lingle, Covey and BowmanBueler and others (2005). Test F is a steady-state solution and test G is a time-dependent perturbation of that steady state. The following solution is described as time-dependent, that is, we describe test G. The derivation of test F is a strict simplification.

We also restrict ourselves, in the remainder of this section, to the angularly symmetric case; cylindrical coordinates r = and z are used. The numerical computations in Results are all on a Cartesian grid, however.

To be reasonably physical, we choose H = H(t, r) based on an angularly symmetric steady-state solution of the isothermal Equation (13), namely the test D thickness function from Reference Bueler, Lingle, Covey and BowmanBueler and others (2005):

where s = r/L. Here r = L is the (constant) location of the margin. This function solves a steady-state version of Equation (13) as shown in Reference Bueler, Lingle, Covey and BowmanBueler and others (2005).

The profile given by Equation (18) is a smooth function of r. The smoothness of H _s(r) is important because of the derivative-intensive exact calculations which (eventually) yield the compensatory accumulation M and heat source Σ_c.

Next, we define a perturbed thickness function

where, concretely,

if 0.3L < r < 0.9L while φ(r) = 0 otherwise, and

The term +φ(r)γ(t) is an annular perturbation with amplitude A _p and period t _p. It also appears in isothermal test D (Reference Bueler, Lingle, Covey and BowmanBueler and others, 2005).

For the steady-state solution test F, we will choose A _p = 0 so the perturbation vanishes. For test G we set A _p = 200 m and t _p = 2000 years. In particular, our choice of A _p in test G is small enough so that ∂H/∂r is always negative for all r and t; a consistent sign for ∂H/∂r is convenient in the computation of compensatory sources.

Next, we define the temperature field for the exact solutions test F and G by

where ν is found from H and T _s by

Then T satisfies boundary conditions (15) and (16), that is, T = T _s at the surface and ∂T /∂z = –G/k at the base.

Illustrations and comments

Figure 1 shows a radial section of the test F ice sheet defined by Equation (19) with temperature field defined by Equation (22). Figure 2 shows the same quantities for test G at t = 500 years, that is, at one-quarter of the perturbation period tp when the perturbation has maximum magnitude.

Fig. 1. Thickness and temperature field (K) in test F.

Fig. 2. Thickness and temperature field (K) in test G, evaluated at t = 500 years.

Equation (22) defines a vertical temperature profile which (unfortunately) is not very realistic. Note that the ice warms as one goes from the surface to the bed, but the magnitude of ∂T/∂z is relatively constant in each column of ice and there is no horizontal advection bulge in the temperature contours.

The compensatory heat source Σ_c, derived from H and T via the full thermomechanically coupled system, is non-zero and thus also non-physical. In fact, Σ_c is larger in magnitude than Σ, though the maximum magnitudes of these functions are within a factor of five.

Because our particular continuum model has a simple Arrhenius term in the constitutive relation (4), and because we do not impose T ≤ T _pmp, where T _pmp is the pressuremelting temperature, we restrict ourselves to cold ice. This is achieved by using a relatively cold surface temperature (the EISMINT II experiment F surface temperature condition, with a range 223 K < T _s < 235 K) and a relatively thin sheet (with dome height H ₀ = 3000 m compared to H ₀ = 3600m in EISMINT II). It follows from Equations (22) and (23) that 223 K < T < 273 K for all points in the ice; see Figures 1 and 2.

A crucial technical point is that T(t , r , z) in Equation (22) satisfies the boundary conditions (15) and (16) and possesses the property that the integrals for horizontal velocity, flux and vertical velocity can be computed analytically. Specifically, because integrals of the form ∫ e_z zⁿ dz can be expressed in terms of elementary functions, and because of the form of the Arrhenius term, we see that functions T with z dependence of the form 1/(C+z) are convenient for purposes of analytical computation. In fact, the integrand in the expression for horizontal velocity (see Equations (4) and (5)) takes the form e^z zⁿ multiplied by a complicated z-independent factor. Thus, the integrals (5), (6) and (7) can each be computed exactly. These admittedly technical issues explain choices (22) and (23).

Interestingly, T defined by Equation (22) satisfies

and thus it can be regarded as having a steady-state profile generated by vertical conduction and downward (vertical) advection with a hypothetical vertical velocity w _hyp = . This may suggest some measure of physical reasonableness, but W _hyp is not the vertical velocity w in the solution which follows. Indeed, T satisfies Equation (17) by design while it satisfies Equation (24) incidentally.

The substitution of H and T into Equations (3) and (5), respectively, and use of the constitutive relation (4), gives analytical formulae for effective stress σ and horizontal velocity U. Taking the horizontal divergence of U and vertically integrating the result (see Equations (6) and (7)) gives analytical formulae for the vertical velocity w and for the term ∇ · Q in Equation (1). Substitution of H and ∇ · Q into Equation (1) gives an analytical formula for the compensatory accumulation M. Substitution of T, U and w into Equations (8) and (17) gives an analytical formula for the compensatory heat source Σ_c.

The formulae derived for M and Σ_c have the following property. If Equations (1), (17) and coupling relations (3–8) are solved using the initial values found from the exact H and T, and if the derived formulae for M, Σ_c are used as accumulation and heat source for the whole coupled system during the solution, then H in Equation (19) and T in Equation (22) solve the whole initial/boundary problem. The previous sentence is at the heart of the matter: we have used the continuum model to find, from chosen H and T, the functions U, w, M, Σ and Σ_c so that all of the coupled equations are simultaneously true. As long as we believe in the uniqueness of solutions of the continuum problem, then we have the exact solution to that problem. We have manufactured a continuum truth against which to compare numerical approximations. In the following section, we will perform such a comparison against the result of a finite-difference scheme.

The calculations summarized above are both tedious and non-trivial, but their detailed form appears in Reference Bueler and BrownBueler and Brown (2006). We can, however, communicate many of the essential features of the exact solutions by illustrating the various coupling quantities in the thermomechanically coupled system: Figures 3–7 show test F, and Figures 7–9 show test G at time t = 500 years.

Fig. 3. Magnitude |U| (m a^-1) of horizontal velocity in test F; the vector field U points in the positive r direction.

Fig. 4. Vertical velocity w (ma^–1) in test F.

Fig. 5. Strain heating Σ(10^–3 Ka^–1) in test F.

Fig. 6. Compensatory heating term Σc (10^–3 Ka^–1) in test F.

Fig. 7. Compensatory accumulation M(m a^-1) in test G (solid line) at time t = 500 years and in test F (dashed line). The equilibrium radius for test F is r = 403.6 km.

Fig. 8. Magnitude |U| of horizontal velocity (m a^–1) in test G at time t = 500 years.

Fig. 9. Strain heating Σ(10^–3 Ka^–1) in test G at time t = 500 years.

In some of these figures the margin is not shown. This is in order to avoid division by zero. All quantities H, T, U, w, Σ, Σ_c and M are bounded all the way to the margin, but in some cases the formulae have 0/0 limits at the margin (Reference Bueler and BrownBueler and Brown, 2006). Numerical errors which might arise by numerically evaluating the exact formulae near these 0/0 limits are avoided in the illustrations here. Similarly, evaluating the exact formulae very close to the (exactly known) location of the margin can and should be avoided when the exact solutions are used for verification. In particular, one should check that a gridpoint is at least a few meters inside the exact margin r = L before evaluating the exact solutions at that gridpoint, in order to determine the numerical error there.

Verification

The exact solutions described above were specifically designed for the purpose of verification. The verification procedure described next is permanently built into our comprehensive ice-sheet model PISM (https://gna.org/projects/ pism/). When we make code changes we redo verification to determine whether numerical results are still close to continuum values.

The data presented here involved substantial grid refinement and represent many thousands of processor hours, but use of exact solutions tests F and G for basic verification should be within the resources of all ice-sheet modelling groups.

We briefly summarize our scheme in Appendix A. The reader not interested in the details of the scheme or in reproducing these results need not read Appendix A but should note that we use an explicit finite-difference scheme for Equation (1), with Reference MahaffyMahaffy (1976)-type evaluation of the diffusivity (see also method 2 in Hindmarsh and Payne, 1996). Furthermore, we use a semi-implicit and first-order upwinded scheme (Reference Morton and MayersMorton and Mayers, 2005) with equal spacing in the vertical, for Equations (2) and (17). Roughly speaking, our scheme uses the simplest consistent choices. The local truncation error (Reference Morton and MayersMorton and Mayers, 2005) of our scheme for Equation (1) is O(Δx², Δy², Δt) while the local truncation error for our temperature scheme is O(Δx, Δy, Δz, Δt); see Appendix B. A major purpose of this paper is to quantify the ways in which actual numerical errors for the coupled free boundary problem relate to these local errors.

To verify our thermomechanically coupled SIA scheme using tests F and G, we initialize with the exact continuum values of all quantities and run the code using the exact values of the compensatory accumulation M and heat source Σ_c at every step. The quantities H, T, u, v, w and Σ are computed numerically at each time-step, and no reference is made to their exact values after initialization except to measure errors at the end of the run. Boundary conditions are held to their exact values at every step. Of course, the numerical solution gradually diverges from the exact values; the results below describe this quantitatively.

The runs reported here were for a fixed time of 25 kyr. This was sufficient time for the scheme to diverge measurably from the exact values. As the exact values are known, there is no reason to run to numerical equilibrium, by whatever standard might be determined. Instead, the actual errors were computed at the final time. Our verification results would not have been significantly different if we had run all computations for 200 kyr.

All runs were completed with adaptive time-stepping, with maximum time-step ▽t limited to 10 years. In particular, Δt satisfied both constraints (25) and (27) at each step; see Analysis below. (Inequality (25) was the active constraint because the maximum velocities in tests F and G are relatively modest. By contrast, the computations of EISMINT II experiment F reported below were affected by both constraints.)

We use a computational ‘box’ with dimensions 1800 km × 1800 km × 4000 m. For example, a spatial grid with spacing Δx = Δy = 10 km and Δz = 40 m corresponds to a threedimensional grid with 181 × 181 × 201 ≈ 6.6 × 10⁶ points. Note that our simple scheme has equal spacing in the vertical (Appendix A).

Our first use of test F for verification is to refine the horizontal grid, namely to use Δx = Δy = 30, 15, 10, 7.5 and 5 km, and to consider thickness errors. Figure 10 shows that there is reasonable decay, under horizontal grid refinement only, of average and dome thickness errors. (Various Δz in the range 6.7–80 m were used to produce Figure 10. The role of Δz is addressed next.) This experiment gives an apparently disappointing result for maximum thickness error, however, as refinement produces maximum thickness errors which are all in the range 25–29 m. These maximum errors do not significantly reduce with refinement (not shown).

Fig. 10. Average and dome thickness errors for test F. Several values of Δz were used for each horizontal refinement level; see text.

The maximum thickness and temperature errors always occur in the vicinity of the inevitably poorly approximated margin (Reference Bueler, Lingle, Covey and BowmanBueler and others, 2005). As discovered by Reference Calvo, Dſaz, Durany, Schiavi and VasquezCalvo and others (2002b), however, the quantity η = H^(2n+2)/n = H ^8/3 is much better approximated than the thickness H itself because η has a bounded gradient at the margin. As shown in Figure 11, the maximum error in η converges at a reasonable rate under purely horizontal grid refinement. For clarity we report the relative maximum of η, that is, max |η∏um ^_ ηexact\/η exact.

Fig. 11. Maximum relative error of η = H ^8/3 in test F. Convergence under grid refinement is clear and occurs at a reasonable rate. Several values of Δz were used for each horizontal refinement level.

The grid refinement study reported in Figures 10 and 11 was carried out with several choices of Δz in the range 6.7–80 m. Thickness errors depend weakly on Δz. We can expect temperature errors to depend on the size of Δz, however. If Δx = Δy = 10 km are fixed and Δz is reduced steadily from 80 m to 6.7 m then the average temperature error in test F falls steadily from 0.12 to 0.02 K (not shown).

Clearly grid refinement is a 3-D matter. Instead of reassessing test F in this light, we turn to test G which is time-dependent and therefore of greater interest anyway.

For test G we chose a fixed 3-D refinement path Δx = Δy = 1800/N km and Δz = 4000/N m with N = 60, 90, 120, 180, 240 and 360. In particular, Δx ranges from 30 km down to 5 km while Δz ranges from 66.7 m down to 11.1m. The time-step Δt is limited to 10 years in all runs, but this constraint is only active in the Δx = 20, 30 km cases since for finer grids the time-step was further reduced by adaptive time-stepping (constraints (25) and (27); see Analysis).

We observed reasonable convergence of all the quantities measured. The convergence is predictable for the part of the refinement path with Δx ≤ 15 km. The observed rates O(Δx¹’⁹⁰) and O(Δx^1.98) in Figures 12 and 13, respectively, for the average thickness and relative maximum η errors, are essentially optimal for our O(Δx²) local truncation error scheme for the mass conservation equation (Reference Morton and MayersMorton and Mayers, 2005). In fact, given that there is a free boundary problem for thickness, with margin location not determined in advance, these rates are remarkably good.

Fig. 12. Average and dome thickness errors for test G. In this and the following two figures the horizontal and vertical grids are simultaneously refined. Rate of convergence given for Δx ≤ 15 km runs.

Fig. 13. Maximum relative error of n = H^8/3 in test G.

The horizontal location of the margin is at best approximated up to O(Δx) error (Reference Calvo, Dſaz, Durany, Schiavi and VasquezCalvo and others, 2002b; Reference Bueler, Lingle, Covey and BowmanBueler and others, 2005), but this rate is hard to observe because the location of the margin is discrete. There either is or is not ice at a given gridpoint.

Figure 14 shows the numerical errors in temperature using test G. Note that the thickness error made at a single point, namely the dome, and the temperature error at a single point, under the dome at the base, can be somewhat unpredictable. Reports of average and maximum errors are more meaningful.

Fig. 14. Temperature errors for test G.

It is believed that Figures 12–14 are the strongest possible indication that we have chosen and correctly coded an effective thermomechanically coupled scheme. In the continuum circumstances of tests F and G, grid refinement produces good agreement between the numerical result and the continuum solution, as well as good rates of convergence.

The angular symmetry of the basal temperature map is of obvious interest. It is shown for test G at the end of a 200 kyr model run with Δx = 30 km in Figure 15. There is no significant angular asymmetry even on this rough grid, and grid refinement clearly reduces the asymmetry. Figure 14 already addresses angular symmetry because the maximum temperature error is less than 2 K at every point in the 3-D domain of ice (for all levels of refinement). This shows that spokes of greater angular variation than 2 K cannot exist, as the reported error compares the numerical result to the exact and angularly symmetric solution to the continuum problem.

Fig. 15. Basal temperatures in test G at 200 kyr with Δx = Δy = 30 km and Δz = 66.7m. Grey shaded region is at or above pressure-melting temperature. Contour interval is 2 K, and the first contour outside the grey shading is 270 K.

The non-existence of spokes in tests F and G is not asserted to be physically interesting. Tests F and G are, as noted, not based on the physical conservation of energy but rather on the more general Equation (17). We emphasize, however, that because an angularly symmetric exact solution is known, any spokes would appear as numerical errors in a figure such as Figure 14.

Note that test G is as cold as EISMINT II experiment F. Figure 15 recalls another fact, however, namely that the temperature regimes of tests F and G are not physical. Tests F and G solve Equation (17) while EISMINT II experiment F describes a temperature field solving Equation (2).

Note that the temperature in tests F and G is allowed to rise above the pressure-melting point, although the absolute temperature remains below 273 K, and thus the shaded region in Figure 15 is merely for illustrating the contrast of basal temperature regime.

On the spokes in EISMINT II experiment F

We now have a ‘verified’ numerical scheme, verified to the extent allowed by known exact solutions. We can apply our numerical scheme to EISMINT II experiment F. When comparing 11 models, experiment F produced spokes for all schemes which included horizontal advection (Reference PaynePayne and others, 2000). For this numerical experiment, one solves the thermomechanically coupled system including Equation (2). We also return to the Reference Paterson and BuddPaterson and Budd (1982) flow law and use a 1500 km × 1500 km × 5000 m computational box.

Our verified scheme produces spokes for EISMINT II experiment F. Furthermore, they do not disappear under, nor are they reduced by, grid refinement. Our results confirm those of Reference Payne and BaldwinPayne and Baldwin (2000) and Reference Saito, Abe-Ouchi and BlatterSaito and others (2006). Figure 17 shows their configuration on a modestly refined Δx = Δy = 12.5 km grid (a factor of 2 finer than the EISMINT II choice). Note that our code produces results with angular symmetry for EISMINT II experiment A (not shown).

Fig. 16. Computed thickness and temperature field (K) in EISMINT experiment F. This is a slice along a flowline of a full 3-D computation with Δx = Δy = 12.5km and Δz = 25m. Compare with Figure 1.

Fig. 17. Basal temperatures in EISMINT experiment F with Δx = Δy = 12.5 km and Δz = 25m. Shaded region is at pressure melting temperature. Contour interval is 2K and innermost, coldest contour is at 244K.

The temperature distribution at depth for experiment F is significantly different from that of tests F and G; (cf. Figs 1 and 16). It may also be the case that the margin in experiment F is more singular than that in tests F and G (Reference Saito, Abe-Ouchi and BlatterSaito and others, 2006), but this is difficult to evaluate as we have only gridded approximate margin shapes for experiment F.

We will analyze our EISMINT II experiment F numerical results as solutions close to continuum solutions of importance, in the Discussion section.

Margin approximation

It is well known that at a stationary grounded ice-sheet margin in the SIA approximation, |∇ H| is unbounded. It follows that all known numerical methods are subject to large numerical thickness errors at the margin. In fact, these errors completely dominate the numerical thickness errors made anywhere else in an ice-sheet computation (Reference Bueler, Lingle, Covey and BowmanBueler and others, 2005). The only solution to this difficulty is to recognize that the (isothermal) SIA is well posed for transformed thickness η = H^(2n+2)/n (Reference Calvo, Dſaz, Durany, Schiavi and VasquezCalvo and others, 2002b). The quantity η has bounded gradient at the margin and is well approximated by many numerical schemes.

We are describing here the issue of agreement between numerical and continuum solutions to our model (the SIA). A distinct issue is whether the SIA is adequate for particular modelling goals, in other words whether margins must be handled by the full, non-shallow Stokes equations in order to adequately model real ice sheets. The results of Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson (2004) suggest that even in surprisingly non-shallow circumstances the isothermal SIA produces results close to those of the isothermal Stokes model.

The margins of interest to users of the thermomechanically coupled SIA have strongly varying temperatures within the ice near the margin. This affects their shape in a poorly understood manner. For the exact solutions described in the previous section, however, the margin is chosen to have the same profile as the isothermal, stationary, ablation-zone margin. For such isothermal margins, if × is the distance from the margin then the thickness satisfies H ~ x ^1/2 for small × (Reference FowlerFowler, 1997).

Based only on numerical computation we believe that realistic thermomechanically coupled, stationary, ablation-zone margins are more singular than isothermal margins (cf. Figs 1 and 16). It is, however, very difficult to diagnose the shape of a margin based purely on numerical evidence. Better numerical techniques may clarify the situation (Reference Saito, Abe-Ouchi and BlatterSaito and others, 2007). No complete asymptotic analysis of thermomechanically coupled margin shape, for the SIA, has been found in the literature, although Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler (2001) analyses the shape of the margin in the case that advection is ignored.

Diffusivity of the mass conservation equation

Because we are considering the SIA, Equation (1) and its isothermal version (13) can each be written in the form of a non-linear diffusion equation

with diffusivity coefficient

(Reference Hindmarsh, Payne, Hooke and HutterHindmarsh and Payne, 1996; Reference Saito and Abe-OuchiSaito and Abe-Ouchi, 2005). The thermomechanically coupled formula combines Equations (1), (5) and (6). It is correct for arbitrary flow laws F(T, σ). This scalar quantity D, the diffusivity coefficient, is known numerically at every mass conservation time-step because it is the scalar factor involved in computing the flux Q. The identification Q = –D∇H is the ice flow analogue of Fourier’s law for heat, although the expression where is the vertically averaged horizontal velocity, also holds. As diffusion equations, (1) and (13) are simultaneously highly non-linear and highly degenerate. (‘Degenerate’ applies because D → 0 at margins and divides, where H → 0 and |∇H| → 0, respectively.)

The maximum value of D anywhere on an ice sheet controls the largest stable time-step for explicit approximations of Equation (1) (Reference Hindmarsh, Payne, Hooke and HutterHindmarsh and Payne, 1996; Reference Hindmarsh, Straughan, Greve, Ehrentraut and WangHindmarsh, 2001). Indeed, in the thermomechanically coupled case using the explicit numerical scheme described in Appendix A, we observe the stability restriction

The maximum is taken over values of D computed at every horizontal grid location (x_i, y_j ) on the simulated ice sheet at the current time-step t_l.

Restriction (25) is reliable enough to use in an adaptive time-stepping scheme if advection in the temperature equation is also well approximated; see the comments on numerical approximation of the advection below. Our observed restriction (25) is only slightly more restrictive than that in the isothermal case, where Reference Hindmarsh, Straughan, Greve, Ehrentraut and WangHindmarsh (2001) predicts a critical value of the righthand side of (25) of 1/(2n) = 1/6. When the time-step violates inequality (25) the resulting stability problems (wiggles) emerge at locations on the ice sheet where the function H ⁿ⁺²|∇H| ^n–1 is largest. This will generically occur away from divides and margins.

Implicit schemes for Equation (1) avoid stability restrictions such as (25), but fully coupled implicit schemes are significantly harder to implement. In any case, there is an additional restriction on thermomechanically coupled timestepping, described next.

Domination of the temperature problem by vertical advection

In simulating real ice sheets, rough surfaces sometimes produce large vertical velocities in columns of ice below the roughness. Similarly, sliding laws which allow an abrupt onset of sliding adjacent to frozen bed areas can produce unreasonably large vertical velocities (through incompressibility). These large velocities can produce a poor approximation of the advection part of the temperature problem if the numerical time-step is not restricted.

To explain why large vertical velocities appear, supposewe combine Equations (5) and (7), and suppose for completeness there is a possibly non-zero basal horizontal velocity field U _b and a possibly non-flat bed z = b. (We now need to distinguish between the surface elevation z = h and the thickness H.) Then

where F = F(T, σ) and ∇² h = ∇·(∇h).One sees that the vertical velocity is affected by both the first and second derivatives of the surface h, by the divergence of U _b and by the gradient of b. The integrals here are not the point;we are simply identifying dependencies on horizontal derivatives and thus on possible numerically induced inaccuracies.

Because the vertical velocity depends on the second derivative of the surface, it is more sensitive to surface roughness than are the horizontal velocity components (which depend only on the first derivatives of the surface elevation). In particular, a gridded (rectangular or triangular mesh) approximation of the surface elevation which is acceptable for computing horizontal velocity may be rough enough to cause unreasonable vertical velocities in some columns of ice. We observe, for instance, that digital elevation maps for Antarctica drive this effect in many locations when these maps are used as initial conditions for ice-sheet simulation.

The dependence of w upon the divergence of the basal horizontal velocity field has been observed before. For instance, EISMINT II experiment H, which involved pressuremelting temperature-activated basal sliding, has been criticized because it used basal velocities U _b proportional to the effective shear stress and involved sudden transitions of the basal velocity as one moves from frozen base to melted base (Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler, 2001). The resulting jump in U _b produces unbounded vertical velocity within the ice in the continuum limit. Such discontinuities in basal velocity are presumably avoidable by using more complete physical models of the basal conditions.

Flow velocities are of importance here because they advect heat, and this advection must be numerically approximated. Numerical schemes for Equation (2) may have conditions on the length of the time-step in terms of the magnitude of the advection velocities. Upwinded schemes of any finite order require some version of the CFL condition (Reference Morton and MayersMorton and Mayers, 2005) for stability.

For example, we use a first-order upwind for temperature Equation (2); see Equation (A2) in Appendix A. This method sacrifices accuracy and it produces some artificial diffusivity but it is essentially as stable as possible (Reference Morton and MayersMorton and Mayers, 2005). For our scheme a precise form of the CFL condition for pure advection is the inequality

See Appendix B and Equation (B3) in particular. Here the numerators of the terms within parentheses are the gridded values of the components of velocity. The maximum is over the whole 3-D gridded region of ice at the time t_l . Condition (27) is combined with condition (25) in our adaptive timestepping scheme, and the shorter of the two time-steps is used.

We observe in practice that condition (27) is more likely to be activated by large vertical velocities than horizontal velocities. Of course, vertical speeds in an ice sheet are typically lower in absolute terms than horizontal speeds, but the vertical mesh is always much finer. As described above, numerically computed vertical velocities can be of surprisingly large magnitude because of surface roughness and/or large gradients in the basal velocity field.

Payne and Dongelmans (1997) introduced a scheme with unequally spaced five-point upwinding which maintains stability with higher-order local truncation error.We believe this scheme, which is subject to some unstated CFL condition, is used in several current ice-sheet models. Such a scheme is definitely worthwhile if it can be shown that the global approximation error – the actual error – is reduced. The exact solutions of this paper are well suited to demonstrating such improved performance.

The strain-heating term: emergence of spokes

A standard finite-difference convergence analysis (Appendix B) of our scheme for Equation (2) reveals that if condition (27) holds for each time-step then the growth of the numerical error is controlled by two quantities. One quantity is the smoothness of the continuum solution to Equation (2), as reflected in the local truncation error of the scheme (Reference Morton and MayersMorton and Mayers, 2005). The other quantity, however, is the size of a derivative of the strain-heating term ∑. Recall that

Let E_l be the maximum of the numerical temperature errors over the entire spatial grid at time-step t_l . That is, consider the maximum of the difference between the computed approximate temperature T_ijkl and the exact temperature T(x, y, z, t) predicted by the continuum equations but evaluated at the gridpoint (x_i , y_j , z_k ) at time tl; this is an instance of the definition of numerical error above. One can show (Appendix B) that

where τ_l is the local truncation error at time t_l . Note that τ_l depends on the size of higher partial derivatives of the exact solution to the temperature equation, that is, it depends on the smoothness of the continuum solution.

By inequality (29) the local truncation error can, at worst, cause arithmetic growth in the actual error E_l . By contrast, in estimate (29) the error E_l multiplies ∂Σ/∂T. The worst-case error can therefore grow exponentially (geometrically) if ∂Σ/∂T is large. The quantity ∂Σ/∂T acts like an interest rate for growth of error. A refined form of inequality (29) is inequality (B7), which shows that the exponential growth can be greatest at the locations where ∂Σ/∂T is maximum.

Consider the spokes which appear in EISMINT II experiment F, and are illustrated in Results above. Note that the EISMINT II experiment F continuum set-up has perfect angular symmetry. The exact solution of this initial/boundary value problem, which is unknown, also has this angular symmetry. Thus the magnitude of the spokes, as measured by the deviation from a given temperature value as one follows a curve of fixed z and r in a circle within the ice around the z axis, is a lower bound for the numerical temperature error. (In this literal sense, the spokes unequivocally are a numerical error.)

Figure 18 shows that ∂Σ/∂T has strong maxima along the warm spokes in EISMINT II experiment F. (In this figure, large values of ∂Σ/∂T for points where the thickness is less than 1000m are set to zero. We assert that these values are corrupted by margin mis-approximation, but that the slow flow very close to the margin (e.g. within 20 km of the margin) is unable to feed back to affect the surface slopes and near-basal strain rates at the upstream locations of the onset of spokes. Better numerical schemes for approximation of the fields near the margin may remove this issue (Reference Saito, Abe-Ouchi and BlatterSaito and others, 2007).)

Fig. 18. Basal values of ∂∑/∂T (10^–12s^–1) from EISMINT experiment F, from the same computation which provided Figure 17. To give clear detail, only one-quarter of the computational region is shown.

Thus the spokes emerge under a version of the well-known thermocoupling feedback mechanism, with warmer spokes locally drawing down the surface and steepening the gradient nearby, thus increasing the strain heating upstream of the warm spokes, and generally accelerating flow to produce more strain heating, and so on. This describes ‘thermal runaway’, which is also observed in fluid systems other than ice (Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler, 2001). The spokes themselves are errors, however, because the continuum problem and its unknown (exact) solution have angular symmetry. Thus the quantity which is ‘running away’ in this analysis is the size of the numbers the numerical temperature errors.

The analysis here usefully predicts that spokes emerge in locations of peak values of the diagnostically computable quantity ∂Σ/∂T. On the other hand, the derivative of Σ with respect to T, using the second form of Σ in Equation (28), is:

This quantity is always positive. From Equation (30) it follows that if a realistic Arrhenius-type term in the flow law is used with the Reference Paterson and BuddPaterson and Budd (1982) law as the representative case, then as the temperature rises from cold to within 10 K of the melting point the quantity ∂Σ/∂T will jump up discontinuously. The jump occurs at a temperature of 263.15 K in the Paterson–Budd relation, where Q changes from 6.0 to 13.9 (in units of 10⁴ J mol^–1): see Figure 19. Note Σ is a continuous function of temperature but ∂Σ/∂T is not in the Paterson-Budd case. The Hooke (1981) flow law, which has been demonstrated to produce spokes for experiment F in practice (Reference Payne and BaldwinPayne and Baldwin, 2000), apparently behaves in an even worse manner because the value of ∂Σ/∂T near melting is even larger, although there is no actual jump. This seems to cause even stronger error growth in the immediate vicinity of the margin.

Fig. 19. Graphs A(T) for the Reference Paterson and BuddPaterson and Budd (1982) and Hooke (1981) flow laws, and of the cold and warm simple Arrhenius parts of the Paterson–Budd law extended to all temperatures. For the Paterson–Budd and Hooke laws the slope of A(T) jumps or changes rapidly as one approaches the melting temperature. Therefore ∂Σ/∂T jumps or changes rapidly; see Equation (28).

However, if the contribution of Σ to the temperature field as in Equation (2) is non-physically ‘smeared’ then the spokes go away. In particular, suppose we smooth the strain-heating term Σ in horizontal directions, at each temperature time-step, by replacing gridded values Σ_ijk,l with an average of horizontally neighboring values. An example of the result for the basal temperature map is shown in Figure 20; the spokes are gone. In this case, the smoothing was by a finite-support averaging kernel which was approximately Gaussian and had a standard deviation of about two grid spaces (25 km) and a maximum range of four grid spaces (50 km) on this 12.5 km grid. From runs on various grids we believe the smoothing is effective in eliminating the spokes if it has a range (standard deviation) of the order of 25 km for a 1 year time-step.

Fig. 20. Basal temperatures from EISMINT II experiment F, exactly as in Figure 17, but recomputed with smoothing of the strain-heating term. Compare also with Figure 15.

The average was globally energy-conserving for the ice, that is, all of the heat generated by strain heating (and no more heat) was put into the ice. The location of that heating was spread in the horizontal by a mechanism not explained by the missing horizontal conduction terms in the conservation-of-energy equation, however. Easy estimates of the thermal diffusivity of ice show that our ad hoc smoothing is significantly greater than would be the effect of reintroducing horizontal conduction of heat through the ice or bed. (Neither of these reintroductions was numerically tested, however.) An interesting possibility is that the inclusion of longitudinal stresses into the coupled system would spread the strain heating itself in the horizontal (Reference HindmarshHindmarsh, 2006) and provide a physical mechanism which smears the spokes.

The ‘thermal runaway’ of the numerical temperature errors is rather slow in that it takes at least thousands of model years to become significant. It can only occur if the conditions causing the runaway remain in place for a significant time and apply over significant spatial extent. In particular, EISMINT II experiment F results show spokes in its flat base and smoothly varying accumulation conditions but they take more than 10 kyr to stabilize after the geometry is roughly in place. It is therefore likely that in conditions of significantly varying bed elevation, or of varying surface temperature and accumulation boundary conditions, the necessary coherence for thermomechanically coupled runaway of error does not occur. Thus the spokes phenomenon might be minimally influential for the simulation of real ice sheets.

The ‘thermal runaway’ of the error, though it may grow exponentially in some cases, has a maximal exponential rate if the geometry is fixed. In fact, one can turn inequality (29) into a proof of convergence of our numerical scheme for temperature under the assumption of fixed geometry and also stringent smoothness assumptions about the velocity fields (Reference Bueler and BrownBueler and Brown, 2006). This proof yields a maximum exponential rate based on the largest value of ∂Σ/∂T as described above.

We see no strong reason to believe that some version of the spokes will not appear when numerical 3-D thermomechanically coupled Stokes simulations are applied to circumstances similiar to EISMINT II experiment F. This is because the temperature equation keeps its same strain-heating term, and thus an analysis of the numerical error in the temperature equation will reveal the same possibility (at least) of ‘thermal runaway’ of the error. Indeed it appears that in at least one higher-order stress balance scheme there are still spokes for the EISMINT II experiment F boundary values (Reference Saito, Abe-Ouchi and BlatterSaito and others, 2006).