1. The Three Thermal Regimes of a Cold Ice Sheet

The theoretical study of temperatures and strains in a cold ice sheet has been developed by Reference LliboutryLliboutry (1963, Reference Lliboutry1964–1965, Reference Lliboutry1966), with the following assumptions:

• a. The flow Iines lie in vertical parallel planes (the plane problem). This assumption is not absolutely necessary. It is sufficient to assume that transverse strains are small.

• b. The bedrock topography is smooth and receives everywhere the same geothermal heat flux, say KG ₀ (K is the thermal conductivity of ice = 70·3 MJ/m deg year, and G ₀ the geothermal gradient in a motionless ice sheet, when the melting point is not reached at the bottom).

• c. The ice flow law, in the restricted temperature range of concern, may be written:

  (1)

  

  denoting the effective strain-rate, τ the effective shear stress, θ the temperature difference from the melting point (this latter reaches −2 deg at the bottom of a 3 000 m thick ice sheet), and k a constant (k ≈ 0.23 deg⁻¹). The bottom temperature will be called T. (It must be noted that in author’s previous publications θ ₀ was used instead of T, and denoted what will be called here exp (kT).)

Glen’s law with n = 3 is generally adopted, and I shall do the same. Nevertheless, as any non-linear and elaborate theory ends in numerical computations, more accurate numerical values may be used, should we have them at hand.

It has been shown that three temperature regimes may exist in a steady state.

Regime I: The ice sheet is cold throughout, including the ice-bedrock interface (T < 0). The entire geothermal flux KG ₀ crosses this interface.

Regime II: The ice sheet is cold throughout, apart from the ice-bedrock interface (T = 0). Only part of the geothermal flux, say KG, enters into the ice, the remainder K(G ₀−G) (some centimetres per year) being lost in melting ice.

Regime III: There exists in the lower part of the ice sheet a temperate ice layer, without temperature gradient (the very small temperature gradient arising from the pressure-induced variation of the melting point is neglected). No geothermal flux enters into the ice (G = 0) It seems that in some regions this temperate bottom layer may be several hundred metres thick (Reference LliboutryLliboutry, 1966), but in many cases the thickness may be only a few metres.

It is only when such a temperate ice layer exists that the ice sheet can slide over its bed (Reference LliboutryLliboutry, 1967).

Thus a temperate bottom layer may exist under every definite outlet of the polar ice sheets, as well as under all parts of the catchment area where the velocities increase strongly and very long arcuate transverse crevasses develop. Ice dynamics, with this third regime, is a matter of sliding theory and does not differ essentially from temperate glacier dynamics (Reference LliboutryLliboutry, 1968, in press).

2. Temperature Profile at a Given Site

Temperatures are computed following the ice in its movement (Lagrange coordinates). In the steady state, the heat withdrawn to the surface through conduction and solid convection must balance the heat coming from the geothermal flux and the heat generated by the plastic deformation of ice. Each of these four processes gives rise to a term in the heat equation.

Another phenomenon, which has been numerically computed by Reference Jenssen and RadokJenssen and Radok (1963), occurs when the flow carries the ice to lower altitudes, where the air temperature is higher. This warming is imperfectly transmitted downwards, and consequently the temperature profile with depth shows a minimum. The calculations in this paper concern the centre of ice sheets where the surface warming remains relatively small, and its effect has therefore been neglected by the artifice of keeping the surface temperature constant at the lowest value, which it had at the ice divide. With this approximation, the problem at a given site is one of linear conduction, the only variable being the distance z from the bedrock. Moreover, it will be assumed that the parameters (thickness Z, superficial slope a, yearly balance A, etc.) change sufficiently slowly to allow us to take them as constants in the computation of the steady state. This will be called a “quasi-permanent regime”. Eulerian variables and Lagrangian variables can then become indistinguishable.

The problem can be solved if we consider in the ice sheet two superimposed layers (in the thermal regimes I and II). Above there is a colder and thus more rigid layer, where the deviatoric stresses are lower; in this layer the heat generation is negligible. Below we have a thin bottom shear layer, where heat generation is important, but where solid convection becomes negligible.

Let us prove this statement. In the upper layer, the rate of vertical strain is −aU/Z, where U denotes the horizontal velocity. The rate of horizontal strain is about +aU/Z and the effective shear strain The energy dissipated per unit volume and unit time is . With the values given later, this is of the order of 10–7 J/m³ year. In the bottom shear layer, the effective shear stress τ differs very little from f, the friction per unit area over the bedrock. The energy dissipated per unit volume and unit time is . With the values given later, it amounts to about 10⁻¹ J/m³ year, a million times more.

In the bottom shear layer, the heatFootnote * equation reduces to

This equation is wrong by a factor 2 in Llihoutry (1963, 1964–65).

The solution of Equation (2) is

where

In the upper part of the bottom shear layer, the curve z(θ) may be replaced by its asymptote:

On this asymptote the value θ = T is reached for z ₁ = z where

z ₁ is a rough estimate of the thickness of this shear layer where the temperatures are perturbed by the heat generation. Gz ₁ appears as a function of (G/G ₁), which reaches a maximum, 1·16/G when G/G ₁ = 0.460 In regime 1 (G = G ₀ ≈ 1 deg/44m), z ₁ remains below 50 m. In regime II, z ₁ may become large, a fact which would make the present theory inaccurate.

Let us examine now the upper part of the ice sheet. Replacing the firn by an ice layer of the same weight, the thickness of the ice sheet would be Z, and the vertical velocity of the ice at the surface A (positive downwards). When the ice sheet is “stationary” (i.e. Z is constant for given geographical coordinates), A equals the balance ( = specific budgetFootnote †), but this equality is not necessary for the moment.

Summer melting is assumed to be very small (or the firn layer very thin), so that melt water does not disturb the temperatures under the layer subjected to seasonal changes. This constant temperature close to the surface will be referred to as “surface temperature” and denoted S. When melting is absent, S is more or less equal to the mean air temperature.

In order to compute the temperature profile upwards, the profile within the shear layer will be replaced by its asymptote. As the bottom shear layer is very thin, the computation is then the same as if all the heat of deformation were generated at the very ice-bedrock interface, the temperature being there θ _c = T+(2/k) ln (2G ₁/(G ₁+G)) and the thermal gradient G ₁ In this simpler case, the thermal gradient at a distance z from the bedrock is (Reference RobinRobin, 1955):

where h denotes the thermal diffusivity of ice.

Integration of Equation (7) gives

As AZ > 2h, erf (AZ/2h)^1/2 ≈ 1 and we obtain

(In previous papers, the last term was omitted. It amounts to about 1 deg in the computed examples.)

Equations (4) and (8) hold in both thermal regimes I (G = G ₀) and II (T = 0, 0 < G < G ₀). Eliminating G ₁, they allow us to compute the adjustable parameter, T in regime I or G in regime II, from measurable or easily computable quantities. In the metrebar-year-degree system:

For G ₀, the value 1 deg/44 m may be adopted, at least for Greenland and eastern Antarctica which are old stable shields.

The friction f is given by the classical equation, deduced from the stresses operating in the body of the ice sheet as

where ρg is the weight of unit volume of ice = 1 bar/11·5 m and a is the mean surface slope within a length of some Z.

Lastly may be deduced from f if the flow law for basal ice were well known quantitatively. Unfortunately, this is not the case. We shall adopt the form

B remains unknown at present. It will be shown later that it may be estimated from the surface profile near the ice divide.

3. Stability of the Basal Temperature in Regime I

The possibility that the bottom temperature may be unstable was suggested by Reference LliboutryLliboutry (1964–65), according to the following heuristic considerations.

The solid convection increases when the ice sheet becomes more active. This is the case when the temperatures increase. Ice flow becomes easier and a slightly smaller thickness is sufficient for a given discharge. This produces a cooling of the bottom. Thus there exists a negative feed-back. Solid convection is a factor of stability. On the other hand ice deformation produces heat, and, as the temperature increases, the deformation is enhanced. This means a positive feedback. Ice deformation is therefore a factor of instability which can become more important than the stabilizing factor.

For a quantitative study of the phenomenon, let us rewrite Equations (4) and (8), using Equation (10) as

In regime I, G = G ₀. Numerically, with the assumed values

Equation (11) gives the thermal gradient which, for a given bottom temperature, results from the heat of deformation. Any change in the bottom temperature T causes an immediate change in the strain-rate in the vicinity, which extends upwards and rapidly involves all the bottom shear layer. A change in the heat generation and in G ₁ ensues. Thus any change in the bottom temperature causes an almost immediate change in R according to Equation (11).

Equation (12) gives the change of R which follows a change in T when a new steady temperature profile is established throughout the entire ice sheet, that is a long time after.

Thus, any local change in T (due for instance to a trough in the bedrock) affects R first according to Equation (11). The physical loop is therefore

described in this direction, and not in the opposite one. This loop is described slowly, and, if this physical process diverges, T would not rise instantaneously up to the melting point. What can be said is simply that Equation (12) is no longer valid; the thermal regime is now unsteady and in the heat equation the term Cρ dθ/dt must be taken into account.

Now a mathematical procedure for computing T and R may be provided by the same loop (although described at computer velocity) and in the same direction. (The starting values may be R = 1 and T ₁ = S+0.1757(Z/A)^1/2.) We assume here that the existence of a steady state corresponds to the mathematical stability of the iteration procedure. When mathematical instability appears, we may restore the stability by reversing the direction in the loop:

While this is a purely mathematical artifice for the example quoted it could have physical significance in cases of dominating surface temperature effects, of the kind mentioned in Section 2 but neglected in this paper.Footnote *

The mathematical convergence condition, for the iterative process here considered, is

or

Let us put

(a = 0·965 in the metre-bar-year system). Because of the modulus signs, Equation (15) gives

In the quite distinct and realistic examples which have been worked out, R/(R ²−1) ≫ 1/R(1+R) In this case it is only the inequality on the right which may be violated, and which limits stability in practice.

We can illuminate the mathematical process in the following way. Let us define the function

the roots of which are possible values for T. Its derivative is

According to Equation (17), dy/dT must be negative in order that the root may be stable. Now all the possible cases are represented in Figure 1. In both cases a₁, and a₂, there exists one stable solution; in case b two solutions, a stable one and an unstable one; in both cases c₁, and c₂, no solution. For regime I, there may moreover be five limiting cases, shown with self-explanatory notations in Figure 1.

Fig. 1. Possible shapes of the curve y(T) (cf. Equation (18)).

Now at the ice divide (R = 1), y(T) reduces to a straight line with negative slope. This means the case a, (if we assume a negative bottom temperature). While flowing off from the ice divide, the curve y(T) deforms in a continuous way, and two solutions are possible:

4. Setting the Problem

During field traverses, measurements which can now be easily donc are:

• the distance from the ice divide X (with tellurometers),

• the surface slope a (by precise levelling),

• the thickness of the ice sheet Z (with high-frequency depth sounders),

• the mean balance during last decades and the present surface temperature S (by coring the firn).

In order to deduce the bottom temperature, the following conditions must hold:

• i. The ice sheet must be in a stationary state; specifically its thickness and temperature may have short-term fluctuations, but no long-term ones. (For polar ice sheets, a fluctuation extending over several centuries remains a short-term fluctuation.)

• ii. In spite of these short-term fluctuations, the measured mean balance and mean surface temperature, which are both mean values for the last decades, must equal the means for the last millennia. In this case, the vertical velocity of ice at the surface A of the theory equals the measured balance (in ice height equivalent).

With the assumption of a steady state, Equations (11) and (12) allow us to compute the bottom temperatures from measurable quantities, provided the constants G ₀ and B are well known.

Unfortunately B (defined by Equation (10)) remains an unknown quantity. Laboratory experiments give B ≈ 0.17 bar⁻³ year⁻¹ for temperate ice, and B ≈ 1.0 bar⁻³ year⁻¹ for cold ice (Reference LliboutryLliboutry, 1964–65, p. 87–89). With cold glacier ice the same value B ≈ 1.0 was found by Reference LandauerLandauer (1959), and by Reference Hansen and LandauerHansen and Landauer (1958). Using Soviet I.G.Y. data, I found B ≈ 1.0 too for the Kupol Dzheksona (Jackson ice cap), in Zemlya Frantsa-losifa (unpublished). Nevertheless at the very bottom of a great ice cap, B may reach bigger values, owing to a very strong ice fabric, as suggested by Landauer’s data at Red Rock (see Reference LliboutryLliboutry, 1964–65, p. 86).

Thus B must be inferred, from another independent equation. This equation comes, in the steady state, from equating the discharge of ice to the nourishment over the distance from the ice divide. To do this, the flow lines of the ice sheet must be well known. It is generally assumed that they follow the direction of the maximum surface slope, but in most cases the slope of the surface as a whole is poorly known. It is a pity that during most ice cap traverses the slope of the surface in every azimuth was not measured at each station.

I shall write first the whole set of equations from which both the temperature profile at a given site and the surface profile could be computed, if B were thoroughly known. Next I shall take an approximate ice-sheet model (plane horizontal bedrock, uniform balance, parallel flow lines in the XY coordinates) for which the computation is simpler. In this case the temperature profile becomes a function of the thickness Z independent of the value of B, and the surface profile allows a very easy estimation of the constant B. It happens that this model fits pretty well the case of the 1959 E.G.I.G. profile in Greenland, and thus a pertinent value of B can be found.

5. Velocity and Discharge of an Ice Sheet

We assume that the bedrock is not too rough at the scale of some kilometres. In this case, the horizontal shear stress approximately equals the effective shear stress, and the horizontal velocity U(z) at a distance z from the bedrock is:

As θ(z) is given by Equation (3), a numerical computation of U(z) would be feasible. As θ and diminish very much above the bottom shear layer, U(z) increases only very little above the top of the shear layer. A quite good approximation of the mean value of U(z) through the whole ice sheet is then U(z ₁) which will he called U. Now from Equation (2)

whence the first integral:

Equation (21) means that the Newtonian energy which is dissipated per unit time and per unit area in the bottom shear layer, (ρgZ) (aU) = fU, equals the increase in the thermal flux as it crosses the shear layer.

By comparing Equations (4) and (21), we find:

Whence it follows that

These equations may be interpreted in the following way. The heat which appears in the bottom shear layer produces at its two limits two opposite thermal gradients ±fU/2K. On these gradients is superimposed a uniform gradient determined by the ice conductivity and the solid convection. In the steady state, the total thermal gradient must be directed upwards everywhere since the bedrock cannot remove heat.

By introducing the ratio R = G ₁/G, already given by Equation (11), Equation (21) may be written

In the thermal regime I, G = G ₀. Next, Equations (11) or (22) lead to the following equation, which may be used to compute U at a given site, whether the balance is stationary or not:

With the assumed law for (Equation (10)) and G ₀ = 1 deg/44m, this may be written numerically:

Let Q be the discharge per unit width, then

In Equations (24) and (25), fU may be substituted by its value given by Equation (27). Then R and T appear as functions of the discharge Q. However, assuming that the ice is in a stationary state, Q may be computed by a quite different way.

Let Y(X) be the distance between two neighbouring lines of flow, as drawn on a map from the maximum surface slopes. Let us define ϵ(X) by:

The mass balance gives

or

If for instance ϵ and A are constants:

If the flow lines are straight and converge at a distance L from the ice divide, and if A is a constant, ϵ = 1/(L−X) and

In the general case, ϵ and A are known functions of X, and Q may be computed step by step.

6. Simultaneous Computation of all the Variables of X

Let ζ(X) be the known altitude of the actual bedrock under an existing ice sheet (otherwise isostasy must be taken into account). Then the longitudinal profile Z(X) can be deduced. This, and Equations (29), (24), (12), (11) and (9) give the following set:

Equation (32-5) may be replaced by the following one, deduced from (32-3), (32-5) and (32–6):

This set may be solved numerically step by step, with the following starting values:

The procedure stops either when T = 0, or when dy,/dT = 0 which means that the temperature gets unstable. Thus at each step the following quantity, given by Equation (19) must be computed:

In order to watch for the appearance of the unstable root of y(T) = 0, the following quantity may also be computed:

The computation has been done for the following model for the ice sheet:

• ζ = 0 (a plane horizontal bedrock)

• ϵ = 0 (parallel flow lines in the (X, Y) plane)

In this case, Equations (32-1) and (32-2) reduce to:

Next it can be proved (see Appendix) that T, Z and dy/dT become functions of R independent of the poorly known constant B; the same applies to (B/X ⁴),(Bf ⁴) and aX.

Numerically, the set (32) becomes:

At each step, the 5 variables Z, R, T, f and a must be computed by an iterative process. Now, given fixed values for R and T, it can easily be seen that the equations giving, Z, f and a would lead to numerical instability. This property persists for the whole set of five equations, but the computation has shown that a simple procedure allows stability to be restored: each new approximation off, R and T is averaged with the old value.

Let us write X[N], Z[N], f[X], …, for the values at the step number N, and Z[N,J] f [N,J], … for their Jth approximation. The following quantities must be computed, in a cyclical way:

with the starting values

As Equation (32-6) (or Equation (37-5)) is only valid for a mean slope, a step ΔX equal to H has been taken. In spite of such large steps a more refined formula than (38–1) is unnecessary, as a varies very slowly with X. The accurate critical values for which y(0) = 0 or for which dy/dT = 0 were obtained afterwards by linear interpolation.

7. Numerical Results and Discussion

The computation has been done with the I.B.M. 7044 computer of the Centre de Mathématiques Appliquées of the Université de Grenoble, for three ice sheet models, which correspond more or less to the Greenland ice sheet, to east Antarctica and to Vatnajökull (Iceland). For each model, the computation has been done with three plausible values of B. The constancy of the critical values of Z, R, T, Bf ⁴, B/X ⁴ and aX for the different values of B was used to check the computation.

With the 3 models, the evolution of the temperature was the one referred to as “evolution II” in the third section. The characteristics of the models, i.e. the assumed values of H, A, S, as well as the critical values independent of B are given in Table I. The subscript 2 refers to the first values for which y(0) = 0, with T ≪ 0 (case a₂ b of Figure 1: a second unstable root for T appears), and the subscript c to the second values for which dy/dT = 0, with T ≪ 0 (case bc₂ the temperature instability appears). The critical values which depend on B are given in Table II.

In Figure 2, R, T and −dy/dT = 0 (which may be called the “stability index”) are given as functions of Z for the three models. In Figure 3 the profiles computed with the “Greenland model” for B = 1 and B = 4 are compared with the profile Station Centrale–Station Jarl-Joset surveyed in 1959 by the International Glaciological Expedition to Greenland (E.G.I.G.) (Reference HofmannHofmann, 1964; Reference MälzerMälzer, 1964). This profile, which follows approximately the line of maximum slope, is very symmetrical and fits quite well between the two theoretical profiles. This would not have been the case with the profile Station Centrale–Terme Nevière (Reference Joset and HoltzschererJoset and Holtzscherer, 1954; Reference TschaenTschaen, 1959), which does not follow the maximum slope on the western side, and is perturbed by an underlying mountain 500 m high on the eastern side. Figure 4 gives the location of these traverses.

Fig. 2. Bottom temperature (T), reduced thermal gradient at the top of the shear layer (R), stability index for the bottom temperature (−dy/dT), reduced abscissa(XB^−1/4) and reduced friction (fB^1/4) as functions of Z.

The curve (Z, XB^−1/4) is the surface profile of the ice sheet, Z being magnified 250 times when B = 1.

• Model I: H = 3 000 m, A = 0.07 m/year, S = −58°C.

• Model II: H = 3 200 m, A = 0.32 m/year, S = −28°C.

• Model III: H = 800 m, A = 2.60 m/year, S = −10°C.

Fig. 3. Surface profiles for model II, and actual E.G.I.G. profile. Vertical exaggeration 500 tires.

Fig. 4. Map of the Greenland traverse area.

The computation with model II and B = 1 gives X _c = 74.39 km for H−Z _c = 149m. On the E.G.I.G. 1959 profile, the same value of (H−Z _c) is reached at two points 180·6 km apart. In order to have X _c = 90·3 km, as B/X _c ⁴ remains constant, it follows that B = 2.18 bar⁻³ year⁻¹. This value is very credible, as it refers to very old bottom ice.

According to these calculations the steady-state bottom temperature becomes inapplicable 29 km before Station Centrale, the velocity being there about 9·5 m/year. It would take at least one thousand years for the ice to travel this distance, and it seems plausible that the bottom temperature would not have been raised to melting point during this time.

The problem of what happens when the bottom temperature is no longer in a steady state remains unsolved. I suppose that the warming is not regular over the whole width of the ice sheet. Even with a horizontal plane bedrock, the ice flow may divide into definite currents of ice with the bottom at melting point, and stagnant masses with a cold bottom between these currents. This “kinematic” instability will be enhanced by the bedrock topography. In every valley of the bedrock, the ice is thicker, the bottom warmer, the instability reached earlier, and this process is a self-increasing one. (The valley may even be deepened by the erosion which takes place.) The surface temperature effect referred to at the beginning of Section 2 may also become important and counteract the basal warming.

All these processes show how crude were the theories which tried to give a single formula for the whole profile of an ice sheet, from edge to edge. The present more elaborate theory must be substituted for these earlier simple models and may serve as a guide for future research in the field.