Introduction

People have interacted with glaciers since historical times; either the glaciers posed a threat to settlements or they hindered mobility in mountainous terrain. Although glaciers seem to be static phenomena at first glance, their waxing and waning over long time periods was soon recognized. The understanding of the physics of glaciers changed from intuitive interpretations of individual observations on glaciers in the late 18th century to peculiar ideas about extrusion flow as late as the mid-20th century.

Shortly after the inception of the Journal of Glaciology in January 1947, physicists started to enter the field of glacier mechanics and thermodynamics. These activities initiated a new chapter in physical and theoretical glaciology. With the availability of increasing computational power, the large body of knowledge that has developed since those pioneering days can now be condensed into ever more sophisticated numerical models of the many aspects of the thermomechanics, dynamics and stability of glaciers and ice sheets.

The Beginning

Horace-Benedict de Saussure (1740-99) observed sliding of glaciers and related the sliding motion to the action of water flowing at the bottom (Reference ClarkeClarke, 1987, and references therein). With this, de Saussure may have been the first to hypothesize a connection between glacier motion (sliding) and the occurrence of water at its base.

Bernhard Friedrich Kuhn (1762-1825) was among the early researchers who tried to understand the working of glaciers (Rothlisberger, 1987). He was probably the first to combine thermodynamics and mechanics to explain the motion of glaciers. Although our present knowledge indicates it is not very realistic, it is interesting to look at Reference KuhnKuhn’s (1787) theory in some detail, partly in his own words, albeit translated from his old-fashioned German. He was already aware of some seasonal variation of the mechanical activities in a glacier:

• At the onset of the cold season it becomes slowly morequiet; the vaults close themselves and the ice sinks back to the ground. Also the thundering and banging, which is so striking for the visitors, becomes more seldom, and is no longer heard in the middle of the winter.

Kuhn experienced the ‘thundering’ and ‘banging’ in the summer on the glaciers of Grindelwald in the Bernese Alps, and realized that

• from 1770 through 1778 both Grindelwald glaciers increased with fast paces; and each time at the end of the beautiful season, that they approached the surrounding scattered objects by 20-30 and more Klafter (˜1.8m), and most astonishing, their expansion was larger during summer than during winter.

He related the motion of the glacier to the heat collected in summer by the rocks surrounding the glaciers and subsequently flowing underneath the ice:

• Much more it seems that its true cause lies in the heating by the sun of the ground in the neighbourhood of the glaciers. The fire material, set on motion by the approaching (of the sun), spreads to all sides, and flows through the surface of the ground towards the glaciers …

• The heat acts as long on the base of the glaciers until the entire weight of the ice masses only rests on pillars here and there. Even these will be destroyed by the continuously working heat flux. If in the end they are no longer able to carry the load, then necessarily the masses must collapse. Now, the glaciers rest always more or less on slopes of the mountains or on steep sloping steps of the valleys; if then the supports give way to the weight, the mass, together with the sinking, moves forwards as far down the slope until an equilibrium between force and resistance is reached again.

Kuhn’s understanding of the motion of the ice in glaciers may have been the weakest (Rothlisberger, 1987) part of his understanding of glaciers in general. Even at that time, he realized that the glaciers had once been much larger and had reached far down into the valleys. He may have been one of the discoverers of ice-age glaciations in the Alps.

Johann Jakob Scheuchzer (1672-1733), Johann von Charpentier (1786-1855) and Louis Agassiz (1807-73) put forward the ‘dilatation theory’ to explain the motion of glacier ice by the forces exerted by the expansion of water freezing in fissures and capillary tubes in the ice. James David Forbes (1809-68) opposed this theory, based on his own observations on Mer de Glace, France, where he observed the ice to flow fastest in the centre of the glacier. He concluded that, if the dilatation theory were true, then the flow would be greatest at sunset and near the margins. He was perhaps the first to propose viscous flow for glacier ice. (For more details see Reference ClarkeClarke, 1987, and references therein.) Systematic observations of the motion of glaciers were conducted during the 19th century on Unteraar- gletscher and Rhonegletscher, Switzerland. Several coloured stone rows across Rhonegletscher were surveyed several times after 1894 (Reference MercantonMercanton, 1916), and visibly demonstrated the differential motion of the ice across the glacier.

Reference DemorestDemorest (1942) recognized several modes of motion of the ice by distinguishing

1. … discontinuous fracture and shearing in which differential movements occur along discrete planes, whereas at depth the ice is plastic and motion is by plastic flow of a streamline nature, …

and further divided streamline flow into pressure-controlled extrusion and drainage-controlled gravity flow (Reference WaddingtonWaddington, 2010).

Physical Glaciology

In the mid-20th century, investigations began to discover the flow properties of the ice. The laboratory experiments of Reference GlenGlen (1952) resulted in the flow law, already anticipated by Reference PerutzPerutz (1950a, Reference Perutzb),

where is the shear rate, t the shear stress and B and n are constants. Measurements of the closure of boreholes and tunnels and subsequent theoretical analyses calibrated and confirmed the power-law formulation. Reference NyeNye (1953) applied the flow law to field observations in the Jungfraufirn, Swiss Alps (Reference Gerrard, Perutz and RochGerrard and others, 1952), and the observed contraction rates of ice tunnels on Skauthore glacier, Norway (Reference McCallMcCall, 1952), and Z’Mutt and Arolla glaciers, Swiss Alps (Reference HaefeliHaefeli, 1951, Reference Haefeli1952; Reference Haefeli and KasserHaefeli and Kasser, 1951).

Reference NyeNye (1953) computed an exact solution for the stress/ strain-rate relation for the closure of cylindrical and spherical holes in glacier ice. Applying the measurements of the tunnels he obtained values of

These studies were the turning point for the increasing acceptance of the idea that glacier ice behaves like a quasiviscous fluid, much like metals and other polycrystalline materials. Thus, the understanding of glacier flow became a problem in fluid mechanics with a non-Newtonian material.

In the 1950s, the rigorous mathematical descriptions of physics started to enter glaciology. This can be seen from the ‘Bader peak’, the first ‘Nye peak’ and the ‘Robin and Rothlisberger peak’, presented by Garry Reference ClarkeClarke (1987) in a memorable talk, ‘A short history of scientific investigations on glaciers’, given in the special session on the occasion of the 50th anniversary of the International Glaciological Society during the Second Symposium on Remote Sensing in Glaciology in 1986, held at the University of Cambridge, UK. The peaks mark large peaks in the number of equations per 100 pages of the Journal of Glaciology, compared to the relatively steady number of maps (Fig. 1). The 1950s also saw the onset of papers in the Journal of Glaciology combining maps and equations (e.g. Reference HaefeliHaefeli, 1952; Reference RöthlisbergerRcithlis-berger, 1955).

Fig. 1. Maps and equations in the Journal of Glaciology. The peaks are labelled to indicate the main contributors to each peak: B (Reference BaderBader, 1950), N1 (Reference NyeNye, 1952a,b), RR (Reference RobinRobin, 1955; Reference RöthlisbergerRo¨thlisberger, 1955), N2 (Reference NyeNye, 1959a,c) and HL (Reference Hutter, Legerer and SpringHutter and others, 1981; Reference Lliboutry and ReynaudLliboutry and Reynaud, 1981). (From Reference ClarkeClarke, 1987.)

The early attempts to describe glaciers in terms of mathematical physics were essentially guided by the need for simplifications. The style of these early works was influenced by the highly developed analytical methods used in contemporary fluid mechanics. The mapping of the complexities of the natural systems considered was a necessary step of the preceding research. However, only the simplest versions of the systems lead to mathematical problems that can be treated analytically or that yield simple numerical solutions. It is only in recent times that computing power and, triggered by it, experience in numerical mathematics allows us to attack more and more complex problems in glacier dynamics. In this sense, the 1950s saw the pioneers who converted glaciology into a physical science, in particular in the treatment of glacier flow (Reference NyeNye, 1951, Reference Nye1952b, Reference Nye1959c) and thermodynamics (Reference RobinRobin, 1955).

Reference NyeNye (1951, Reference Nye1952b) first described the fundamental patterns of the flow of ice in glaciers and ice sheets. He exploited symmetries in specific examples of geometries to obtain solutions of the corresponding differential equations. However, these examples already illustrated the most important patterns of ice flow that occur in valley glaciers and ice sheets. The parallel-sided slab (thin-layer) solution for the difference between the ice velocities u₀ at the ice surface and u at depth d, measured perpendicular to the surface of the slab, can be given analytically,

where (in Nye’s terminology) K = (ρg/B)ⁿ , B and n are constants of the flow relation, and ρ, g and α are the ice density, the acceleration due to gravity and the inclination angle of the slab, respectively. Exploiting the symmetries of an infinitely deep channel of width d yields the same basic equations as for ice flow in a slab, except that the coordinate axes are exchanged.

Another example that approximates a valley glacier even more realistically is a channel defined by a half-circular cylinder. Here again, the corresponding solution of the equations corresponds to the solution for the parallel slab, Equation (3), except for an additional factor (1 /2)ⁿ on the right-hand side, and d now denotes the radial distance from the circular axis (Reference NyeNye, 1952b). However, the more detailed explanation was only given in Nye’s (1965a) paper on the flow of ice in channels of various shapes. This solution exploits the circular symmetry by showing that the Hagen- Poiseuille flow in a circular tube is described by the same equations, if the pressure gradient along the tube is given by ρg sin α.

For some more general and more realistic cases, the plasticity assumption n → ∞ was used to allow simpler mathematical treatments (Reference NyeNye, 1951). The stress and velocity solution for slab flow of uniform slope was used as a first approximation for slab flow with a slowly varying slope. For a steady geometry, the vertical component of the velocity, V, the radius of curvature of the bed, R, and the longitudinal gradient of the flow rate, 0, are closely related by (Reference NyeNye, 1951)

where ϕ^′ denotes the longitudinal derivative of ϕ and α is the local inclination of the slab. With these considerations, Reference NyeNye (1951) derived the connections between active (compressive) and passive (extending) flow and the trajectories of maximum shear stress, along which a possible shear fracture or fault could occur. Based on these considerations, Nye (1952b) distinguished the three different possible crevassing patterns in a glacier flowing in a channel (Fig. 2), depending on whether the flow is compressive, neutral or tensile.

Fig. 2. The unbroken lines show the theoretical positions and directions of crevasses in three possible cases. The diagrams at the top indicate the stresses acting near the margin. (From Reference NyeNye, 1952b.)

In the discussion of the mechanical properties of ice, Nye (1952b) mentioned that

• the values of the constants depend rather sensitively on the amount of bubbliness in the ice and on the temperature.

Thus it became necessary to study the englacial temperatures to fully treat the dynamics of glaciers and ice sheets. Reference RobinRobin (1955) first attempted to estimate the temperature distribution in an ice sheet. He began by solving the Fourier equation for a vertical profile in the centre of an ice sheet, assuming the vertical velocity component grew proportionally to the distance to the bed, matching the accumulation rate at the surface. To estimate temperature profiles at some distance from the centre, several simplifying assumptions, similar to the modern shallow-ice assumptions, were made (e.g. that the horizontal heat diffusion is negligible). The influence of the increasing surface temperature conditions with distance to the centre on such a vertical profile was transformed into a vertical temperature gradient at the surface. Furthermore, no horizontal advection or diffusion was considered explicitly, but horizontal advection and diffusion were mimicked by a constant heat source along the vertical profiles (see Equation (15) of Reference RobinRobin, 1955). Although the match with observed data from Antarctica was not particularly good, the principal patterns of the temperature distribution in an ice sheet were correctly reproduced (Fig. 3), in particular, the core of cold ice advection and the possibility of temperate ice at the bed in the marginal area.

Fig. 3. Sketch of the temperature distribution in a simplified ice sheet similar to central Greenland. (From Reference RobinRobin, 1955.)

The simplification that the flow of glaciers can be described as quasi-stationary was justified by Reference NyeNye (1951):

1. The rate of deformation in a glacier or an ice-sheet is so slow, and the times involved are so long, that the deformation will be almost entirely of the type which is called in metals ‘quasi-viscous creep’ or K-flow’ (the transient component of creep, or β-flow, will be negligible).

Thus, the problem reduces to a Stokes problem (Equations (1) and (2) of Reference NyeNye, 1951). He also assumed incompressibility (his Equation (4)) and isotropy (Equation (5)). Reference NyeNye (1953) additionally presented the tensorial formulation of the isotropic flow law of the Glen type. Together with the need to know the temperature field to obtain the temperature- dependent viscosity, the basic formulation of the thermomechanically coupled flow of the ice in glaciers was established at that time. In modern notation the equations for the flow and stress fields are

and the thermodynamic equation for the evolving interior temperature field is

where T, T ^d and D are the stress, deviatoric-stress and strain-rate tensors, respectively, u is the velocity vector, A the temperature-dependent rate factor, the function f(τ) defines the stress/strain-rate relation, O is the total temporal derivative of the temperature O, k is the heat conductivity, c is the specific heat and Q is the density of the heat production rate. The above equations for interior fields must be complemented by an equation for the evolution of the ice thickness, h, in a given climate (surface mass balance, a _s, and basal mass balance, a _b),

where t is the time and q the horizontal volume flux vector.

Analytical (exact) solutions were the preferred solutions for mathematical problems in the pre-computer era, but even today they are sometimes considered superior to numerical solutions. However, the above set of equations (Equations (5-7)) allows for exact solutions only for a small set of simple geometries and physical conditions. In the appendix of Nye’s (1952b) paper an analytical solution of the shape of a circular ice sheet resting on a horizontal bed for perfect plastic ice was given,

where h is the ice thickness at the position x and h ₀ is the maximum thickness at the centre of the ice sheet. Using this solution, Nye (1952c) tried to estimate the thickness of the Greenland ice sheet. Reference VialovVialov (1958) presented an exact isothermal steady-state solution for constant accumulation and for ice with Glen-type flow properties, f(t)=t ^n-1. Solutions with more realistic accumulation functions and time-dependent solutions, both isothermal and thermomechanically coupled, were proposed by Reference HalfarHalfar (1981, Reference Halfar1983) and Reference Bueler, Lingle, Kallen-Brown, Covey and BowmanBueler and others (2005, Reference Bueler, Brown and Lingle2007). The early motivation to find solutions of the ice-sheet equation was the wish to estimate the thickness of existing ice sheets (Reference NyeNye, 1952c; Reference WeertmanWeertman, 1961), whereas later the motivation was the need for exact solutions to verify the numerical methods of ice-sheet models.

The topic of exact solutions of the ice-sheet equation demonstrates a pattern in mathematical glaciology: if one searches the literature for some specific topic in glacier physics, one finds a pioneering paper at the very beginning written by J.F. Nye. A few selected examples are (1) the response of glaciers to climatic changes (Reference NyeNye, 1958, Reference Nye1960), where the kinematic wave theory was also introduced, (2) papers related to the interpretation of field measurements, such as the determination of strain-rate tensors at the surface of glaciers (Reference NyeNye, 1959b) and (3) the correction of measurement of past firn accumulation in ice sheets from ice cores (Reference NyeNye, 1963). Furthermore, the first calculation of ice flow in a circular channel (Reference NyeNye, 1952b) was extended to the flow in channels of various shapes in Nye’s seminal paper (Reference NyeNye, 1965a), the influence of longitudinal stress on glacier flow was briefly mentioned by Reference NyeNye (1952b), then extended (Reference NyeNye, 1969b), and its influence on sliding over wavy beds was also discussed by (1969 Reference NyeNye, 1970).

Reference Jenssen and RadokJenssen and Radok (1961) used a CSIRAC electronic computer to compute transient vertical temperature profiles in ice sheets. In Reference NyeNye’s (1965a) paper on channel flow and his paper on the frequency response of glaciers (Reference NyeNye, 1965b), results were obtained by numerical solution using a digital computer (IBM 1620). Despite the limited speed of computation at that time, finite-difference solutions with damping of the oscillations and acceleration of the convergence with the alternate direction method and over-relaxation could be obtained within a reasonable time (hours to days). These papers mark the beginning of the new field of computational glaciology, addressed in the following section.

Computational Glaciology

With the dawn of digital computing, scientists quickly started to use the new tools to obtain numerical solutions to previously intractable problems. The basic problem of ice- sheet and glacier modelling is to compute the ice flow in a moving domain with additional internal evolving fields (Reference Hindmarsh and HutterHindmarsh and Hutter, 1988), where the domain itself is a result of the ice flow and the surface mass balance.

Owing to the limited computing power available, the first numerical models of ice sheets and glaciers were restricted to far-reaching approximations. Reference Campbell and RasmussenCampbell and Rasmussen (1969, Reference Campbell and Rasmussen1970) and Reference Rasmussen and CampbellRasmussen and Campbell (1973) presented a vertically integrated, hydrostatic glacier model with various flow laws that parameterize basal stresses and velocities. Reference Budd and JenssenBudd and Jenssen (1975) formulated the threedimensional (3-D) model equations of the dynamics of glacier systems and studied transient evolution of the dynamic into equilibrium states in two-dimensional (2-D) simplifications. The first 3-D ice-sheet model was developed by Reference MahaffyMahaffy (1976) and applied to Barnes Ice Cap (located in central Baffin Island, Nunavut, Canada). The model was based on a Glen-type rheology and the shallow-ice approximation (SIA) of the force balance, a consistent approximation that was rigorously derived later by Reference HutterHutter (1983) and Reference MorlandMorland (1984). A major simplification was that the strong temperature dependence of the ice viscosity (contained in the rate factor A; Equation (5)) was not considered, which allows the governing equations to be vertically integrated and the temperature-evolution equation (Equation (6)) to be ignored.

Reference JenssenJenssen (1977) applied the first thermomechanically coupled model to the Greenland ice sheet. However, the resolution of only 12 × 12 × 10 (vertical) gridpoints was too coarse to resolve the detailed features of the geometry of the ice sheet, and computations had to be interrupted after 1000years integration time, due to numerical problems. In this work, he introduced the terrain-following coordinates that are used in most subsequently developed ice-sheet models. Similar models by the same working group were applied to the North American and Antarctic ice sheets (Reference Budd and JenssenBudd and Smith, 1981, Reference Budd and Smith1982;Reference Budd, Jenssen and SmithBudd and others, 1984 and later publications). A different approach was taken in a series of studies in the early 1980s by (1982a,b,Reference OerlemansOerlemans c, 1983). The numerical solution of the thermomechanically coupled problem of ice-sheet dynamics was made more efficient by using a spectral approach, with polynomial basis functions for the vertical temperature profile and the velocity profile (Reference OerlemansOerlemans, 1982a). Bedrock adjustment due to glacial isostasy was accounted for, and the Antarctica version of the model even included a simple treatment of ice shelves.

Based on Oerlemans’ work, the first comprehensive, 3-D, thermomechanically coupled ice-sheet model was developed by Reference Huybrechts and OerlemansHuybrechts and Oerlemans (1988) and Reference HuybrechtsHuybrechts (1990a, Reference Huybrechtsb). The version for the Antarctic ice sheet comprised coupled ice shelves in the shallow-shelf approximation (SSA; Reference Morland, Van der Veen and OerlemansMorland, 1987), a sheet/shelf transition model including membrane stress gradients, basal sliding and isostatic bed adjustment. A relatively fine horizontal resolution of 40 km and 10 layers for the scaled vertical coordinate was applied, which required the model to run on contemporary supercomputers. A Greenland version of the model was first presented by Letreguilly and others (1991a,b).

At approximately the same time, Reference HerterichHerterich (1988, Reference Herterich1990) and Bohmer and Herterich (1990) came up with a model of similar complexity. Like the Huybrechts model, the Antarctica version of the Herterich model included coupled ice- sheet/ice-shelf dynamics based on the SIA and SSA, respectively, and had an explicit treatment of the transition zone. Major differences were that the Herterich model did not use a scaled vertical coordinate, did not account for basal sliding and ran at a coarser horizontal resolution (100 km). This model was also used to investigate the possible existence of an extended glaciation in Tibet during the last glacial period (Reference Kuhle, Herterich and CalovKuhle and others, 1989) and, later and complemented by basal sliding, to simulate the response of the Greenland ice sheet to various climate scenarios (Reference Calov and HutterCalov and Hutter, 1996).

During the 1990s, several other 3-D, thermomechanically coupled ice-sheet models entered the scene. In the first intercomparison topic of the European Ice-Sheet Modelling Initiative (EISMINT; Reference Huybrechts and PayneHuybrechts and others, 1996) five models of that kind (plus several simpler ones) participated, while in the second topic (Reference PaynePayne and others, 2000, and references therein) the number was already up to 10. The latter topic focused, in particular, on the effects of thermomechanical coupling, and revealed that the radial symmetry implied in the experimental design can be broken by the formation of radial spokes of the velocity, ice- thickness and temperature distributions. While this feature appeared to be common to all participating models, details varied between them, and the preferential alignment of the spokes with the lines of symmetry of the numerical grid pointed, at least partly, to a numerical cause for the phenomenon. Nevertheless, it was also concluded that the models showed considerable agreement in their predictions of global-scale response to imposed climate change.

None of the ice-sheet models mentioned above used a force balance beyond the SIA. However, for realistic simulations of smaller glaciers, which are less shallow than ice sheets and have a more complex topography, the SIA is inadequate, and thus during the past ˜15 years 3-D glacier models that solve the full-Stokes problem of ice flow have been developed (Reference GudmundssonGudmundsson, 1999; Reference Luthi and FunkLuthi and Funk, 2000;Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007;Reference Jouvet, Picasso, Rappaz and BlatterJouvet and others, 2008, Reference Jouvet, Huss, Blatter, Picasso and Rappaz2009). This is computationally much more demanding than solving the SIA and requires modern, powerful computers.

Integrating higher-order or full-Stokes dynamics in ice- sheet models is even more of a challenge due to the size of ice sheets. Reference Saito, Abe-Ouchi and BlatterSaito and others (2003, Reference Saito, A. and Blatter2006) implemented the first-order approximation (FOA) by Reference BlatterBlatter (1995) (see also Reference Greve and BlatterGreve and Blatter, 2009) and investigated the impact (compared with the SIA) on the solutions of the above- mentioned EISMINT experiments by Reference PaynePayne and others (2000). A similar model was devised by Reference PaynePattyn (2003), later upgraded to full-Stokes dynamics and applied to problems of ice flow over subglacial lakes (Reference PattynPattyn, 2008). However, a fully operational, 3-D, thermomechanically coupled, prognostic ice-sheet model with higher-order or full-Stokes dynamics is not yet available.

Efforts towards higher-order and full-Stokes modelling culminated in the Higher-Order Model (HOM) intercomparison topic of the Ice-Sheet Model Intercomparison Project (ISMIP) (Reference PattynPattyn and others, 2008), in which 27 numerical models of different complexity and one analytical model from 20 contributors participated. These models were compared and verified in a series of experiments applied to 2- and 3-D geometries of both steady-state (diagnostic) and time-dependent (prognostic) type. It was found that all models produced results that are in approximate agreement, even for high aspect ratios when the SIA is clearly violated. The full-Stokes models were most consistent with each other and an analytical solution that exists for one of the experiments, while the spread among the different flavours of higher-order models was somewhat larger.

Most recently, in the Heinrich Event INtercOmparison (HEINO) topic of ISMIP, internal large-scale ice-sheet instabilities in different contemporary ice-sheet models were explored (Reference CalovCalov and others, 2010). For the experiments a simplified geometry that reproduces the main characteristics of the Laurentide ice sheet (including the sedimented region over Hudson Bay and Hudson Strait) and a temporally constant glacial climate were employed. It was found that all participating models (eight SIA models and one combined SIA/SSA model) are capable of producing Heinrich-type free oscillations if the boundary conditions are sufficiently favourable. However, the large differences between the results of different models (some are much more prone to produce oscillations than others) clearly showed that further model improvements are crucial for adequate, robust and reliable simulations of ice-sheet instabilities.

Three operational ice-sheet models, which all took part in ISMIP HEINO, are currently available on the internet as free software:

• Glimmer-CISM (Community Ice-Sheet Model), http://glimmer-cism.berlios.de/ (Reference Rutt, Hagdorn, Hulton and PayneRutt and others, 2009);

• PISM (Parallel Ice-Sheet Model), http://www.pism-docs.org/ (Reference Bueler, Brown and LingleBueler and others, 2007);

• SICOPOLIS (SImulation COde for POLythermal Ice Sheets), http://sicopolis.greveweb.net/ (Reference GreveGreve, 1997).

These models include 3-D coupling of temperature and velocity fields in the SIA. They differ in the numerical grids and some numerical schemes, but also allow for various physical processes such as polythermal conditions and sliding parameterization. SICOPOLIS (Reference GreveGreve, 1997) includes a polythermal scheme that solves the thermodynamics equation for temperature in the cold and for moisture content in the temperate domain. The cold/temperate transition surface (CTS) is tracked by matching the thermal and kinematic conditions along the boundary. PISM offers a polythermal scheme that solves the thermodynamics equation for enthalpy, and the CTS is given by a defined value of the enthalpy (Reference Aschwanden and BlatterAschwanden and Blatter, 2005). PISM also offers a novel scheme for sliding including the membrane stresses at the bed by balancing the driving forces at the bed of the SIA ice sheet and the resistive forces of a subglacial deforming layer in a shallow-shelf-type membrane in between (Reference Bueler and BrownBueler and Brown, 2009). Glimmer-CISM is particularly designed to be interfaced to a range of global climate models. All three models are currently under rapid development, and at different stages towards the implementation of higher-order dynamics, ice-shelf/stream dynamics and basal hydrology. Figure 4 illustrates schematically the various components of the models.

Fig. 4. Scheme of a modern, dynamic/thermodynamic ice-sheet model. The rectangular boxes correspond to prognostic model components, the oval boxes to input quantities.

Prospects

Despite the general opinion that waiting for faster computers is not good advice in numerical model development, the fast progress in computing power has helped in glacier modelling. The development of more and more efficient and stable numerical algorithms, combined with the ever-growing availability of storage space and advances in computation speed, enabled us to model ever increasingly complex glacial systems at increasing spatial and temporal resolution. In modelling geophysical systems (and perhaps any type of physical system) the wishes always aim beyond the capacity of the computers.

Computation power and efficiency of numerical methods have progressed to the point where the models are likely to be able to include full-Stokes solutions for entire ice sheets with adaptive grids to obtain high spatial resolution where required. Currently, at least three projects have started to develop such next-generation ice-sheet models,

• http://www.elmerfem.org/wiki/index.php/Elmer_Ice_Sheet_modeling

• http://issm.jpl.nasa.gov/

• http://trac.mcs.anl.gov/projects/sisiphus

The development of this type of ice-sheet model lies beyond the problem of thermomechanics of ice masses and requires close collaboration between specialist physicists, mathematicians and computer scientists.

In computational glaciology, increasing expectations of solutions to pressing problems, such as the future of ice sheets and glaciers in a warming climate and especially the possibility of catastrophic disintegration of large ice shelves in Antarctica, have led to the present surge in ice-sheet and glacier model development. Petascale computation power, adaptive finite-element meshes and efficient preconditioning may enable the computation of the systems behaviour to reasonable resolution and accuracy. However, true predictive power also rests on sufficient knowledge of the boundary conditions and the physical mechanisms. Thus, a combination of precisely targeted (even then expensive) field observations, advances in knowledge of the physical processes involved and advancing all aspects of numerical modelling technology such as algorithms, coding and computation power will remain the tasks of glaciology in the future.