Basic equation set

A thermodynamic ice-sheet model typically consists of two key prognostic equations: one for ice-thickness evolution (which includes ice flow) and one for the evolution of ice temperature. In addition, there are several ancillary diagnostic equations for variables such as vertical velocity, and snow accumulation and ablation. All of the models included in the intercomparison were time-(t-)dependent and threedimensional (with horizontal dimensions x and y, and vertical dimension z which is positive upwards from sea level).

Ice-thickness evolution is given by

where ∇ · is the two-dimensional horizontal divergence operator, M is the local accumulation/ablation rate on the upper surface (basal melt rates are several orders of magnitude smaller and are ignored) and is the vertically averaged horizontal velocity vector.

Ice flows by both internal deformation and basal slip (referring to any differential movement between the ice base and the fixed underlying bedrock). Ice flow by deformation is assumed to be driven solely by horizontal shear stresses which, at the spatial scales over which ice sheets operate, are a purely local function of ice thickness and surface slope (density ρ and acceleration due to gravity g being assumed constant). This assumption, in combination with the non-linear flow law for ice (Reference GlenGlen, 1955), allows horizontal velocity (u_i (z)) to be calculated from

where u_i (h) is the slip velocity, s is ice-surface elevation, h is bedrock elevation, n is usually taken as 3, T* is absolute temperature corrected for the dependence of melting point on pressure (see Equation (6)), and the flow factor A introduces the temperature dependence of ice deformation. Basal slip is not generally included in the experiments described here, and its parameterization is left to the description of the relevant experiments below.

When Equation (2) is averaged vertically from bedrock to ice surface and substituted into Equation (1), a non-linear parabolic equation results. This equation is usually referred to as the “ice-sheet equation” (Reference Hindmarsh and PayneHindmarsh and Payne, 1996), and its solution is at the core of most ice-sheet models (although a small number of models do solve Equations (1) and (2) separately). Horizontal ice velocities are then found diagnostically from the ice-sheet equation, and vertical velocity (w) is, in turn, found diagnostically using the divergence of this horizontal velocity field.

Before outlining ice-temperature evolution, the key link between ice-thickness/flow evolution and ice temperature is introduced. The flow factor A in Equation (2) is found from ice temperature using an Arrhenius relation

where a is a constant of proportionality, Q is the activation energy for ice creep and R is the universal gas constant. Ice temperature therefore controls the flow of ice locally via its non-linear relationship to the ice-flow factor. No additional flow-enhancement factor is used in the experiments.

The ice temperature (T) evolves according to

where k is ice conductivity and c _p is its specific heat capacity. These and the other thermal parameters used in the model are assumed to be independent of ice temperature. The terms in Equation (4) represent vertical diffusion, horizontal and vertical advection, and dissipation, respectively. The presence of the dissipation term in Equation (4) represents a further link between ice flow and ice temperature.

Two boundary conditions are required for the integration of Equation (4) forward through time. The boundary condition at the ice-sheet base is

where G is the geothermal heat flux (assumed constant). The evolving temperature is constrained so that it cannot exceed the melting point (T′):

where β is the Clausius–Clapeyron gradient. A melt rate may be calculated diagnostically once basal temperatures have reached melting point. The upper boundary condition is a specified ice-sheet surface temperature (T _surface), the details of which are discussed below.

Table 2 provides a list of the constants introduced above with the values assigned to them in the SGEs.

Table 2.

Constants used in the experiments

Experimental design

The experimental design is similar to the “moving margin” benchmark of Reference Huybrechts and PayneHuybrechts and others (1996). As in the original benchmark, the ice-accumulation/ablation rate (M in m a⁻¹) is a function of geographical position (x and y in km) alone:

where M _max is the maximum accumulation rate and S _b is the gradient of accumulation-rate change with horizontal distance. The accumulation rate becomes zero at a prescribed radial distance, R _el, from the summit (x _summit, y _summit). This parameterization results in a large accumulation area, over most of which accumulation is held constant at M _max. Around this area, M falls linearly with distance from (x _summit, y _summit) and rapidly becomes negative (ablation).

The ice-surface air temperature (T _surface in K) is also made a function of geographical position (in contrast to the benchmark where it was solely a function of ice thickness):

where T _min is the minimum surface air temperature and S _b is the gradient of air-temperature change with horizontal distance. Both Equations (7) and (8) are dependent solely on position and not on ice-surface elevation (as was the case in Reference Huybrechts and PayneHuybrechts and others (1996)). Although this situation is unrealistic, it does simplify the interpretation of the results which are now caused solely by internal thermomechanical interactions.

It is important to stress that Equations (7) and (8) are radially symmetric, and that all other external forcings are constant or radially symmetric over the model domain. The glaciological variables predicted by the numerical models should also therefore have radial symmetry.

A further change from the benchmark was the use of a 25 × 25 km² horizontal grid (formerly 50 × 50 km²). The model domain remains 1500 × 1500 km² (61 × 61 grid-cells). The ice-sheet divide is therefore located at the point (750 km, 750 km).

The choice of time-step was left to the individual modeller. All experiments lasted for 200 kyr. Most of the experiments employed a flat bedrock, and the effects of isostasy were ignored throughout. The similarity of these experiments to those of the “moving margin” benchmark allowed participants to check for coding errors before embarking on the SGEs.

A brief outline of the 12 SGEs is given in Table 3. Each SGE uses the specification for SGE A (below) and changes a single input parameter to investigate what effect this will have. This paper concentrates on seven of these experiments, which produced particularly interesting results. The remaining experiments produced results which were either difficult to interpret (experiments I–L) or provided trivial checks on model consistency (experiment E). The seven specific experiments are summarized below in greater detail.

Table 3.

Brief summary of the full SGE set

Experiment A

This experiment ran the coupled thermodynamic equations to equilibrium from an initial condition of zero ice on a flat bedrock topography. The following constants were used in the boundary-condition Equations (4) and (7)

Ice temperature and rheology (Equation (3)) were coupled using (based on Reference Paterson and BuddPaterson and Budd (1982))

Experiment B

This experiment used the final, steady-state ice sheet of experiment A (t = 200 kyr) as an initial condition, and applied the following altered boundary-condition constants:

All other constants had the same value as in experiment A.

Experiment C

This experiment also used the final, steady-state ice sheet of experiment A as an initial condition, and applied the following altered boundary-condition constants:

This experiment assessed ice-sheet response to a stepped accumulation-rate change. The change reduces the “plateau” accumulation rate from 0.5 to 0.25 m a⁻¹ and also the area over which this maximum value operates. It will generate two types of change: change caused by the reduced ice-sheet span and change caused by reduced accumulation rates.

Experiment D

This experiment also used the final, steady-state ice sheet of experiment A as an initial condition, and applied the following altered boundary-condition constant

This experiment assessed ice-sheet response to a stepped accumulation-rate change. The change reduces only the area over which the maximum accumulation rate operates. It will therefore lead only to the changes that are caused by reduced ice-sheet span, which contrasts with the two-fold changes of experiment C.

Experiment F

This experiment repeated experiment A (zero ice initial conditions) with the following alteration to the boundary-condition constants:

Experiment G

This experiment aimed at assessing the response of ice-sheet models to the incorporation of basal slip. It repeated experiment A using a linear sliding law

where B is a free parameter whose value was set to 1 × 10^–3 ma^–1 Pa^–1 that of basal melt.

Experiment H

This experiment repeated experiment G using

where T _pmp is the melting point for ice. This experiment assessed the full effect of the slip and its interaction with the distribution of basal temperature.

Experiment A

An initial comparison of global quantities such as volume and area, as well as divide thickness and basal temperature, is given in Table 4. These values represent the equilibrium state, which was normally reached after approximately 50 kyr.

Table 4.

Results for basic glaciological quantities after 200 kyr in experiment A

In this experiment, the ice-sheet area is fixed by the spatial distribution of accumulation and ablation. Volume and divide thickness therefore show a fairly good degree of comparability between models. However, there are large differences between the thermal characteristics (in particular, the fraction of the bed at pressure-melting point) predicted by individual models.

Figure 1 shows the predicted steady-state patterns of basal temperature for each model in the intercomparison. It is clear that the models can be divided into two groups. The first group predicts radial symmetry in basal temperature (models W, Y and Z fall into this group), while this symmetry breaks down in the second group (this occurs with increasing severity in models X, U, R, T, Q, S and V). In the latter group, the border between basal ice which is frozen and melting is broken by spokes of cold ice extending into the melt zone. This pattern is also present in the distributions of velocity and the ice-flow factor (A in Equation (3)), and introduces slight irregularities into the ice-thickness distribution.

Fig. 1.

Predicted steady-state basal temperatures in experiment A for each model in the intercomparison. Temperatures are in K and uncorrected for melting-point variation. The ice-covered area is shaded grey. The units of the x and y axes are kilometers.

Figure 2 shows the spatial distribution of steady-state thickness, ice-flux magnitude and flow factor for model W. This set of results was chosen as being typical of the group of models which exhibited radial symmetry.

Fig. 2.

Predicted steady-state distributions of (a) ice thickness in m, (b) flux magnitudes in m² a⁻¹ (calculated as the 2-norm of the ice-flux vectors ) and (c) flow factor A in 10⁻²⁵ Pa⁻³a⁻¹ for model W in experiment A. Only the lower left quadrant of the domain is shown, and the ice-covered area is shaded grey. The units of the x and y axes are kilometers.

Experiment F

Experiment F forms a useful comparison to A because the only difference is that F has a surface air-temperature boundary condition which is 15 K cooler. Figure 3 shows the predicted steady-state patterns of basal temperature for each model in experiment F. The development of cold-ice spokes within the basal melt zone is now evident in the majority of models. The form of these spokes differs from 12 to 16 broad, well-defined zones (e.g. models Q, U, V, W, X and Z) to numerous narrower fingers (e.g. S and T). Models R and Y lose their radial symmetry but do not develop the spoked pattern. The reduction in surface air temperatures appears to have enhanced the instability which was present in some models during experiment A, to the extent that nearly all models are now affected.

Fig. 3.

Predicted steady-state basal temperatures in experiment F for each model in the intercomparison.

Figure 4 indicates that the patterning evident in basal temperature is reflected in the flow and form of the ice sheet predicted by model W. Cold-ice spokes are associated with reduced rates of ice flow, reduced flow-factor values and inflections in the ice-sheet surface. These relationships are also found in the output of other models with well-defined spokes.

Fig. 4.

Predicted steady-state distributions of (a) ice thickness, (b) flux magnitudes and (c) flow factor for model W in experiment F.

Experiment B

In Experiment B, the steady-state ice sheet of experiment A experiences a 5 K warming in surface air temperature. Table 5 summarizes the equilibrium changes that this warming invokes in the ice sheet. The vast majority of the warming (89–99%) is transferred to the bed. This leads to an expansion of the area of basal melt (total area remaining constant) and a general thinning of the ice sheet (5% at the divide and 3% globally). Time series of the changes in ice thickness and basal temperature at the divide are shown in Figure 5. They generally require 40 kyr to take full effect.

Table 5.

Differences between experiments B and A

Fig. 5.

Time series of thickness (a) and basal-temperature (b) change at the divide (750 km, 750 km) during the first 40 kyr of experiment B for each model in the intercomparison.

Figure 6 shows the predicted patterns of basal temperature at the end of experiment B. The group of models (X, U, R, T, Q, S and V) which exhibited the development of cold spokes of ice in experiment A are now all virtually radially symmetric. This implies that the 5 K warming was sufficient to remove the instability generating these irregularities in flow, and marries nicely with the results from experiment F in which cooling exacerbated the instability.

Fig. 6.

Predicted steady-state basal temperatures in experiment B for each model in the intercomparison.

Experiment C

In experiment C, a uniform shift in the accumulation/ablation field is imposed so that the local mass balance is reduced by 0.25 m a⁻¹ everywhere. Table 6 summarizes the equilibrium changes that this initiates. These changes arise through both a reduction in local accumulation/ablation and a slight reduction in the size of the accumulation area.

Table 6.

Differences between experiments C and A

The reduction in accumulation leads to a general thinning (13% at the divide) and areal contraction of the ice sheet. The basal temperature at the divide warms by 4 K because of the reduction in vertical advection of cold ice consequent on reduced accumulation (the reduction in thermal insulation consequent on reduced ice thickness appears to play a secondary role). The reduced fraction of the ice bed which is melting (−27%) is a consequence of the reduced spatial extent of the ice sheet (−19%), which appears to outweigh the local warming occasioned by decreased vertical advection. The steady-state patterns of basal temperature in experiment C (not shown here), although not perfectly symmetrical, showed a considerable reduction in the extent of their cold-ice spokes. Time series of the changes in ice thickness and basal temperature at the divide are shown in Figure 7.

Fig. 7.

Time series of thickness (a) and basal-temperature (b) change at the divide during the first 80 kyr of experiment C for each model in the intercomparison.

Experiment D

The accumulation/ablation field in experiment D is altered by reducing the size of the accumulation area to the same extent as in experiment C but without changing the accumulation rate. This experiment therefore investigates the effect of changes in ice-sheet size in isolation from changes in accumulation rate. Table 7 summarizes the equilibrium changes in this experiment.

Table 7

Differences between experiments D and A

The changes generated by the reduction in accumulation-area size are small compared to those in experiment C. Volume and area fall by 12% and 10%, respectively. The thickness at the divide is reduced by 2%, which leads to a reduction in thermal insulation and a consequent basal temperature drop of 0.2 K. There is no clear pattern to the predicted changes in the fraction of the bed at melting point. The steady-state patterns of basal temperature in experiment D (not shown here) remain similar to those of experiment A. Time series of the changes in ice thickness and basal temperature at the divide are shown in Figure 8.

Fig. 8.

Time series of thickness (a) and basal-temperature (b) change at the divide during the first 80 kyr of experiment D for each model in the intercomparison.

Experiment G

This experiment incorporated basal slip using a parameterization based solely on local gravitational driving stresses (Equation (9)). The constants used in the parameterization ensured that ice flow was predominantly (approximately 90%) by slip and that it occurred under the entire ice sheet irrespective of basal temperature.

Table 8 summarizes this experiment. The steady-state values in the table were reached after approximately 10 kyr. The table indicates a far greater degree of agreement between models than was evident in experiment A. Figure 9 shows the predicted steady-state patterns of basal temperature for each model in experiment G. The vast majority of models now show radial symmetry.

Table 8.

Results for basic glaciological quantities after 200 kyr in experiment G

Fig. 9.

Predicted steady-state basal temperatures in experiment G for each model in the intercomparison.

Model U is an outlier in this experiment and is generally too thick.

Experiment H

This experiment changed the parameterization of basal slip such that it occurred only where the basal ice was at melting point (Equation (10)).

Table 9 summarizes this experiment. It differs from the other experiments reported here in that a time-independent equilibrium was not produced by the majority of models. Variability was concentrated in the basal thermal regime and is illustrated in Figure 10. This variability could be seen to a far lesser extent in the ice-thickness and volume time series.

Table 9.

Results for basic glaciological quantities after 200 kyr in experiment H

Fig. 10.

Time series of changing basal-melt area fraction during the course of experiment H for each model in the intercomparison.

Figure 11 shows the basal temperature patterns at the end of this experiment, which are a snapshot of the models’ behaviour in this case rather than a long-term steady state. The breakdown of radial symmetry is virtually complete for all models. The degree to which experiment A and F’s coherent spokes of cold basal ice are developed varies between models. Models V, W and Z have a relatively small number of well-defined spokes, while in other models this pattern is complicated by small-scale variability (increasingly so in models S, T, X, Y, Q and R). Model U does not display this instability.

Fig. 11.

Predicted basal temperatures in experiment H at 200 kyr for each model in the intercomparison.

Loss of radial symmetry: indications of local thermomechanical instability

The SGEs represent the first model intercomparison exercise which explicitly investigates the effects of thermomechanical coupling. The main consequence of this is that several potential feedbacks are introduced which could not occur in the uncoupled experiments reported by Reference Huybrechts and PayneHuybrechts and others (1996). The most important result is that thermomechanical coupling appears to generate instabilities in modelled ice flow. Experiments A, B and F can be used to investigate the effects of this instability because they differ only in the temperature boundary condition applied to the ice surface. Divide air temperature is 238 K in experiment A, 223 K in experiment F and 243 K in experiment B (which differs from the other two experiments in that it takes the results of A, and not zero ice cover, as an initial condition).

The instability appears as a series of spokes of cold ice which break the expected radial symmetry of the experiments. The details of the patterns produced differ between models: the number of spokes differs as does the character of the boundary between an individual spoke and the surrounding “warm” ice. In addition, the instability manifests itself in a varying number of models for each experiment. With the exception of model Y, however, all models exhibit spoked patterning at some point. Although basal temperature fields illustrate the effect most clearly, the instability does lead to spoked patterning in the ice-velocity and flow-factor fields. It also causes local changes in ice thickness (and hence surface elevation).

Reference Payne and DongelmansPayne and Dongelmans (1997) identified similar behaviour in a series of experiments utilizing a rectangular ice sheet. They attribute the instability to the coupling of ice temperature and flow. This type of feedback was first investigated by Clarke and others (1977) who introduced the term “creep instability” which links ice flow in Equation (2) with ice-temperature evolution in Equation (4). This is a consequence of the dissipation term in Equation (4), which relates to the heat generated by ice flow, and the temperature dependence of the ice-flow factor in Equation (2) found using Equation (3). There is the possibility of an explosive increase in ice velocity as small perturbations in flow trigger enhanced warming, higher ice-flow factors and faster ice velocities. This feedback is constrained only by the eventual onset of melting.

It is believed that a localized form of this effect is responsible for the cold-ice spokes described above. This relies on the process of creep instability exaggerating initially small differences in temperature between adjacent gridpoints. However, a further effect may also be involved, which is the link between ice thickness and flow. The enhanced ice velocities discussed above imply that the same flux of ice can now be discharged through a reduced ice thickness. Areas affected by creep-induced velocity increase may therefore suffer from drawdown and local ice-surface lowering. This could, in turn, lead to a further concentration of flow because ice flow following the local ice-surface gradient would become channelled towards the area of enhanced ice flow. The thermomechanical coupling inherent in the current experiments can therefore be used to explain the development of the spoked patterns. Further work is needed, however, in order to establish whether this proposed mechanism is correct.

Experiments A, B and F indicate that there is a relationship between the severity of patterning and overall ice-sheet temperature as controlled by the air-temperature boundary condition. It is clear that the formation of cold-ice spokes is strongest in the coldest experiment and progressively disappears in the warmer experiments. The inverse nature of this temperature dependence is counter-intuitive: one would expect a flow instability to be more important in warmer, faster-flowing ice (especially if it has a partially numerical origin). Accumulation rate can also affect spoke formation through its control on the vertical advection term in the temperature-evolution equation (Equation (4)) as shown in experiment C. The patterning was not found in earlier EISMINT intercomparisons (Reference Huybrechts and PayneHuybrechts and others, 1996) simply because it relies on thermomechanical coupling which the earlier work omitted.

An explanation of this dependence awaits further work. However, one possibility is that the spatial location of the boundary separating cold basal ice from that at melting point is determined by overall ice-sheet temperature, and when this boundary lies closer to the steep surface slopes of the ice-sheet margin the effects of creep instability on ice flow and its deflection are important. In contrast, when the boundary is located nearer the flat surface slopes of the divide, the instability cannot produce the same changes in flow direction.

The results of experiment B are particularly interesting because they show that the instability can be reversed. The initial conditions for this experiment come from experiment A, and in some cases show patterning. However, through the course of experiment B the patterning disappears. This suggests that the cold-ice spokes are not a permanent, numerical feature but arise through the time-dependent interaction of ice flow and temperature.

The present study shows that thermomechanical instabilities can develop within ice-sheet models and affect their predictions. However, there remains the question whether the processes discussed here rely entirely on problems associated with numerics (or the inapplicability of the reduced model) or reflect processes which occur in the real world. If the instability does reflect processes which are not entirely numerical then there exists an intriguing link with the formation of ice streams.

Further work is obviously needed to answer this question. A few ideas can, however, be suggested here. The structure of the cold-ice spokes produced by some models implies that numerical details play a role. Examples are the prevalence of spokes running along the x and y axes in some models, and high-frequency noise (non-linear instability) in others. This intercomparison focuses attention on the details of how thermomechanical coupling is implemented within individual models, in particular, the discretization used to represent dissipation in Equation (4) and the effect of the commonly used parameterization of Equation (3) from Reference Paterson and BuddPaterson and Budd (1982). Longitudinal stresses are likely to be important at the boundaries of the the cold-ice spokes and the surrounding faster-flowing, warm ice. The consequences of their omission in the assumed stress regime on which the models are based are therefore another area needing inquiry.

Response to stepped changes in climate

The second set of results to come from the intercomparison stems from the experiments in which a stepped change in climate was imposed (experiments B–D). These experiments address one of the key uses of ice-sheet models: providing context for ice-core palaeoenvironmental reconstructions.

Experiment B imposed a uniform 5 K warming on the equilibrium ice sheet produced in experiment A. This led to a general thinning of the ice sheet, which amounted to 2–3% on average and 4–5% at the divide. The overall area of the ice sheet remained unchanged because it was determined solely by the unaltered snow-accumulation parameterization (Equation (7)). The thinning occurs as a consequence of globally warmer ice which flows faster (Equation (3)) and builds up reduced thicknesses. At the divide, the majority of the 5 K warming is felt by the basal ice. Reduced ice thickness leads to reduced thermal insulation which partially cools the base and accounts for the difference between imposed warming and basal warming. There is also a large increase in the area of the base which suffers melt. The temporal response requires approximately 40 kyr and has a tanh-type profile with a maximum rate of change at 20 kyr.

Experiment C imposed a uniform reduction in accumulation/ablation rate of 0.25 m a⁻¹. This led to the expected thinning and retreat of the ice sheet. At the divide, basal temperatures warm by 3–4 K. This is a consequence of the reduced vertical ice advection caused by the reduction in accumulation rate. However, the ice sheet as a whole cools and suffers a reduction in the area experiencing basal melt. This occurs because the smaller ice sheet both has a reduced insulating effect and generates less dissipational heat (with the reduced ice fluxes now involved). In this and the other climate-change experiments, there appears to be no relationship between the existence of spoked patterning and a model’s response to climate change. This implies that, while the patterning is locally important, it does not significantly affect the global behaviour of an ice sheet.

In experiment D, the size of the accumulation area is reduced without changing maximum accumulation rate. This change also occurred in experiment C, but its effects were overwhelmed by those of the accumulation-rate reduction. There is a slight thinning and retreat of the ice sheet, as expected. This is accompanied by a slight cooling of the ice sheet, which is felt at the divide (0.1–0.2 K) and as a reduction in basal melt area (9–10%). These temperature changes are brought about by the decreased insulating effect of the thinner ice mass.

The range of conflicting responses generated in these experiments illustrate the complexity introduced with thermomechanical coupling. They therefore emphasize the need for numerical ice-sheet modelling in interpreting both ice-sheet response to climate change and the information obtained from ice cores.

Incorporation of basal slip

Experiments G and H investigate the effects of incorporating basal slip. In the former, the whole ice sheet experiences basal slip which amounts to approximately 90% of total velocity. This essentially simplifies the ice-flow problem by replacing the non-linear flow law of Equation (2) with a linear flow law (Equation (9)). In addition, the interaction between ice temperature and flow is lessened considerably. Most of the models cope with this problem well and there is a very close clustering of their results. This experiment also tests how well a model’s temperature calculation copes with plug flow. Again there is very close agreement between most of the models.

Experiment H allows sliding only where the ice-sheet base reaches melting point. This is perhaps more realistic than the assumption used in experiment G, because basal meltwater is likely to lubricate the ice-sheet bed. However, this parameterization does introduce a severe discontinuity between ice which is frozen to its bed and does not experience slip, and adjacent ice at melting point which slips. The thermomechanical effects of the discontinuity have been discussed by Reference PaynePayne (1995) and Reference Payne and DongelmansPayne and Dongelmans (1997). It essentially introduces a second means by which ice temperature and flow can interact. This experiment produces spoked patterning similar to that described for experiment F above. The principal difference between the two experiments is that experiment H produces time-dependent behaviour, which arises because the basal-slip discontinuity migrates as a consequence of a “nick point” developing in the ice-sheet surface profile (Reference PaynePayne, 1995).

Basal slip parameterizations similar to Equation (9) are widely used and rely on an assumed local balance between gravitational driving stress and basal shear stress. However, the known importance of longitudinal stresses in areas of fast ice flow means that a fuller, non-local momentum balance should be employed. Examples of this more realistic approach are Reference MacAyealMacAyeal (1989) and Reference Marshall and ClarkeMarshall and Glarke (1997).

Model intercomparison

The experiments reported here reveal points of both difference and similarity between the ten models involved. The main difference is associated with the the prediction of spoked patterns for quantities such as basal temperature, ice flux and ice-flow factor. There are differences in both the exact form taken by the patterning and the experiments in which it appears. This emphasizes the need to focus future model-development work on both the numerical discretization of these thermomechanical processes within ice-sheet models and the effect of omitting longitudinal stresses from such models. This work should also address issues of numerical stability and their possible role in generating the observed patterning.

It should, however, be stressed that these experiments have identified a major qualitative similarity between models, in that nine out of ten models do show some form of patterning. The climate-change experiments also reveal strong similarities between models. The magnitudes and directions of response to imposed climate change are consistently reproduced between models. This is often in spite of the existence of several competing processes. Further, the presence or absence of spoked patterning does not appear to have a significant impact on the global response of the ice-sheet models to climate change. These findings should encourage the use of ice-sheet models in applications such as response to future climate change and the interpretation of ice-core data.