Skip to main content Accessibility help


  • Access



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Polythermal three-dimensional modelling of the Greenland ice sheet with varied geothermal heat flux
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Polythermal three-dimensional modelling of the Greenland ice sheet with varied geothermal heat flux
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Polythermal three-dimensional modelling of the Greenland ice sheet with varied geothermal heat flux
        Available formats
Export citation


Computations over 50 000 years into steady state with Greve’s polythermal ice-sheet model and its numerical code are performed for the Greenland ice sheet with today’s climatological input (surface temperature and accumulation function) and three values of the geothermal heat flux: (42, 54.6, 29.4) mW m−2. It is shown that through the thermomechanical coupling the geometry as well as the thermal regime, in particular that close to the bed, respond surprisingly strongly to the basal thermal heat input. The most sensitive variable is the basal temperature field, but the maximum height of the summit also varies by more than ±100m. Furthermore, some intercomparison of the model outputs with the real ice sheet is carried out, showing that the model provides reasonable results for the ice-sheet geometry as well as for the englacial temperatures.


  • A i,b Ice-covered basal area

  • A t,b Basal area covered by temperate ice

  • A(T′) Rate factor for cold ice

  • A t(ω) Rate factor for temperate ice

  • Volume flux through the CTS

  • c Specific beat of ice

  • D (ω) Water drainage function

  • E Enhancement factor

  • g Gravity acceleration

  • h z coordinate of the free surface

  • h max Maximum h of the entire ice sheet

  • H t,max Maximum thickness of temperate ice layer

  • L Latent heat of ice

  • Q Geothermal heat flux

  • T Temperature

  • T′ Homologous temperature

  • T M Pressure-melting point

  • T 0 Melting point at zero pressure

  • T CC Basal temperature at Camp Century ice core

  • T Dye3 Basal temperature at Dye 3 ice core

  • T GRIP Basal temperature at GRIP ice core

  • V tot Total ice volume

  • V temp Volume of temperate ice

  • x,y Horizontal cartesian coordinates

  • z Vertical cartesian coordinate (height a.s.l.)

  • β Clausius—Clapeyron gradient

  • κ Heal conductivity of ice

  • ρ Density of ice

  • σ Effective shear stress

  • ω Water content (mass fraction of water in temperate ice)

  • ω max Threshold water content in the drainage function

1. Introduction

In this paper we present a summary of computations that were performed with the polythermal ice-sheet model of Greve dissertation in preparation) using the numerical code SICOPOLIS that was developed for it. The ice within the Greenland ice sheet is allowed to become temperate; in the regions where it is cold the energy equation is used as an evolution equation for the temperature field; in the temperate region it serves as an equation governing the production and advection of the water content of the ice. The two regions are separated by the so-called cold temperate transition surface (CTS) whose motion is described by the kinematic equation and the dynamic conditions peculiar to this type of Stefan problem. The theory is essentially known (Fowler and Larson, 1978: Hutter, 1982: Hutter and others, 1988; Blatter, 1991; Blatter and Hutter, 1991; Hutter, 1993) and is used here in the form given by Greve (dissertation in preparation). Its implementation to a numerical code using finite-difference techniques has been performed by Greve; it makes efficient use of the σ-transformation see model-feature compilation in section 3) and stretched coordinates in the vertical, to achieve the utmost accuracy in the vicinity of the base where the ice may become temperate in rather thin layers. We shall not present any details, but in the following section will set out briefly the boundary value problem that describes the problem.

Computations will be performed for the climate-driving forces as follows. The surface temperature distribution will be that of Ohmura (1987) as parameterized by Calov (1994). Similarly, the snow balance and melting rate at the surface will be implemented according to Ohmura and Reeh’s (1991) and Reeh’s (1991) suggestions as implemented computationally by Calov (1994). Both these input quantities will be held constant. The basal surface will be held fixed (using the database of Letréguilly, and others (1990)), as only steady-state conditions are studied, and a constant vertical geothermal heat flux will be applied. This heal flux will be varied, i.e., three computations will be performed for Q = (42, 54.6, 29.4) mW m−2 with variations of ±30% about the mean.*

We shall demonstrate that the thermal state of we ice close to the base depends strongly on this variation of the geothermal heat. Due to the thermomechanical coupling, however, both the surface topography and the total volume of the Greenland ice sheet are equally affected. Such results are important findings, in particular in relation to their impact on future anthropogenic mass-balance predictions.

2. Theoretical Background

The theory of polythermal ice presented here has been formulated by Greve (dissertation in preparation) on the basis of earlier approaches provided by the authors cited in the introductory paragraph. We do not wish to give the derivation here, but mention that it is essentially based on the shallow-ice approximation (Hutter, 1983; Morland, 1984), i.e., the model equations are scaled with respect to the aspect ratio (typical height to length) of the ice sheet, and only zero-order terms are kept. Since we carry out only steady-state experiments in this paper, neither bedrock sinking nor thermal inertia of the lithosphere is taken into account. As for the stress-stretching relation. Glen’s flow law with an exponent n = 3 is used (see e.g. Paterson, 1981) as is common in ice-sheet modelling.

Below, we list those model equations that describe the temperate regions within the ice sheet. The complete set will be compiled be Greve (dissertation in preparation). For the meaning of the different quantities see the notation list.

Water content in the temperate regions:



. (2.2)

Transition conditions at the CTS:

  • (i) (melting condition):


  • (ii) (freezing condition):


3. The Three-Dimensional Polythermal Ice-Sheet Model, Sicopolis

SICOPOLIS (Simulation Code for Polythermal Ice Sheets) is a three-dimensional ice-sheet model developed by Greve (dissertation in preparation). It allows time-dependent integration of the thermomechanical model for polythermal ice sheets as described in section 2. Hereby, the ice-sheet flow is obtained as the response lo a given climatic input, namely surface temperature, accumulation and ablation. A further boundary condition is provided by the geothermal heat flux, whose influence on the results is investigated in this work.

The main new feature of this model as compared to previous three-dimensional ice-sheet models is that it accounts for polythermal conditions within the ice sheet, i.e. the possible presence of cold as well as temperate regions, in a physically adequate way. That is to say that the water content in the temperate regions is computed by solving Equation (2.1); the dependence of the rate factor A in the flow law on the water content is considered; and the positioning of the CTS is carried out by fulfilling Equations (2.3)(2.5). This entails the possibility of discontinuities of the temperature gradient and the water content at the CTS in case of freezing conditions.

The model essentially assumes that water in temperate regions travels with the same velocity as ice. However, preliminary computations have revealed that, without any water-drainage mechanism included, the water content can exceed 100% for several grid points, a physically meaningless result, of course. Therefore, it is necessary to prescribe a parameterization for water drainage. On a higher level of model complexity, this could be achieved by using something like an extended Darcy law (Hutter, 1993) for the interaction force between water and ice travelling with different vetocities; however, our formulation with one single momentum balance for the mixture ice plus water does not allow this approach. Thus, we deal with the problem by introducing a drainage function D (ω) in the water-content Equation (2.1). Based on measurements of typical water-content values in glaciers, we choose


This means that any water surplus exceeding the threshold value ωmax is assumed to be instantaneously drained into the ground. We are aware of the fact that our water-drainage model violates the local mass balance owing in water annihilation out of the interior of the temperate ice. However, global mass balance is satisfied by accounting for the amount of drained water in the computation of the evolution of the ice surface.

Summary of further model features

Three-dimensional model, based on finite differences Polythermal calculation (as described above) Horizontal resolution: 10 km

Time step: 10 a for the calculation of the velocities and the topography, 100 a for the calculation of the temperature and water content

In the vertical: σ-coordinates for each of the three regions (lithosphere: 11 equidistant grid points; temperate ice: 11 equidistant grid points: cold ice: 51 grid points with densification towards the bottom), i.e., vertical columns are mapped on [0.1| intervals

Temperature and water-content equation: implicit discretization of z derivatives, explicit discretization of x and y derivatives (upwind scheme tor the horizontal advection terms)

Height evolution equation: ADI scheme (i.e., one of the two directions is discretized implicitly in an alternating way)

Glen’s flow law with n = 3 and an enhancement factor E = 3 (accounting for the reduced stiffness of Wisconsinan ice in the near-basal regions where the main shearing takes place) used; rate factor A depending on homologous temperature in case of cold ice (common exponential law with two different activation energies Q for T′ < −10°C and T′ > −10°C, respectively; see Paterson, 1981) and on water content in case of temperate ice (At(ω) = A(T′ = 0°C) × (1 + 184ω), according to Lliboutry and Duval (1985)

Parameterization of surface temperature, accumulation and ablation as used by Calov (1994), based on Ohmura (1987), Ohmura and Reeh (1991) and Reeh (1991)

Basal sliding in case of a temperate base (Weertman-type sliding law as used by Calov (1994)): no basal sliding in ease of a cold base.

The values for the physical quantities occurring in the model are listed in Table 1.

Table 1. Values for the distinct physical quantities as used in the model calculations

4. Three Model Runs with Varied Geothermal Heat Flux

In this section we want to investigate the influence of the geothermal heat flux on modelled steady states of the Greenland ice sheet under present climate conditions. To this end, three experiments have been conducted with SICOPOLIS. In experiment rc001 the standard value = 42 mW m−2 is used, experiment rc004 is carried out with a 30% higher value ( = 54.6 mW m−2), and experiment rc005 with a 30% lower value ( = 29.4 mW m−2). These variations cover roughly the range of uncertainly of (for a detailed discussion of measurements of geothermal heat fluxes see Lee (1970)). The iterations start with the present topography as obtained from measurements by Letréguilly and others (1990) and an isothermal state with T =−10°C in the entire ice sheet; they are conducted for 50 000 model years, being sufficient to reach approximately the steady state.

In Figure (1a−c), the results of the runs rc001, rc004 and rc005 for the topography of the free surface are depicted. Figure (2) shows the model output for the homologous temperature at the ice base and the distribution of temperate ice above the base, and Table 2 gives the results for h max, V tot, V temp, Ai,b, H t.max, T GRIP, T CC and T Dye3 (see notation list), supplemented by the observed values at hand (topography data from Letréguilly and others (1990); ice-core data from Weertman (1968), Gundestrup and Hansen (1984), Johnsen and others (1992) and Dahl-Jensen (personal communication, 1994)).

Fig. 1. Steady-state free surface topography for the model runs described in the text with differing geothermal heat flux. Panels (a,b,c) are for = (42,54.6,29.4) m W m−2, all other conditions being the same. Equidistance of the level lines is 200 m.

Fig. 2. Isotherms of the steady-state homologous basal temperatures shown as solid lines, the homologous temperature being indicated in centrigrade. Open diamond symbols indicate positions when the basal ice is at the pressure-melting point, yet with no temperate layer above; full diamonds (full circles) indicate positions where there is a basal layer of temperate ice with a melting (freezing) CTS, Panels (a,b,c) are for (42,54.6,29.4) m W m−2, all other conditions being the same.

Table 2. Model output for low (rc005), standard (rc001) and high (rc004) , respectively, and observed values (obs), given for hmax, Vtot, Vtemp, Ai,b, At,b, Ht,max, TGRIP, TCC and TDye3 as defined in the notation list

Figure (1a—c) indicates that the free surface geometry depends surprisingly strongly on the amount of the geothermal heat flux. Not only do the heights of the two domes differ from each other by roughly 100 m, the location and some details in the form are also distinct. Thus, an error in the estimate of the geothermal heat flux of 30% (which is realistic) affects the total mass balance of the ice sheet in a non-negligible way. This is also evident from row 2 in Table 2 where the total ice volume V tot is listed. These results, incidentally, corroborate the opinion often expressed by theoreticians that the thermomechanical coupling is significant in quantitative ice-sheet dynamics. They show that the thermal effects feed considerably on the mechanics, and the results are qualitatively to be expected. In run rc005 (Fig.1c) the ice is coldest and thus the stiffest of the three cases, implying that the ice sheet must be the thickest. In run rc004 (Fig. 1b) the heat input into the ice from below is largest, the bottom-most ice weakest, horizontal advection largest and height arid total volume therefore smallest. Note that the total ice-covered surface is nearly insensitive to variations of the geothermal heat flux.

Even more affected by the variation of the geothermal heat flux is the thermal regime of the ice sheet, in particular its response close to the base, as presented in Figure (2a) c and in rows 3 and 5 of Table 2. In all three cases of our model run, quite a broad belt of the marginal basal ice is temperate. This basal ring of temperate ice is interrupted at a few broad stripes in the east, north and west of the Greenland ice sheet. In the three cases of Figure (2a—c) the temperate basal surfaces amount to 39%, 51% and 33% of the entire ice-covered basal surface and are therefore surprisingly large. This result qualitatively verifies earlier results by (Calov and Hutter 1994a,b) who demonstrated using a much simpler model that the thermal conditions of the basal ice depend critically on the thermal input from the surrounding boundaries.

Figure (2) also shows the positions in the model runs where there occurs a non-vanishing basal layer of temperate ice with either melting or freezing CTS. It can be seen that a freezing CTS with the possibility of discontinuities of the temperature gradient and water content always occurs at positions near the ice margin, in agreement with what is expected. The associated volume of the temperate ice is listed in Table 2 and is of order 3000 km3 (0.1% of the total ice volume with varations of up to 10%.

Table 2 gives some more results of the simulations as an overview, supplemented by observations. Obviously, the thermal regime at the base, namely the basal area covered by temperate ice, A t,b, and the basal temperatures at the ice-core sites, T GRIP/CC/Dye3, respond most critically in the different values of the geothermal heat flux. However, there is a notable sensitivity of h max and V tot as well, i.e., the overall geometry of the ice sheet is affected by the different values in a non-negligible way, as stated above. The ice-covered basal area, A i,b, appears to be almost insensitive.

The agreement between the modelled ice sheets and the real world is remarkably good. With the exception of T Dye3, the observed values of the quantities listed in Table 2 lie within or at least very close to the range of the modelled ones; the rather big deviation for T Dye3 is supposedly due to the rather small number of grid points covering the southern part of the ice sheet. However, one must bear in mind that the actual ice sheet is certainly not in steady state due to its experienced climate history. This affects much less the mass balance, whose response time is rather small, than it does the thermal regime. A transient SICOPOLIS run driven with the standard and a climate history as suggested bt the GRIP ice core provides an ice sheet of comparable geometry (somewhat thicker), yet with notably lower englacial temperatures, and thus with a V temp reduced by about 18% and an A t,b reduced by about 15% compared to the steady-state results (Greve, dissertation in preparation), due to the memory of the Wisconsinan ice age. In agreement with these considerations the modelled steady-state basal ice-core temperatures tend to be higher than the observed ones (except T GRIP from runs rc005 and rc001).

Unfortunately it is difficult to check the model output for the quantities describing temperate ice, namely V temp, A t,b and H t,max, against the real world because of data sparsity. First, the three deep ice cores mentioned above are all situated in regions where no temperate ice is found in the entire vertical column, in agreement with the model results for all three runs presented here. Secondly, boreholes drilled at Jakobshavns Isbræ (West Greenland, south of Disco Island) have revealed a temperate ice base and suggested the probable existence of an overlying temperate layer (Iken and others, 1993). This is confirmed by the three model runs (see Fig. 2), so that the model output in terms of polythermal features is at least promising, even though quantitative checks cannot be carried out at the moment. Remote-sensing techniques based on the measurement of radio-wave velocities (Macheret and others, 1993) may remedy ibis lack of data in the future.


While performing this work R. Greve was supported by the Studienstiftung des Deutschen Volkes and the Department of Geophysical Sciences. University of Chicago.

* Amazingly, the numerical value 42 of the standard represents exactly the Answer to the Ultimate Question of Life, the Universe and Everything (Adams, 1979). Further study will be required to figure out whether this correspondence can lead to a more profound understanding of all these things or is just an accident.

These temperatures are calculated by averaging the modelled values at the four adjacent grid points, weighed by the inverse distances of the grid-point positions to the actual borehole posilions.


Adams, D. 1979. The hitch hiker’s guide to galaxy. London, Pan Books Ltd.
Blatter, H. 1991. Effect of climate on the cryosphere. Climate conditions and the potythermal structure of glaciers. Zürcher Geogr. Schr., 41, 1-98.
Blatter, H. and Hutter, K., 1991. Polythermal conditions in Arctic glaciers. J. Glaciol., 37 (126), 261-269.
Calov, R. 1994. Das thermomechanische Verhalten des grönländischen Eisschildes unter der Wirkung verschiedener Elimaszenarien Antworten eines theoretisch-numerisichen Modells. (Dissertation, Technische Hochschule, Darmstadt.)
Calov, K. and Hutter., K. In press a. The thermomechanical response of the Greenland ice sheet to various climate scenarios. Part I: Analysis of the temperature regime. Climate Dynamics.
Calov, K, and Hutter., K. In press b. The thermomechanical response of the Greenland ice sheet to various climate scenarios. Part II: Analysis of flow field and geometry. Climate Dynamics.
Fowler, A. C. and Larson, D.A., 1978. On the flow of polythermal glaciers. I. Model and preliminary analysis. Proc. R. Soc. London, Ser. A, 363 (1713), 217-242,
Gundestrup, N. S. and Hansen, B.L., 1984. Bore-hole survey at Dye 3, south Greenland. J. Glaciol., 30 (106), 282-288.
Hutter, K. 1982. A mathematical model of polythermal glaciers and ice sheets, Geophys. Astrophys, Fluid Dyn., 21 (3-4), 201-224.
Hutter, K. 1983. Theoretical glaciology; material science of ice and the mechanics of glaciers and ice sheets. Dordrecht, etc., D. Reidel publishing Co./Tokyo, Terra Scientific Publishing Co.
Hutter, K. 1993. Thermo-mechanically coupled ice sheet response—cold, polythermal, temperate. J. Glaciol., 39 (131), 65-86.
Hutter, K., Blatter, H. and Funk, M., 1988. A model compulation of moisture content in polythermal glaciers. J. Geophys. Res., 93 (B10), 12,205-12.214.
Iken., A., Echelmeyer, K., Harrison, W. and Funk., M. 1993, Mechanisms of fast flow in Jocobshavns Isbrae, West Greenland: Part I. Measurements of temperature and water level in deep boreholes. J. Glaciol., 39(131), 15-25.
Johnsen, S.J. and 9 others, 1992. Irregular glacial Interstadials recorded in a new Greenland ice core. Nature, 359(6383), 311-313.
Lee, W. H. K. 1970. On the global variations of terrestrial heat flow. Phys. Earth Planet. Inter., 2(5), 332-341.
Letréguilly, A., Reeh, M. and Huybrechts., P. 1990, Topographical date for Greenland, Report. Bremerhaven, Alfred-Wegener-Institut für Polarforschung.
Lliboutry, L. and Duval., P. 1985. Various isotropic and anisotropic ices found in glaciers and polar ice caps and their corresponding rheologies. Annales Geophysicae, 3(2), 207-224.
Macheret, Yu. Ya., Moskalevsky, M. Yu. and Vasilenko, E.V. 1993. Velocity of radio waves in glaciers as an indicator of their hydrothermal state, structure and regime. J. Glaciol., 39 (132), 373-384.
Morland, L. W. 1984. Thermo-mechanical balances of ice sheet flow, Geophys. Astrophys. Fluid Dyn., 29, 237 266.
Ohmura, A. 1987. New temperature distribution maps for Greenland, Z. Gletscherkd. Glazialgeol., 23(1), 1-45.
Ohmura, A. and Reeh., N. 1991. New predpitation and accumulation maps for Greenland. J. Glaciol., 37(125), 141-148.
Paterson, W.S.B. 1981. The physics of glaciers. Second edition. Oxford, etc., Pergamon Press.
Reeh, N. 1991. Parameterization of melt rate and surface temperature on the Greenland ice sheet, Polarforschung, 59(3), 113-128.
Weertman, J, 1968. Comparison between measured and theoretical temperature profiles of the Camp Century, Greenland, borehole. J. Geophys, Res., 79(8), 2691-2700.