Two-Dimensional Thermo-Mechanical Modelling of Flow and Depth-Age Profiles Near the Ice Divide in Central Greenland

Abstract

The ice divide on an ice sheet is a special place on the ice sheet, where the velocity and stress distributions are different from those found just two ice thicknesses from the ice divide. A two-dimensional thermo-mechanical model has been used to model the flow near the ice divide in central Greenland. The results show that the surface strain-rates are increased by 50% of the values found 15 km down-stream and that the basal temperatures rise by 3°C. Estimates of a depth–age distribution in the vicinity of Crête show that the isochrones rise at the dome. The model calculations predict 1600 m of Holocene ice, 1000 m of Wisconsin ice, and 400 m of ice older than 115kyear. 15 km down-stream from the dome, the thickness of ice older than 115 kyear will be reduced by 50%. In these calculations it has been assumed that the ice thickness at the dome is time-invariant. The existence of the old ice depends on the basal temperature. The estimates of the basal temperature depend strongly on past accumulation rates and on the geothermal heat flux but the ice has most likely been below the pressure-melting point throughout the last glacial period.

Introduction

The ice divide on an ice sheet has been shown to be a location where the temperature, velocity, and stress distributions in the ice are special. Many ice-sheet flow models have been made to describe the flow of ice along an entire flow line or locally at some distance from the ice divide. In order to simplify the very complex set of thermo-mechanical equations, those special requirements essential to a good ice-divide solution have often been omitted. Since the ice-divide region is narrow, poor solutions here will not seriously affect model calculations of surface elevation, temperature distribution, stress, and velocity distribution away from the ice divide.

However, the ice-divide region, though it is a complicated region to model, is of great interest as a deep ice-core drilling site. There is no horizontal ice motion, so there is no need for corrections for advective transport. If the position of the ice divide has been unchanged with time, the annual layer thickness, depth–age relationship, isotopic composition, concentrations of impurities, etc. in ice cores must be explained by the temporal changes of precipitation rates, surface temperatures, and atmospheric conditions at the ice divide. Hence, we do not have to make up-flow corrections along an entire flow line from the ice divide to the drill location.

In recent years, ice-divide conditions have been modelled by Raymond (1983), Reeh (1988), and Dahl-Jensen (in press). In this paper, Dahl-Jensen's fully coupled thermo-mechanical flow model is used to discuss the situation at Crête, central Greenland, where an ice core will hopefully be drilled to bedrock in the near future.

The Ice-Divide Model

The model used to calculate an ice-divide solution has been described by Dahl-Jensen (in press). It is a steady-state two-dimensional thermo-mechanical model in which the free surface, the velocity, stress, and temperature distributions are calculated along a flow line from the ice divide to the ice-sheet terminus. The model calculations have been done by finite differences, and an iteration procedure has been used in order to de-couple the energy-balance equation from the rest of the field equations. This is essential in order to obtain a good ice-divide solution. The shallow-ice approximation (Hutter, 1982, 1983; Morland, 1984) is used to obtain a set of lead-order equations.

The inputs to the model are: the bedrock topography and the geothermal heat flux along the flow line, the ice thickness at the ice divide, and as functions of the unknown surface elevation, the surface temperature and accumulation/ablation rates. Glen's (1955) flow law is used as the constitutive relationship with a flow-law exponent n = 3.

where

are the strain-rates,

are the stress deviators, and T_θ is the second invariant of the stress-deviator tensor

Fig. 1. Site-location map of central Greenland in the vicinity of Crête (CR). Sites A–Η were established during the 1984–85 post-GISP campaign. The crosses in a west–east direction around 71°N show the E.G.I.G. sites in the area. The solid line on the E.G.I.G. line is the flow line used in the model calculations. The numbers on the surface-contour lines show the elevation in m a.s.l. The solid lines give the areal distribution of accumulation rates (in mm of ice equivalent per year). (From Clausen and others, 1988.)

Table 1

(Paterson, 1981, p. 31). B(T) is the temperature-dependent part of the flow-law parameter, and , the enhancement factor, depends on other ice properties. B(T) is an Arrhenius function of temperature and with E = 1, (−20°C) is 5.5 × 10⁻³ bar⁻³ year⁻¹ (Paterson, 1981, p. 39).

The model is used for the conditions found at Crête (71.12°N, 322.68°W) in central Greenland. Crête is located on the main north–south ice divide of the Greenland ice sheet and is close to the summit point of the east–west E.G.I.G. line, which in the west terminates north of Jakobshavns Isbræ (Malzer, 1964; Reeh, 1983; Clausen and others, 1988). The velocity data from this area (Hofmann, 1986) suggest a northward velocity which must be doubted, because the surface slopes towards the south (Fig. 1) (Clausen and others, 1988, fig. 1). Table I lists the data input used in the model calculations in order to determine the flow between Crête and 30 km down-stream along the E.G.I.G. line.

The topography along the E.G.I.G. flow line is shown in Figure 2a. The observations show a flat bedrock in the Crête area, so the mean slope of the bedrock for the first 30 km is set to 0, and this mean bed is used in the calculations.

The model results for E = 1 are shown in Figure 2b and c. It is seen that the surface longitudinal strain-rate at the ice divide increases by 50% of the strain-rate 30 km away. The basal temperature rises in a narrow area near the ice divide, with 3°C creating a “hot spot” here (Paterson and Waddington, 1986). Figure 3 shows the stress, velocity, and strain-rate profiles over depth at the ice divide and five ice thicknesses away. The horizontal velocity and the shear stress are 0 at the ice divide. The profiles of the vertical velocity and the longitudinal stress deviator at the divide are different from the similar profiles five ice thicknesses away. Compared to off-ice-divide locations, the vertical velocities found near the base at the divide are smaller than those found further down-stream. The longitudinal stress deviator at the ice divide has a very strong vertical gradient close to bedrock and the corresponding horizontal strain-rate is concave with an inflection point 800 m above the bedrock (Dahl-Jensen, in press). These ice-divide results are similar to the results obtained by Raymond (1983), Reeh (1988), and Reeh and Paterson (1988).

Fig. 2. Steady-state model calculations for the first 30 km of the flow line from Crête (x = 0 km) along the E.G.I.G. line. The dashed lines in (a) are the bedrock and surface elevations. The solid line at the bed is the mean bedrock for the first 30 km of the E.G.I.G. line. This is the bedrock used in the model calculations. The solid surface elevation is the calculated surface. The surface strain-rates in (b) and the basal temperatures in (c) demonstrate the special conditions found near an ice divide.

Fig. 3. Depth profiles of the shear stress τ _xz , the longitudinal stress deviator σ′_x, horizontal velocity , vertical velocity w, and horizontal strain-rate ∂u/∂x at Crête (solid curves, x = 0 km) and at 15 km distance from Crête along the E.G.I.G. flow line (dashed curves). At the ice divide, the shear stress τ _xz and the horizontal velocity are zero.

Depth–Age Profile At Créte

In order to discuss the advantages/disadvantages of drilling a deep ice core on an ice divide, calculations of a depth–age profile at Crête have been made.

Results from existing ice cores and bore holes in Greenland and northern Canada show that the large ice sheets are not in a steady state. There are several reasons for this:

• (a) The ice deposited on the ice sheet during the Wisconsin period, the last glaciation, deforms more easily than ice from the Holocene (Reeh, 1985; Fisher and Koerner, 1986; Dahl-Jensen and Gundestrup, 1987). The boundary between the Holocene and the Wisconsin ice moves downward in the ice with time and this changes the rheology of the ice sheet. As a first-order approximation, this phenomenon can be included in Glen's flow law by an enhancement-factor value of 1 in Holocene ice and of 3 in Wisconsin ice.

• (b) The temperatures in the central parts of the large ice sheets are still cooled by the cold climate of the Wisconsin period. At Dye 3, this cooling is calculated to be 5°C ( Dahl-Jensen and Johnsen,1986).

• (c) The ice sheets are not in mass balance. It has not yet been possible to measure whether the central parts are thinning or thickening (Paterson, 1981, p. 56–57; Reeh, 1985).

I have investigated what effect (a) and (b) would have on the model solution in the vicinity of the ice divide, with the assumption that the ice thickness at Crête is kept constant at its present value. The values from Table I are used unless otherwise noted. Figure 4 shows the vertical velocity at the ice divide and the surface slope in the vicinity of the ice divide has been calculated for each of the following cases:

• (a) Enhancement factor E = 1 for 0 ≤ ≤ H. This is a steady-state Holocene situation.

• (b) Enhancement factor Ε = 3 for 0 ≤ ≤ H. This is still climatically a Holocene situation but the enhancement factor has been raised to 3 as for the Wisconsin ice.

• (c) Enhancement factor E = 1 for 0 ≤ ≤ 0.5. 3 for 0.5 ≤ ≤ . This enhancement-factor distribution is close to, that expected at Crête at present, where the transition between Holocene and Wisconsin ice is believed to be 1500 m above the bedrock.

• (d) Enhancement factor Ε = 1, as in (a), but the geothermal heat flux has been lowered from 40 mW/m², as used in (a)–(c), to 30 mW/m². This lowers the basal temperatures by 4.5°C.

• (e) Enhancement factor E = 1, as in (a), geothermal heat flux 50 mW/m². This increases basal temperatures by 4.5°C.

• (f) Enhancement factor Ε = 1, and geothermal heat flux 40 mW/m² as in (a), but the surface temperature has been lowered by 12°C from −30.4°C, as used in (a)–(e), to −42.4 °C. This is the cooling expected during the glaciation.

• (g) Enhancement factor Ε = 3, geothermal heat flux, and surface temperature as in (f), but the accumulation rate has been changed from 0.29 m/year, as used in (a)–(f), to 0.10m/year. These are the conditions which could be found in a Steady-state Wisconsin ice sheet.

It can be seen from Figure 4a that the vertical velocity profiles (a)–(f), in which the accumulation rate has been kept constant, are nearly identical. They are not seriously affected by changes of the enhancement-factor distribution ((a)–(c)), nor by changes of the temperature distribution ((a), (d), (e)). However, the increase of the surface slope with distance from the ice divide does change in order to maintain mass balance (Fig. 4b).

Close to the ice divide, where the horizontal velocity is small, the depth–age profile of the ice is controlled by the vertical velocity. It is therefore seen that the depth–age profile close to an ice divide cannot be seriously affected by the distribution of the enhancement factor or by the temperature through the ice. This is an astonishing result which is due to the special situation at the ice divide. The stress configuration here is dominated by the longitudinal stress deviator, which can be interpreted as the drag of the ice sheet on the ice divide. The shear stress is small in this area. Away from the ice divide, where the shear stress dominates, the flow will depend on the temperature and the enhancement-factor distribution.

Fig. 4. Vertical velocity–depth profiles for cases (a)–(g) in Figure 4a. In Figure 4b, the surface slope from the ice divide (x = 0 km) and 5 km down-stream from cases (a)–(g) are shown.

During the Wisconsin period, the surface-accumulation rate on the Greenland ice sheet was lower than at present. Most derivations of the average Wisconsin accumulation rate suggest values of one-half to one-third of the present (Nye, 1963; Hammer and others, 1978; Beer and others, 1984; Paterson and Waddington,1984; Dahl-Jensen and Johnsen, 1986). Although temperature and enhancement-factor distributions at the ice divide have little influence on the vertical velocity profile, the accumulation rate has a large influence when the ice thickness at the ice divide is kept constant. In Figure 4a, the vertical velocity profile with an accumulation rate of 0.10 m/year, which is one-third of the present, is shown. The velocity profile (g) is seen to have the same shape as (a)–(f) but scaled with 0.33 as the accumulation rate.

In the area close to the ice divide, it is possible to derive a “non-steady-state” model by combining a steady-state Holocene model and a steady-state Wisconsin model in order to calculate a depth–age profile here. In the Wisconsin model, the surface temperature is −12°C colder than the present temperature (−42.4°C), the enhancement factor is 3, and the accumulation rate is reduced to 33% (λ_gl = 0.15m/year) or to 50% (λ_gl = 0.15 m/year) of the present accumulation rate. In the Holocene model, the present surface temperature (−30.4°C), the accumulation rate (0.29 m/year), and the enhancement factor Ε = 1 are used.

The fundamental and critical assumption is that the ice thickness at the ice divide has been unchanged with time. How the ice thickness at Crête has actually changed is unknown (Reeh, 1985). Lowering of the surface temperature will produce a cold wave that will reach the bedrock after 10–20 kyear. This will make the ice more difficult to deform and thereby increase the ice thickness at Crête. The ice deposited during the Wisconsin period, however, deforms more easily than the ice from the Holocene period. After 50 000 years of glaciation, most of the ice will be of the Wisconsin type and this softer ice will result in a decrease of ice thicknesses at Crête. Decreasing the accumulation rate will both change the mass balance of the ice sheet and increase the basal temperatures. To obtain a complete answer to this question involves non-steady-state modelling (Budd and Young, 1983). In this reconstruction the ice thickness at the dome is kept constant.

The “non-steady-state” model is created by repeating a glacial cycle with 100 kyear of Wisconsin conditions (the steady-state Wisconsin model) and 10 kyear of Holocene conditions (the steady-state Holocene model), assuming constant ice thickness at the ice divide. 2.5 kyear are allowed between each step in which the temperature and accumulation rate are shifted. This attempt to make a “non-steady” model near an ice divide is justified by the fact that the changing temperatures and enhancement factor in the ice has little effect on the vertical velocities. The vertical velocity, however, is sensitive to the mass balance and some time must be expected before the ice sheet has adjusted to “new” kinematic conditions, so the steady-state models give reasonable predictions of the velocities. Isochrones are constructed by following the particle paths with time. In Figure 5, the isochrones marking the transitions between Wisconsin and Holocene conditions are presented. At the ice divide it is seen that the isochrones rise, creating an isochrone bump. At the ice divide there are 1600 m of Holocene ice, 1000 m of Wisconsin ice, and 500 m of ice older than 115 kyear when λ_g1 = 0.10m/year. For λ_g1 = 0.15 m/year, there are 1600 m of Holocene ice, 1100 m of Wisconsin ice, and 400 m of ice older than 115 kyear. For both examples, it can be seen that five ice thicknesses from the ice divide (15 km) the thickness of ice older than 115 kyear is reduced to half that at the ice divide. The model predicts 60–80 m of Eemian ice at the ice divide and 30–50 m of Eemian ice at a 15 km distance from the ice divide.

Fig. 5. Model calculations of the depth–age distribution in the vicinity of an ice divide. The climatic inputs used in the calculations are the surface temperatures and accumulation rates through a 115 kyear glacial–interglacial cycle as shown at the top of the figure. Two cases are shown, λ_g1 = 0.10m/year and λ_g1 = 0.15 m/year. The isochrones are those from the start and the end of glacial periods, i.e. 10 kyear B.P., 115 kyear B.P., 125 kyear B.P., 235 kyear B. P. …, The shaded areas represent intermediate glacial periods. The isochrone bump at the ice divide will result in a tilted annual layer, which is shown in the last frame of Figure 5.

The isochrone bump at the ice divide is only one to two ice thicknesses wide, so the isochrones will have a slope because of this bump. Annual layers in a core drilled at an ice divide will therefore be tilted. The slope of the annual layers between the ice divide and one-half ice thickness away are shown in Figure 5. The maximum slope is 5–7° 600 m above the bedrock. Increasing layer tilting from 5° at 20 m depth to a maximum of 13° at 105 m depth have been observed in two surface to bedrock cores (total ice thickness 130 m) taken within 30 m of the top of the Agassiz ice cap, Ellesmere Island, Canada. These have been attributed to the proximity to the top of the ice cap (personal communication from R.M. Koerner, D.A. Fisher, and N. Reeh).

Basal Temperature at Créte

The calculations of the isochrones at Crête are based on the assumption that there has never been melting at hte bedrock, i.e. the bedrock temperature has never exceeded the pressure-melting point of −2.4° C. Paterson and Waddington (1986) have discussed the range of geothermal heat fluxes that would give non-melting conditions at the bedrock. It was concluded that, with a heat flux of less than 48 mW/m², the bedrock temperature is lower than the pressure-melting point. If the bedrock is Precambrian rock with a heat flux of 38.7 ± 8 mW/m² (Lee and Uyeda, 1965), the ice has most likely been frozen to the bedrock throughout the last glacial–interglacial period.

The non-steady-state temperature model used to calculate the temperature profile at Dye 3 (Dahl-Jensen aned Johnsen, 1986) is applied to Crête. The velocity and stresses used in the temperature calculations are those derived from the steady-state flow model for Crête. The results agree reasonably well with those of Paterson and Waddington (1986).

In order to avoid repeating the results presented by Paterson and Waddington, we shall discuss here the basal temperature-dependence on the glacial accumulation rate. The geothermal heat flux is set to the mean value for Precambrian bedrock, i.e. 40 mW/m², and the climatic input; in this case, the surface temperature and accumulation rate are kept simple as shown in Figure 5. It is known that the temperatures and the precipitation rate during the Wisconsinan stage was not constant but, for this purpose, it is sufficient to use average values. The glacial surface temperature is −42.4 °C, 12°C colder than the present surface temperature at Crête. It is seen in Figure 6 that the basal temperature changes 8°C when the accumulation rate during the glaciation is one-third to two-thirds of its present value. When the accumulation rate is two-thirds of the present, the basal temperature is close to steady state at the end of the glaciation (15 kyear B.P.). When the warm period begins, the basal temperature first decreases because the increased accumulation rate increases the advection of cold from near the surface ice to the bedrock. It takes 15 kyear after the termination of the warm period for the “warm” wave from this warm period to reach the bedrock and for the basal temperatures to reach maximum values. When the glacial accumulation rate is reduced to one-third of its present value, the temperature decrease under the warm period is more pronounced, while the temperature increase after the termination of the warm period is smoothed out more and a stable basal temperature at the termination of the glacial period (12.5 kyear B P.) is not reached.

Fig. 6. Calculations of the basal temperatures through the 115 kyear glacial–interglacial cycle shown at the top of Figure 5. In all of the three cases presented, the mean glacial temperature is −42.4° C.

In general, a decreased glacial accumulation rate results in increased temperatures at the base. For each case, the basal temperature changes 2° C throughout the glacial–interglacial cycle. The basal temperatures depend strongly on the unknown geothermal heat flux and the past precipitation rate at Crête.

Discussion

In order to discusss the assumptions for the “non-steady-state” ice-divide model and the possibility of observing rises in the isochrones, basal temperatures, and surface strain-rates at an ice divide, we must investigate the time-scales needed to reach steady-state conditions.

The temperature calculations made through a glacial cycle show that the basal temperature reaches a constant value after 100 000 year or more, which must be the time-scale needed for the temperature to reach static conditions at an ice divide.

Outside the ice-divide area, steady state of the ice thickness, i.e. mass balance of the ice sheet, depends on the temperature of the ice and on ice rheology. Since there are internal layers in the ice consisting of intermediate ice with an enhancement factor of 1 and of glacial ice with an enhancement factor of 3 that moves with ice flow, the rheology changes with time. The depth–age profiles in Figure 5 show that several glacial periods must be expected in the ice, all slowly sinking down through the ice. The time-scale needed for the ice sheet to reach mass balance is therefore more than a glacial-cycle length of 115 000 year.

Finally, the velocity and stress distributions through the ice sheet are believed to adjust themselves within some thousands of years to the actual temperature distribution, accumulation, and geometry of the ice sheet because they are mostly determined by the momentum balance of the ice.

At the ice divide, the ice thickness is assumed to be constant, and ten ice thicknesses away from the ice divide the ice thickness changes less than 1% between the situations (a)–(g) in Figure 4. The area considered is therefore always close to mass balance. The vertical velocities in this area are somewhat independent of the temperature distribution and the depth variation of the enhancement factor, so the positions of the isochrones found by the “non-steady-state” model are believed to be reasonable.

The isochrone and temperature results are very dependent on the temporal changes in precipitation rate, as well as the assumption of constant ice thickness at the ice divide. The precipitation rates during the last glacial period and earlier are poorly known, so the isochrones can only be considered reliable as far back as the Wisconsin/Holocene transition at = 1500 m, 10 000 year B.P. (Reeh, in press). The very simplified climatic input and the approach to a non-steady-state model by combining two steady-state models is therefore believed to give a view of what can be expected at Crête. The phenomenon of an isochrone bump, which implies that the amount of ice-core material older than 115 000 year from a core drilled at the divide is twice as large as in an ice core drilled 15 km away, is interesting. It is deduced from the model calculations with any climatic history, as long as no bottom melting has occurred.

The presence of an increased amount of old ice at Crête depends on how stable the position of the ice divide has been. The isochrone bump is very narrow and is only found in a 12 km area around the ice divide. Model calculations on the glacial ice cover of Greenland with a perfectly plastic ice-sheet model (Reeh, 1982, 1984) estimate an ice-divide movement of less than 50 km to the west. Earlier positions of the ice divide are unknown but even small movement would make the isochrone bump and the temperature bump less pronounced. The increase in the surface strain-rate, however, is expected to adjust within some thousands of years to the present position of the ice divide. E.G.I.G. data do suggest surface strain-rates as large as 1.3 × 10⁻³ year⁻¹ at the ice divide (Hofmann, 1986; personal communication from S.J. Johnsen).

The situation of Crête on the north-south-trending ice divide (Fig. 1) and the south-trending surface slope suggests a small south-trending flow along the ice divide. A study of how this affects the west-east-trending bump cannot be done by the applied two-dimensional model but this probably also tends to smooth the bumps of temperatures and isochrones.

Internal layers seen from radar echo-sounding in this area do not show a rise near Crête. However, the quality and cover of the area are insufficient to disprove that the internal layers have a local bump (Gudmandsen, 1973). New radio echo-soundings made in 1987 could change this picture.

Conclusions

The model calculations presented here show that a deep core at an ice divide has the advantage of yielding more old ice. It is predicted that there would be up to twice the amount of ice-core material older than 100 kyear in an ice core from the ice divide than 15 km away.

The isochrone rise at the ice divide at Crête is probably less pronounced because of slow flow along the ice divide from the site Summit, 100 km north of Crête, towards Crête, and a possible movement of the ice divide with time.

The estimates of the basal temperatures at Crête show a wide range of varying values, mainly due to the unknown geothermal heat flux and past accumulation rates. Most of the estimates, however, suggest that the basal ice has been frozen throughout the last glacial period, so the old ice is preserved at Crête and along the first part of the flow line.

Since the data material from Crête is so sparse, model predictions are bound to be uncertain. Surface velocities and surface strain-rates in the area would be very helpful. And eventually, one day, perhaps data from an ice core and a bore hole at Crête will increase our knowledge of past climatic conditions and of the ice flow close to an ice divide.

Acknowledgements

This work has been funded by the Danish Science Research Council, the European Economic Communities, XII Directorate General (contract CLI.067.DK), and the Danish Commission for Scientific Research In Greenland.