Notation

• 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 ₀ 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 (Reference Fowler and LarsonFowler and Larson, 1978: Reference HutterHutter, 1982: Reference Hutter, Blatter and FunkHutter and others, 1988; Reference BlatterBlatter, 1991; Reference Blatter and HutterBlatter and Hutter, 1991; Reference HutterHutter, 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 Reference OhmuraOhmura (1987) as parameterized by Reference CalovCalov (1994). Similarly, the snow balance and melting rate at the surface will be implemented according to Reference Ohmura and Reeh.Ohmura and Reeh’s (1991) and Reference ReehReeh’s (1991) suggestions as implemented computationally by Reference CalovCalov (1994). Both these input quantities will be held constant. The basal surface will be held fixed (using the database of Reference Letréguilly, Reeh and Huybrechts.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⁻² with variations of ±30% about the mean.Footnote *

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 (Reference HutterHutter, 1983; Reference MorlandMorland, 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. Reference PatersonPaterson, 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.1)

with

. (2.2)

Transition conditions at the CTS:

(2.3)

• (i) (melting condition):

  

  (2.4)

• (ii) (freezing condition):

  

  (2.5)

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 (Reference HutterHutter, 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

(3.6)

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 Reference PatersonPaterson, 1981) and on water content in case of temperate ice (A_t(ω) = A(T′ = 0°C) × (1 + 184ω), according to Reference Lliboutry and Duval.Lliboutry and Duval (1985)

Parameterization of surface temperature, accumulation and ablation as used by Reference CalovCalov (1994), based on Reference OhmuraOhmura (1987), Reference Ohmura and Reeh.Ohmura and Reeh (1991) and Reference ReehReeh (1991)

Basal sliding in case of a temperate base (Weertman-type sliding law as used by Reference CalovCalov (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⁻² is used, experiment rc004 is carried out with a 30% higher value ( = 54.6 mW m⁻²), and experiment rc005 with a 30% lower value ( = 29.4 mW m⁻²). These variations cover roughly the range of uncertainly of (for a detailed discussion of measurements of geothermal heat fluxes see Reference LeeLee (1970)). The iterations start with the present topography as obtained from measurements by Reference Letréguilly, Reeh and Huybrechts.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, A_i,b, H _t.max, T _GRIP, T _CC and T _Dye3 † (see notation list), supplemented by the observed values at hand (topography data from Reference Letréguilly, Reeh and Huybrechts.Letréguilly and others (1990); ice-core data from Reference WeertmanWeertman (1968), Reference Gundestrup and HansenGundestrup and Hansen (1984), Reference JohnsenJohnsen 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⁻², 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⁻², 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 h_max, V_tot, V_temp, A_i,b, A_t,b, H_t,max, T_GRIP, T_CC and T_Dye3 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 (Reference Calov and Hutter.Calov and Hutter 1994a,Reference Calov and Hutter.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 km³ (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 (Reference Iken., Echelmeyer, Harrison and Funk.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 (Reference Macheret, Moskalevsky and VasilenkoMacheret and others, 1993) may remedy ibis lack of data in the future.