# Results from the Ice-Sheet Model Intercomparison Project-Heinrich Event INtercOmparison (ISMIP HEINO)

R. Calov, R. Greve, A. Abe-Ouchi, E. Bueler, P. Huybrechts, J. V. Johnson, F. Pattyn, D. Pollard, C. Ritz, F. Saito, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

R. Calov, R. Greve, A. Abe-Ouchi, E. Bueler, P. Huybrechts, et al. Results from the Ice-Sheet Model Intercomparison Project-Heinrich Event INtercOmparison (ISMIP HEINO). Journal of Glaciology, 2010, 56 (197), pp.371-383. 10.3189/002214310792447789. insu-00554220

HAL Id: insu-00554220

https://insu.hal.science/insu-00554220

Submitted on 24 Mar 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Results from the Ice-Sheet Model Intercomparison Project–Heinrich Event INtercOmparison (ISMIP HEINO)

Reinhard CALOV,¹ Ralf GREVE,² Ayako ABE-OUCHI,³ Ed BUELER,⁴

Philippe HUYBRECHTS,⁵ Jesse V. JOHNSON,⁶ Frank PATTYN,⁷ David POLLARD,⁸

Catherine RITZ,⁹ Fuyuki SAITO,¹⁰ Lev TARASOV¹¹

¹Potsdam Institute for Climate Impact Research, PO Box 601203, D-14412 Potsdam, Germany
E-mail: calov@pik-potsdam.de

²Institute of Low Temperature Science, Hokkaido University, Kita-19, Nishi-8, Kita-ku, Sapporo 060-0819, Japan

³Centre for Climate System Research, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8568, Japan

⁴Department of Mathematics and Statistics, University of Alaska Fairbanks, Fairbanks, Alaska 99775-6660, USA

⁵Earth System Sciences & Departement Geografie, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussels, Belgium

⁶Department of Computer Science, Room 417, Social Science Building, University of Montana, Missoula, Montana 59812-5256, USA

⁷Laboratory of Glaciology, CP 160/03, Université Libre de Bruxelles, Avenue F.D. Roosevelt 50, B-1050 Brussels, Belgium

⁸EMS Earth and Environmental Systems Institute, The Pennsylvania State University, University Park, Pennsylvania 16802-2711, USA

⁹Laboratoire de Glaciologie et Géophysique de l'Environnement, CNRS/Université Joseph Fourier - Grenoble 1, 54 Rue Molière, BP 96, 38402 Saint-Martin-d'Hères Cedex, France

¹⁰Japan Agency for Marine-Earth Science and Technology, 3173-25 Showamachi, Kanazawa-ku, Yokohama, Kanagawa 236-0001, Japan

¹¹Department of Physics and Physical Oceanography, Memorial University of Newfoundland, St Johns, Newfoundland A1B 3X7, Canada

ABSTRACT. Results from the Heinrich Event INtercOmparison (HEINO) topic of the Ice-Sheet Model Intercomparison Project (ISMIP) are presented. ISMIP HEINO was designed to explore internal large-scale ice-sheet instabilities in different contemporary ice-sheet models. These instabilities are of interest because they are a possible cause of Heinrich events. A simplified geometry experiment reproduces the main characteristics of the Laurentide ice sheet, including the sedimented region over Hudson Bay and Hudson Strait. The model experiments include a standard run plus seven variations. Nine dynamic/thermodynamic ice-sheet models were investigated; one of these models contains a combination of the shallow-shelf (SSA) and shallow-ice approximation (SIA), while the remaining eight models are of SIA type only. Seven models, including the SIA-SSA model, exhibit oscillatory surges with a period of ~1000 years for a broad range of parameters, while two models remain in a permanent state of streaming for most parameter settings. In a number of models, the oscillations disappear for high surface temperatures, strong snowfall and small sediment sliding parameters. In turn, low surface temperatures and low snowfall are favourable for the ice-surge cycles. We conclude that further improvement of ice-sheet models is crucial for adequate, robust simulations of cyclic large-scale instabilities.

## 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 (Heinrich, 1988; Bond and others, 1992; Broecker, 1994; Andrews, 1998; Clarke 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, McManus 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 (MacAyeal, 1993), ice-shelf break-up (Hulbe, 1997; Hulbe and others, 2004) and tidewater instability (Meier and Post, 1987), or a combination of these, have been considered. Based on a critical discussion of these possibilities and some computer modelling, Clarke 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 Bond and others (1993) and Alley 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
<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

synchronization has been demonstrated by Calov and others (2002).

MacAyeal (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 Payne (1995), Greve and MacAyeal (1996) and Hindmarsh and Le Meur (2001), multi-millennial internal oscillations of ice sheets were simulated with two-dimensional (one vertical and one horizontal direction) models. Marshall and Clarke (1997a,b) 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.

Calov 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'; Fowler 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, Papa and others (2005) used a simplified version of the model by Marshall and Clarke (1997a,b) with a parameterization of the basal sliding similar to Calov 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) (Huybrechts and others, 1996; Payne 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) (Pattyn 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 (Payne and others, 2000). However, based loosely on Payne (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. Greve 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.pik-potsdam.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).

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

$$b = b_{\min} + \frac{b_{\max} - b_{\min}}{R} \times d. \quad (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$:

$$T_s = T_{\min} + S_T d^3, \quad (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’).

For the flow law of ice, the usual power law for incompressible, polycrystalline, isotropic ice (Glen’s flow law) in the form given by Greve 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

$$T' = T + \beta p, \quad (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*<sub>b</sub>, is then computed using the linear sliding law:

$$v_b = -\frac{C_S}{\rho g} \tau_b, \quad \text{for } T'_b = T_0 \text{ and soft sediment,} \quad (4)$$

which relates *v*<sub>b</sub> to the basal drag, *τ*<sub>b</sub> (both quantities are vectors). Here *ρ* is the density of ice, *g* is the gravity acceleration and *C*<sub>S</sub> 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

$$v_b = -\frac{C_R}{\rho g} \frac{|\tau_b|^2}{N_b^2} \tau_b, \quad \text{for } T'_b = T_0 \text{ and hard rock,} \quad (5)$$

where *C*<sub>R</sub> is the rock sliding parameter (Table 1) and *N*<sub>b</sub> the basal normal stress. For both cases, no-slip conditions are assumed when the basal temperature is below the pressure-melting point,

$$v_b = 0, \quad \text{for } T'_b < T_0. \quad (6)$$

Since in the given model set-up the ice-sheet bed is flat, in Equations (4) and (5) we have simply $\tau_b = -(t_{xz}^b, t_{yz}^b, 0)$ and $N_b = -t_{zz}^b$, where $t_{xz}^b, t_{yz}^b$ and $t_{zz}^b$ denote components of the stress tensor, $T$, at the base.

In those models using the shallow-ice approximation (SIA) (Hutter, 1983; Morland, 1984), $\tau_b$ is equal to the negative of the basal value of the classical driving stress, so $\tau_b = \rho g H$ grad $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. Hindmarsh, 2004; Greve and Blatter, 2009) and the shallow-shelf approximation with basal drag (SSA) (MacAyeal, 1989; Weis 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 \cdot \tau_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 Greve and others (2006). The hard-rock sliding law (Equation (5)) and the standard value of $C_R$ are those of Greve 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 $|\tau_b|^2/N_b^2$ in Equation (5) is $\sim 10^{-6}$, so the sliding velocities resulting from Equation (4) are typically three orders of magnitude larger than those resulting from Equation (5) (Greve and others, 2006).

<table><caption>**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)</caption><thead><tr><th>Quantity</th><th>Value</th></tr></thead><tbody><tr><td>Gravity acceleration, g</td><td>9.81 m s<sup>-2</sup></td></tr><tr><td>Density of ice, ρ</td><td>910 kg m<sup>-3</sup></td></tr><tr><td>Power law exponent, n</td><td>3</td></tr><tr><td>Flow enhancement factor, E</td><td>3</td></tr><tr><td>Rate factor, A(T')</td><td>A<sub>0</sub> e<sup>-Q/(RT')</sup></td></tr><tr><td>Pre-exponential constant, A<sub>0</sub></td><td>3.61 × 10<sup>-13</sup> s<sup>-1</sup> Pa<sup>-3</sup><br>(T' ≤ 263.15 K)</td></tr><tr><td>Activation energy, Q</td><td>1.73 × 10<sup>3</sup> s<sup>-1</sup> Pa<sup>-3</sup><br>(T' ≥ 263.15 K)</td></tr><tr><td>Melting point at atmospheric pressure, T<sub>0</sub></td><td>60 kJ mol<sup>-1</sup><br>(T' ≤ 263.15 K)</td></tr><tr><td>Clausius–Clapeyron constant, β</td><td>8.314 J mol<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Universal gas constant, R</td><td>2.1 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Heat conductivity of ice, κ</td><td>2009 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Specific heat of ice, c</td><td>335 kJ kg<sup>-1</sup></td></tr><tr><td>Minimum surface mass balance, b<sub>min</sub></td><td>0.15 m ice equiv. a<sup>-1</sup></td></tr><tr><td>Maximum surface mass balance, b<sub>max</sub></td><td>0.3 m ice equiv. a<sup>-1</sup></td></tr><tr><td>Domain radius, R</td><td>2000 km</td></tr><tr><td>Minimum surface temperature, T<sub>min</sub></td><td>233.15 K</td></tr><tr><td>Surface temperature gradient, S<sub>T</sub></td><td>2.5 × 10<sup>-9</sup> K km<sup>-3</sup></td></tr><tr><td>Sliding parameter</td><td></td></tr><tr><td>– for hard rock, C<sub>R</sub></td><td>10<sup>5</sup> a<sup>-1</sup></td></tr><tr><td>– for soft sediment, C<sub>S</sub></td><td>500 a<sup>-1</sup></td></tr><tr><td>Geothermal heat flux, q<sub>geo</sub></td><td>42 mW m<sup>-2</sup></td></tr><tr><td>Length of year</td><td>31 556 926 s</td></tr></tbody></table>

## 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. The standard run, ST, is defined by the settings and parameters given in section 2 and Table 1.

2. The minimum surface temperature, $T_{min}$, in Equation (2) is varied by $\pm 10$ K, resulting in the settings $T_{min} = 223.15$ K for run T1 and $T_{min} = 243.15$ K for run T2.

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<sup>-1</sup>, $b_{max} = 0.15$ m ice equiv. a<sup>-1</sup> for run B1 and $b_{min} = 0.3$ m ice equiv. a<sup>-1</sup>, $b_{max} = 0.6$ m ice equiv. a<sup>-1</sup> for run B2 (m ice equiv. indicates metres ice equivalent).

4. The sediment sliding parameter, $C_S$, in Equation (4) is varied according to $C_S = 100$ a<sup>-1</sup> (run S1), $C_S = 200$ a<sup>-1</sup> (run S2) and $C_S = 1000$ a<sup>-1</sup> (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><caption>**Table 2.** Participating models</caption><thead><tr><th>Author(s)</th><th>Model</th><th>URL</th><th>Source</th></tr></thead><tbody><tr><td>E. Bueler</td><td>PISM* in SIA mode</td><td>http://www.pism-docs.org/</td><td>Bueler and others (2007)</td></tr><tr><td>R. Greve and R. Calov</td><td>SICOPOLIS†</td><td>http://sicopolis.greveweb.net/</td><td>Greve (1997); http://hdl.handle.net/2115/34755</td></tr><tr><td>P. Huybrechts</td><td>Huybrechts model</td><td>—</td><td>Huybrechts (1990, 2002); Huybrechts and de Wolde (1999)</td></tr><tr><td>J.V. Johnson</td><td>Glimmer</td><td>http://glimmer-cism.berlios.de/</td><td>Rutt and others (2009)</td></tr><tr><td>F. Pattyn</td><td>Pattyn model in SIA mode</td><td>—</td><td>Pattyn (2003)</td></tr><tr><td>D. Pollard</td><td>Pollard model</td><td>—</td><td>Pollard and DeConto (2007)</td></tr><tr><td>C. Ritz</td><td>GRISL‡ in SIA mode</td><td>—</td><td>Ritz and others (2001)</td></tr><tr><td>F. Saito and A. Abe-Ouchi</td><td>IcIES§</td><td>—</td><td>Saito and Abe-Ouchi (2005)</td></tr><tr><td>L. Tarasov</td><td>MUN/UofT GSM¶</td><td>—</td><td>Tarasov and Peltier (1997, 1999, 2002)</td></tr></tbody></table>

*Parallel Ice Sheet Model.

†Simulation COde for POLythermal Ice Sheets.

‡GRenoble model for Ice Shelves and Land Ice.

§Ice sheet model for Integrated Earth system Studies.

¶Memorial University of Newfoundland/University of Toronto Glacial Systems Model.

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.

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.

Figure 4 shows the fraction of warm-based ice over the sediment area, *A/A*<sub>sed</sub>, 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*<sub>sed</sub> returns to ~0 during the recovery phases, whereas maximum values during

<table><caption>**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 (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 (Hindmarsh and Payne, 1996, only applicable to the SIA); *σ*: sigma transformation in the vertical</caption><thead><tr><th>Model ID</th><th>Type</th><th>Method (H eq./T eq.)</th><th>Grid (H eq./T eq.)</th><th>Vertical resolution</th><th>Time-step</th></tr></thead><tbody><tr><td>(a)</td><td>SIA</td><td>FD/FD</td><td>ACH2/AA, σ</td><td>26 layers</td><td>0.1 years</td></tr><tr><td>(b)</td><td>SIA</td><td>FD/FD</td><td>ACH3/AC, σ</td><td>91 layers</td><td>0.25 years</td></tr><tr><td>(c)</td><td>SIA/SSA</td><td>FD/FD</td><td>ACH/ACH, σ</td><td>10 layers</td><td>1 year</td></tr><tr><td>(d)</td><td>SIA</td><td>FD/FV</td><td>ACH3/AC, σ</td><td>65 layers</td><td>0.25 years</td></tr><tr><td>(e)</td><td>SIA</td><td>FD/FD</td><td>ACH3/AA, σ</td><td>41 layers</td><td>0.1 years</td></tr><tr><td>(f)</td><td>SIA</td><td>FD/FD</td><td>ACH2/ACH, σ</td><td>20 layers</td><td>0.1 years</td></tr><tr><td>(g)</td><td>SIA</td><td>FD/FD</td><td>ACH3/AA, σ</td><td>31 layers</td><td>0.2 years</td></tr><tr><td>(h)</td><td>SIA</td><td>FD/FD</td><td>ACH2/AA</td><td>25 m</td><td>0.05–0.25 years</td></tr><tr><td>(i)</td><td>SIA</td><td>FD/FD</td><td>ACH1/ABH, σ</td><td>21 layers</td><td>0.125 years</td></tr><tr><td></td><td></td><td></td><td>Remark: Concentration of layers towards the bottom.</td><td></td><td></td></tr><tr><td></td><td></td><td></td><td>Remark: Adaptive time-stepping. All runs except ST include a thermal bedrock model.</td><td></td><td></td></tr><tr><td></td><td></td><td></td><td>Remark: Concentration of layers towards the bottom. Larger time-step of 0.25 years for the T eq.</td><td></td><td></td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

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.

Four particular time slices are inspected in Figures 5–10, namely for the times $t_1$ (maximum average ice thickness over the sediment area during the last 50 ka), $t_2$ (minimum average ice thickness over the sediment area during the last 50 ka), $t_3$ (minimum average basal temperature over the sediment area during the last 50 ka) and $t_4$ (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.pik-potsdam.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.

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_1$, whereas models (e)–(i) have distinct surface depressions over the sediment region. By contrast, the geometries at time $t_2$ 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_1$ and $t_2$, always exhibiting ice thicknesses between the minima and maxima of the other models.

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

The basal temperatures shown in Figures 7 and 8 for times $t_3$ and $t_4$, 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_3$ (in the middle of the recovery phase), the detailed temperature pattern varying significantly between the models. At time $t_4$ (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

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Ice thickness at time $t_1$ (maximum average ice thickness over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Ice thickness at time $t_2$ (minimum average ice thickness over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Basal temperature (relative to pressure melting) at time $t_3$ (minimum average basal temperature over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Basal temperature (relative to pressure melting) at time $t_4$ (maximum extent of warm-based ice over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.</figcaption></figure>

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_3$ and $t_4$ (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_3$ (recovery phase), whereas the basal sliding velocities exceed $1000 \text{ m a}^{-1}$ at time $t_4$ (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 \text{ 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. usage of non-symmetrically discretized, but mathematically correct, numerical schemes, whose small asymmetries can grow to macroscopic size;

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

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
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Basal sliding velocity at time $t_3$ (minimum average basal temperature over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Basal sliding velocity at time $t_4$ (maximum extent of warm-based ice over the sediment area during the last 50 ka) in run ST for each model in the intercomparison.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

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.

### 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.

<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

**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.

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. Bueler 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 (Bueler 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.

<figure><img src="image_16.png" /><figcaption>**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.</figcaption></figure>
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, Greve 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\gamma$. 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 $\gamma = 0.2°C$ also yields strong oscillations, while sub-melt sliding with $\gamma \ge 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. Zwinger and others, 2007; Seddik 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 (Pollard and DeConto, 2007; Bueler 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.

## ACKNOWLEDGEMENTS

ISMIP (Ice-Sheet Model Intercomparison Project) arose from the Numerical Experimentation Group of CliC (Climate and Cryosphere, a core project of the World Climate Research Programme co-sponsored by the Scientific Committee on Antarctic Research). Constructive reviews by O. Souček and an anonymous reviewer resulted in improvements to the paper. The authors thank S. Dolaptchiev for producing a huge number of raw plots and S. Lubrich for assistance with creating plone-powered password-protected group space accounts for the HEINO results. R.C. thanks A. Ganopolski for many fruitful discussions on Heinrich events. R.G. thanks T. Dunse for clarifying remarks about surging and streaming. F.P. thanks C. Delcourt for her help with the model runs.

## REFERENCES

* Alley, R.B. and 6 others. 2006. Outburst flooding and surge initiation in response to climatic cooling: an hypothesis. *Geomorphology*, **75**(1–2), 76–89.

* Andrews, J.T. 1998. Abrupt changes (Heinrich events) in late Quaternary North Atlantic marine environments: a history and review of data and concepts. *J. Quat. Sci.*, **13**(1), 3–16.

* Arakawa, A. and V.R. Lamb. 1977. Computational design of the basic dynamical processes of the UCLA general circulation model. In Chang, J., ed. *General circulation models of the atmosphere*. New York, Academic Press, 173–265.

* Bond, G. and 13 others. 1992. Evidence for massive discharges of icebergs into the North Atlantic Ocean during the last glacial period. *Nature*, **360**(6401), 245–249.

* Bond, G. and 6 others. 1993. Correlations between climate records from North Atlantic sediments and Greenland ice. *Nature*, **365**(6442), 143–147.

* Broecker, W.S. 1994. Massive iceberg discharges as triggers for global climate change. *Nature*, **372**(6505), 421–424.

* Bueler, E. and J. Brown. 2009. Shallow shelf approximation as a 'sliding law' in a thermomechanically coupled ice sheet model. *J. Geophys. Res.*, **114**(F3), F03008. (10.1029/2008JF001179.)

* Bueler, E., J. Brown and C. Lingle. 2007. Exact solutions to the thermomechanically coupled shallow-ice approximation: effective tools for verification. *J. Glaciol.*, **53**(182), 499–516.

* Calov, R., A. Ganopolski, V. Petoukhov, M. Claussen and R. Greve. 2002. Large-scale instabilities of the Laurentide ice sheet simulated in a fully coupled climate-system model. *Geophys. Res. Lett.*, **29**(24), 2216. (10.1029/2002GL016078.)

* Clarke, G.K.C., S.J. Marshall, C. Hillaire-Marcel, G. Bilodeau and C. Veiga-Pires. 1999. A glaciological perspective on Heinrich events. In Clark, P.U., R.S. Webb and L.D. Keigwin, eds. *Mechanisms of global climate change at millennial time scales*. Washington, DC, American Geophysical Union, 243–262.

* Fowler, A.C. and E. Schiavi. 1998. A theory of ice-sheet surges. *J. Glaciol.*, **44**(146), 104–118.

* 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), 901–918.

* Greve, R. and H. Blatter. 2009. *Dynamics of ice sheets and glaciers*. Berlin, etc., Springer.

* Greve, R. and D.R. MacAyeal. 1996. Dynamic/thermodynamic simulations of Laurentide ice-sheet instability. *Ann. Glaciol.*, **23**, 328–335.

* Greve, R., M. Weis and K. Hutter. 1998. Palaeoclimatic evolution and present conditions of the Greenland ice sheet in the vicinity of Summit: an approach by large-scale modelling. *Palaeoclimates*, **2**(2–3), 133–161.

* Greve, R., R. Takahama and R. Calov. 2006. Simulation of large-scale ice-sheet surges: the ISMIP HEINO experiments. *Polar Meteorol. Glaciol.*, **20**, 1–15.
Heinrich, H. 1988. Origin and consequences of cyclic ice rafting in the northeast Atlantic Ocean during the past 130 000 years. *Quat. Res.*, **29**(2), 142–152.

Hindmarsh, R.C.A. 2004. A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling. *J. Geophys. Res.*, **109**(F1), F01012. (10.1029/2003JF000065.)

Hindmarsh, R.C.A. and E. Le Meur. 2001. Dynamical processes involved in the retreat of marine ice sheets. *J. Glaciol.*, **47**(157), 271–282.

Hindmarsh, R.C.A. and A.J. Payne. 1996. Time-step limits for stable solutions of the ice-sheet equation. *Ann. Glaciol.*, **23**, 74–85.

Hulbe, C.L. 1997. An ice shelf mechanism for Heinrich layer production. *Paleoceanography*, **12**(5), 711–717.

Hulbe, C.L., D.R. MacAyeal, G.H. Denton, J. Kleman and T.V. Lowell. 2004. Catastrophic ice shelf breakup as the source of Heinrich event icebergs. *Paleoceanography*, **19**(1), PA1004. (10.1029/2003PA000890.)

Hutter, K. 1983. *Theoretical glaciology: material science of ice and the mechanics of glaciers and ice sheets*. Dordrecht, etc., D. Reidel Publishing Co./Tokyo, Terra Scientific Publishing Co.

Huybrechts, P. 1990. A 3-D model for the Antarctic ice sheet: a sensitivity study on the glacial–interglacial contrast. *Climate Dyn.*, **5**(2), 79–92.

Huybrechts, P. 2002. Sea-level changes at the LGM from ice-dy-
namic reconstructions of the Greenland and Antarctic ice sheets
during the glacial cycles. *Quat. Sci. Rev.*, **21**(1–3), 203–231.

Huybrechts, P. and J. de Wolde. 1999. The dynamic response of the Greenland and Antarctic ice sheets to multiple-century climatic warming. *J. Climate*, **12**(8), 2169–2188.

Huybrechts, P., T. Payne and the EISMINT Intercomparison Group. 1996. The EISMINT benchmarks for testing ice-sheet models. *Ann. Glaciol.*, **23**, 1–12.

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), 4071–4087.

MacAyeal, D.R. 1993. Binge/purge oscillations of the Laurentide ice sheet as a cause of the North Atlantic's Heinrich events. *Paleoceanography*, **8**(6), 775–784.

Marshall, S.J. and G.K.C. Clarke. 1997a. A continuum mixture model of ice stream thermomechanics in the Laurentide ice sheet. 1. Theory. *J. Geophys. Res.*, **102**(B9), 20,599–20,614.

Marshall, S.J. and G.K.C. Clarke. 1997b. A continuum mixture model of ice stream thermomechanics in the Laurentide ice sheet. 2. Application to the Hudson Strait ice stream. *J. Geophys. Res.*, **102** (B9), 20,615–20,637.

McManus, J.F., D.W. Oppo and J.L. Cullen. 1999. A 0.5-million-year record of millennial-scale climate variability in the North Atlantic. *Science*, **283**(5404), 971–975.

Meier, M.F. and A. Post. 1987. Fast tidewater glaciers. *J. Geophys. Res.*, **92**(B9), 9051–9058.

Morland, L.W. 1984. Thermomechanical balances of ice sheet flows. *Geophys. Astrophys. Fluid Dyn.*, **29**(1–4), 237–266.

Papa, B.D., L.A. Mysak and Z. Wang. 2005. Intermittent ice sheet discharge events in northeastern North America during the last glacial period. *Climate Dyn.*, **26**(2–3), 201–216.

Pattyn, F. 2003. A new three-dimensional higher-order thermomech-
anical ice-sheet model: basic sensitivity, ice stream development,
and ice flow across subglacial lakes. *J. Geophys. Res.*, **108**(B8),
2382. (10.1029/2002JB002329.)

Pattyn, F. and 20 others. 2008. Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM). *Cryosphere*, **2**(1), 95–108.

Payne, A.J. 1995. Limit cycles in the basal thermal regime of ice sheets. *J. Geophys. Res.*, **100**(B3), 4249–4263.

Payne, A.J. and 10 others. 2000. Results from the EISMINT model intercomparison: the effects of thermomechanical coupling. *J. Glaciol.*, **46**(153), 227–238.

Pollard, D. and R.M. DeConto. 2007. A coupled ice-sheet/ice-shelf/sediment model applied to a marine margin flow-line: forced and unforced variations. *In* Hambrey, M.J., P. Christoffersen, N.F. Glasser and B. Hubbard, * eds*. *Glacial sedimentary processes and products*. Malden, MA, Blackwell, 37–52.

Ritz, C., V. Rommelaere and C. Dumas. 2001. Modeling the evolution of Antarctic ice sheet over the last 420,000 years: implications for altitude changes in the Vostok region. *J. Geophys. Res.*, **106**(D23), 31,943–31,964.

Rutt, I.C., M. Hagdorn, N.R.J. Hulton and A.J. Payne. 2009. The Glimmer community ice sheet model. *J. Geophys. Res.*, **114**(F2), F02004. (10.1029/2008JF001015.)

Saito, F. and A. Abe-Ouchi. 2005. Sensitivity of Greenland ice sheet simulation to the numerical procedure employed for ice-sheet dynamics. *Ann. Glaciol.*, **42**, 331–336.

Seddik, H., R. Greve, T. Zwinger and L. Placidi. 2009. A full-Stokes ice flow model for the vicinity of Dome Fuji, Antarctica, with induced anisotropy and fabric evolution. *Cryos. Discuss.*, **3**(1), 1–31.

Tarasov, L. and W.R. Peltier. 1997. A high-resolution model of the 100 ka ice-age cycle. *Ann. Glaciol.*, **25**, 58–65.

Tarasov, L. and W.R. Peltier. 1999. The impact of thermo-mechanical ice sheet coupling on a model of the 100 kyr ice-age cycle. *J. Geophys. Res.*, **104**(D8), 9517–9545.

Tarasov, L. and W.R. Peltier. 2002. Greenland glacial history and local geodynamic consequences. *Geophys. J. Int.*, **150**(1), 198–229.

Weis, M., R. Greve and K. Hutter. 1999. Theory of shallow ice shelves. *Contin. Mech. Thermodyn.*, **11**(1), 15–50.

ZWinger, T., R. Greve, O. Gagliardini, T. Shiraiwa and M. Lyly. 2007. A full Stokes-flow thermo-mechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka. *Ann. Glaciol.*, **45**, 29–37.

MS received 30 November 2009 and accepted in revised form 24 April 2010