We present a series of benchmark experiments designed for testing and comparing numerical ice-sheet models. Following the outcome of two EISMINT workshops organized to intercompare large-scale ice-sheet models currently in operation, model benchmark experiments ate described for ice sheets under fixed and moving margin conditions. These address both steady-state and time-dependent behaviour under schematic boundary conditions and with prescribed physics. A comparison was made of each model’s prediction of basic geophysical variables such as ice thickness, velocity and temperature. Consensus achieved in the model inter-comparison provides reference solutions against which the accuracy and consistency of ice-sheet modelling codes can be assessed.

An ice-shelf model which features efficient numerical techniques is developed to determine the back-force exerted by sides and pinning points, such as islands of an embayed ice shelf. The model is applied to three ideal geometries and shows that the restraint exerted by a small island, even far downstream from the grounding line, can represent about one-half of the total restraint due to the embayment. Our results are further interpreted to suggest several criteria useful for testing any ice-shelf model.

A simple computer scheme developed by Budd and Smith (1985) and modified by D. Jenssen has been further developed to provide a rapid computation of steady-state balance fluxes over arbitrary ice masses, given the surface elevations and net accumulation distribution. The scheme provides a powerful diagnostic tool to examine the flux and state of balance over whole ice masses or limited regions to interpret field observations for dynamics or the state of balance.

In many cases the uncertainty in the state of balance may be much less than the uncertainty in the deformation and sliding properties of the ice and so the flux and velocities derived from balance could provide a useful guide for the dynamics where direct observations are sparse.

The scheme assumes that, on a horizontal scale of many ice thicknesses, the ice-flow direction is approximately down the steepest surface slope. The continuity equation is used to compute steady-state implied downslope fluxes at each grid point from integrations of the net accumulation over the area from the summits to the edges. The algorithm ensures the exact integral balance of the surface net flux over the area with flow through boundaries.

Applications are demonstrated for the whole of Antarctica and for regional areas. Comparisons are made between fluxes computed from observed ice thicknesses and velocities and those computed from balance. The observed ice thicknesses can also be used to compute surface velocities from assumed column-to-surface velocity ratios. The combined fluxes from observations and balance can be used to compute rates of change of elevation with time.

The equation describing the surface evolution of a large ice sheet is examined on the basis of a scale analysis applied to Antarctic conditions. Changes in the surface elevation are mainly driven by mass-balance fluctuations which approximately follow global atmospheric temperature variations. The essential spatial non-uniformity of the accumulation rate and the resultant difference between central and coastal regions in reaction time-scales are taken into account. The dynamic interaction of the time-lagging interior with the quasi-stationary margin is described. As a result, a simplified model is deduced to simulate the surface-elevation variations in the central parts of the Antarctic ice sheet caused by mass-balance perturbations corresponding to the main Milankovich cycles with the periods of 19–100 kyears. Special computational tests are performed to validate the model through intercomparison with the predictions obtained with a two-dimensional thermomechanical model. The sensitivity of the model to physical factors (represented by dimensionless tuning parameters) is discussed. Climatically controlled variations of the ice-sheet thickness in the vicinity of Vostok Station during the past 200 kyears are estimated.

A degree-day model and an energy-balance model for the Greenland ice sheet are compared. The two models are compared at a grid with 20 km spacing. Input for both models is elevation, latitude and accumulation. The models calculate the annual ablation over the entire ice sheet. Although on the whole the two models yield similar results, depending on the tuning of the models, regional discrepancies of up to 45% occur, especially for northern Greenland. The performance of the two types of model is evaluated by comparing the model results with the sparsely available (long-term) mass-balance measurements. Results show that the energy-balance model tends to predict a more accurate mass-balance gradient with elevation than does the degree-day model.

Since so little is known about the present-day climate of the ice sheet, it is more useful to consider the sensitivity of the ablation to various climate elements than to consider the actual present-day ablation. Results show that for a 1 K temperature perturbation, sea-level rise is 0.31 mm year−1 for the energy-balance model and 0.34 mm year−1 for the degree-day model. After tuning the degree-day model to a value of the ablation, equivalent to the ablation calculated by the energy-balance model, sensitivity of the degree-day model increases to 0.37 mm sea-level change per year. This means that the sensitivity of the degree-day model for a 1 K temperature perturbation is about 20% higher than the sensitivity of the energy-balance model. Another set of experiments shows that the sensitivity of the ablation is dependent on the magnitude of the temperature perturbation for the two models. Both models show an increasing sensitivity per degree for larger perturbations. The increase in the sensitivity is larger for the degree-day model than for the energy-balance model. The differences in the sensitivity are mainly concentrated in the southern parts of the ice sheet.

Experiments for the Bellagio temperature scenario. 0.3°C increase in temperature per decade, leads to sea-level rise of 4.4 cm over a period of 100 years for the energy-balance model. The degree-day model predicts for the same forcing a 5.8 cm rise which is about 32% higher than the result of the energy-balance model.

A standard numerical experiment featuring the Ross Ice Shelf, Antarctica, is presented as a test package for the development and intercomparison of ice-shelf models. The emphasis of this package is solution of stress-equilibrium equations for an ice-shelf velocity consistent with present observations. As a demonstration, we compare five independently developed ice-shelf models based on finite-difference and finite-element methods. Our results suggest that there is little difference between finite-element and finite-difference methods in capturing the basic, large-scale flow features of the ice shelf. We additionally show that the fit between model and observed velocity depends strongly on the ice-shelf temperature field, for which there is presently little observational control. The main differences between model results are due to the equations being solved, the boundary conditions at the ice from and the discretization method (finite element vs finite difference).

Ice-sheet modelling is an essential tool for estimating the effect of climate change on the Greenland ice sheet. The large spatial and long-term temporal scales of the ice-sheet model limits the amount of data which can be used to test model results. The geological record is useful because it provides test material on the time-scales typical for the memory of ice sheets (millennia). This paper compares modelled ice-margin positions with a geological scenario of ice-margin positions since the Last Glacial Maximum to the present in West Greenland. Morphological evidence of ice-margin positions is provided by moraines. Moraine systems are dated by 14C-dated marine shells and terrestrial peat. Three Greenland ice-sheet models are compared. There are distinct differences in modelled ice-margin positions between the models and between model results and the geological record. Disagreement between models and the geological record in the near-coastal area is explained by the inadequate treatment of marginal processes in a tide-water environment. A smaller than present ice sheet around the warm period in the Holocene (Holocene climatic optimum) only occurs if such a period appears in the forcing (ice-core record) or used temporal resolution. Smoothing of the GRIP record with a 2000 year average eliminates the climatic signal related to the Holocene climatic optimum. This underlines the importance of short-term and medium-term variations (decades, centuries) in climatic variables in determining ice-margin positions in the past but also in the future.

All parts of a two-dimensional, isothermal, stationary marine glacier (grounded ice sheet, ice shelf and transition zone) with constant viscosity are analysed by perturbation methods. In so doing, all zones of different flow patterns can be considered separately. Correlations between spatial scales for all parts can be expressed in terms of the typical ice-surface slope distant from the ocean, which reflects exterior conditions of the glacier’s existence. In considering the ice-sheet–ice-shelf transition zone, a small parameter characterizing the difference between ice and water densities is used. Such an analysis allows us to find boundary conditions at the grounding line for the grounded ice mass. Glacier-surface profiles are determined by numerical methods. The grounding-line position found by using the boundary conditions derived in this paper differs from that obtained by using Thomas and Bentley’s (1978) boundary conditions by about 10% of the grounded ice-stream length.

A simplified model of ice-sheet behaviour is described. It combines the assumptions of rapid ice flow, high viscous activation energy and realistic sediment-based sliding dynamics to form a non-linear diffusion-type equation which can display relaxation oscillations analogous to those of surging glaciers, and which may be relevant to large-scale surges of the Hudson Strait and Cabot Strait ice streams of the Laurentide ice sheet.

When the physics of this model is applied to a laterally extensive unidirectional ice flow, such as that in the Siple Coast of Antarctica, an appropriate mechanism may exist for the spontaneous generation of ice streams.

Various spatial discretizations for the ice sheet are compared for accuracy against analytical solutions in one and two dimensions. The computational efficiency of various iterated and non-iterated marching schemes is compared.

The stability properties of different marching schemes, with and without iterations on the non-linear equations, are compared. Newton–Raphson techniques permit the largest time steps. A new technique, which is based on the fact that the dynamics of unstable iterated maps contain information about where the unstable root lies, is shown to improve substantially the performance of Picard iteration at a negligible computational cost.

An adaptive-grid finite-volume glacier model is described. The model is an implicit one-dimensional flowline model. The discretized implicit finite-volume equations are solved by an iterative predictor–corrector method. The grid adapts as the terminus moves in response to changes in surface mass balance. Only the terminus grid point and the penultimate grid point are adapted as the glacier-terminus position changes in order to minimize computation. Several modelling experiments are carried out to demonstrate the performance of the model. Comparisons are made with a fixed-grid finite-volume model and a fixed-grid finite-difference model. Comparisons are made on two levels. The differences in methods, finite-volume method versus finite-difference method, arise from differences in accuracy and programming efficiency. The differences in grids, adaptive-grid versus fixed-grid, arise from differences in the numerical smoothness of the motion of the moving terminus. This affects questions of stability and accuracy.

Perturbation of divide position is considered by a linearization about the Vialov–Nye solution and also about related solutions with O(1) relief. Relaxation times of one-sixteenth the fundamental thickness/accumulation-rate time-scale are found for the Vialov–Nye configuration, while substantial basal topography can halve the rate of relaxation. Steady divide position is most sensitive to anti-symmetric accumulation-rate distributions near the divide, but transient divide motion is most strongly excited by anti-symmetric accumulation rate variations halfway between the margin and the divide. Relaxation times for the Antarctic Peninsula divide position are estimated to be around 200 years, while the larger Greenland ice sheet has a divide-position relaxation time of around 600 years.

Modelling accumulation rate as a white-noise process permits analysis of divide perturbation as a (stochastic) Ornstein–Uhlenbeck process, where the standard deviation of the response is proportional to the standard deviation of the forcing. If observed accumulation-rate variability in the Antarctic Peninsula were anti-symmetric about the divide, it would be sufficient to force the divide position to fluctuate with standard deviation 10–20 times the depth of the ice sheet. There appears to be sufficient noise to cause Raymond bumps to be spread significantly. More data on the statistical variation of accumulation with position are needed. Random forcing will increase the complexity of any fold structures created in the divide region and in particular the number of such structures intersecting any borehole.

An analysis of the linear stability of marine ice sheets uncoupled from associated ice shelves is presented. The principal feature is a zero eigenvalue associated with infinitesimal shifts along the line of neutral equilibrium in phase space, termed the “equilibrium manifold”. A finite-difference scheme is constructed which respects this stability properly.

The zero eigenvalue appears to allow modelling errors to accumulate rather than dissipate as occurs in land-based ice sheets. The practical significance of this is that even rather fine spatial grids may allow substantial numerical error to accumulate.

A finite-element model is developed in order to calculate the coupled ice and heat flow and the surface topography in cold, steady-state ice sheets. The model decouples the heat-flow equation and the surface mass-balance condition from the rest of the equations and solves the problem by an iterative method. The model is used to examine the thermomechanics of ice divides. Initial studies of a symmetric, plane ice divide and an axisymmetric ice divide have led to the following conclusions, which are consistent with previous results. The ice-divide zone is a narrow region, only a few ice thicknesses wide, where the surface slope drops to zero and the flow solution changes. The longitudinal strain rate is high, especially in the upper layers, and the vertical velocity is smaller than away from the divide. This causes the basal temperatures to increase and the isochrones to rise. Divergent-flow conditions widen the ice-divide zone, whereas they do not influence the solution at the ice divide.

In this work, we formulate a model for the isothermal flow of a (basal) ice–till mixture that is overlain by a layer of pure ice. Such a model is relevant to the case of a glacier or ice sheet possessing a till at its base. To this end, ice is treated as usual as a constant true-density, very viscous fluid, while till, which is assumed to consist of sediment and bound (i.e. moving with the sediment) interstitial water, is also assumed in a first approximation to behave as such a fluid. Since the mixture is assumed isothermal, only the mass- and momentum-balance relations for till and ice need be considered. To complete the model, no-slip and stress-free boundary conditions are assumed at the base and free surface, respectively. By working with the former conditions, we neglect the process of entrainment of sediment into the basal layer, concentrating rather on its flow behaviour and thickness. The transition from the till–ice mixture layer to the overlaying pure ice layer is idealized in the model as a moving interface representing in the simplest case the till material boundary, at which jump-balance relations for till and ice apply. As in the basal layer, till and ice are assumed to interact mechanically at this interface. In the context of the thin-layer approximation, numerical solutions of the lowest-order form of the model show that it is predominantly the thickness of the basal (mixture) layer that is influenced by the ice–till momentum interaction.

A simple model is developed based on the notion that on active ice streams the resistance to flow is partitioned between basal drag and lateral drag. The relative roles of these sources of resistance is determined by a friction parameter that effectively describes the strength of the bed under the ice stream. Reduction in the basal strength is caused by meltwater production, taken proportional to the product of basal drag and ice speed. The width of the ice stream is governed by the balance between entrainment or erosion of ice from the slow-moving inter-stream ridges and advection from the ridges into the ice stream. Entrainment of ridge ice is parameterized as a function of the shear stress at the lateral margins, in one case proportional to the lateral shear stress and in the second case scaled to ice-stream width. In the first formulation, the model rapidly becomes unstable but, using the second formulation, a steady state is reached with lateral drag providing all or most of the resistance to flow. The results point to the great importance of achieving an understanding of entrainment. With the second model and a wide range of parameter values, there is no cyclic behavior, with rapid flow being followed by a quiescent phase.

A finite-element model for simulating multi-dimensional air flow with heat, mass and chemical species transport through firn is discussed. The model is applied to an investigation of near-surface layering effects on ventilation rates. Field measurements of permeability at Summit, Greenland, are presented that show that permeability varies by at least a factor of 10 over the top 3 m, with the surface windpack having much lower permeability, in general, than the underlying firn. The effect of a lower-permeability surface layer is to decrease the air flow in the underlying firn, yet there is still sufficient air flow in the top meters of the firn so that ventilation must be considered for species transport. Channeling, or increased air flow in a layer overlain by a less-permeable layer, can occur even if the microstructure of each layer is isotropic. Conventional estimates of chemical transport due to diffusion alone are likely to underestimate transport, while estimates of ventilation that consider the firn as a homogeneous half-space may overestimate ventilation effects at the near-surface. Effects of firn layering are important for ventilation and must be considered for accurate assessment of firn–air transport mechanisms.

In Kronprins Christian Land at 80° N in eastern north Greenland, it has been observed that the surface of the Wisconsin ice is significantly darker than the Holocene ice found immediately upstream from a transition located 710 m from the ice margin. δ18O analysis has shown the dark surface is of Wisconsin origin. Deep ice cores from the Greenland ice sheet all indicate that the Wisconsin ice contains orders of magnitude more microparticles which could he the reason for the dark appearance of the Wisconsin surface.

Photographic documentation, spectral surface-albedo measurements and satellite-image analysis all indicate a reduced albedo of Wisconsin ice. The effect of this reduced albedo is not confirmed by the ablation measurements. However, measurements of ablation variability within small test sites has documented that large errors will arise if only one stake per measuring point is used. Energy-balance calculations show ablation rates should be 10–70%, i.e. 1.8–8.4 mm d−1, respectively, less than experienced. Additionally, a satellite-image analysis shows even higher albedo contrasts 20–30 km to the south of our transect locality.

Immediately after the termination of the ice age, most of the surface in the ablation zone consisted of this low-reflectance ice. So, in the early Holocene, the dark ice of Wisconsin origin is likely to have resulted in higher ablation rates than previously considered. This may probably partly explain the fast rate of retreat/disintegration of the ice sheets in the Northern Hemisphere, after the termination of the Wisconsin ice age.

Simple models to calculate melt and refreezing are reviewed. Both degree-day and energy-balance models can give distributed melt inputs to ice-dynamics models but have only been tested extensively in West Greenland, and more data are needed from the remoter parts of Greenland. The energy-balance model is more realistic but needs input data that are not generally available over the whole ice sheet. On the other hand, degree-day factors vary from situation to situation although a value of about 8 kg m−2 d−1 deg−1 for ice ablation is a reasonable first approximation as assumed in recent ice-dynamics models. Meltwater refreezing in the accumulation area can be modelled very simply but more detailed physical models are needed to describe the shifts in accumulation zones under different climates.

A two-dimensional vertically integrated ice-flow model has been used to simulate the current state of the ice cap of King George Island, South Shetland Islands, Antarctica, as well as the sensitivity of this state to climate change. The model was forced by an energy-balance model that generates the specific mass balance from climatological input data of two research stations. It proved difficult to simulate-satisfactorily the entire geometry of the present-day ice cap. Nevertheless, it was possible to simulate a steady-state ice cap whose volume and areal extent approximate the (estimated) current situation. Several experiments have indicated that this state is highly sensitive to climate change. The model predicts that cooling by 1 K will increase the ice volume by 10% and warming by 1 K will decrease it by 36%. A 10% change in precipitation will alter the ice volume by less than 8%. Application of the IPCC-90 Business-as-Usual scenario leads to a 55% reduction in the ice volume by the year AD 2100, compared to the present-day situation. The response of the ice cap to warming is therefore totally different from the response of the main Antarctic ice sheet which is believed to gain mass by increasing temperatures.