Introduction

The GRIP European and the nearby GISP II American deep-drilling projects at the Summit ice divide in central Greenland (lat. 72°18'N., long. 37°55'W.) expect to retrieve ice cores dating back several glacial cycles. The temperatures measured in the bore hole, the crystal growth and fabric development within the cores, and the age of the oldest ice they contain will all depend on the thermal history of the ice around the coring site.

The detailed thermal history is also needed to determine if the basal ice has melted, a matter of great importance, because it would mean the loss of the oldest ice at the drilling sites or down-stream. We first examine some steady models to show that melting is possible. We assume that the vertical velocity pattern v(y) at a steady ice divide can be approximated by (Reference RaymondRaymond, 1983)

where y is the height above the bed, H is the ice-sheet thickness, and v _s is the vertical velocity at the surface (equal and opposite to accumulation). Then we solve the steady heat equation

where к = K/pc is thermal diffusivity, and Κ, ρ, and c are thermal conductivity, density, and specific heat capacity at constant volume, all assumed constant for this steady-state analysis. All body sources and sinks of heat (e.g. strain heating) are contained in the source term f(y). The boundary conditions are temperature T _s at the upper surface, and geothermal heat flux Q = К∂Т/∂у at the bed.

We solve Equation (2) using an integrating factor F (y) where

to get

The first term is a constant, the second term is the shape of the temperature—depth curve in the absence of body sources f (y) and the third term is the correction due to those body sources. Neglecting the third term for now, and using constant thermal parameters for ice at 0°C (Reference PatersonPaterson, 1981, p. 186), we plot in Figure 1a the steady temperature profiles for conditions appropriate to the present ice sheet (H = 3000 m) in central Greenland (q = 41.7 mWm⁻²= 1.0HFU, typical of Precambrian rocks). Curve A shows the temperature under steady interglacial conditions, i.e. T_s = –32.0°C and v_$ = –0.23 m year⁻¹. Curve B shows the temperature under plausible, steady ice-age conditions, i.e. T_s = –40.0°C and v_s = –0.10 m year⁻¹. In both cases the basal temperature falls within 1 °C of the pressure melting-point. To assess sensitivity to the unknown basal heat flux Q, we reduced Q by 20% to get curve C It is not inconceivable that Q could exceed 1.0 HFU by a similar amount, resulting in basal melting for both glacial and interglacial conditions.

Fig. 1. Fig. 1. Simple models suggest basal melting may destroy old ice at Summit, (a) Steady-state temperature profiles for constant K and κ, and linear vertical strain-rate (Reference RaymondRaymond. 1983). H = 3000 m. Curve A: Q = 1.0 HFU (41.7 mW m²). T_s = –32°C. and v_s = –0.23 m year⁻¹, modern interglacial conditions. Curve B: Q = 1.0 HFU (41.7 mW m⁻²). T_s = –40°C and v_s = –0.10 m year⁻¹, typical ice-age conditions. Curve C: same as B, except Q = 0.8 HFU (33.36 mW m⁻²). showing sensitivity to unknown basal heat flux Q. Curve M: the pressure melting-point for ice. (b) Steady-slate time-scales for interglacial (dashed) and ice-age (solid) conditions for various basal melting rates m. Curve A: m = 0.0. Curve B: m = 0.001 m year⁻¹. Curve C: m = 0.005 m year⁻¹. Curve D: m = 0.01 m year⁻¹.

If the base reaches the pressure melting-point for all or part of a glacial cycle, old ice will be lost. We assess the sensitivity of the time-scale to basal melt rate m for the steady interglacial and ice-age sheets. We modify the velocity pattern in Equation (1) to incorporate basal melt rate m < 0 so that

The depth–age relation is then

If all the geothermal heat flux of Q= 1.0 HFU went to melting ice in an isothermal basal boundary layer, |m| would be about 0.005 m year⁻¹. If Tertiary volcanic rocks exist under Summit, a substantially larger melt rate could exist. If the basal temperature reaches melting during only part of a glacial cycle or if some heat is conducted away from the interface, a steady-state model with a melt rate of |m| = 0.001 m year⁻¹ may be a useful analog. Figure 1b shows the time-scales for melting rates of 0.0, 0.001, 0.005, and 0.01 m year⁻¹. In the worst case (|m|=0.01 m year⁻¹), there would be no ice older than 100 kyear left. A core could not even recover ice from the previous interglacial period. Because the Greenland ice sheet is probably never in steady state (Reference WaddingtonWaddington, 1987) and we have made other simplifying assumptions, these analyses are intended only to suggest that the problem of basal melt needs serious study. Similar concerns motivated Reference Paterson and WaddingtonPaterson and Waddington (1986) to examine the basal temperature at Crěte, Greenland, using a one-dimensional time-dependent model.

Two-Dimensional Transient Model

Because of the size and long response time of the Greenland ice sheet, neither steady-state models nor borehole measurements alone can provide an adequate thermal history. A transient thermal model is required to extrapolate those measurements over space and time. More generally, the heat equation for a deforming medium can be written in terms of temperature T, time t, and Cartesian coordinates x _m as (Reference Paterson and ClarkePaterson and Clarke, 1978)

For an ice mass, Equation (7) couples to the momentum-conservation equations (Reference HutterHutter, 1983, chapter 3) through the velocity terms v _m and the variation of ice rheology with temperature.

In this study, we decouple the coupled heat and mass-conservation equations and run an ice-flow model using present-day surface mass balances and temperatures to establish a modern velocity pattern. We assume that the ice-divide geometry has not changed over time. We then scale the velocity pattern to an assumed mass-balance history and run a separate heat-flow model with the specified velocities.

We use a two-dimensional finite-element ice-flow model consisting of an 8 × 15 mesh of four-node quadrilaterals. We specify a stress-free upper surface, a frozen bed, and velocities along the vertical boundaries as inReference Raymond Raymond (1983). For surface and bed topography, we follow data provided by Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others (1990), who also brought to our attention the 1974 value of the 10 m firn temperature. We assume a temperature profile that is near isothermal in its upper part, with the maximum gradient near the bed where the ice is near the pressure melting-point. The model solves the stress-equilibrium equations for velocities which vary linearly between nodes. We make minor adjustments to rheology and to a lateral divergence factor to make the vertical velocity solution match the present-day mass balance. The final lateral divergence factors indicate flow is very near plane strain.

For the heat-flow model, we have adapted a general purpose, non-linear, finite-element program capable of solving steady-state or transient problems in one to three dimensions using three- to nine-node elements, with internal heat sources, surface fluxes, and variable material properties (Reference Taylor and ZienkiewiczTaylor, 1977, chapter 22). We use nine-node elements in two dimensions, 2× 2 Gaussian quadrature, and a fully implicit time-stepping scheme. This formulation incorporates horizontal conduction and advection. By setting f = 0, we neglect internal heat sources (e.g. strain heating).

In the model, we represent the divide by approximately 3 km of ice over 2 km of bedrock using columns of seven ice elements and three rock elements. We assume a continuous heat flux at the ice-bedrock boundary (i.e. no melting) and introduce a standard geothermal flux of 41.7 mWm⁻² at the bottom of the bedrock to include the effect of transient thermal storage in the bedrock on heat flux into the basal ice (Reference Paterson and ClarkePaterson and Clarke, 1978;Reference Budd, Young and Robin Budd and Young, 1983; Reference RitzRitz, 1987). At the surface, we force the temperature to follow its probable historical values (Figure 2a, derived fromReference Paterson and Waddington Paterson and Waddington (1986), adjusted to present 10 m firn temperature at Summit). We use the velocity pattern taken from our two-dimensional ice-flow model with the velocities scaled to a probable mass-balance history (Figure 2b, derived fromReference Paterson and Waddington Paterson and Waddington (1986), scaled to present mass balance at Summit). Within the ice, we match the boundary of each heat-flow element to the boundary of an ice-flow element and assume the nodal velocities of both vary bilinearly. This forces the heat-flow model, whose temperatures and velocities may vary biquadratically within each element, to interpolate the same bilinear velocity field assumed by the ice-flow model. The model geometry is shown in Figure 3.

Fig. 2. Assumed surface temperatures and ice-equivalent precipitation rate.

Fig. 3. Temperature model element pattern, nodal velocities, and thermal boundary conditions. For clarity, only the velocities at the corners of each element are plotted.

We start the model using a steady-state solution corresponding to the temperature and mass-balance conditions at 135 kyear b.p. We step the model through time for six glacial cycles using 5000 year time steps to let the model “forget” its initial conditions. We then step through two additional cycles using 500 year time steps to get and check a final solution.

Results

To check the numerical model and to provide a base line to study the effects of two-dimensional ice flow, we first duplicated results from the one-dimensional finite-difference model used by Reference Paterson and WaddingtonPaterson and Waddington (1986) by using both single and multiple columns of elements with insulated vertical sides, zero horizontal velocities, and vertical velocities that did not vary horizontally. The results for conditions at the divide can be seen in Figure 4. Curve A shows the basal ice temperature over the last glacial cycle for an ice divide with a flow pattern in Equation (1) for isothermal ice (Raymond, 1983) and with ice conductivity and diffusivity fixed at their 0°C values. Curve B shows the basal temperature using the vertical velocity pattern from the central column of nodes in Figure 3 in our two-dimensional velocity model. Results from the one-dimensional finite-difference model (dotted lines) agree with the results from the two-dimensional heat-flow model to within 0.1 °C.

Fig. 4. . Basal temperatures at Summit ice divide from one-dimensional (dashed curves) and two-dimensional models assuming: curve A: ice-flow pattern at an isothermal ice divide (Raymond. 1983). no horizontal heat diffusion or advection, ice conductivity and diffusivity fixed at their 0 ° C values; curve B: as in curve A but flow pattern now taken from the two-dimensional velocity model; curve C: as in curve B but ice conductivity and diffusivity now assume temperature-dependent values; curve D: as in curve C but horizontal heat diffusion and advection now included in full two-dimensional treatment.

We then ran the model under more realistic conditions. In Figure 4, curve C, we allowed the ice conductivity and diffusivity to assume temperature-dependent values (Reference YenYen, 1981, Equation (33), fig. 1, table 2):

We overcame the non-linearity this introduced by using Newton–Raphson iteration at each time step. The colder temperatures relative to curve B result from higher conductivity at low temperatures.

Finally, we modelled the two-dimensional ice divide at Summit using ice and bed topography determined byReference Hodge, Wright, Bradley, Jacobel, Skou and Vaughn Hodge and others (1990): we added columns of elements to both sides of the divide so that the ice and rock extend 30 km to east and west, and we supplied the model with the full two-dimensional flow pattern from the velocity model. Because the flow and isotherm patterns are nearly horizontal far from the divide, we again used a simple insulated boundary condition at the sides without causing significant error in the temperatures near the divide. Curve D shows the resulting temperatures at the base of the divide. We find that the basal temperatures are again colder by an additional 2.4 °C due to horizontal advection in the ice, and to horizontal diffusion in both ice and bedrock. By comparing two steady-state two-dimensional finite-element solutions, one with constant heat flux Q applied at the base of the ice, and the other with Q applied deep in the rock (Fig. 3), we estimate that 0.5 °C of this cooling can be attributed to horizontal conduction in the bedrock. Application of the analysis of Waddington (1987) to a two-dimensional steady-state solution also predicts this.

Our mass-balance and surface-temperature histories have contrasting effects on the basal ice temperature. In Figure 5, curves A–C, we vary the mass-balance history. In curve B, we hold it constant over time. With no temporal variation in advection, the basal temperature curve is a damped replica of the surface-temperature curve, delayed about 20 kyear, and a mirror image of the basal temperature of our standard model. In curve C, we let the mass balance vary but go no higher than today’s value (dashed curve in Figure 2b). With less downward advection of cold surface ice early in the glacial cycle, the basal temperatures are warmer and reach their minimum value 30 kyear later. A damped remnant of the surface-temperature peak also appears at 110kyear b.p.

Fig. 5. Basal temperatures at Summit ice divide for precipitation rates: (A) following Figure 2b; (B) constant at today’s value; (C) never exceeding today’s value (dotted line of Figure 2b).

In Figure 6, we increase the geothermal heat flux and show how it raises the maximum basal temperature reached in a glacial cycle. In curve A, we assume the standard mass-balance history of Figure 2b; in curve B, the reduced mass-balance history of Figure 2b. We find that the basal ice melts if the geothermal heat flux exceeds 1.3–1.4 HFU (56–59 mWm ⁻²).

Fig. 6. Sensitiviiy of maximum basal temperature to geothermal flux Q and to assumed interglacial precipitation rate. Curve A derived with solid precipitation curve in Figure 2b. Curve B derived with dashed precipitation curve in Figure 2b. Dashed line indicates pressure melting-point.

So far, we have concentrated on the temperature at the base of the divide because the oldest ice is likely to be found there and because basal ice down-stream originated at the divide. Several studies have predicted that the temperature under an ice divide should be warmer than under the flanks because the downward advection is less there (Reference Paterson and WaddingtonPaterson and Waddington, 1986;Reference Dahl—Jensen Dahl-Jensen, 1989). Such a “hot spot” appears in our model. In Figure 7, we show its evolution over time. The hot spot is most pronounced at the onset of glaciation when the mass balance is high; this confirms its advective origin.

Fig. 7. Basal temperatures across the Summit ice divide over a glacial cycle.

In Figure 8, we show the vertical temperature profile at the divide over time. For comparison, we also plot the temperature profile used by the ice-flow model to calculate the velocity. Below a depth of about 2.0 km, the temperature gradient is nearly constant, with temperatures at any particular depth varying by about 3°C over time. Above that depth, the temperature variations are large; at 125 kyear b.p. the temperature gradient reverses sign.

Fig. 8. Vertical temperature profiles at the Summit ice divide at 125 kyear b.p.. (curve A), 115 kyear B.P. (curve B). 90 kyear B.P. (curve C). 40 kyear B.P. (curve D). and the present ( curve E). Dotted line is the temperature profile used in the velocity model.

Discussion

The close agreement of our two-dimensional finite-element model simulating one-dimensional heat flow with a one-dimensional finite-difference heat-flow model confirms its proper mathematical behavior. That the results for our two-dimensional model simulating two-dimensional heat flow confirm the predictions made by the one-dimensional model Reference Paterson and Waddington(Paterson and Waddington, 1986) lends further support that our two-dimensional model has a proper formulation. For example, as in the one-dimensional model, we find that a hot spot arises under the divide and that the basal temperatures are dominated by the assumed geothermal heat-flux and mass-balance boundary conditions. We also find little likelihood of basal melting over the last 130 kyear if the geothermal flux falls in the range characteristic of Precambrian shield terrain. (A compilation of observed Precambrian Shield flux values (Reference Sclater, Jaupart and GalsonSclater and others, 1980) suggests roughly a 20% chance of this having happened.) In fact, as a result of the refinements incorporated in our two-dimensional model, we predict basal temperatures lower than those predicted byReference Paterson and Waddington Paterson and Waddington (1986), even after accounting for the different boundary conditions between Crête and Summit.

The GISP II drilling site lies roughly at the western edge of our element grid where we have imposed boundary conditions, whereas the GRIP site lies at the divide proper (Fig. 3). Thus, we expect our computations to be most accurate near the latter and rather less accurate near the former. We have concentrated on the area near the divide because, should the basal ice melt there, there is little chance of finding older ice elsewhere. To provide an accurate analysis for the GISP II site, we should extend the western boundary of our model 10–30 km.

Near the divide, the two-dimensional model results are limited first by the central simplification, i.e. the decoupling of the heat and mass-conservation equations. The temperature pattern at the divide is dominated by advective effects (e.g. see the increased influence of the surface-temperature history in Figure 5, curves A, C, and B, as the variability of advection is reduced) but our model velocity field (driving the advection) does not yet respond to temperature variations, which even at the bed exceed several degrees over the course of a glacial cycle. This must influence the velocity pattern, particularly near the bed where the temperature may at times approach the pressure melting-point. Figure 4, curves A and B, shows that changing from the velocity pattern for isothermal ice (Raymond, 1983) to a pattern calculated for cold ice at Summit changes calculated basal temperatures by about 2.5 C. If velocity and temperature were fully coupled and mutually consistent at each time step, calculated basal temperature would be altered in a complicated way, but by an amount less than 2.5 ° C.

The temperature assumed by our velocity model exceeds the output of the heat-flow model by up to 10 °C (Fig. 8). By assuming warmer, softer ice at depth, we have exaggerated the vertical advection in the heat-flow model and calculated temperatures significantly lower than those assumed in the velocity model. However, the error this introduces in the basal temperature is small. This is because we calculate the velocity by assuming a temperature pattern and adjusting the flow-law parameter, A ₀, (Reference PatersonPaterson, 1981, chapter 3) until the surface velocity matches the known surface mass balance. Were we to cool the ice in the velocity model to make its temperature consistent with the temperature model, we would have to soften A ₀ to remove much of the temperature-induced hardening in order to match this constraint. From tests using a one-dimensional model, we expect that the resulting changes in advection would alter our calculated basal temperatures at Summit by less than 0.5 °C.

Strain heating has been neglected. Its influence can be estimated by an independent steady-state analysis. Differentiating Equation (1) gives the vertical strain-rate –∂u/∂xwhich by continuity equals -ди/дх. At the divide, the shear stresses vanish, and the strain-heating term (f(y) term in Equation (2)) is set to (Reference Paterson and ClarkePaterson and Clarke, 1978)

where n and A are the constants for the Glen flow law for ice (Paterson, 1981, p. 85). Numerical integration of the final term in Equation (4) using Equation (10) for f(y) shows that strain heating would increase the basal temperature by about 0.1 °C.

The two-dimensional results are limited secondly by our assumption that the ice thickness is invariant. This probably exaggerates the influence of mass-balance changes on basal temperatures through advection. The central Greenland ice sheet was perhaps 150 m thicker along the north–south main ice divide during the last ice age (Reference ReehReeh, 1984) and may have been much thinner during the Sangamon interglacial (Reference KoernerKoerner, 1989;Reference Reeh and BleilReeh, in press). These changes would introduce different surface-temperature boundary conditions and alter the stresses controlling the ice-flow pattern. The net effect on the basal temperature is again complicated and unknown.

In order to study ice-flow and temperature patterns at higher resolution required for detailed ice-core interpretation, the coupled ice and heat-flow equations must be solved. The full two-dimensional temperature model advances us toward that goal.

Acknowledgements

This work was supported by the U.S. National Science Foundation under grant NSF-DPP 8613935. We thank S.M. Hodge for providing data on the Summit area and R.L. Taylor for providing his finite-element program. We thank W.S.B. Paterson for discussions that have contributed to this work.