Alfred-Wegener-Institut für Polar- und Meeresforschung, Postfach 120161, D-27515 Bremerhaven, Germany; and Departement Geografie, Vrije Universiteit Brüssel, Pleinlaan 2, B-1050 Brüssel, Belgium
Ice sheets respond dynamically to changes in boundary conditions, such as climate variations, basal thermal conditions, and isostatic adjustments of the underlying bedrock. These cause the ice sheets to evolve towards a new equilibrium. Long response time-scales of up to 104 years are involved, determined by the ratio of ice thickness to yearly mass turnover, physical and thermal processes at the bed, and processes affecting ice viscosity and mantle viscosity. The response of the ice sheets is further complicated by feedback processes which may amplify or mitigate the ice sheet's adjustment to the forcing or by internal instabilities that may cause rapid changes in ice volume due to changes in the dynamic flow regime. A primary motivation for developing numerical models of ice flow is to gain a better understanding of the spatial and temporal behaviour of ice sheets and glaciers and to predict their response to external forcing. Modelling ice-sheet dynamics presents a powerful framework to investigate the complex interactions between the ice sheets and the climate system in a quantitative way, in past as well as future environments. Ice-flow models are commonly based on fundamental physical laws and assumptions thought to describe glacier flow.
At the top end of the class of ice-sheet models are so-called three-dimensional thermomechanical models, which are able to describe the time-dependent flow and shape of real ice sheets. These models are akin to general circulation models developed in other branches of climate science. Their development closely follows technical progress in such fields as computer power, ice-core and sediment drilling, remote sensing and geophysical dating techniques, which are providing both the required calculating means and the necessary data to feed and validate such models. Models of this type have been applied to the existing ice sheets of Greenland and Antarctica, and to those which covered the continents of the Northern Hemisphere during the Quaternary ice ages. Typical studies have concentrated on mechanisms and thresholds of ice-sheet inception during the Tertiary (Huybrechts, 1994a; DeConto & Pollard, 2003), ice-sheet form and extent during glacial-interglacial cycles (Marshall et al., 2000; Ritz et al.,
2001; Charbit et al., 2002; Huybrechts, 2002), and the response of the polar ice sheets to future climatic warming (Huybrechts & de Wolde, 1999; Van de Wal et al., 2001). In this context, the key interactions being investigated are between the effects of a change in climate on the accumulation and ablation fields and the ice sheet's response in terms of changed geometry and flow, including the ice sheet's contribution to the worldwide sea-level stand. Related work has considered the ice sheets as a boundary condition for other components of the Earth's geophysical system, providing changes in surface loading for isostasy and gravity models (Le Meur & Huybrechts, 2001; Tarasov & Peltier, 2004), or providing changes in freshwater fluxes for ocean models, especially to investigate changes of the thermohaline circulation of the North Atlantic Ocean (Schmittner et al., 2002; Fichefet et al., 2003). Three-dimensional thermomechanical ice-sheet models have also been used to investigate the potential for internally generated flow instability (Payne, 1995; Payne & Dongelmans, 1997; Marshall & Clarke, 1997a). In this application, the crucial interactions are between the thermal and flow regimes. In addition, models of the Greenland and Antarctic ice-sheets are in use to assist with the location and dating of ice cores (Greve, 1997b; Huybrechts et al., 2004a), estimating internal distributions of passive tracers such as oxygen isotope ratios (Clarke & Marshall, 2002), yield information about fields that are inaccessible for direct observation such as at the ice-sheet base (Huybrechts, 1996), or assess the component of their present-day evolution due to adjustment to past climate changes (Huybrechts & Le Meur, 1999).
In this chapter the discussion concentrates on three-dimensional whole ice-sheet models applied to the ice sheets of Antarctica and Greenland. This has the advantage that the quality of the model simulations can be assessed against available observations. Also, the range of physical characteristics, climate regimes, and flow mechanisms encountered in both polar ice sheets is probably a good representation of the range of behaviour which occurred in the palaeo-ice-sheets at various times during their evolution. The Antarctic ice sheet is located in a very cold climate, where almost no surface melting occurs and precipitation amounts are limited by low air temperature. Therefore virtually all Antarctic ice is eventually transported into floating ice shelves that experience melting or freezing at their underside and eventually break up to form icebergs. External forcing of such a cold ice sheet is mainly through changes in accumulation rate and basal melting rates below the ice shelves. The ice-sheet extent is limited mainly by the depth of the surrounding ocean and by the capacity of the ocean to float the ice shelves, touching upon the crucial issue of grounding-line dynamics. An important characteristic of the West Antarctic ice sheet is that it is a marine ice sheet, resting on a bed far below sea level. Much of its ice transport towards the coast occurs in ice streams, which are distinct fast-flowing features that rest on smooth sedimentary beds, have very flat surface profiles, and are sharply bordered by relatively stagnant ice at their sides. Their ice flux is dominated by basal flow, mostly through deformation of the underlying water-saturated sediments (Alley et al., 1986a). By contrast, the Greenland ice sheet is situated in a much warmer climate, with a temperature difference of 10-15°C in the annual mean. Summer temperatures are high enough to initiate widespread summer melting. All around the ice-sheet margin, mean annual ablation exceeds the accumulation. A negative surface budget results at elevations below about 1000m in the north and 1600-1800m in the southwest. High coastal temperatures do not favour ice shelves, but there are a few along the north and northeast coast. The Greenland ice sheet loses mass by calving of icebergs, mostly at grounding lines in a tidewater environment, and by meltwater runoff from the surface, in roughly equal shares.
80.2 Building a three-dimensional thermomechanical ice-sheet model
Planform time-dependent modelling of ice sheets largely stems from early work by Mahaffy (1976) and Jenssen (1977), extending on the pioneering Antarctic study of Budd et al. (1971). These papers develop work by Nye (1957) on what has become known as the shallow-ice approximation (Hutter, 1983). This approximation recognizes the disparity between the vertical and horizontal length scales of ice flow, and implies grounded ice flow by simple shear. This means that the gravitational driving stress is balanced by shear stresses and that transverse and longitudinal strain rate components are neglected. Although the assumption is not valid at all places in the ice sheet, such as at the ice divide or near the ice-sheet margin (Baral et al., 2001), it has shown general applicability in large-scale ice-sheet modelling as long as surface slopes are evaluated over horizontal distances an order of magnitude greater than ice thickness. Under the shallow-ice approximation, both components of the horizontal velocity can be represented as algebraic functions of the local ice geometry (surface slope and ice thickness), which greatly simplifies the numerical solution. The model by Mahaffy (1976) was vertically integrated and was developed as a computer program to find the heights of an arbitrary ice sheet on a rectangular grid. It incorporated Glen's flow law (Glen, 1955) for polycrystalline ice deformation by dislocation creep. That is an empirical relation derived from laboratory tests in analogy with the behaviour of metals at temperatures near to their melting point, and is most commonly used in ice flow modelling. It considers ice as a non-Newtonian viscous fluid, relating strain rates to stresses raised mostly to the third power. However, in ice the rate of deformation for a given stress is also a strong function of temperature. For the range of temperatures found in natural ice masses (-50°C to 0°C), the effective viscosity changes by more than three orders of magnitude. The first model that dealt with the flow-temperature coupling in a dynamic fashion was developed by Jenssen (1977). Jenssen introduced a stretched vertical coordinate, transformed the relevant continuity and thermodynamic equations, and presented a framework to solve the system numerically.
At the heart of three-dimensional thermomechanical models is the simultaneous solution of two evolutionary equations for ice thickness and temperature, together with diagnostic representations of the ice velocity components. These express fundamental conservation laws for momentum, mass and heat, supplemented with a constitutive equation (the flow law). Plate 80.1 shows the structure of one such model as it was described in Huybrechts (1992), and further refined in Huybrechts & de Wolde (1999) and Huybrechts (2002). This model was first developed for the Antarctic ice sheet. It solves the thermomechanically coupled equations for ice flow in three subdomains, namely the grounded ice sheet, the floating ice shelf, and a stress transition zone in between at the grounding line. The flow within the three subdomains is coupled through the continuity equation for ice thickness, from which the temporal evolution of ice-sheet elevation and ice-sheet extent can be calculated by applying a flotation criterion. Grounded ice flow is assumed to result both from internal deformation and from basal sliding over the bed in those areas where the basal temperature is at the pressure melting point and a lubricating water layer is present. Ice deformation in the icesheet domain results from vertical shearing, most of which occurs near to the base. For the sliding velocity, a generalized Weertman relation is adopted, taking into account the effect of the subglacial water pressure. Ice shelves are included by iteratively solving a coupled set of elliptic equations for ice-shelf spreading in two dimensions, including the effect of lateral shearing induced by sidewalls and ice rises. At the grounding line, longitudinal stresses are taken into account in the effective stress term of the flow law. These additional stress terms are found by iteratively solving three coupled equations for depth-averaged horizontal stress deviators. Calving of ice shelves is ignored. Instead, the ice shelves extend to the edge of the numerical grid but this has little influence on the position of the grounding line. The temperature dependence of the rate factor in Glen's flow law is represented by an exponential Arrhenius equation.
The distinction between ice-sheet flow and ice-shelf flow was also made in the Antarctic model developed by Ritz et al. (2001), but they introduced an additional subdomain representing a 'dragging ice shelf' to incorporate ice-stream dynamics. In their model, inland ice is differentiated from an ice stream zone by the magnitude of basal drag. This is based on the observation that ice stream zones are characterized by low surface slopes and thus low driving stresses, but yet have fast sliding, as is the case at the Siple Coast in West Antarctica. Ritz et al. (2001) treat these zones as semi-grounded ice shelves by replacing the shallow ice approximation by the set of ice-shelf equations to which basal drag is added. Gross model behaviour turns out to be quite similar to that of the Huybrechts model, except that the West Antarctic ice sheet has a lower surface slope near the grounding line and that the break in the slope now occurs further upstream at the place where the dragging ice shelf joins the inland ice subject to the shallow ice approximation. One consequence is that grounding-line retreat in the Ritz model occurs more readily in response to rising sea levels, and that the West Antarctic sheet contains less additional ice for an expanded grounding line.
Three-dimensional flow models applied to the Greenland ice sheet have been developed along similar lines, except that ice flow is only considered in the grounded ice domain and ice shelves are not dealt with (Ritz et al., 1997; Greve, 1997a; Huybrechts & de Wolde, 1999). The position of the calving front in these models is not predicted from a self-consistent treatment of calving dynamics, as its physics are poorly understood and a convincing calving relation does not exist. Instead, these models prescribe the position of the coastline, beyond which all ice is removed as calf ice. Another characteristic of most Greenland models concerns the incorporation of variable ice fabric in the flow law to account for different ice stiffnesses as established for Holocene and Wisconsin ice in Greenland ice cores. This is achieved by prescribing a variable enhancement factor in the rate factor of the flow law, and necessitates the simultaneous calculation of ice age to track the depth of the Holocene/Wisconsin boundary. The model developed by Greve (1997a) furthermore incorporates polythermal ice and considers the possibility of a temperate basal ice layer, in which the water content and its impact on the ice viscosity are computed.
In whole ice-sheet models it is necessary to take into account the isostatic adjustment of the bedrock to the varying ice load. Early studies considered a damped return to local isostatic equilibrium (e.g. Oerlemans & van der Veen, 1984) but in most recent models the bedrock model consists of a rigid elastic plate (lithosphere) that overlies a viscous asthenosphere. This means that the isostatic compensation not only considers the local load, but integrates the contributions from more remote locations, giving rise to deviations from local isostasy. For an appropriate choice of the viscous relaxation time, this treatment produces results close to those from a sophisticated self-gravitating spherical visco-elastic earth model, while at the same time being much more efficient in terms of computational overhead (Le Meur & Huybrechts, 1996). Another feature common to most thermomechanical models is the inclusion of heat conduction in the bedrock, which gives rise to a variable geothermal heat flux at the ice-sheet base depending on the thermal history of the ice and rock.
Interaction with the atmosphere and the ocean in large-scale ice-sheet models is effectuated by prescribing the climatic input, consisting of the surface mass-balance (accumulation minus ablation), surface temperature, and, if applicable, the basal melting rate below the ice shelves. Changes in these fields are usually heavily parametrized in terms of air temperature. Precipitation rates are often based on their present distribution and perturbed in different climates according to temperature sensitivities derived from ice cores or climate models. This is principally because of a lack of a better method and implies that interaction between the pattern of precipitation and an evolving ice sheet cannot be accounted for properly. Meltwater runoff, if any, is usually obtained from the positive degree-day method (Braith-waite & Olesen, 1989; Reeh, 1991). This is an index method pro viding the bulk melting rate depending on air temperature only, but is very efficient in its use and generally gives very acceptable results (Ohmura, 2001). Models of this type are usually driven by time series of regional temperature changes (available from ice-core studies) and by the eustatic component of sea-level change, relative to present values.
Ice-sheet models are typically implemented using finite-difference techniques on a regular grid of nodes in the two horizontal dimensions, and using a stretched co-ordinate system in the vertical. Horizontal grid resolutions are mostly in the range of 20 to 50 km with between 20 and 100 layers in the vertical, concentrated towards the base where the bulk of the velocity shear takes place. Finite element implementations exist (e.g. Hulbe & MacAyeal, 1999) although these are often performed on a regular grid (Fastook & Prentice, 1994). Recent model applications have used much improved compilations of crucial input data such as bed elevation that became available on high-resolution grids from the BEDMAP (Lythe et al., 2001) and PARCA projects (Bamber et al., 2001a; Gogineni et al., 2001).
Three-dimensional models of the Antarctic and Greenland ice sheets have been used to address two main issues: the expansion and contraction of these ice sheets during the glacial-interglacial cycles, and the likely effects of greenhouse-induced polar warming.
Was this article helpful?