Skip to main content Accessibility help


  • Access
  • Cited by 75



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Results from the EISMINT model intercomparison: the effects of thermomechanical coupling
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Results from the EISMINT model intercomparison: the effects of thermomechanical coupling
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Results from the EISMINT model intercomparison: the effects of thermomechanical coupling
        Available formats
Export citation


This paper discusses results from the second phase of the European Ice Sheet Modelling Initiative (EISMINT). It reports the intercomparison of ten operational ice-sheet models and uses a series of experiments to examine the implications of thermomechanical coupling for model behaviour. A schematic, circular ice sheet is used in the work which investigates both steady states and the response to stepped changes in climate. The major finding is that the radial symmetry implied in the experimental design can, under certain circumstances, break down with the formation of distinct, regularly spaced spokes of cold ice which extended from the interior of the ice sheet outward to the surrounding zone of basal melt. These features also manifest themselves in the thickness and velocity distributions predicted by the models. They appear to be a common feature to all of the models which took part in the intercomparison, and may stem from interactions between ice temperature, flow and surface form. The exact nature of these features varies between models, and their existence appears to be controlled by the overall thermal regime of the ice sheet. A second result is that there is considerable agreement between the models in their predictions of global-scale response to imposed climate change.


The European Ice Sheet Modelling Initiative (EISMINT) is a programme funded by the European Science Foundation, which aims to examine the critical links between global climate change and ice sheets by improving mathematical modelling in a number of key areas. One of these areas is the intercomparison of the thermomechanical ice-sheet models which have been developed by a variety of groups over the last decade or so.

This paper reports findings from the second phase of EISMINT model intercomparison (January 1996–December 1997) which culminated in a workshop in Grindelwald, Switzerland (24–28 September 1997). The simplified geometry experiments (SGEs, in which only very simplified, regular boundary conditions and model domains are employed) formed one part of the intercomparison and were partnered by experiments related to the modelling of the Greenland ice sheet, the Antarctic ice sheet, ice shelves and the grounding line. The results from these experiments will be reported separately.

This paper follows Huybrechts and others (1996) but differs slightly in its aim. The aim of Huybrechts and others (1996) was to provide a benchmark to aid in the development of new ice-sheet model codes, while the aim of the present paper is to highlight some of the problems associated with coupling ice-temperature and flow modelling (as such it is a level one intercomparison by the definition of Huybrechts and others (1996)). The models used in Huybrechts and others (1996) included both ice-temperature and flow evolution, but the two evolved independently and were not allowed to interact. The results discussed in this paper are especially significant in that they were independently reproduced by a group of models which use a wide range of numerical implementations. The findings should therefore reflect general phenomena and not numerical features specific to a certain model.

A total of ten groups took part in the SGEs. Table 1 identifies these groups and the key authors associated with their work.

Table 1. The groups which took part in the EISMINTSGEs

Description of the Experiments

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 (Glen, 1955), allows horizontal velocity (ui (z)) to be calculated from


where ui (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” (Hindmarsh 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 Huybrechts and others (1996). As in the original benchmark, the ice-accumulation/ablation rate (M in m a−1) 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 Huybrechts 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 km2 horizontal grid (formerly 50 × 50 km2). The model domain remains 1500 × 1500 km2 (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 Paterson 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−1 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.


The results of the intercomparison will be discussed anonymously. The ten groups which submitted results will consistently be referred to as groups Q–Z in the tables, figures and discussion below.

All models successfully simulated the divide migration under changing climate forcing, with the time-scales of surface and temperature adjustment robustly reproduced (experiment E). Similarly, experiments with more realistic and complicated bed topography gave straightforward results with little variability between groups (Experiments I–L). We do not discuss these experiments further in this paper.

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 m2 a−1 (calculated as the 2-norm of the ice-flux vectors ) and (c) flow factor A in 10−25 Pa−3a−1 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−1 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.


The discussion of these results can be divided into four themes: the loss of radial symmetry in some experiments; the response to stepped changes in climate; the incorporation of basal slip; and the degree of comparability between models.

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 Huybrechts 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).

Payne 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 (Huybrechts 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 Paterson 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−1. 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 Payne (1995) and Payne 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 (Payne, 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 MacAyeal (1989) and Marshall 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.


The intercomparison has proven useful in two major respects. First, it highlights the possible existence of a thermomechanical instability in the flow of ice sheets. The nature of an intercomparison exercise meant that we could not explore the causes of this instability fully. However, further work in this direction is forthcoming (Payne and Baldwin, 2000). We are unable to say at this juncture whether the instability is purely numerical and therefore has no real-world equivalent, or whether it reflects a physical process occurring in real ice sheets. If the latter is true, then there are links to the formation of ice streams and the interpretation of the geomorphology of formerly glaciated areas (Punkari, 1984; Clark, 1993; Kleman and others, 1997; Payne and Baldwin, 1999). However, irrespective of whether the instability is physical or numerical, these findings cast doubt on the details of the temperature regimes predicted by large-scale ice-sheet models of the type routinely applied to the Greenland and Antarctic ice sheets.

The second outcome is that, despite the surprise discovery of the thermomechanical instability, the majority of models agree closely with one another. This is true of both equilbrium characteristics such as ice-thickness distribution (and hence ice volume), and the response of divide thickness and temperature to stepped changes in climate. This is encouraging, because the typical applications of large-scale ice-sheet models are to predict the response of ice volume to changes in climate (and its effect on global sea levels), and to identify potential ice-core drilling sites or place the results from existing ice cores in context. The response of ice divides to a simple climate change is complicated because several, opposing effects are often triggered. It is therefore particularly encouraging that the models agreed well with one another.


We would like to extend our thanks to the European Science Foundation for its support of EISMINT and to P. Rowe for her work in organizing the initiative. We would also like to acknowledge the work of C. Doake and the other members of the EISMINT Organizing Committee. Finally, we would like to thank D. Dahl-Jensen, C. Schott Hvidberg and an anonymous referee for the many helpful comments made on the original draft of this paper.


Galov, R. and Marsiat, I.. 1998. Simulations of the Northern Hemisphere through the last glacial–interglacial cycle with a vertically integrated and a three-dimensional thermomechanical ice-sheet model coupled to a climate model. Ann. Glaciol., 27, 169176.
Glark, G. D. 1993. Mega-scale glacial lineations and cross-cutting ice-flow landforms. Earth Surf. Processes Landforms, 18(1), 129.
Glarke, G. K. G., Nitsan, U. and Paterson, W. S. B.. 1977. Strain heating and creep instability in glaciers and ice sheets. Rev. Geophys. Space Phys., 15(2), 235247.
Fastook, J. L. and Prentice, M.. 1994. A finite-element model of Antarctica: sensitivity test for meteorological mass-balance relationship. J. Glaciol., 40(134), 167175.
Glen, J. W. 1955. The creep of polycrystalline ice. Proc. R. Soc. London, Ser. A, 228(1175), 519538.
Greve, R. 1997. Application of a polythermal three-dimensional ice sheet model to the Greenland ice sheet: response to steady-state and transient climate scenarios. J. Climate, 10(5), 901918.
Hindmarsh, R. G. A. and Payne, A. J.. 1996. Time-step limits for stable solutions of the ice-sheet equation. Ann. Glaciol., 23, 7485.
Huybrechts, P. 1990. A 3-D model for the Antarctic ice sheet: a sensitivity study on the glacial–interglacial contrast. Climate Dyn., 5(2), 7992.
Huybrechts, P., Payne, T. and The EISMINT Intercomparison Group. 1996. The EISMINT benchmarks for testing ice-sheet models. Ann. Glaciol., 23, 112.
Kleman, J., Hättestrand, C., Borgström, I. and Stroeven, A.. 1997. Fennoscandian palaeoglaciology reconstructed using a glacial geological inversion model. J. Glaciol., 43(144), 283299.
MacAyeal, D. R. 1989. Large-scale ice flow over a viscous basal sediment: theory and application to Ice Stream B, Antarctica. J. Geophys. Res., 94(B4), 40714087.
Marshall, S. J. and Clarke, G. K. C.. 1997. A continuum mixture model of ice stream thermomechanics in the Laurentide ice sheet. 1. Theory. J. Geophys. Res., 102(B9), 20,59920,614.
Paterson, W. S. B. and Budd, W. F.. 1982. Flow parameters for ice sheet modelling. Cold Reg. Sci. Technol., 6(2), 175177.
Payne, A. J. 1995. Limit cycles in the basal thermal regime of ice sheets. J. Geophys. Res., 100(B3), 42494263.
Payne, A. J. and Baldwin, D. J.. 1999. Thermomechanical modelling of the Scandinavian ice sheet: implications for ice-stream formation. Ann. Glaciol., 28, 8389.
Payne, A. J. and Baldwin, D. J.. 2000. Analysis of ice-flow instabilities identified in the EISMINT intercomparison exercise. Ann. Glaciol., 30, 204219.
Payne, A. J. and Dongelmans, P. W.. 1997. Self-organization in the thermomechanical flow of ice sheets. J. Geophys. Res., 102(B6), 12,21912,233.
Punkari, M. 1984. The relationship between glacial dynamics and tills in the eastern part of the Baltic Shield. Striae, 20, 4954.
Ritz, C., Fabre, A. and Letréguilly, A.. 1997. Sensitivity of a Greenland ice sheet model to ice flow and ablation parameters: consequences for the evolution through the last glacial cycle. Climate Dyn., 13(1), 1124.
Tarasov, L. and Peltier, W. R.. 1999. The impact of thermo-mechanical ice sheet coupling on a model of the 100 kyr ice-age cycle. J. Geophys. Res., 104(D8), 95179545.