1. Introduction

Heinrich events (HEs), which have been discovered in North Atlantic sediments as layers of ice-rafted debris (IRD), are associated with quasi-periodic episodes of massive iceberg discharge from the Laurentide ice sheet through Hudson Bay and Hudson Strait into the Atlantic Ocean (Reference HeinrichHeinrich, 1988; Reference BondBond and others, 1992; Reference BroeckerBroecker, 1994; and rews, 1998; Reference Clarke, Marshall, Hillaire-Marcel, Bilodeau, Veiga-Pires, Clark, Webb and KeigwinClarke and others, 1999). Six major events, labelled H1–H6, have been identified for the Wisconsinan, with a recurrence interval of 7–13 ka. HEs appear during the last glacial cycle, and have also been detected during the other glacials through massive IRD. Inspecting marine records for the last five glacial cycles, Reference McManus, Oppo and CullenMcManus and others (1999) found a dramatic increase in IRD when (relative) benthic oxygen isotopes, which are an indicator for ice volume, exceed a certain threshold.

HEs are regarded as profound and catastrophic events, but are still poorly understood. Different glaciological mechanisms for ice-sheet instability, such as large-scale ice-stream surging (Reference MacAyealMacAyeal, 1993), ice-shelf break-up (Reference HulbeHulbe, 1997; Reference Hulbe, MacAyeal, Denton, Kleman and LowellHulbe and others, 2004) and tidewater instability (Reference Meier and PostMeier and Post, 1987), or a combination of these, have been considered. Based on a critical discussion of these possibilities and some computer modelling, Reference Clarke, Marshall, Hillaire-Marcel, Bilodeau, Veiga-Pires, Clark, Webb and KeigwinClarke and others (1999) concluded that episodic surging of a large ice stream in Hudson Strait is the most plausible mechanism. In addition to their connection to ice dynamics, HEs also show a clear connection to past climate change. They tend to occur at the culmination of a longer-term cooling cycle (Bond cycle) and are followed by a rapid warming. However, a number of possible triggers, which can give the ice sheet the final push when it is ready to discharge, are proposed by Reference BondBond and others (1993) and Reference AlleyAlley and others (2006). Such triggering does not contradict the idea that HEs are essentially an intrinsic ice-sheet instability. Triggering can also explain the observed synchrony of IRD allocated to the Laurentide ice sheet and to other ice sheets. The possibility of such a synchronization has been demonstrated by Reference Calov, Ganopolski, Petoukhov, Claussen and GreveCalov and others (2002).

Reference MacAyealMacAyeal (1993) proposed a ‘binge/purge’ free oscillatory mechanism, which explains HEs as transitions between two modes of operation of ice sheets: slow movement of ice over a frozen base vs a fast sliding mode when the ice bed is molten. In subsequent studies by Reference PaynePayne (1995), Reference Greve and MacAyealGreve and MacAyeal (1996) and Reference Hindmarsh and Le MeurHindmarsh and Le Meur (2001), multi-millennial internal oscillations of ice sheets were simulated with two-dimensional (one vertical and one horizontal direction) models. Reference Marshall and ClarkeMarshall and Clarke (1997a, Reference Marshall and Clarkeb) developed a three-dimensional (3-D) model in order to investigate the role of ice streams in Hudson Strait in Heinrich-type events. However, they found only small-scale instabilities, restricted to the area of the mouth of Hudson Strait.

Reference Calov, Ganopolski, Petoukhov, Claussen and GreveCalov and others (2002) demonstrated that large-scale instabilities of the Laurentide ice sheet resembling HEs in periodicity, amplitude, spatial extent and discharge rate can be simulated with a 3-D dynamic/thermodynamic ice-sheet model (SICOPOLIS). In their simulations, a full Heinrich cycle is characterized by four distinct phases. During the recovery phase after a previous surge, the temperature at the ice bed in Hudson Bay is well below the pressure-melting point, and the ice flows by slow deformation movement. After the ice sheet over Hudson Bay has become sufficiently thick, the temperate basal area expands rapidly from the mouth of Hudson Strait towards Hudson Bay (‘activation wave’; Reference Fowler and SchiaviFowler and Schiavi, 1998). Consequently, fast basal sliding sets in, and the ice sheet enters the surge phase, during which it thins rapidly and its slope towards Hudson Strait decreases. Eventually, melting conditions at the ice base can no longer be sustained, the temperate basal area retreats downstream (‘deactivation wave’) and the surge dies out, thus terminating the Heinrich cycle.

Later, Reference Papa, Mysak and WangPapa and others (2005) used a simplified version of the model by Reference Marshall and ClarkeMarshall and Clarke (1997a,Reference Marshall and Clarkeb) with a parameterization of the basal sliding similar to Reference Calov, Ganopolski, Petoukhov, Claussen and GreveCalov and others (2002), and also succeeded in simulating millennial-scale oscillations of the Laurentide ice sheet under steady external forcing. These advances in modelling have further supported the plausibility of the internal oscillation theory.

This study joins a series of ice-sheet model intercomparison exercises within the European Ice-Sheet Modelling Initiative (EISMINT) (Reference Huybrechts and PayneHuybrechts and others, 1996; Reference PaynePayne and others, 2000; Huybrechts, http://homepages.vub.ac.be/∼phuybrec/eismint/antarctica.html; Ritz, http://homepages.vub.ac.be/∼phuybrec/eismint/greenland.html) and the Ice-Sheet Model Intercomparison Project (ISMIP) (Reference PattynPattyn and others, 2008; further ongoing topics). It is referred to as ISMIP HEINO (Heinrich Event INtercOmparison). The model domain resembles that of the EISMINT Phase 2 Simplified Geometry Experiments (Reference PaynePayne and others, 2000). However, based loosely on Reference PaynePayne (1995), the geometry and boundary conditions have been modified such that they are more similar to the Laurentide ice sheet during the last glacial period, and soft-sediment areas have been introduced which allow rapid basal sliding. The dependence of large-scale ice-sheet instabilities on different atmospheric and basal conditions is investigated. Reference Greve, Takahama and CalovGreve and others (2006) carried out the ISMIP HEINO exercise with SICOPOLIS, and found strong internal oscillations for the standard run and all seven variations specified in the ISMIP HEINO description (http://www.pikpotsdam.de/∼calov/heino/he_setup_2006_11_02.pdf). The objectives of ISMIP HEINO are to demonstrate how far the contemporary large-scale ice-sheet models can reproduce such internal oscillations, to find out under which boundary conditions they appear and to compare the strength and frequency of the oscillations in the different models.

2. Ismip Heino Set-Up

For the ISMIP HEINO experiments, a simplified geometry is employed (Calov and Greve, http://www.pik-potsdam.de/∼calov/heino/he_setup_2006_11_02.pdf). It consists of a horizontal square with a 4000 km side length, in which a circle of 2000 km radius defines the land area prone to glaciation. A soft-sediment area has been chosen to resemble Hudson Bay and Hudson Strait, and is surrounded by hard rock (Fig. 1). The basal topography is flat and rigid (no isostasy).

Fig. 1. Model domain of ISMIP HEINO (Calov and Greve, http://www.pik-potsdam.de/∼calov/heino/he_setup_2006_11_02.pdf).The land area is shown in white, the ocean is shaded grey. The areas inside the square ABCD (‘Hudson Bay’) and the rectangle EFGH (‘Hudson Strait’) correspond to soft sediment. The remaining land area is hard rock.

Since our intercomparison is of internal ice-sheet dynamics, a steady (temporally constant) glacial climate is assumed. The surface mass balance over the land area increases linearly with distance, d, from the centre of the model domain, i.e. from the centre value, b _min, over the domain radius, R, to the margin value, b _max. This yields

(1)

Ablation (surface melting) is not considered, because it is only important for the southern margin of the glacial Laurentide ice sheet, while at the mouth of Hudson Strait calving is the dominant mass-loss process. Calving is modelled implicitly by assuming that the coastward ice mass flux is discharged into the surrounding ocean, where the ice thickness is set to zero. The surface temperature, T _s, is assumed to increase with the third power of distance, d:

(2)

with T _min the minimum temperature and S_T the horizontal gradient.

The geothermal heat flux entering the ice body from below is assumed to be spatially and temporally constant. All parameters and their values are listed in Table 1 (standard set-up ‘ST’).

Table 1. Physical parameters of the standard ISMIP HEINO set-up ‘ST’ (Calov and Greve, http://www.pik-potsdam.de/∼calov/heino/he_setup_2006_11_02.pdf)

For the flow law of ice, the usual power law for incompressible, polycrystalline, isotropic ice (Glen’s flow law) in the form given by Reference Greve and BlatterGreve and Blatter (2009, section 4.3.1) is used. The flow enhancement factor, E, the power exponent, n, and the rate factor, A(T ^′), are listed in Table 1. The temperature, T ^′, relative to pressure melting is given by

(3)

where T is the absolute temperature and β the Clausius–Clapeyron constant (Table 1). Rapid basal sliding is assumed for the sediment area (‘Hudson Bay’ and ‘Hudson Strait’) if the basal temperature reaches the pressure-melting point. The sliding velocity, v_b, is then computed using the linear sliding law:

(4)

which relates v_b to the basal drag,τ_b (both quantities are vectors). Here ρ is the density of ice, g is the gravity acceleration and C _S is the sediment sliding parameter (Table 1). By contrast, slow hard-bed sliding is assumed for the rock area, given by the Weertman-type sliding law

(5)

where C _R is the rock sliding parameter (Table 1) and N _b the basal normal stress. For both cases, no-slip conditions are assumed when the basal temperature is below the pressure-melting point,

(6)

Since in the given model set-up the ice-sheet bed is flat, in Equations (4) and (5) we have simply and , where , and denote components of the stress tensor, T, at the base.

In those models using the shallow-ice approximation (SIA) (Reference HutterHutter, 1983; Reference MorlandMorland, 1984),τ_b is equal to the negative of the basal value of the classical driving stress, so τ_b = ρgHgrad h (where H is the ice thickness and h the free surface elevation; here H = h due to the flat and rigid base). In other force balances, namely full-Stokes flow, higher-order ice-sheet approximations (e.g. Reference HindmarshHindmarsh, 2004; Reference Greve and BlatterGreve and Blatter, 2009) and the shallow-shelf approximation with basal drag (SSA) (Reference MacAyealMacAyeal, 1989; Reference Weis, Greve and HutterWeis and others, 1999), the basal shear stress also depends on the flow field, which provides a physically more adequate model of basal sliding than the SIA with purely geometrically controlled basal drag. Since basal heat production is given by |v_b τ_b|, this also affects the thermal conditions at the ice base.

A justification of the sediment sliding law (Equation (4)) in terms of shear deformation of a linear-viscous sediment layer is given by Reference Greve, Takahama and CalovGreve and others (2006). The hard-rock sliding law (Equation (5)) and the standard value of C _R are those of Reference Greve, Weis and HutterGreve and others (1998). Note also that, even though the numerical value of C _S is smaller than that of C _R, sediment sliding is much stronger than hard rock sliding. This is because the additional factor in Equation (5) is ∼10⁻⁶, so the sliding velocities resulting from Equation (4) are typically three orders of magnitude larger than those resulting from Equation (5) (Reference Greve, Takahama and CalovGreve and others, 2006).

3. Set of Runs

For all ISMIP HEINO runs, the horizontal resolution is 50 km, which leads to 81 × 81 gridpoints in the horizontal plane. The choice of vertical resolution and the time-step are left open. The ice is built up from ice-free initial conditions (zero ice thickness) over a time t =0–200 ka. The set of model runs is as follows:

1. 1. The standard run, ST, is defined by the settings and parameters given in section 2 and Table 1.

2. 2. The minimum surface temperature, T _min, in Equation (2) is varied by ±10 K, resulting in the settings T _min = 223.15 K for run T1 and T _min = 243.15 K for run T2.

3. 3. The surface accumulation in Equation (1) is changed by a factor of 0.5 and 2, giving values b _min = 0.075 m ice equiv. a⁻¹, b _max = 0.15 m ice equiv. a⁻¹ for run B1 and b _min = 0.3 m ice equiv. a⁻¹, b _max = 0.6 m ice equiv. a⁻¹ for run B2 (m ice equiv. indicates metres ice equivalent).

4. 4. The sediment sliding parameter, C _S, in Equation (4) is varied according to C _S = 100 a⁻¹ (run S1), C _S =200 a⁻¹ (run S2) and C _S = 1000 a⁻¹ (run S3).

4. Participating Models

In total, nine models from eleven contributors participated in the intercomparison exercise (Table 2). In the following, the models are referred to in an anonymous fashion by ID letters, which are given along with the main model features in Table 3.

Table 2. Participating models

Table 3. Main features of the participating models (made anonymous). SIA: shallow ice approximation; SSA: shallow shelf approximation with basal drag; H eq.: ice thickness equation; T eq.: temperature equation; FD: finite-difference method; FV: finite-volume method; AA: Arakawa A grid in 3-D (Reference Arakawa and LambArakawa and Lamb, 1977); ABH: Arakawa B grid in the horizontal plane; AC: Arakawa C in 3-D grid; ACH: Arakawa C grid in the horizontal plane; ACH1/2/3: Arakawa C grid in the horizontal plane with method 1/2/3 (Reference Hindmarsh and PayneHindmarsh and Payne, 1996, only applicable to the SIA); σ: sigma transformation in the vertical

Eight of the nine models employ the SIA, whereas in model (c) a combination of the SIA and SSA is used. The discretization of the model equations is carried out by the finite-difference (FD) method in eight of the nine models; model (d) partly uses the finite-volume (FV) method. In addition, the models differ in the numerical grid, the vertical resolution and the applied time-steps.

5. Results

5.1. Standard run ST

A graphical representation of the average ice thickness over the sediment area, H, as a function of time in run ST for each model is given in Figure 2. In order to eliminate spin-up effects, the models were run for 200 ka, and only the last 50 ka of the total simulation times are shown. Models (a)–(g) show clear oscillations, and the surging and recovering phases are clearly visible. However, both phases differ in duration from model to model, and the recurrence times vary between ∼5 and 17 ka. Model (g) shows an interruption of oscillatory behaviour that lasts ∼20 ka, and models (h) and (i) show no noticeable oscillations at all.

Fig. 2. Average ice thickness over the sediment area, H, as a function of time in run ST for each model in the intercomparison. Only the last 50 ka are shown.

The power spectrum of H for each model in the intercomparison is depicted in Figure 3 (normalized to unity for each model separately). The plots confirm the above-mentioned range in periodicity. Further, the spectra vary from very clear, single-peak types to those with one or more secondary peaks. It is interesting to note that even models (h) and (i) have peaks (though with very low spectral power) at ∼7 ka, which is in agreement with the other models that show clearly visible oscillations.

Fig. 3. Power spectrum of the average ice thickness over the sediment area in run ST (see Fig. 2) for each model in the intercomparison. The maximum power has been normalized to unity for each model separately.

Figure 4 shows the fraction of warm-based ice over the sediment area, A/A _sed, as a function of time. This quantity is a measure of the expansion of streaming over the sediment area. For all oscillating models (a)–(g), A/A _sed returns to ∼0 during the recovery phases, whereas maximum values during surges vary from model to model between ∼0.5 and ∼1. The shapes of the curves are highly varied. In the case of the non-oscillating models (h) and (i), A/A _sed remains nearly constant at ∼0.2 and ∼0.15, respectively, so these models remain in a state of slight streaming throughout time.

Fig. 4. Fraction of warm-based ice over the sediment area, A/A _sed, as a function of time in run ST for each model in the intercomparison. Only the last 50 ka are shown.

Four particular time slices are inspected in Figures 5–10, namely for the times t ₁ (maximum average ice thickness over the sediment area during the last 50 ka), t ₂ (minimum average ice thickness over the sediment area during the last 50 ka), t ₃ (minimum average basal temperature over the sediment area during the last 50 ka) and t ₄ (maximum extent of warm-based ice over the sediment area during the last 50 ka), defined in the reference document (Calov and Greve, http://www.pikpotsdam.de/∼calov/heino/he_setup_2006_11_02.pdf). Note that these times differ from model to model because of the irregular appearance of the surges.

Fig. 5. Ice thickness at time t ₁ (maximum average ice thickness over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.

Fig. 6. Ice thickness at time t ₂ (minimum average ice thickness over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.

Fig. 7. Basal temperature (relative to pressure melting) at time t ₃ (minimum average basal temperature over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.

Fig. 8. Basal temperature (relative to pressure melting) at time t ₄ (maximum extent of warm-based ice over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.

Fig. 9. Basal sliding velocity at time t ₃ (minimum average basal temperature over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.

Fig. 10. Basal sliding velocity at time t ₄ (maximum extent of warm-based ice over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.

Comparison between respective panels in Figures 5 and 6 illustrates the amplitude of the oscillations. Models (a)–(d) show almost circularly symmetric ice-sheet geometries at time t ₁, whereas models (e)–(i) have distinct surface depressions over the sediment region. By contrast, the geometries at time t ₂ are rather similar for all oscillating models (a)–(g), with strong surface depressions over the sediment region due to the previous surge. As a consequence of the permanently maintained slight streaming, the respective geometries of models (h) and (i) look rather similar at times t ₁ and t ₂, always exhibiting ice thicknesses between the minima and maxima of the other models.

The basal temperatures shown in Figures 7 and 8 for times t ₃ and t ₄, respectively, both indicate that, away from the sediment region, a ring of warm-based ice occurs close to the margin, and basal temperatures decrease towards the interior of the ice sheet. This is the expected distribution. By contrast, in the sediment region, cold-based conditions prevail for the oscillating models (a)–(g) at time t ₃ (in the middle of the recovery phase), the detailed temperature pattern varying significantly between the models. At time t ₄ (during a surge), warm-based conditions dominate in the sediment region, ranging from ∼50% to ∼100% of the area (as mentioned above). For the non-oscillating models (h) and (i), a channel along ‘Hudson Bay’ and ‘Hudson Strait’ towards the ice margin is constantly at the pressure-melting point, thereby sustaining the weak, but continuous, streaming state.

The distributions of the basal sliding velocity at times t ₃ and t ₄ (Figs 9 and 10) closely reflect that of the basal temperature. Away from the sediment region, the basal sliding velocity increases monotonically from the interior towards the margin. In the sediment region, the oscillating models (a)–(g) show little or no basal sliding at time t ₃ (recovery phase), whereas the basal sliding velocities exceed 1000 m a⁻¹ at time t ₄ (surge). The detailed structure of the surge differs from model to model. The sustained streaming conditions of the non-oscillating models (h) and (i) are also clearly visible.

It is striking that the fields shown in Figures 5–10 are not fully symmetric with respect to the line y = 2000 km, whereas the geometry and boundary conditions are. The strength of symmetry breaking clearly differs from model to model. There are three possible reasons for this phenomenon, namely:

1. 1. usage of non-symmetrically discretized, but mathematically correct, numerical schemes, whose small asymmetries can grow to macroscopic size;

2. 2. usage of the usual symmetrically discretized numerical schemes allows asymmetrical order of floating-point evaluation, and these slight asymmetries can grow to macroscopic size;

3. 3. real bugs in the coding. Since the participating models have been tested in previous model intercomparison exercises, severe coding errors are unlikely, so either of the first two reasons (or a combination of them) is probably the cause of the broken symmetry.

5.2. Variations of the surface boundary conditions (runs T1, T2, B1, B2)

We now discuss the effect of varied surface boundary conditions (temperature, accumulation) on the different model results. For runs T1 (10°C colder), T2 (10°C warmer) and ST (reference for comparison), the average ice thickness over the sediment area, H, is depicted as a function of time for each model in Figure 11. The mean value of H over the shown 50 ka increases with decreasing surface temperature, which is a direct consequence of the temperature dependence of the rate factor (colder ice is stiffer; Table 1). More interestingly, it becomes clear that lower temperatures favour the occurrence of oscillations, whereas higher temperatures hamper them. While seven out of the nine models show oscillations in run ST, oscillations occur for all nine models in the colder run, T1, but only for three models in the warmer run, T2. By contrast, there is no evident correlation between the period of oscillations (if existing) and the surface temperature.

Fig. 11. Average ice thickness over the sediment area, H, as a function of time in runs ST (solid), T1 (dashed) and T2 (dotted) for each model in the intercomparison. Only the last 50 ka are shown.

For the runs with varied accumulation rate, B1 (half accumulation) and B2 (double accumulation), H is shown in Figure 12 (along with ST for comparison). An immediate observation is that the mean value of H over the shown 50 ka increases in the order B1, ST, B2, i.e. with increasing accumulation. Similar to the dependence on surface temperature, the occurrence of oscillations is favoured by lower accumulation rates: seven out of eight models (no results from model (f) for runs B1 and B2) produce oscillations in run B1, but only four models in run B2. Further, the period of oscillations (if existing) is clearly correlated with the accumulation rate. This is because the accumulation rate directly affects the growth time required to build up the ice sheet to the critical thickness at which a surge can be released.

Fig. 12. Average ice thickness over the sediment area, H, as a function of time in runs ST (solid), B1 (dashed) and B2 (dotted) for each model in the intercomparison. Only the last 50 ka are shown. The straight line in (e) belongs to run B2; it should be dotted but appears solid due to the used plot utility. There are no results from (f) for B1 or B2.

5.3. Variations of the basal sliding (runs S1, S2, S3)

The effect of varied basal sliding conditions over the sediment area becomes evident in Figure 13, which shows the average ice thickness over the sediment area, H, as a function of time for each model in runs S1 (20% basal sliding), S2 (40% basal sliding), S3 (200% basal sliding) and ST (reference for comparison). Since sediment sliding is the crucial process for the occurrence of large-scale surges, it is not surprising that a significant influence of the sliding parameter for soft sediment, C _S, on the oscillatory behaviour becomes evident: two out of eight models (no results from model (f) for runs S1–S3) produce noticeable oscillations in run S1, three models in run S2 and six models in run S3. As long as oscillations are present, their amplitudes tend to increase in the order S1, S2, ST, S3 (with increasing sliding). However, the period of oscillations is only affected to a small extent and in a non-systematic fashion.

Fig. 13. Average ice thickness over the sediment area, H, as a function of time in runs ST (solid), S1 (dashed), S2 (dotted) and S3 (dash–dotted) for each model in the intercomparison. Only the last 50 ka are shown. There are no results from (f) for S1, S2 or S3.

6. Discussion and Conclusion

Section 5 has shown that Heinrich-type oscillations occur for a broad range of parameters and for all participating models. Low surface temperatures, low accumulation rates and high sliding velocities over sediment favour oscillatory surges, whereas high surface temperatures, high accumulation rates and low sliding velocities hamper them. However, there are significant differences between the models. At one end of the range, oscillations are produced by models (b) and (d) for all eight runs, whereas, at the other end, model (h) produces oscillations only for run T1 (Fig. 14). Also, strong variabilities in amplitudes, recurrence times and detailed structures of the surge cycles occur. These results are consistent with findings of two of the participating modellers who experimented with different schemes. They observed a dependence of the oscillations on discretization, ranging from strong oscillations to almost non-oscillatory behaviour with permanent ice streaming.

Fig. 14. Illustration of the parameter space for each model in the intercomparison. Upper panel: variation of the surface temperature; middle panel: variation of the surface mass balance; lower panel: variation of the sliding parameter. Black boxes indicate that oscillations occur (criterion: period 5–20 ka, power >10⁷ km² a², distinct peak visible in the power spectrum); white boxes show that there are no oscillations; grey boxes denote no data.

We are aware that simulation of the abrupt transition between slow and fast ice movement at the position where the basal ice changes from cold to warm is problematic in the SIA, which yields a strong localization in the ice velocities. Reference Bueler and BrownBueler and Brown (2009) analysed the situation considering ice flow in an inclined ice slab, and demonstrated that in the case of SIA the vertical velocity becomes infinite in the ice column (or rather along the line in the ice slab) that is situated at the transition between cold- and warm-based ice. It is from this high-frequency transition in basal velocity that some numerical schemes may produce macroscopic asymmetric behaviours, and other schemes do not. To some extent this sharp transition is hidden by the rather coarse horizontal grid used in the computations (Reference Bueler and BrownBueler and Brown, 2009). It is apparent that for non-localized equations, like full-Stokes or SSA, the situation will improve, and the velocities in the column over the transition from cold to warm basal ice will behave more smoothly. However, note that SIA model (d) features explicit smoothing of the abrupt transition between cold-based no-slip conditions and warm-based sliding by linearly increasing the sliding parameters, C _S and C _R, from zero to their regular values (Table 1) as the basal temperature increases from 0.1 K below the melting point to the melting point (Table 3). Despite this smoothing, model (d) produces oscillations for all eight runs, indicating that a sharp transition between no-slip and basal sliding is not essential for oscillating behaviour of SIA models. Further, the combined SIA/SSA model (c) is one of the more oscillation-prone models (Fig. 14), which demonstrates that oscillations are not limited to pure SIA models.

In addition to the standard run and the seven variations specified in the ISMIP HEINO description (Calov and Greve, http://www.pik-potsdam.de/∼calov/heino/he_setup_2006_11_02.pdf) and discussed here, Reference Greve, Takahama and CalovGreve and others (2006) carried out further simulations with the ice-sheet model SICOPOLIS. A doubled horizontal resolution of 25 km was applied, different time-steps of 2, 1, 0.5 and 0.25 years were tested, the sediment mask was rotated by 5, 10, …, 45° around the centre of the quadratic model domain, and a continuous formulation for sub-melt sliding was employed as an alternative to the abrupt switch between no-slip conditions and fully developed basal sliding. The latter involves a temperature parameter, γ. Its meaning is that basal sliding is 1/e ^m times fully developed, warm-based sliding if the basal temperature is below the pressure-melting point by mγ. We found that the oscillations of the simulated ice sheet are very robust with respect to resolution, time-step and grid rotation. Sub-melt sliding with γ = 0.2°C also yields strong oscillations, while sub-melt sliding with γ ≥ 0.5°C produces continuous states of streaming without oscillations, similar to models (h) and (i) for the standard run ST. This gives rise to the conjecture that in these models the oscillations are suppressed by too strong implicit dampening of the transition from warm-based sliding to cold-based no-slip due to the applied numerical scheme.

A solution to the problems of SIA models in adequately capturing the physics in the ice column where the flow changes from slow to fast could be to treat this column as a boundary-layer problem. Alternatively, one could make use of the full-Stokes solution for the entire ice sheet. However, it should be kept in mind that full-Stokes solvers for ice sheets are computationally very expensive. An attempt to carry out the ISMIP HEINO exercise with the full-Stokes model Elmer/Ice (e.g. Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007; Reference Seddik, Greve, Zwinger and PlacidiSeddik and others, 2009) failed because of overwhelming computational demand. Other already existing approaches use a combination of SIA and SSA, like model (c) of this study (Reference Pollard, DeConto, Hambrey, Christoffersen, Glasser and HubbardPollard and DeConto, 2007; Reference Bueler and BrownBueler and Brown, 2009), and are more applicable because of computational efficiency, and because of delocalization due to membrane stresses in these models. Nevertheless, the traditional, pure-SIA models should not be dismissed. Pragmatically speaking, one could try a simple parameterization at the transition from slow to fast flow and regard such modified SIA models as a workaround. The capturing of observed frequencies and recurrence times of Heinrich-type oscillations with such a workaround would be of some advantage. In the end, the applicability depends on the desired application of the model. Heinrich events play an important role in palaeoclimatic research, and a model such as that presented here offers a tool to simulate some of their characteristics.

There are two important lessons to learn from our intercomparison exercise. On the one hand, all models which took part (including the combined SIA/SSA model (c) and the SIA model (d) with a smoothed transition between cold-based no-slip conditions and warm-based sliding) are capable of producing Heinrich-type free oscillations if the boundary conditions are sufficiently favourable. However, the large differences between the results of different models demonstrate that it is a great asset to maintain a plurality of ice-sheet models in the community. This allows us to establish, using model intercomparison, which aspects of simulated ice-sheet dynamics are robust, which are not and where there is a need for improvement.