# Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM)

F. Pattyn, L. Perichon, A. Aschwanden, B. Breuer, B. de Smedt, Olivier

Gagliardini, G. H. Gudmundsson, R.C.A. Hindmarsh, A. Hubbard, J. V.

Johnson, et al.

► To cite this version:

F. Pattyn, L. Perichon, A. Aschwanden, B. Breuer, B. de Smedt, et al., Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM). The Cryosphere, 2008, 259, pp.95 à 108. insu-00379774

HAL Id: insu-00379774

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

Submitted on 16 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM)*

F. Pattyn¹, L. Perichon¹, A. Aschwanden², B. Breuer³, B. de Smedt⁴, O. Gagliardini⁵, G. H. Gudmundsson⁶, R. C. A. Hindmarsh⁶, A. Hubbard⁷, J. V. Johnson⁸, T. Kleiner³, Y. Konovalov⁹, C. Martin⁶, A. J. Payne¹⁰, D. Pollard¹¹, S. Price¹⁰, M. Rückamp³, F. Saito¹², O. Souček¹³, S. Sugiyama¹⁴, and T. Zwinger¹⁵

¹Laboratoire de Glaciologie, Université Libre de Bruxelles, CP160/03, Av. F. Roosevelt 50, 1050 Brussels, Belgium

²Institute for Atmospheric and Climate Science, ETH Zurich, Universitaetstrasse 16, 8092 Zurich, Switzerland

³Institute for Geophysics, University of Muenster, Corrensstrasse 24, 48149 Muenster, Germany

⁴Vakgroep Geografie, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium

⁵Laboratoire de Glaciologie et de Géophysique de l'Environnement (LGGE), CNRS, UJF-Grenoble I, BP 96, 38402 Saint Martin d'Hères Cedex, France

⁶Physical Science Division, British Antarctic Survey, Natural Environment Research Council, High Cross, Madingley Road, Cambridge CB3 0ET, UK

⁷Centre for Glaciology, Institute of Geography and Earth Sciences, Aberystwyth University, Ceredigion, SY23 3DP, UK

⁸Department of Computer Science, Social Science Building Room 417, Univ. of Montana, Missoula MT, 59812-5256, USA

⁹Moscow Engineering Physics Institute, Moscow, Russia

¹⁰Bristol Glaciology Centre, School of Geographical Sciences, University Road, University of Bristol, Bristol BS8 1SS, UK

¹¹Earth and Environmental Systems Institute, College of Earth and Mineral Sciences, 2217 Earth-Engineering Sciences Bldg., Pennsylvania State University, University Park, PA 16802, USA

¹²Frontier Res. Center for Global Change, 3173-25 Showamachi, Kanazawa-ku, Yokohama City, Kanagawa 236-0001, Japan

¹³Department of Geophysics, Charles University Prague, V. Holešovičkách 2, 18000 Praha 8, Czech Republic

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

¹⁵CSC-Scientific Computing Ltd., Keilaranta 14, P.O. Box 405, 02101 Espoo, Finland

*Ice Sheet Model Intercomparison Project for Higher-Order Models; http://homepages.ulb.ac.be/~fpattyn/ismip

Received: 20 December 2007 – Published in The Cryosphere Discuss.: 19 February 2008

Revised: 7 July 2008 – Accepted: 8 August 2008 – Published: 26 August 2008

**Abstract.** We present the results of the first ice sheet model intercomparison project for higher-order and full-Stokes ice sheet models. These models are compared and verified in a series of six experiments of which one has an analytical solution obtained from a perturbation analysis. The experiments are applied to both 2-D and 3-D geometries; five experiments are steady-state diagnostic, and one has a time-dependent prognostic solution. All participating models give results that are in close agreement. A clear distinction can be made between higher-order models and those that solve the full system of equations. The full-Stokes models show a much smaller spread, hence are in better agreement with one another and with the analytical solution.

## 1 Introduction

According to the recent IPCC report (IPCC, 2007), dynamical ice-flow processes not included in current models but suggested by recent observations could increase the vulnerability of the ice sheets to warming, increasing future sea level rise. Understanding of these processes is limited, and there is no consensus on their magnitude. It was also stressed that a net ice mass loss could occur if dynamical ice discharge dominates the ice sheet mass balance (IPCC, 2007). Although the viscous flow of ice is rather well understood on a theoretical level, confidence in models is low because processes at the ice base and the seaward margin are poorly understood and have not been represented in models. Several stress components come into play in regions of high variability in basal topography and/or basal slipperiness.

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

Correspondence to: F. Pattyn
(fpattyn@ulb.ac.be)
Despite the lack of comprehensive predictive ice sheet modeling, the ice sheet modeling community has evolved considerably over the last decade. Increasing computational power has led to the development of more complex ice sheet models, with varying degrees of approximations to the Stokes equations. However, progress has been hampered by the lack of a universal verification framework and in particular by a lack of full-Stokes analytical solutions. While methods exist for constructing exact solutions to the Newtonian full-Stokes equations in the flowline case (e.g. Ladyzhenskaya, 1969), obtaining analytical solutions for the 3-D case is less straightforward. Nevertheless, solutions based on perturbation analysis exist (e.g. Gudmundsson, 2003), and were used for this study. Benchmark validation exercises were carried out on large-scale ice-sheet and ice-shelf models in the 1990s (Huybrechts et al., 1996; MacAyeal et al., 1996; Payne et al., 2000), but these tests were largely restricted to zero-order (SIA) models and solutions. Here we present an intercomparison exercise which involves 28 higher-order flow models of varying complexity. The experiments described in this paper are designed to evaluate the conditions under which different higher-order solutions are viable and to determine whether numerical issues affect the result.

During the first and second EISMINT<fnref n="1" /> model intercomparison exercises, a number of benchmarks were proposed specifically for ice sheet models (Huybrechts et al., 1996, 1998; Payne et al., 2000) and ice shelf models (MacAyeal et al., 1996). These ice sheet models were based on the zeroth-order shallow-ice approximation (SIA; Hutter, 1983), incorporating only vertical shear stresses in the force balance. The ISMIP-HOM exercise focuses on so-called higher-order models, i.e. models that incorporate further mechanical effects, principally longitudinal stress gradients, as well as those that solve the full system of equations of the Stokes problem.

The six experiments which comprise this benchmark exercise are designed to be universally accessible to many different types of models, i.e. flowline (2-D), vertically integrated planform (2.5-D) and full 3-D models. Furthermore, the experiments are defined as well-posed continuum problems so that their application is not limited to any specific numerical methodology. The equations have been solved using well established finite-difference and finite-element methods, in addition to more esoteric spectral techniques. The latter hold particular promise for providing high-quality results in the absence of analytical solutions, but since only one participant provided spectral model results, they were not isolated for comparison.

With the exception of Exp. F, all experiments are diagnostic; i.e. time evolution is not considered. This means that for a given ice geometry, a Glen-type flow law, and appropriate boundary conditions, the stress and velocity fields

can be calculated. Exp. F considers the time-dependent response (the experiment is run until the free surface and velocity field reach a steady state) for a constant viscosity (linear flow law). Constant viscosity is assumed because in this case there exist analytic solutions derived from a first-order perturbation analysis of flow down an inclined plane (Gudmundsson, 2003). In all experiments, thermomechanical effects are neglected and ice is considered to be isothermal and isotropic.

## 2 General model setup

### 2.1 Model physics, parameters and constants

Higher-order models are ice-sheet or glacier models that incorporate effects not present in the shallow-ice approximation. In most cases this implies the inclusion of longitudinal stress gradients apart from the two horizontal-plane shear components (Hindmarsh, 2004). Longitudinal stresses have recently been termed “membrane stresses” when considered in three dimensions (Hindmarsh, 2006). The suite of models is based on conservation laws of mass and momentum, i.e.

$$ \nabla \cdot \mathbf{v} = 0, \qquad (1) $$

$$ \rho \frac{dv}{dt} = \nabla \cdot \mathbf{T} + \rho \mathbf{g}, \qquad (2) $$

where $\rho$ is the ice density, $\mathbf{g}$ is gravitational acceleration, $\mathbf{v}$ is the velocity vector, and $\mathbf{T}$ is the stress tensor. Values for parameters and constants are given in Table 1. Generally, acceleration terms in Eq. (2) are neglected. Ice incompressibility is more easily described if the stress tensor is split into a deviatoric part and an isotropic pressure $P$,

$$ \mathbf{T} = \mathbf{T}' - P\mathbf{I}, \qquad (3) $$

where $P = -\frac{1}{3}\text{tr}(\mathbf{T})$. The constitutive equation for ice then links deviatoric stresses to strain rates:

$$ \mathbf{T}' = 2\eta \dot{\mathbf{e}}', \qquad (4) $$

where $\mathbf{T}'$ and $\dot{\mathbf{e}}$ are the deviatoric stress and strain-rate tensor, respectively, and $\eta$ is the effective viscosity. Both linear and nonlinear ice rheologies are considered. In the latter case (Glen's flow law), $\eta$ is strain-rate dependent and defined by

$$ \eta = \frac{1}{2} A^{-1/n} \dot{\varepsilon}_e^{(1-n)/n}, \qquad (5) $$

where $\dot{\varepsilon}_e$ is the effective strain rate. For the case of linear rheology, Eq. (5) reduces to $\eta=(2A)^{-1}$, where $A$ is spatially uniform (i.e., the ice is isothermal). Neglecting acceleration terms, the momentum balance is written as:

$$ \mathrm{div} \mathbf{T} + \rho \mathbf{g} = \mathrm{div} \mathbf{T}' - \mathrm{grad} P + \rho \mathbf{g} = 0, \qquad (6) $$

<footnote n="1">¹EISMINT: European Ice Sheet Model INTercomparison; http://homepages.vub.ac.be/~phuybrec/eismint.html</footnote>
<table><caption>**Table 1.** Constants for the numerical model.</caption><thead><tr><th>Constant</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>A</td><td>Ice-flow parameter</td><td>10<sup>-16</sup></td><td>Pa<sup>-n</sup> a<sup>-1</sup></td></tr><tr><td>ρ</td><td>Ice density</td><td>910</td><td>kg m<sup>-3</sup></td></tr><tr><td>g</td><td>Gravitational constant</td><td>9.81</td><td>m s<sup>-2</sup></td></tr><tr><td>n</td><td>Exponent in Glen's flow law</td><td>3</td><td></td></tr><tr><td></td><td>Seconds per year</td><td>31 556 926</td><td>s a<sup>-1</sup></td></tr></tbody></table>

Since acceleration due to gravity acts only in the vertical, this leads to

$$ \frac{\partial \tau'_{xx}}{\partial x} - \frac{\partial P}{\partial x} + \frac{\partial \tau'_{xy}}{\partial y} + \frac{\partial \tau'_{xz}}{\partial z} = 0, \qquad (7) $$

$$ \frac{\partial \tau'_{yx}}{\partial x} + \frac{\partial \tau'_{yy}}{\partial y} - \frac{\partial P}{\partial y} + \frac{\partial \tau'_{yz}}{\partial z} = 0, \qquad (8) $$

$$ \frac{\partial \tau'_{zx}}{\partial x} + \frac{\partial \tau'_{zy}}{\partial y} + \frac{\partial \tau'_{zz}}{\partial z} - \frac{\partial P}{\partial z} = \rho g. \qquad (9) $$

Solving Eqs. (7)–(9) leads to the full-Stokes solution. Simplifications of these equations lead to the various higher-order approximations discussed below.

## 2.2 Boundary conditions

In Exps. A, B, E1 and F1 the ice is frozen to the bed ($v_b=0$). For the other experiments, basal sliding is introduced through a friction law, characterized by a friction coefficient $\beta^2$. This friction law has the form of

$$ \beta^2 \mathbf{t} \cdot \mathbf{v} = \mathbf{t} \cdot (\mathbf{T}\mathbf{n}_b) = \tau_b, \qquad (10) $$

where $n_b$ is the unit normal vector pointing into the bedrock, $\mathbf{t}$ is the unit tangent vector, and $\beta^2$ (Pa a m$^{-1}$) is a positive scalar (MacAyeal, 1993). Basal shear stress $\tau_b$ is not equal to the driving stress but is part of the solution. The stress is negligible at the upper ice surface, implying that $\mathbf{n}_s \cdot (\mathbf{T}\mathbf{n}_b) = P_{\text{atm}} \approx 0$.

Kinematic boundary conditions apply at the upper and lower surfaces of the ice mass, i.e.

$$ \frac{\partial z_i}{\partial t} + v_x(z_i) \frac{\partial z_i}{\partial x} + v_y(z_i) \frac{\partial z_i}{\partial y} - v_z(z_i) = 0, \qquad (11) $$

for $i=(s, b)$. Since the vertical velocity field must obey the incompressibility condition Eq. (1), and the surface accumulation/ablation is zero, the vertical velocity at the surface contains the local imbalance as well and becomes a model output.

## 2.3 Model domain

For the 3-D experiments the model domain is a square of side L, and for 2-D experiments the domain is a flowline of length L in the x-z plane. The minimum number of grid points is

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Basal topography $z_b$ (m) for Exp. A according to Eq. (13) for $L=80$ km. Ice flow is from left to right.</figcaption></figure>

not predefined, and any type of discretization scheme can be used. The number of grid points in the horizontal and vertical directions can be chosen freely. The basic parameter for the experiments is the length scale L of the domain. Exps. A–D are carried out for $L=160, 80, 40, 20, 10$ and $5$ km, respectively, which results in aspect ratios $\epsilon=H/L$ varying from 0.006 to 0.2. Periodic boundary conditions are applied in the horizontal, so that the domain is surrounded by an infinite number of copies of itself.

# 3 Experiment description

## 3.1 Exp. A: Ice flow over a bumpy bed

Exp. A considers a parallel-sided slab of ice with a mean ice thickness $H=1000$ m lying on a sloping bed with a mean slope $\alpha=0.5^\circ$. This slope is maximum in x and zero in y. The basal topography is then defined as a series of sinusoidal oscillations with an amplitude of 500 m. The surface elevation is defined as

$$ z_s(x, y) = -x \cdot \tan \alpha. \qquad (12) $$

The basal topography is then given by

$$ z_b(x, y) = z_s(x, y) - 1000 + 500 \sin(\omega x) \cdot \sin(\omega y), \quad (13) $$

where $x \in [0, L]$ and $L=160, 80, 40, 20, 10$ and $5$ km, respectively. The basal bumps have a frequency of $\omega=2\pi/L$. The bed topography is shown in Fig. 1.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Basal friction coefficient $\beta^2$ for Exp. C.</figcaption></figure>

## 3.2 Exp. B: Ice flow over a rippled bed

The only difference with Exp. A is that the basal topography does not vary with $y$, so that the experiment is suitable for 2-D flowline models as well. The basal topography is thus formed by a series of ripples with an amplitude of 500 m:

$$z_s(x, y) = -x \cdot \tan \alpha \quad (14)$$

$$z_b(x, y) = z_s(x, y) - 1000 + 500 \sin(\omega x). \quad (15)$$

## 3.3 Exp. C: Ice stream flow I

The experiment setup is similar to Exp. A, except that the bedrock topography is flat, so that the ice thickness is spatially uniform ($H=1000$ m):

$$z_s(x, y) = -x \cdot \tan \alpha \quad (16)$$

$$z_b(x, y) = z_s(x, y) - 1000, \quad (17)$$

where $x \in [0, L]$ and $L=160, 80, 40, 20, 10$ and $5$ km, respectively, and where $\alpha=0.1^\circ$. The basal friction coefficient is prescribed as

$$\beta^2(x, y) = 1000 + 1000 \sin(\omega x) \cdot \sin(\omega y). \quad (18)$$

The $\beta^2$-field is shown in Fig. 2. The basal friction oscillations have a frequency $\omega=2\pi/L$.

## 3.4 Exp. D: Ice stream flow II

The only difference with Exp. C is that the basal friction coefficient does not vary with $y$, so that the experiment is suited

<figure><img src="image_6.png" /><figcaption>Fig. 3. Surface and bedrock profile of Haut Glacier d’Arolla.</figcaption></figure>

for 2-D flowline models as well. The basal friction field is thus formed by a series of ripples defined as

$$\beta^2(x, y) = 1000 + 1000 \sin(\omega x). \quad (19)$$

Note that in Exps. C and D the basal friction coefficient $\beta^2$ goes to zero within the domain.

## 3.5 Exp. E: Haut Glacier d’Arolla

Exp. E is a diagnostic experiment along the central flowline of a temperate glacier in the European Alps (Haut Glacier d’Arolla), based on earlier experiments by Blatter et al. (1998) and Pattyn (2002). Model input consists of the longitudinal surface and bedrock profiles of Haut Glacier d’Arolla, Switzerland, according to the Little Ice Age geometry (Fig. 3). The longitudinal profile of this glacier has a very simple geometry, hence the resulting stress field is not influenced by geometrical perturbations such as the presence of a steep ice fall. In a first experiment (E1), a zero basal velocity is considered ($\beta^2=\infty$), and the width of the drainage basin is kept equal to 1 along the entire flowline. The flow-law rate factor $A=10^{-16} Pa^{-n} a^{-1}$ is assumed constant. Upstream and downstream boundary conditions imply zero ice thickness and velocity. The input file has a resolution $\Delta x=100$ m, but the authors were free to choose any grid resolution.

A second experiment (E2) considers a narrow zone of zero traction, similar to the experiment described in Blatter et al. (1998):

$$\beta^2 = 0 \quad \text{for } 2200 \le x \le 2500 \text{ m}$$

$$\mathbf{v}_b = 0 \quad \text{otherwise}$$
## 3.6 Exp. F: Prognostic experiment

Exp. F is a prognostic experiment for which the free surface is allowed to relax until a steady state is reached for zero surface mass balance:

$$ \lim_{t \to \infty} \frac{\partial H}{\partial t} = \lim_{t \to \infty} \left[ -\nabla_h \cdot \int_{z_b}^{z_s} \mathbf{v}_h \, dz \right] = 0, \quad (20) $$

where $\mathbf{v}_h$ is the horizontal velocity vector (m a⁻¹). Basic model setup differs from the setup in Exps. A and C. A slab of ice with mean ice thickness $H^{(0)}=1000$ m is considered, resting on a sloping bed with a mean slope $\alpha=3.0°$ (Fig. 4). This slope is maximum in $x$ and zero in $y$. The bedrock plane is parallel to the surface plane and is perturbed by a Gaussian bump. Initial bedrock $B^{(0)}$ and unperturbed surface $S^{(0)}$ elevation are thus governed by

$$ S^{(0)}(x, y) = 0 \quad (21) $$

$$ B^{(0)}(x, y) = -H^{(0)} + a_0 \left( \exp \left[ \frac{-(x^2 + y^2)}{\sigma^2} \right] \right), \quad (22) $$

where $\sigma=10000$ m = $10H^{(0)}$ and where $x, y$ (m) are the horizontal coordinates with respect to the center of the Gaussian bump. The basal perturbation has a maximum height of one-tenth of the mean ice thickness, i.e. $a_0=100=0.1H^{(0)}$ (Fig. 4). The domain size $L$ is taken to be $100 H^{(0)}$ in $x$ and $y$. The horizontal coordinates for output are scaled by

$$ \hat{x} = \frac{x}{H^{(0)}} \qquad \hat{y} = \frac{y}{H^{(0)}}. \quad (23) $$

Periodic boundary conditions are applied in the horizontal. The major difference with the previous experiments is that $n=1$ in Eq. (5), so that the effective viscosity is constant and becomes $\eta=(2A)^{-1}$. Therefore, the unperturbed velocity field at the surface is defined by

$$ U^{(0)} = AH^{(0)} \tau_b^{(0)} = \rho g A \left[ H^{(0)} \right]^2 \sin \alpha, \quad (24) $$

where $\tau_b^{(0)}=\rho g H^{(0)} \sin \alpha$ is the unperturbed basal shear stress, and $A=2.140373 \times 10^{-7} \text{ Pa}^{-1} \text{ a}^{-1}$, so that $U^{(0)}=100 \text{ m a}^{-1}$.

Experiments are carried out for different values of the slip ratio $c$, which determines the relation between the basal velocity and basal drag. The basal velocity is written in terms of a basal friction coefficient $\beta^2$, or

$$ U_b = \frac{\tau_b}{\beta^2}. \quad (25) $$

Following the scalings given by Gudmundsson (2003), the basal friction coefficient is related to the slip ratio $c$ by

$$ \beta^2 = (cAH^{(0)})^{-1}. \quad (26) $$

<figure><img src="image_7.png" /><figcaption>Fig. 4. Tilted coordinate system used for Exp. F.</figcaption></figure>

<table><caption>Table 2. Constants for the model setup according to Exp. F.</caption><thead><tr><th>Constant</th><th>Value</th></tr></thead><tbody><tr><td>A</td><td>Ice-flow parameter</td></tr><tr><td>n</td><td>Flow law exponent</td></tr><tr><td>α</td><td>Mean surface slope</td></tr><tr><td>a<sub>0</sub></td><td>Amplitude Gaussian bump</td></tr><tr><td>σ</td><td>Width Gaussian bump</td></tr></tbody></table>2.140373 × 10<sup>-7</sup> Pa<sup>-1</sup> a<sup>-1</sup><br>1<br>3°<br>100 m<br>10 000 m

Experiments are run for slip ratios $c=0$ and 1 (F1 and F2, respectively). It is easily demonstrated that $U_b^{(0)}=cU^{(0)}$. Table 2 lists the main constants used for Exp. F. Using these settings, the model should run until a steady state of the free surface is reached.

# 4 Model classification

A total of 27 numerical models and one analytical model from 20 contributors participated in the intercomparison exercise. Table 3 summarizes model characteristics and contributions. The different Stokes approximants extend the shallow-ice approximation in various ways (Hutter, 1983; Hindmarsh, 2004). We will follow here the classification scheme for higher-order models by Hindmarsh (2004), who gives a detailed description of the different “longitudinal stress schemes” widely used in ice sheet modeling. The most common longitudinal stress approximations introduce the two horizontal velocity components as field variables. This leads to an elliptic system with two rather than four variables of the full system at points in three-dimensional space (Pattyn, 2003; Hindmarsh, 2004), and the resulting linear systems are generally better conditioned than those resulting from the numerical analysis of the full system (Hindmarsh, 2004). These models are termed “multilayer models”. A number of these models solve an elliptic system at one
<table><caption>**Table 3.** List with the 28 participating models. Model: model acronym based on the initials of each author; Type: the model type (see text for description); Dims: model dimensions; Method: numerical method (FE = finite elements, FD = finite differences, Sp = spectral method, FV = finite volume, An = analytical); A–F participation in the different experiments is marked with an x.</caption><thead><tr><th>Model</th><th>Type</th><th>Dim</th><th>Method</th><th>A</th><th>B</th><th>C</th><th>D</th><th>E</th><th>F</th><th>Reference</th></tr></thead><tbody><tr><td>aas1</td><td>full-Stokes</td><td>2-D</td><td>FE</td><td></td><td>x</td><td></td><td>x</td><td>x</td><td></td><td>unpublished</td></tr><tr><td>aas2</td><td>full-Stokes</td><td>3-D</td><td>FE</td><td>x</td><td>x</td><td>x<sup>a</sup></td><td>x</td><td></td><td></td><td>unpublished</td></tr><tr><td>ahu1</td><td>LMLa</td><td>3-D</td><td>FD</td><td>x<sup>a</sup></td><td>x<sup>a</sup></td><td>x<sup>a</sup></td><td>x<sup>a</sup></td><td>x</td><td></td><td>Hubbard et al. (1998); Nienow et al. (2005)</td></tr><tr><td>ahu2</td><td>LMLa</td><td>2-D</td><td>FD</td><td></td><td>x</td><td></td><td>x<sup>c</sup></td><td></td><td></td><td>Hubbard et al. (1998); Nienow et al. (2005)</td></tr><tr><td>bds1</td><td>LMLa</td><td>2-D</td><td>FE</td><td></td><td>x</td><td></td><td>x</td><td>x</td><td></td><td>unpublished</td></tr><tr><td>cma1</td><td>full-Stokes</td><td>3-D</td><td>FE</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>Martin et al. (2003)</td></tr><tr><td>cma2</td><td>LMLa</td><td>3-D</td><td>FE</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td>x</td><td>unpublished</td></tr><tr><td>dp01</td><td>L1L2</td><td>2-D</td><td>FD</td><td></td><td></td><td>x</td><td>x</td><td></td><td></td><td>Pollard and DeConto (2007)</td></tr><tr><td>fpa1</td><td>LMLa</td><td>3-D</td><td>FD</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>Pattyn (2003)</td></tr><tr><td>fpa2</td><td>full-Stokes</td><td>3-D</td><td>FD</td><td>x<sup>c</sup></td><td></td><td>x<sup>b</sup></td><td></td><td></td><td></td><td>Pattyn (2008)</td></tr><tr><td>fsa1</td><td>LMLa</td><td>3-D</td><td>FD</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td>x</td><td>based on Pattyn (2003); Colinge and Rappaz (1999)</td></tr><tr><td>ghg1</td><td>full-Stokes</td><td>3-D</td><td>An</td><td></td><td></td><td></td><td></td><td></td><td>x<sup>d</sup></td><td>Gudmundsson (2003)</td></tr><tr><td>jjv1</td><td>full-Stokes</td><td>3-D</td><td>FE</td><td>x<sup>b</sup></td><td>x</td><td></td><td></td><td></td><td></td><td>Johnson and Staiger (2007)</td></tr><tr><td>lpe1</td><td>L1L1</td><td>2-D</td><td>FD</td><td></td><td></td><td></td><td>x</td><td></td><td></td><td>based on MacAyeal (1989); Pattyn (2003)</td></tr><tr><td>mbr1</td><td>LMLa</td><td>3-D</td><td>FD</td><td>x<sup>b</sup></td><td>x<sup>b</sup></td><td>x<sup>c</sup></td><td>x<sup>c</sup></td><td>x</td><td>x</td><td>Breuer et al. (2006)</td></tr><tr><td>mmr1</td><td>full-Stokes</td><td>3-D</td><td>FE</td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>unpublished</td></tr><tr><td>mtk1</td><td>LTSML</td><td>3-D</td><td>FD</td><td>x</td><td>x</td><td>x<sup>a</sup></td><td></td><td></td><td>x</td><td>based on Pattyn (2003); Hindmarsh (2004)</td></tr><tr><td>oga1</td><td>full-Stokes</td><td>3-D</td><td>FE</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>Zwinger et al. (2007); Gagliardini and Zwinger (2008)</td></tr><tr><td>oso1</td><td>SIA'</td><td>3-D</td><td></td><td>x<sup>a</sup></td><td>x<sup>a</sup></td><td></td><td></td><td></td><td></td><td>unpublished</td></tr><tr><td>rhi1</td><td>full-Stokes</td><td>3-D</td><td>Sp</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td>Hindmarsh (2004)</td></tr><tr><td>rhi2</td><td>LMLa</td><td>3-D</td><td>Sp</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td>Hindmarsh (2004)</td></tr><tr><td>rhi3</td><td>full-Stokes</td><td>3-D</td><td>Sp</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td>Hindmarsh (2004)</td></tr><tr><td>rhi4</td><td>L1L2</td><td>3-D</td><td>Sp</td><td>x<sup>a</sup></td><td>x<sup>a</sup></td><td>x</td><td>x</td><td></td><td></td><td>Hindmarsh (2004)</td></tr><tr><td>rhi5</td><td>L1L1</td><td>3-D</td><td>Sp</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td>Hindmarsh (2004)</td></tr><tr><td>spr1</td><td>full-Stokes</td><td>2-D</td><td>FV</td><td></td><td>x<sup>a</sup></td><td></td><td></td><td>x<sup>d</sup></td><td></td><td>Price et al. (2007)</td></tr><tr><td>ssu1</td><td>full-Stokes</td><td>2-D</td><td>FE</td><td></td><td>x</td><td></td><td>x</td><td>x</td><td></td><td>Sugiyama et al. (2003)</td></tr><tr><td>tpa1</td><td>LMLa</td><td>3-D</td><td>FD</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>based on Pattyn (2003); Hindmarsh (2004)</td></tr><tr><td>yko1</td><td>full-Stokes</td><td>3-D</td><td>FD</td><td>x</td><td></td><td></td><td></td><td>x<sup>d</sup></td><td></td><td>unpublished</td></tr></tbody></table>

elevation only (generally the upper surface), giving a computationally two-dimensional problem. These higher-order approximations are labeled L1L1, L1L2, LMLa and LTSML (Hindmarsh, 2004).

The L1L1 approximation is a one-layer longitudinal stress scheme using $\tau'_{xx}$ at the surface computed by solving elliptic equations and is identical to the approximation used by MacAyeal (1989). An alternative approximation is the L1L2 approach, or one-layer longitudinal stress scheme, using $\dot{\varepsilon}_{xx}$ at the surface computed by solving elliptical equations with a vertical correction of $\tau'_{xx}$. Here, the surface velocities used in computing the non-horizontal plane stresses are computed using the vertical shear stresses in the shear strain relationship and in the sliding relationship.

The most common approximation is the LMLa or multi-layer longitudinal stress scheme. This is the classic longitu-

dinal stress scheme used by Blatter (1995) and Pattyn (2003). Compared with L1L2, the longitudinal stresses use the velocity at the corresponding elevations rather than at the surface, and the stress-invariant calculations are self-consistent rather than using the SIA stress (Hindmarsh, 2004). Finally, there is the LTSML or multilayer longitudinal stresses scheme with horizontal shear stress gradients approximated by SIA. Here, horizontal gradients of the vertical velocity are neglected. Horizontal plane shear stresses, when needed to calculate the horizontal gradient of such shear stresses, are approximated by SIA values. This approach is similar to LMLa, but with inclusion of the vertical resistive stress, as in Van der Veen and Whillans (1989).

<footnote><sup>a</sup> not for L=5 km</footnote>

<footnote><sup>b</sup> not for L=5 and 10 km</footnote>

<footnote><sup>c</sup> not for L=5, 10 and 20 km</footnote>

<footnote><sup>d</sup> only no-sliding case</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Results for Exp. A: norm of the surface velocity across the bump at $y=L/4$ for different length scales $L$. The mean value and standard deviation are shown for both types of models.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Results for Exp. B: norm of the surface velocity for different length scales $L$. The mean value and standard deviation are shown for both types of models.</figcaption></figure>

# 5 Results

A graphical representation of all the results for each of the contributing authors as well as the submitted data files are found in the supplemental files http://www.the-cryosphere.net/2/95/2008/tc-2-95-2008-supplement.zip. The detailed

description of the experimental setup is given in http://www.the-cryosphere.net/2/95/2008/tc-2-95-2008-supplement.zip. An analysis on the CPU performance of each of the experiments is presented in an accompanying paper (Gagliardini and Zwinger, 2008).
<figure><img src="image_10.png" /><figcaption>Fig. 7. Surface velocity profiles along y=L/4 according to the analytical SIA solution for Exps. A–D, based on Eqs. (27) and (28). Results are independent of L.</figcaption></figure>

The numerical results are displayed in both numerical and graphical format. A distinction is made between full-Stokes (FS) and non-full-Stokes (NFS) models, such as LMLa, LTSML, L1L2 and L1L1. All parameters refer to horizontal surface velocities, which are determined as the norm of the horizontal components of the velocity field, defined by $||v_s||=\sqrt{v_x^2+v_y^2}$. Velocities are displayed along a section y=L/4 for Exps. A and C and along a central line y=L/2 for Exp. F. For the 2-D experiments we show results along the flowline. The graphs show the mean velocity and its standard deviation along each section or flowline for both FS and NFS models. The tables present two parameters for each flowline or section: the maximum value of the norm of the surface velocity and the mean value of the velocity. For both parameters the mean and standard deviations are given for both FS and NFS models.

## 5.1 Experiments A and B

The results for Exps. A and B are shown in Figs. 5 and 6, respectively, and the numerical results are given in Tables 4 and 5. Each graph displays the norm of the surface velocity across the bumps at y=L/4 (for Exp. A) and along the central flowline (for Exp. B), for the different length scales and model groups (FS and NFS). The experiments were set up such that for this longitudinal profile the SIA gives a solution independent of L, which is not the case for higher-order models. The surface velocity according to the SIA is given by

$$v_x(z_s) = v_x(z_b) + \frac{2A}{n+1} (\rho g \tan \alpha)^n H^{n+1}, \quad (27)$$

where $v_x(z_b)=0$ is the basal velocity (Fig. 7). The maximum surface velocity according to the SIA remains constant for all length scales (119.69 m a⁻¹). However, whenever topographic differences occur, longitudinal stress gradients must

develop which tend to smooth out the velocity field. For high aspect ratios $ϵ=H/L$ (hence low values of L) this leads to a more or less constant surface velocity field as the ice sheet does not “feel” the individual bedrock undulations. Rather, it feels the fast sequence of large bed undulations as a viscous drag. The aspect ratio $ϵ$ determines the amplitude of the horizontal surface velocity field, and the surface velocity decreases from ~100 to ~10 m a⁻¹.

Full-Stokes models closely agree with each other when calculating the velocity field for different length scales, compared to the larger spread of solutions for the higher-order approximations (Tables 4 and 5). Several factors could be responsible for the larger spread among NFS models: (i) More models are participating; (ii) These models are solving different continuum equations (LMLa, LTSML, L1L1, and L1L2); (iii) At the highest aspect ratios, the different approximations are not valid, so that the full-stress model needs to be solved; and (iv) There are numerical errors relative to the unknown exact solutions of the continuum equations. The FS results are only subject to a spread from the latter cause. For the smallest length scales the standard deviation for the full-Stokes models reduces to <0.2 m a⁻¹ (Tables 4 and 5). The coarsest grid used by the participating models had dimensions 40×40, leading to a numerical mismatch of the order of 2.5%, which is far smaller than the standard deviation shown in Fig. 5 for the higher-order models. It is therefore very likely that the large spread associated with the higher-order models is due to the invalidity of the approximations compared to the full-Stokes solution.

The flowline experiments (Exp. B) show similar results compared to the 3-D experiments (Exp. A). There are differences, however, associated with the absence of transverse stress gradients. A peculiarity is observed in the surface horizontal velocity for the smallest length scale L=5 km. The surface velocity according to all full-Stokes models is larger over the bump than over the trough, hence anti-correlated with ice thickness (Fig. 6). All the other approximations (LMLa, LTSML, L1L2 and L1L1) predict that ice velocity and thickness are positively correlated for all length scales. The marked difference can be attributed to mass conservation, as at such high aspect ratios the horizontal ice flux cannot be balanced by the vertical flux at the free surface since the vertical velocity would be too large for the given depth. The horizontal ice flux is therefore more or less constant, inducing larger velocities for smaller depth and vice versa. This behavior is noticeable only for the flowline experiments, because in 3-D the ice can flow around the bumps. The flow inversion is an artifact stemming from the diagnostic nature of the experiments and would disappear if the free surface were allowed to respond to the applied stress field. Higher-order models fail to produce the velocity inversion, since the stress field is determined solely from horizontal strain rate components and vertical velocities are an a posteriori model result (e.g. Pattyn, 2003).
<table><caption>**Table 4.** Mean values ($\mu$), standard deviation ($\sigma$) and number of participating models ($n$) of the maximum horizontal ice velocity at the surface in the direction of the flow. Results are listed for Exp. A–D and for each length scale. Units are $\text{m a}^{-1}$.</caption><thead><tr><th rowspan="2">L (km)</th><th colspan="2">5</th><th colspan="2">10</th><th colspan="2">20</th><th colspan="2">40</th><th colspan="2">80</th><th colspan="2">160</th></tr><tr><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th></tr></thead><tbody><tr><td><strong>Exp. A</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>15.33 (1.76) 7</td><td>26.04 (4.26) 10</td><td>41.03 (3.92) 11</td><td>64.73 (4.83) 11</td><td>88.33 (5.15) 11</td><td>104.77 (4.99) 11</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>14.56 (0.19) 6</td><td>24.36 (0.47) 6</td><td>39.73 (1.30) 7</td><td>63.89 (2.35) 8</td><td>87.10 (4.04) 8</td><td>102.63 (6.51) 8</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>Exp. B</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>10.87 (1.40) 8</td><td>23.51 (4.29) 12</td><td>47.85 (4.14) 13</td><td>74.88 (5.13) 13</td><td>96.43 (5.76) 13</td><td>109.36 (4.52) 13</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>11.76 (0.20) 9</td><td>22.82 (0.72) 10</td><td>46.91 (0.91) 10</td><td>73.77 (0.89) 10</td><td>95.12 (0.89) 10</td><td>108.33 (0.75) 10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>Exp. C</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>12.14 (5.23) 9</td><td>15.39 (1.50) 11</td><td>18.31 (1.02) 11</td><td>28.48 (1.77) 11</td><td>60.99 (6.05) 11</td><td>141.38 (22.40) 12</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>16.00 (0.01) 5</td><td>16.53 (0.31) 6</td><td>19.05 (0.57) 7</td><td>29.32 (1.09) 7</td><td>60.15 (2.25) 7</td><td>143.75 (11.81) 7</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>Exp. D</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>12.86 (4.88) 9</td><td>16.55 (1.08) 10</td><td>21.48 (2.14) 10</td><td>41.30 (4.29) 11</td><td>103.77 (28.82) 12</td><td>244.97 (31.18) 12</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>16.48 (0.16) 8</td><td>17.11 (0.30) 8</td><td>21.33 (0.56) 8</td><td>41.51 (0.94) 8</td><td>97.64 (1.50) 8</td><td>238.44 (1.21) 8</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 5.** Mean values ($\mu$), standard deviation ($\sigma$) and number of participating models ($n$) of the mean horizontal ice velocity at the surface in the direction of the flow. Results are listed for Exp. A–D and for each length scale. Units are $\text{m a}^{-1}$.</caption><thead><tr><th rowspan="2">L (km)</th><th colspan="2">5</th><th colspan="2">10</th><th colspan="2">20</th><th colspan="2">40</th><th colspan="2">80</th><th colspan="2">160</th></tr><tr><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th><th>μ</th><th>(σ) n</th></tr></thead><tbody><tr><td><strong>Exp. A</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>14.61 (1.79) 7</td><td>20.62 (3.23) 10</td><td>24.93 (2.20) 11</td><td>31.99 (2.07) 11</td><td>37.54 (1.58) 11</td><td>40.36 (1.07) 11</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>14.20 (0.18) 6</td><td>20.02 (0.36) 6</td><td>24.74 (0.79) 7</td><td>31.89 (1.07) 8</td><td>37.31 (1.30) 8</td><td>39.98 (1.52) 8</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>Exp. B</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>10.54 (1.36) 8</td><td>18.24 (2.67) 12</td><td>27.80 (2.08) 13</td><td>35.55 (1.86) 13</td><td>39.76 (1.38) 13</td><td>41.38 (0.76) 13</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>11.04 (0.17) 9</td><td>19.09 (0.56) 10</td><td>28.28 (0.60) 10</td><td>35.75 (0.48) 10</td><td>39.76 (0.28) 10</td><td>41.40 (0.24) 10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>Exp. C</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>12.12 (5.22) 9</td><td>15.19 (1.46) 11</td><td>16.36 (0.94) 11</td><td>19.25 (0.95) 11</td><td>27.24 (1.44) 11</td><td>40.83 (3.79) 12</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>15.99 (0.00) 5</td><td>16.24 (0.16) 6</td><td>16.84 (0.24) 7</td><td>19.76 (0.38) 7</td><td>27.38 (0.55) 7</td><td>41.68 (1.83) 7</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>Exp. D</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NFS</td><td>12.85 (4.88) 9</td><td>16.28 (0.97) 10</td><td>18.33 (1.13) 10</td><td>24.20 (1.09) 11</td><td>38.46 (6.48) 12</td><td>58.30 (4.86) 12</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FS</td><td>16.43 (0.16) 8</td><td>16.81 (0.18) 8</td><td>18.40 (0.24) 8</td><td>24.63 (0.29) 8</td><td>37.00 (0.33) 8</td><td>57.17 (0.35) 8</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

## 5.2 Experiments C and D

In this series of experiments, variations in basal conditions (slipperiness) determine where longitudinal stress gradients must develop. Due to the importance of basal sliding, the ice behaves as in an ice stream, in which vertical shearing is present, though minimal. Ice flow in this experiment can be considered as ice shelf flow with minimal basal traction. The invalidity of the SIA solution is shown in Fig. 7, where the analytical SIA solution is plotted for a simplified basal sliding relationship in which the basal shear stress is supposed to balance the driving stress without longitudinal stress gradients, so that

$$v_x(z_b) = (\rho g H \tan \alpha) \beta^{-2}, \qquad (28)$$

in Eq. (27). However, a singularity in the velocity occurs in Eq. (28) for $\beta^2=0$. Not all velocities are therefore plotted in Fig. 7. As for Exps. A and B, the SIA solution is independent of L. Again, full-Stokes models show a smaller spread than the higher-order approximations (Figs. 8–9 and Table 4). The SIA is definitely not suited for simulating this type of ice flow where longitudinal stresses dominate over vertical shear stresses (compare with Fig. 7). As in Exps. A–B, the amplitude of the surface velocity field decreases with increasing aspect ratio $\epsilon$. At the smallest length scales the surface velocity field is almost constant as the transition between high and low friction regions is smoothed out by longitudinal stress transmission. Similar to Exp. B, a surface velocity field anti-correlated to basal friction is observed for full-Stokes models at L=5 km (Exp. D). For the other higher-order approxima-
<table><caption>**Table 6.** Mean values (μ), standard deviation (σ) and number of participating models (n) of the maximum and mean horizontal ice velocity at the surface in the direction of the flow for Exp. E. Units are m a⁻¹.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Sliding</th><th colspan="2">No sliding</th></tr><tr><th>μ</th><th>(σ)</th><th>μ</th><th>(σ)</th></tr></thead><tbody><tr><th colspan="5" style="text-align: left;">Maximum velocity</th></tr><tr><td>NFS</td><td>67.01</td><td>(3.03)</td><td>6</td><td>122.44</td><td>(41.90)</td><td>5</td></tr><tr><td>FS</td><td>65.95</td><td>(0.63)</td><td>5</td><td>110.62</td><td>(32.86)</td><td>5</td></tr><tr><th colspan="5" style="text-align: left;">Mean velocity</th></tr><tr><td>NFS</td><td>32.00</td><td>(0.94)</td><td>6</td><td>47.15</td><td>(9.24)</td><td>5</td></tr><tr><td>FS</td><td>32.21</td><td>(0.13)</td><td>5</td><td>44.61</td><td>(8.98)</td><td>5</td></tr></tbody></table>

tions this is not observed (but due to the larger disparity in solutions, this effect is unnoticeable in Fig. 9).

In general, the spread in results of the modeled velocity field is higher than for the experiments over the bedrock bumps. The smallest spread is obtained with full-Stokes models, and this spread reduces with increasing ε, contrary to the results from Exps. A–B. L1L1 and L1L2 models have larger spreads than the LMLa models, despite the fact that they were designed for coping with such type of ice flow in the first place.

## 5.3 Experiment E: Haut Glacier d'Arolla

Although the input file lists the bedrock and surface data along the flowline of Haut Glacier d'Arolla with a fixed grid spacing of Δx=100 m, most participants interpolated this dataset at a higher resolution (Fig. 10). The effect of subsampling is captured in Fig. 11, where the oscillations in the basal shear stress along the flowline are either jagged when undersampled or smoother when a sufficiently small grid size is chosen. Again, the spread in surface velocity for full-Stokes models is smaller than for the other approximants, albeit that for the no-slip case, the standard deviation is small for all models.

Inclusion of a sliding zone (an area of zero basal friction) leads to larger differences among the different participating models. Also the full-Stokes models show a much larger spread of solutions. Here, increasing resolution results in other complexities compared to the no-sliding case, such as the occurrence of oscillations in the basal shear stress. The slip/no-slip boundaries are very sensitive to model resolution, as they can be regarded as singularities where the friction parameter β² suddenly jumps from zero to infinity and vice versa. In particular, the linear interpolation leads to break points in basal and surface topography that influence the result. The results of the sliding experiment underline the difficulty of simulating end-member behavior in basal sliding (slip/no-slip).

<table><caption>**Table 7.** Mean values (μ), standard deviation (σ) and number of participating models (n) of the maximum and mean horizontal ice velocity at the surface in the direction of the flow for Exp. F. Units are m a⁻¹.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Sliding</th><th colspan="2">No sliding</th></tr><tr><th>μ</th><th>(σ)</th><th>μ</th><th>(σ)</th></tr></thead><tbody><tr><th colspan="5" style="text-align: left;">Maximum velocity</th></tr><tr><td>NFS</td><td>98.14</td><td>(0.35)</td><td>5</td><td>197.55</td><td>(0.48)</td><td>5</td></tr><tr><td>FS</td><td>98.64</td><td>(0.16)</td><td>2</td><td>197.85</td><td>(0.01)</td><td>2</td></tr><tr><th colspan="5" style="text-align: left;">Mean velocity</th></tr><tr><td>NFS</td><td>96.11</td><td>(0.40)</td><td>5</td><td>193.05</td><td>(1.16)</td><td>5</td></tr><tr><td>FS</td><td>96.42</td><td>(0.05)</td><td>2</td><td>194.67</td><td>(0.04)</td><td>2</td></tr></tbody></table>

## 5.4 Experiment F: prognostic run

Benchmarking of numerical ice sheet models is possible when analytical solutions exist for a particular problem. The analytical solutions used here are derived from first-order perturbation analysis of flow down a uniformly inclined plane (Gudmundsson, 2003). It is inherent in this type of analysis that the resulting flow perturbations are linear functions of basal amplitudes. Numerical solutions are usually not limited by this assumption and can therefore, for any finite amplitude perturbation in basal properties, be better approximations to the Stokes equations than the analytical solutions given by Gudmundsson (2003). For an accurate comparison with the analytical solutions, numerical solutions must be calculated for a number of different amplitudes and then scaled by forming the ratio between each solution and the respective basal amplitude. If this ratio is found to be independent of amplitude for small amplitudes, the scaled numerical solutions can be compared to the analytical ones. This kind of test was done by Raymond and Gudmundsson (2005). The exact error estimate depends on wavelength, amplitude, and slip ratio. (The reader is referred to Raymond and Gudmundsson (2005) and Gudmundsson (2008) for a detailed discussion.) For example, the analytical solutions were found to be generally accurate to within ~1% for sinusoidal perturbations with wavelength larger than ~10H and amplitudes less than ~0.1H. For Exp. F, we expect a similar degree of agreement between the analytical solutions and exact Stokes solvers.

Only a few models participated in this experiment, since not all of them treat an evolving free surface. The objective of the test was to run the models forward in time until a steady state was reached. The definition of steady state was left to the interpretation of each participant. The resulting steady-state surface elevation and velocity are shown in Figs. 12–13.

We present the analytical solution here for the case without basal sliding. The two full-Stokes models show a very good agreement with the analytical solution. Higher-order model solutions also fit also well but show a slightly larger vari-
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Results for Exp. C: norm of the surface velocity at $y=L/4$ for different length scales $L$. The mean value and standard deviation are shown for both types of models.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Results for Exp. D: norm of the surface velocity for different length scales $L$. The mean value and standard deviation are shown for both types of models.</figcaption></figure>

ability, especially for the velocity. This variability increases when basal sliding is introduced. As mentioned before, the intercomparison exercise is based on numerical solutions of shallow and less-shallow continuum approximations of full-Stokes and numerical solutions of the full-Stokes problem it-

self, compared with an asymptotic analysis of the full-Stokes system. The full-Stokes model solutions for the surface elevation lie closer to the analytical solution than do the higher-order solutions, suggesting that the analytical model and the numerical full-Stokes models are both adequate descriptions
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Surface velocity in the direction of the ice flow for Exp. E for the no-sliding (top) and sliding (bottom) experiment.</figcaption></figure>

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

<figure><img src="image_16.png" /><figcaption>**Fig. 11.** Basal shear stress in the direction of the ice flow for Exp. E for the no-sliding (top) and sliding (bottom) experiment.</figcaption></figure>

of the full-Stokes equations. However, the analytical solution for the velocity is well away from the tightly-clustered FS numerical results. Here, the FS solution might be more accurate, as the analytical solution remains a solution based on a perturbation expansion and therefore is an approximate method for solving the full-Stokes equations.

# 6 Conclusions

In this paper we present the results of the first intercomparison exercise of higher-order and full-Stokes ice sheet models. A total of 27 different numerical models participated in this benchmarking effort. Six experiments were designed to evaluate complex ice flow with high spatial variability in basal topography and slipperiness. All experiments were designed in such a way that the Shallow-Ice Approximation (SIA) is not valid, especially at high aspect ratios. Although the SIA is valid for large parts of ice sheets, higher-order models are necessary for describing ice flow in areas of high basal topographic variability and slipperiness, which are generally the most dynamic regions of ice sheets.

Compared to previous benchmark experiments (Huybrechts et al., 1996; MacAyeal et al., 1996; Payne et al., 2000), a significantly higher number of ice sheet models participated in this benchmark. Despite the greater complexity

of the problem, all models produce results that are in approximate agreement, even for high aspect ratios. This shows that numerical ice sheet models have improved considerably during the past decade, and are capable of simulating ice flow in regions where longitudinal stress gradients are important.

Full-Stokes models produce reliable results in the sense that (i) their spread of results is very low ($<0.2 \text{ m a}^{-1}$) and (ii) they give a result in concordance with analytical solutions based on perturbation theory. Results of the higher-order approximants show a significant larger spread that cannot be attributed solely to numerical issues, such as the discretization error. The greatest spread is found for high aspect ratios where all stress components (not only membrane/longitudinal stresses) are equally important, and higher-order approximations are insufficient. Also, coding errors could be present, since all higher-order models are coded by the authors (which is less true for the full-Stokes models).

Only a limited number of L1L1 and L1L2 models participated. They show the same spread in solutions for basal topographic perturbations (Exps. A and B) as the LMLa models. However, they produce a larger spread for the ice stream simulations at relatively low aspect ratios (Exps. C and D) compared to LMLa models, although they were designed for such flows. The LMLa results differ from the full-Stokes results
<figure><img src="image_17.png" /><figcaption>**Fig. 12.** Steady state surface elevation along the central flowline for Exp. F for the no sliding (top) and sliding (bottom) experiment. The black line indicates the analytical solution.</figcaption></figure>

for high aspect ratios and in regions of variable basal sliding. All models (including full-Stokes) agree poorly when sudden variations in basal friction are considered, such as the slip/no-slip jumps in Exp.E.

Finally, the full-Stokes models give the most consistent results. Results from different models show a very small spread and are in good agreement with the analytical solution based on perturbation theory. For most experiments a clear distinction can be made between results from full-Stokes models and from higher-order approximants.

Future model intercomparisons should definitely focus on transient problems, as preliminarily explored in Exp. F, and slip/no-slip boundary problems, such as Exp. E. Analytical solutions to specific (simplified) problems are something to look into as well, which would lead to more reliable benchmarks and verification.

*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). The authors are indebted to FRIA (Fonds pour la formation à la Recherche dans l'Industrie et dans l'Agriculture). This paper forms a contribution to the Belgian Research Programme on the Antarctic (Belgian Federal Science Policy Office), Project nr. SD/CA/02. The authors are deeply

<figure><img src="image_18.png" /><figcaption>**Fig. 13.** Norm of the steady state surface velocity along the central flowline for Exp. F for the no sliding (top) and sliding (bottom) experiment. The black line indicates the analytical solution.</figcaption></figure>

indebted to Ed Bueler for his insightful remarks and suggestions that improved the initially submitted manuscript substantially, and to the Editor Bill Lipscomb for improving the readability of the manuscript. We should like to thank James Fishbaugh for writing the scripts that produced the graphs.

Edited by: W. Lipscomb

## References

* Blatter, H.: Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344, 1995.

* Blatter, H., Clarke, G., and Colinge, J.: Stress and velocity fields in glaciers: Part II. sliding and basal stress distribution, J. Glaciol., 44, 457–466, 1998.

* Breuer, B., Lange, M., and Bl window, N.: Sensitivity studies on model modifications to assess the dynamics of a temperate ice cap, such as that on King George Island, Antarctica, J. Glaciol., 52, 235–247, 2006.

* Colinge, J. and Rappaz, J.: A strongly nonlinear problem arising in glaciology, Mathem. Modelling Numer. Analysis, 33, 395–406, 1999.

* Gagliardini, O. and Zwinger, T.: The ISMIP–HOM benchmark experiments performed using the Finite-Element code Elmer, The Cryosphere, 2, 67–76, 2008.
Gudmundsson, G.: Transmission of basal variability to a glacier surface, J. Geophys. Res., 108(B5), 2253, doi:10.1029/2002JB002107, 2003.

Gudmundsson, G.: Analytical solutions for the surface response to small amplitude perturbations in boundary data in the shallow-ice-stream approximation, The Cryosphere, 2, 77–93, 2008.

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

Hindmarsh, R. C. A.: The role of membrane-like stresses in determining the stability and sensitivity of the antarctic ice sheets: back pressure and grounding line motion, Phil. Trans. R. Soc., A, 1733–1767, 2006.

Hubbard, A., Blatter, H., Nienow, P., Mair, D., and Hubbard, B.: Comparison of a three-dimensional model for glacier flow with field data from Haut Glacier d'Arolla, Switzerland, J. Glaciol., 44, 368–378, 1998.

Hutter, K.: Theoretical Glaciology, Dordrecht, Kluwer Academic Publishers, 1983.

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

Huybrechts, P., Abe-Ouchi, A., Marsiat, I., Pattyn, F., Payne, T., Ritz, C., and Rommelaere, V.: Report of the Third EISMINT Workshop on Model Intercomparison, European Science Foundation (Strasbourg), 1998.

IPCC: Contribution of working group i to the fourth assessment report of the intergovernmental panel on climate change, in: Climate Change 2007: The Physical Science Basis., edited by: Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Avryekt, K., Tignor, M., and Miller, H., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, p. 996, 2007.

Johnson, J. and Staiger, J.: Modeling long-term stability of the Ferrar Glacier, East Antarctica: Implications for interpreting cosmogenic nuclide inheritance, J. Geophys. Res., 112, F03S30, doi:10.1029/2006JF000599, 2007.

Ladyzhenskaya, O.: The Mathematical Theory of Viscous Incompressible Flow, New York, Gordon and Breach, 1969.

MacAyeal, D.: Large-scale ice flow over a viscous basal sediment: Theory and application to Ice Stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, 1989.

MacAyeal, D.: A tutorial on the use of control methods in ice-sheet modeling, J. Glaciol., 39, 91–98, 1993.

MacAyeal, D., Rommelaere, V., Huybrechts, P., Hulbe, C., Determann, J., and Ritz, C.: An ice-shelf model test based on the Ross ice shelf, Ann. Glaciol., 23, 46–51, 1996.

Martin, C., Navarro, F., Otero, J., Cuadrado, M., and Corcuera, M.: Three-dimensional modelling of the dynamics of Johnson Glacier (Livingstone Island, Antarctica), Ann. Glaciol., 39, 1–8, 2003.

Nienow, P., Hubbard, A., Hubbard, B., Chandler, D., Mair, D., Sharp, M., and Willis, I.: Hydrological controls on diurnal ice flow variability in valley glaciers, J. Geophys. Res., 110, F04002, doi:10.1029/2003JF000112, 2005.

Pattyn, F.: Transient glacier response with a higher-order numerical ice-flow model, J. Glaciol., 48, 467–477, 2002.

Pattyn, F.: A new 3D higher-order thermomechanical ice-sheet model: Basic sensitivity, ice-stream development and ice flow across subglacial lakes, J. Geophys. Res., 108(B8), 2382, doi:10.1029/2002JB002329, 2003.

Pattyn, F.: Investigating the stability of subglacial lakes with a full Stokes ice-sheet model, J. Glaciol., 54, 353–361, 2008.

Payne, A., Huybrechts, P., Abe-Ouchi, A., Calov, R., Fastook, J., Greve, R., Marshall, S., Marsiat, I., Ritz, C., Tarasov, L., and Thomassen, M.: Results from the EISMINT model intercomparison: the effects of thermomechanical coupling, J. Glaciol., 46, 227–238, 2000.

Pollard, D. and DeConto, R.: A coupled ice-sheet/ice-shelf/sediment model applied to a marine-margin flowline: Forced and unforced variations, in: International Association of Sedimentologists Special Publication, 39, 37–52, 2007.

Price, S., Waddington, E., and Conway, H.: A full-stress, thermomechanical flowband model using the finite volume method, J. Geophys. Res., 112, F03020, doi:10.1029/2006JF000724, 2007.

Raymond, M. and Gudmundsson, G.: On the relationship between surface and basal properties on glaciers, ice sheets, and ice streams, J. Geophys. Res., 110, B08411, doi:10.1029/2005JB003681, 2005.

Sugiyama, S., Gudmundsson, G., and Ebling, J.: Numerical investigation of the effects of temporal variations in basal lubrications on englacial strain-rate distribution, Ann. Glaciol., 37, 49–54, 2003.

Van der Veen, C. and Whillans, I.: Force budget: I. theory and numerical methods, J. Glaciol., 35, 53–60, 1989.

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