Introduction

This analysis was prompted by the need to estimate basal temperatures beneath the ice caps which covered Britain in the last Devensian period, particularly to determine the areas in which melting was likely to occur. Reference Budd,, Budd,, Jenssen and Radok.Budd and others ([1970]) tackled a similar problem for the Antarctic ice sheet numerically, but this procedure is not always convenient. Reference Robin,Robin (1955) gives a formula for the temperature profile in the absence of horizontal velocity variation, which Reference Weertman,Weertman (1968) modified in an attempt to fit it to the measured temperature profile at Camp Century.

In this paper the problem of heat transfer in moving ice sheets is reconsidered ab initio. In common with other authors the following assumptions are made:

• 1. The ice sheet is in a steady state. That is, accumulation rates, mean surface temperatures, and surface geometry of the ice sheet do not vary with time.

• 2. The horizontal ice velocity is uniform with depth. An immediate consequence of this assumption is that the vertical ice velocity is proportional to the distance from the bed.

• 3. There is no internal heat generation in the ice sheet. At the bottom of the ice sheet we have a constant geothermal heat flux and a frictional heating term proportional to the horizontal velocity.

Two basic models are considered. In the first the ice is assumed to be spreading from a central ridge line with stream-lines which are parallel in plan (see Fig. 1.), while in the second (Fig. 2.) the flow diverges uniformly from a central origin.

Fig. 1. Parallel flow model: streamlines.

Fig. 2. Radial flow model: streamlines.

The first model will be treated in some detail, while the essential features only will be given in the second case. Since these latter results involve confluent hypergeometric functions, which are not normally tabulated, a short table of values is included.

Parallel Flow Model

Using X as the horizontal variable in the direction of ice flow and Z as the vertical variable, the heat-convection-diffusion equation is

where K is the thermal diffusivity (36.3 m²/year for ice), T is the temperature, and U and V are the horizontal and vertical velocities respectively.

If the ice thickness is H and the net surface accumulation rate is À, both of which are initially assumed to be independent of X, it follows that U = ÀX|H and V= —ÀZ|H. If we substitute these values into Equation (1) and scale the variables by a factor β so that we use non-dimensional co-ordinates

where β = A/2KH, we obtain

for which a general solution can be written in the form (Reference Pillow,Pillow, 1964)

where H _n (x) is the nth Hermite polynomial, and ierfc (z) is the nth iterated integral of the complementary error function.Footnote *

The constants  _n and B _n are found from the boundary conditions. At the base, the ice is heated by the geothermal heat flux and by frictional heating. This yields

where r_g is the geothermal heat-flux gradient, rb is the basal shear stress (constant), and k is the ice conductivity. In non-dimensional co-ordinates this becomes

For the upper surface z = β? ( = h, say) we assume initially a constant temperature T_o.

Applying these boundary conditions to Equation (3) we obtain

The first two terms are the Robin equation, and the sum represents the correction to this to account for the movement of the ice.

While this equation is theoretically valid, it is difficult to use computationally. Some information can, however, be obtained. At the base of the ice sheet the temperature is

The second of these sums is O (exp (_ h2)) , while the first reduces to

which takes the value 1/ τ at x = o and is asymptotic to |X|/ τ1/2 as |X| tends to infinity. This suggests that we can usefully assume that the temperature field is linear in x. This assumption is equivalent to ignoring the effects of horizontal diffusion of heat compared with the horizontal convection of heat. As the above results indicate, this assumption is likely to break down for small x and small h, that is when the convection velocities involved are themselves small. Even so, the difference between the two solutions is negligible. For x = o, h = 2.0, the sum in Equation (5) takes its maximum value of 0.64 when z = 0, so that the discrepancy is at most 0.64 πb K/k≈ 0.015 deg.

In addition to neglecting horizontal diffusion the upper-surface boundary condition will be varied to include a variable surface temperature. In keeping with the linear approximations made above, we assume that the surface temperature is given by T = T ₀ +T _gX for X > 0, i.e. that T = T₀ +Tgx/β for x > 0. We now have the equation

with boundary condition

whose solution is

This solution can be reduced to tabulated functions (Reference Abramowitz, and Stegun,Abramowitz and Stegun, 1964) by using the formulae

so that Equation (6) becomes

where the first term is again the Robin equation, the second term is the contribution from the velocity field, and the final term arises from the variation in the surface temperature.

Graphs of the expressions in square brackets in Equation (7) for selected values of h are given in Figures 3a. and 4a.

Fig. 3a. Velocity-dependent température gradients, parallel flow.

Fig. 4a. Surface température variation effect, parallel flow.

This solution can now be modified to account for variations in net accumulation, surface temperature, and ice thickness.

If the local accumulation rate is A(x), then the vertical velocity at the surface is —A(x) = —A(x) + U(x) dH|dx, so that A is replaced by A(x) and the parameter β by β(χ) = (Â(x)/2KH(x))_1/2.

The factor 2 τ_bKx/k which comes from the basal temperature gradient is equal to τbU(x)/k β,so that, finally, we have for the temperature field

Note that all these quantities, with the exception of T(0, H(0)), the surface temperature at the centre of the ice field, are local. Hence the temperature profile can be estimated without knowing values for intermediate sites.

Radial Flow Model

If we denote the radial variable by R and the vertical variable by Z and again scale with respect to β =A/2KH, the solution corresponding to Equation (3) is

where γ= βR and L_n is the _nth Laguerre polynomial, while that corresponding to Equation (7) is

where ϕ(Z) is the confluent hypergeometric function

and

A short table of values of φ and ψ is given in the Appendix.

Finally, in terms of the physical quantities we have

Comparisons and Comments

Figures 3b. and 4b present graphs of (ψ (h) φ(z)/φ(h)—ψ(z)) and (ι—φ(z)/φ(h)) for the same values of h as used in Figures 3a. and 4a. The graphs in Figure 3 represent the effects of heating due to the motion of the ice. It will be seen that there is more heating in the case of radial flow than in that of parallel flow and that the heat penetrates further into the ice. This is because, for a fixed accumulation rate, the ice must travel twice as far in the radial case as in the parallel case to reach the same velocity. Hence both frictional heating and penetration are enhanced in this case.

Fig. 3b. Velocity-dependent temperature gradients, radial flow.

Fig. 4b. Surface temperature variation effect, radial flow.

On the other hand, Figure 4 shows that the radial flow is less sensitive to temperature variations than the parallel flow. This is to be expected since the amount of ice at the local temperature is proportional to the distance from the origin in the radial case but is constant in the parallel case.

Figures 5a. and 5b show the effects of the velocity and temperature difference on a typical Robin profile for the parallel and radial flows respectively. The values have been chosen arbitrarily as β = o.oo1 m⁻¹, h = 2.0, T₀—T = 10 deg and u = 10 m/year. In both cases the temperature initially decreases with depth before rising again as the bottom is approached. The temperature inversion is more pronounced in the case of parallel flow and also penetrates further. It is immediately seen that the temperature inversion is a consequence of the variation in surface temperature. In Figure 5c, we have the case of parallel flow for β =?.??? 5, h = 1.0, T₀—T= 10 deg, and u = 10 m/year. This corresponds to Figure 5a with the accumulation rate cut to a quarter of its previous value. The temperature variation is now monotonic, showing that the inversion will normally only occur in regions with high accumulation rates, or where the surface temperature differences are large. Also, the temperature inversion is more likely to occur in parallel flow than in radial flow.

Fig. 5. Comparative graphs. 1: Robin's equation; 2: Robin’s equation with velocity dependent correction; 3: Full solution (a) parallel flow, β = 0.001 m¹, (b) radial flow, β = 0.001 m,⁻¹, (c) parallel flow, β = 0.0005m^-1τ.

Applying the formula for the case of parallel flow to the “Byrd” bore hole, very good agreement between the theoretical and actual temperature profiles is obtained using À = 0.30 m/year, u = 29 m/year, and T(o, H(o)) = — 3O.7C, as can be seen from Figure 6. It should be noted however that these values for the parameters were arrived at in an ad hoc fashion since the author has been unable to obtain the correct figures.

Fig. 6. Comparison of measured and computed temperature profiles for the "Byrd" bore hole. (Measured values are denoted O.)

Conclusions

While the approximations derived in this paper are not universally applicable, they should prove useful in many cases where an estimate of the temperature field of an ice sheet is required.

The implicit requirement that A be positive restricts the application of the method to regions above the equilibrium line, and the surface temperature should be a monotonic function of distance from the centre. Since the surface temperature usually increases with height above sea-level, this latter requirement is normally satisfied.

The formulae have already been used by Reference Boulton,, Boulton,, Jones, Clayton, Kenning and Shotton,Boulton and others (1977) to estimate bedrock temperatures and regions of basal melting under the Quaternary British ice sheet, and these estimates seem to agree with the features measured in the field.