Introduction
The aim of this paper is to address several issues arising from a recent intercomparison exercise supported by the European Ice Sheet Modelling Initiative (EISMINT) and reported by Reference PaynePayne and others (in press). The exercise concerned the effects of introducing thermomechanical coupling into icesheet models. Timedependent models of icesheet evolution are typically based on two prognostic equations: one for the evolution of ice thickness and one for icetemperature evolution. The former comprises a statement of icemass continuity and a diagnostic relationship between icesheet form and velocity. The two prognostic equations are linked in a number of ways so the evolution of the coupled system is likely to be complex and nonlinear.
The original experiments used an idealized geometry in which all boundary conditions had radial symmetry. One of the key findings was that, in some cases, this symmetry was broken by the development of spokes of cold ice extending outwards from the icesheet divide towards its margin. Interactions between ice form (and flow) and temperature were identified as the cause of this pattern development. There was, however, some doubt about the extent to which numerical instability affected these results. This arose principally because some elements of the patterning followed the geometry of the underlying numerical grid. The original experiments were performed by ten different groups and the primary aim was model intercomparison. They did not therefore provide a convenient means of addressing these issues, which are therefore addressed here.
The patterning is believed to arise through the interaction of the evolving ice thickness, temperature and velocity fields. It is essentially a localized form of the "creep instability" discussed by Reference Clarke, Nitsan and PatersonClarke and others (1977). A slight positive anomaly in ice temperature softens ice locally. This leads to increased horizontal ice velocities because softer ice will flow faster for a given gravitational driving stress. This fasterflowing ice then generates greater quantities of heat via dissipation, which enhances the original temperature anomaly. Eventually, ice thickness (and icesurface elevation) is locally reduced because of the increased ice flow. This then favours further warming because the horizontal flow of ice is always downhill (over the typical length scales associated with ice sheets) and so increasing quantities of ice become channelled through the site of the original temperature anomaly. The result is the differentiation of an initially uniform iceflow field into a pattern where zones of warm, fast flow alternate with ones of cold, slow flow.
The influence of model numerics on this instability arises through the nature of the initial icetemperature anomaly. Finitedifference models, such as the one used in this paper, employ a regular, rectangular mesh on which discrete versions the model’s equations are solved. When the geometry of the modelled ice sheet does not match this rectangular mesh (as is the case with the circular ice sheet modelled here), slight differences in modelled quantities between adjacent gridpoints can occur. These differences are then sufficient to trigger the instability described above. Model numerics therefore play a role in the initiation of the instability A related issue is whether the actual discretization employed by the numerical model is influencing the patterning.
The specific questions addressed in this paper follow directly from the discussion of Reference PaynePayne and others (in press) and are:

(1) how the patternforming instability is initiated;

(2) whether the process of pattern development is reversible;

(3) whether the particular relationship between ice temperature and viscosity used in the intercomparison played a role in the pattern formation; and

(4) whether pattern formation has any dependency on horizontal grid size or timestep.
Reference PaynePayne and others (in press) describe three key experiments, differing only in their parameterization of air temperature which provides the upper boundary condition for the threedimensional, advectiondiffusion equation governing ice temperature. There was a 20 K difference in air temperatures between these experiments and this difference effectively controlled whether the patterning developed or not. In the experiment with colder air temperatures (223 K at the divide), all models exhibited the patterning. While in the warmer experiments (243 K at the divide), the models showed no sign of pattern formation.
These observations prompted the first two questions to be investigated. The model behaviour which occurs in the phase space between the cold, patterned and warm, symmetric end members was not investigated. A series of five experiments which address this range in more detail are therefore used to investigate states which are intermediate between no and fully developed patterning. Further, the initial experiments gave some indication that the patterning could be reversed if air temperatures were warmed sufficiently Here we investigate whether an airtemperature cooling (from a divide air temperature of 243 K to one of 223 K) is sufficient to induce pattern formation and whether this patterning disappears once the cooling is reversed (divide temperature rising from 223 K back to 243 K).
Reference PaynePayne and others (in press) identify one potential cause of the pattern development as the form of the relationship used to calculate ice viscosity from temperature. This relationship was proposed by Reference Paterson and BuddPaterson and Budd (1982) and is exponential. The form of this relationship is, however, different for ice above and below 263 K. There is in fact a discontinuity at 263 K and the softening of ice with increasing temperature is far more dramatic above this point. It is possible that the observed patterning is triggered as ice crosses this discontinuity An alternative form of the relationship between ice temperature and viscosity, continuous over its entire range, was proposed by Reference HookeHooke (1981). The idea that patterning is caused by the specific form of the Reference Paterson and BuddPaterson and Budd (1982) relationship can therefore be tested.
Finally, the hypothesis that the patterning owes its origin to numerical instability can be addressed by repeating the original experiments employing a range of timesteps and horizontal grid spacing.
Model Description
This section summarizes the thermomechanical icesheet model. Reference Huybrechts and PayneHuybrechts and others (1996) and Reference PaynePayne and others (in press) provide a more detailed discussion. The model is time, t, dependent and is quasithree dimensional (with horizontal dimensions × and y and a vertical dimension z which is positive upwards from sea level). Isostatic deflection of the underlying bedrock is not incorporated into the model.
The evolution of ice thickness, H, is given by the continuity equation
where ▽. is the twodimensional horizontal divergence operator, b is the local snowaccumulation/ablation rate, and is the vertically averaged horizontalvelocity vector. Ice deformation is assumed to be driven solely by horizontal shear stresses, τ_{xz} and τ_{yz}, which, at the spatial scales over which ice sheets operate, can be approximated by
where s is icesurface elevation, p is ice density and g is acceleration due to gravity. The values used in this study for all constants are given in Table 1. Assuming a nonlinear flow law for ice (Reference GlenGlen, 1955), one can find horizontal velocity as
where h is bedrock elevation, n is usually taken as 3 and the parameter A introduces the temperature dependence of ice viscosity Similar expressions are used to find τ_{yz}(z) and U_{y}(z), and vertical velocity, w(z), is found from the horizontalvelocity field using the incompressibility condition. Basal slip is not included in the experiments described here.
The ice accumulation/ablation rate, b, is specified as a radially symmetric function of horizontal position, × and y,
where b_{max} is the maximum accumulation rate, S_{b} the gradient of accumulationrate change with horizontal distance and E is the distance from the centre of the model domain at which accumulation rate is zero. This parameterization results in a large accumulation area, over most of which accumulation is held constant at 6_{max}. Around this area, b falls linearly with distance from and rapidly becomes negative (ablation).
The temperature dependence of ice viscosity is incorporated initially using an Arrhenius relation, choosing constants based on work by Reference Paterson and BuddPaterson and Budd (1982)
where T* (in K) is absolute temperature corrected for the dependence of melting point on pressure, a is a constant of proportionality, Q is the activation energy for ice creep, and R is the universal gas constant. It should be noted that the values of a and Q are different above and below 263 K.
This relationship is replaced in some experiments by that due to Reference HookeHooke (1981)
In Equations (5) and (6), T* is defined as
where T_{pmp} is the pressuremeltingpoint temperature and T_{0} is the triple point of water.
The ice temperature, T, evolves according to
where U is the horizontalvelocity vector, k is the conductivity of ice and c is its specificheat capacity. These and the other thermal parameters used in the model are assumed to be independent of ice temperature. Equation (8) forms the second prognostic equation on which the model is based. It contains terms representing vertical diffusion, horizontal advection, vertical advection, and dissipational heat generation, respectively.
Two boundary conditions are required for the integration of Equation (8) forward through time. The upper boundary condition is a specified surface air temperature (T_{a})
where T_{min} is the minimum surface air temperature and ST the gradient of airtemperature change with horizontal distance. Air temperatures therefore increase in a linear fashion away from the centre of the model domain and are radially symmetric. Several of the experiments reported below employ changes to the value of T_{min}, which is a control on the overall temperature of the ice sheet and on the presence of patterning.
The boundary condition at the icesheet base is
where G is the geothermal heat flux (assumed constant).
The evolving temperature field is constrained so that it cannot exceed the pressuremeltingpoint temperature
An equivalent melt rate is determined where calculated temperatures exceed the pressuremelting point.
The model summarized above is solved using finite differences on a regular horizontal and a stretched vertical grid using mostly implicit methods. A discussion of the numerical methods employed in the model can be found in Reference Hindmarsh and PayneHindmarsh and Payne (1996) and Reference Payne and DongelmansPayne and Dongelmans (1997).
Experiments
The results from a series of experiments designed to address the issues raised in the Introduction will now be discussed. All experiments lasted for 200 kyr by which time a longterm equilibrium had developed. A horizontal grid spacing of 25 km and a timestep of 50 years were employed except in the final set of experiments. Unless otherwise stated, the initial conditions are zero ice on a flat bedrock plane.
The nature of the patterning is such that regions of higher temperatures are always associated with softer ice (increased values of A), faster horizontal velocities and reduced ice thicknesses (Reference PaynePayne and others, in press). In order to simplify the discussion of the results, we will therefore concentrate on the predicted patterns of basalice temperature.
The onset of patterning
The onset of patterning was investigated in series of five experiments in which T_{min} in Equation (9) was given the values 238, 235.5, 233, 228 and 223 K (they will be referred to as experiments AE). These and all other experiments reported here are summarized in Table 2. The simulated basalice temperature field in A is radially symmetric and is dominated by patterning in E (Fig. 1). The extensions of cold, interior ice toward the margin will be referred to as "cold spokes", while "warm spokes" are the intervening divideward extensions of warm, margin ice.
Much of the information in Figure 1 is redundant. The comparison of a variety of experiments is therefore eased by using simplified representations of this information. This is done in Figure 2 for experiment E. The basal temperature at each model gridpoint is plotted as a function of its distance from the domain centre or ice divide (the black dots in Fig. 2). These temperatures show very little scatter near the divide and towards the margin. This is indicative of radial symmetry. In contrast, points falling 150 to 500 km from the divide show much scatter, which is indicative of pattern formation (see Fig. lb).
This information can be further summarized by dividing the abscissa into 25 km sections and recording the maximum and minimum within each section (as indicated by the bars in Fig. 2). The difference between the maximum and minimum for each section provides a useful index of the degree of patterning.
The degree of patterning in all five experiments is summarized in Figure 3. The patterning develops between experiments with T_{min} = 233 K and T_{min} = 235.5 K. The development of patterning is therefore triggered fairly abruptly over a small temperature range. The contrast between warm and cold spokes increases as T_{min} decreases. The warm spokes are invariably at pressuremelting point. This difference can therefore be interpreted entirely in terms of cooling within the cold spokes and the majority of this cooling is simply a reflection of the cooler upper boundary condition (T_{a}).
The reversibility of patterning
Two further experiments (each lasting 200 kyr) were undertaken to test whether the patterning shown in Figure 1 is reversible. Experiment F used the final equilibrium state reached by experiment A (which shows no patterning) as an initial condition and imposed a stepped change of T_{min} from 238 K to 223 K. The results, again in terms of radial temperature anomalies, are shown in Figure 4a. The final equilibrium state of experiment F is directly comparable to that of experiment E (which also employed T_{min} = 223 K). The relevant lines in Figures 3 and 4a indicate a close but not perfect correspondence. Investigation of the actual basal temperature field indicated that rather than the 12 coldice spokes shown in Figure 1b, a total of 20 had formed. The diagonal spokes shown in Figure 1b each split into one central, major spoke and two smaller neighbours.
Reversibility was tested in experiment G by using the final state produced in experiment Fas an initial condition and reversing the change of T_{min} back from 223 K to 238 K. The results can be seen in Figure 4b. The patterning rapidly breaks down and radial symmetry is restored after approximately 50 kyr. The final state reached in experiment G is essentially identical to that of experiment A (basal temperature at the divides are within 0.0004 K of each other).
The iceviscosity relationship used
Experiments H and I are identical to experiments A and E except that the Reference Paterson and BuddPaterson and Budd (1982) (Equation (5)) iceviscosity relationship is substituted with that of Reference HookeHooke (1981) (Equation (6)). A comparison of the values of A predicted by the two relationships over a temperature range from 223 K to 273 K is shown in Figure 5. The Reference Paterson and BuddPaterson and Budd (1982) relationship predicts higher values of A at temperatures below 250 K and in a narrow range between approximately 268 K and 272 K.
The basaltemperature distribution at the end of the two experiments is shown in Figure 6, which is directly comparable to Figure 1. The patterning is much stronger in the experiments employing the Reference HookeHooke (1981) relationship. A Tmm of 238.0 K was previously too warm for patterning to develop when employing Reference Paterson and BuddPaterson and Budd (1982), however the patterning is very clearly developed with the Reference HookeHooke (1981) relationship. Patterning in the experiment with T_{min} = 223.0 K is now so pronounced that the concentration of flow leads to the extension of the icesheet margin downstream of warmice spokes. This intensification of the patterning is probably because of the step (but continuous) rise in the value of A close to the melting point, which is a feature of the Reference HookeHooke (1981) relationship (Fig. 6).
Changing timestep and grid spacing
The timestep (dt) and horizontal grid spacing (da;) used in all of the experiments reported so far are 50 years and 25 km, respectively. It should be noted that the maximum horizontal velocity in these experiments is approximately 80 ma^{−1}, which with dx = 25 km would equate to a maximum timestep (in terms of Courant stability) of 312.5 years. A useful test of the role that numerics plays in the development of patterning is to change the numerical timestep and grid spacing.
A series of three further experiments repeated experiment E with dt = 10, 25 and 100 years (experiments J, K, and L). The results, again presented as radial temperature anomalies, are summarized in Figure 7a. Although the degree of patterning is consistent in all four experiments, there are differences in the geometry of the spoked pattern produced in each experiment. In particular, the pattern produced in experiment J (dt = 10 years) features a total of 20 spokes and was very similar to that described for experiment F above. The results for experiments E and K (dt = 50 and 25 years) were virtually identical.
Experiments with dx = 50 km (M) and dx = 12.5 km (N), in addition to the standard dx = 25 km experiment E are shown in Figure 7b. The radial temperature anomalies reveal very different patterning associated with each horizontalgrid resolution. Inspection of the basaltemperature field produced in each experiment indicates that the number of individual coldice spokes rose to 32 in the 12.5 km experiment, while this number fell to 8 in the 50 km experiment (from the 12 of experiment E).
Discussion
The four specific questions raised in the Introduction can now be answered (at least partially).
Patterning develops fairly abruptly between two experiments separated by a 5 K difference in air temperature. Within this 5 K range there is an intermediate form (experiment B) which shows minor undulations in its basaltemperature field. This abrupt onset is thought to be related to the extension of relatively cold ice towards steeper icesurface slopes near the margin. The mechanism discussed in the Introduction partially relies on the deflection of horizontal ice flow towards warmer, fasterflowing areas of the ice sheet. This type of flow adjustment can occur more easily in areas of steeper surface slope near the margin, rather than in the flatter areas near the ice divide. One of the conditions necessary for pattern formation may therefore be that relatively cold ice (with the potential for warming and hence flow differentiation) occurs in close proximity to relatively steep icesurface slopes near the margin.
The development of patterning has been shown to be reversible. The details of the final form of patterning do, however, show some dependence on the pathway taken. The number of spokes developed in experiments starting from zeroice and radialicesheet initial conditions differs. This is an indication that the form of the patterning instability is very sensitive to precise numerical details (as discussed further below).
The discontinuity present in the Reference Paterson and BuddPaterson and Budd (1982) relationship between ice temperature and the value of A is not a necessary condition for pattern development. The relationship derived by Reference HookeHooke (1981) contains no such discontinuity but produces similar behaviour. The exponential nature of these relationships is probably a necessary condition for pattern formation (this could be tested by constructing artificial, linear relationships between temperature and A).
Finally, the details of the patterns developed show dependence on both the numerical timestep and horizontal grid spacing employed. Reference Hindmarsh and PayneHindmarsh and Payne (1996) show that gridspacing effects predicted icesheet geometry because of varying truncation errors. In uncoupled experiments, this affect can amount to as much as 100 m near the margin. These changes in icesheet geometry, exaggerated by thermomechanical coupling, may therefore be responsible for the slight differences in the form of the patterning observed here.
It should be stressed that the type of behaviour discussed in this paper is generic to the majority of the largescale models currently used to simulate the Antarctic and Greenland ice sheets as well as the former, midlatitude ice sheets. This is clear from the intercomparison reported by Reference PaynePayne and others (in press) and implies that the basal temperatures predicted by these models should be treated with some caution.
There is a strong coincidence between spoke orientation and the underlying finitedifference grid (or diagonals to it) which implies an underlying, numerical control. It is unlikely that this arises simply because we model a circular ice sheet using a rectangular grid. This is because Reference Payne and DongelmansPayne and Dongelmans (1997) found similar behaviour when modelling both square ice sheets and ones whose divide ran parallel to the rectangular finitedifference grid. In addition, icesheet models should be able to cope with arbitrary geometries because their main application is to irregular, realworld ice sheets. Future work must address why this grid dependence arises and must seek ways of removing it. One likely area is in the discretization of the dissipation term in the temperatureevolution equation (Equation (8)).
The model discussed in this paper does not incorporate longitudinal stresses. This omission is unlikely to affect the observed patterning because the flow is relatively slow (always ≤ 100 ma^{−1}). The warmice spokes are therefore not equivalent to ice streams. Longitudinal stresses may be important at the head of the warmice spokes where there are relatively steep horizontalvelocity gradients. However, the incorporation of the dissipation^ heating produced by these longitudinal stresses is likely to exaggerate creep instability and the patterning rather than lessen the effect.
Acknowledgements
This work is supported by a grant from the Natural Environment Research Council (GR3/11532). The authors are grateful to A. AbeOuchi, D. DahlJensen and F. Pattyn for the many helpful comments they made on the original draft of this paper.