I. Introduction
Concerning the movement of large ice sheets, two limiting cases have been calculated, pure gliding over the bedrock (Reference NyeNye, 1952) and the behaviour of the ice like a Newtonian liquid (Reference NyePhilberth, 1956). These can both be considered special cases of Glen’s law, with the exponent n = ∞ and n = 1 respectively. As has been shown by Reference HaefeliHaefeli (1961[a],[b]) and by Reference PhilberthPhilberth and Federer (1970), an excellent agreement with the real surface profile of the Greenland ice sheet is obtained if the exponent is taken as n = 3.5.
The temperature profile depends on the horizontal and the vertical velocity profiles (Reference RamseierRobin, 1955). These velocities are functions of the shear strain-rate, which is itself a function of the shear stress and the temperature. An exact calculation of the mutual dependence of temperature and velocity profiles would lead to very complicated expressions. Therefore one has to rely on simplified models. In the following we shall present such a model, which is sufficiently accurate and relatively simple.
Since the fundamental calculation of the temperature profile by Reference RamseierRobin (1955), a number of different refinements have been published (Reference Robin.deWeertman, 1961,1968; Reference Hansen and LangwayLliboutry, 1968). But so far the significance of the vertical velocity for the temperature profile has not been taken into account sufficiently. Generally one still uses the simplified assumption that the vertical velocity decreases in proportion to the distance from the bedrock. This assumption leads to a rather imprecise temperature profile. In a recent paper by Reference Dansgaard and JohnsenDansgaard and Johnsen (1969[b]) the vertical velocity Vh has been calculated from a simplified vx -profile by use of the continuity equation for incompressible media. The simplification consists in the assumption, that v x increases linearly up to a certain height from the bedrock and then remains constant up to the surface. If this height is chosen correctly, the derived vh is a good approximation to the true vh. It is difficult, however, to determine the value of this height, the vx -profile is not known from measurement or from the theory. Reference Dansgaard and JohnsenDansgaard and Johnsen (1969[a]) also calculate the temperature profile with their improved function for vh , and obtain good agreement with measurement. In the vicinity of the bedrock the vertical temperature gradient is approximately given by the sum of geothermal heat and the heat of friction, divided by the thermal conductivity of the ice (Reference LliboutryNye, 1951; Reference Hansen and LangwayLliboutry, 1968). Above this lower region the temperature profile has a monotonie curvature. The curvature lies in the region where the conflict between heat conduction from the lower parts and transport of cold ice from above is most pronounced, i.e. where the product of height h above the bedrock and vertical velocity vh of the ice has an absolute value of the order of the diffusivity K of the ice (cf. Appendix A). The vertical velocity vh in this region has a large influence on the temperature profile. Thus the calculations in this paper aim at a more accurate estimate of the vertical velocity in this region.
II. Symbols Used in this Paper
Symbol Units Description
x km Horizontal coordinate, distance from ice divide.
h m Vertical coordinate, height above bedrock.
A ma-1 Long-time average of total accumulation (in ice thickness).
—A m a-1 Vertical downward velocity, measured from the surface.
vx m a-1 Horizontal velocity.
vxm m a-1 Mean horizontal velocity.
v h m a-1 Vertical velocity.
—VH m m a-1 Vertical downward velocity in the immediate vicinity of the surface for a profile which moves with vxm-
y = kG(H-h) Running dimensionless depth parameter.
T = kGH Full dimensionless depth parameter.
H m Total height of ice sheet, H = H{x).
H m Standard total height (2 500 m).
hc m Critical height, where curvature of temperature profile has its maximum.
ox— cry bar Longitudinal stress.
T,T xh bar Shear stresses.
é é0 a-1 Shear strain-rates.
e Base of natural logarithm.
g m s 2 Acceleration due to gravity.
p Mg nr3 Density of ice.
G deg m_1 Real thermal gradient near the bottom.
G0 deg m-1 Standard value of G (1/44 deg m_I).
Gg deg m-1 Geothermal gradient.
Gt deg m-1 Thermal gradient due to heat generation in shear layer.
n Stress exponent.
K m2 a-1 Thermal diffusivity of ice (38 m2 a"1).
t a Age of the ice.
ps bar Hydrostatic pressure.
Q* J mo]-1 Activation energy of creep.
Symbol Units Description
R J mol-1 deg-1 Universal gas content.
k deg-1 Temperature coefficient, (o. 1—0.25 deg-1).
T K Temperature.
TB K Temperature at the bottom.
ΔT deg Real temperature difference between bedrock and a point vertically above.
ΔT0 deg Temperature difference as read From Table I,
S K Pressure melting point.
θ deg Difference between actual temperature and pressure melting point.
α 0 Surface slope relative to horizontal plane.
β 0 Slope of bedrock relative to horizontal plane.
φ(y,γ) Profile function for horizontal How.
ψ(y,γ) Profile function for vertical How.
III. Assumptions
-
The surface slope a and its horizontal gradient ∂ α/∂x are small.
-
The bedrock is horizontal (β = o).
-
The ice sheet docs not glide over the bedrock (vx = o at h = o).
-
The density p is constant throughout the ice sheet.
-
∂T/∂x and ∂G/∂x are very small.
-
The horizontal gradient of the longitudinal stress
is very small. -
The horizontal gradient of the accumulation ∂A/∂x is small.
-
Only the two-dimensional case is considered.
-
All the values are stationary.
is very small.Assumption (2) is made to simplify the calculations, although diese will be approximately valid for small ß and very small ∂ß∂x. If β ≠ o the x-coordinate is parallel and the h-coordinate orthogonal to the bedrock. Tor the case of a circular ice sheet the same values for ∂T/∂h,T and t are obtained if the linear distance x is changed into the radial distance r.
The validity of assumption (6) is a matter of discussion in the vicinity of the ice divide, and also in the outer regions of ice sheets where ice flows (Reference HaefeliHaefeli, 1968) and other types of spatial instabilities (Reference Hansen and LangwayLliboutry, 1968) may occur.
IV. Calculations
We start from the generalized Glen’s law (Reference WeertmanWeertman, 1968):
In the region of interest the temperature interval and the gradient of the pressure melting point are relatively small, so we can use (Reference BuddBudd, 1968; Reference Hansen and LangwayLliboutry, 1968):
Now for T in equation (1) we use the linear expression
The pressure melting point depends on H—h, This dependence being very small, however, we shall neglect it, so that S depends only on x. For the factor G we put numerically the geothermal gradient Gg or, for the case of an additional heat of friction, Gg+Gt. If heat of friction is present the linear form (2) differs slightly from the true temperature profile in the immediate vicinity of the bedrock, because the heat of friction is not formed at the interface between bedrock and ice, but in the lowermost layers; the deviation from the linear profile, however, is small (Reference Hansen and LangwayLliboutry, 1968).
Above a certain height he the linear form {2) ceases to be a valid description of the true temperature profile. Nevertheless we can use Equation (2) for a rather precise calculation of vx and vh in the whole range of h. This is proved in Appendix A.
Using assumption (6) and Equation (B1) (see Appendix B) one obtains, upon integration
As shown in Appendix B, Equation (3) can be inserted into Equation (1) which becomes
In the following, the movement of the ice will not be calculated from the slope x but by means of the mass budget. The first three factors of Equation (4) are merely a function f (x) which is of no further interest. With n = 3 Equation (4) becomes:
Upon substitution of Y = kGH and y = kG(H-h), integration of Equation (5) gives
where v xm is the mean value of vχ and
With
and
Differentiation of Equation (6) with respect to x yields:
or
The third term in Equation (7) is smaller than the second term for all values ofy, Y and x, but in the central region both terms are small compared with the first one. In many cases, especially if H > 2 000 m, the ∂φ/∂H term is only a small (positive) correction to the φ/H term so that we shall neglect it in further calculations.
Equation (7) now reads
We define
or:
where vxm ∂H/∂x is the change of ice thickness per unit time considered for a profile which moves with the mean horizontal velocity vxm, and —v is the vertical downward velocity of ice particles in the immediate vicinity of the surface considered for a profile which moves with the mean horizontal velocity vzm. —vR is measured from a horizontal plane and in the case of β ≠ o from a plane parallel to the bedrock.
Under normal conditions — vxm ∂H/∂x increases a little more than linearly with x, while A can be assumed to be almost independent of x. Therefore < ivH is nearly constant as long as A in Equation (9) is the predominant term, i.e., according to Equation (9a), as long as
Because of assumption (4) the continuity equation can be written
Integration of Equation (8) and using Equation (9) yields:
Where
The equation of stationary heat transport reads:
Under the condition (10) and with assumption (5), 
and 
can be neglected. On the E.G.I.G. profile in Greenland, for example, the gradient ∂TB/∂x is very small (Reference PhilberthPhilberth and Federer, 1970). Equation (12) now reads:
which yields, upon integration
where
From Equation (13) the temperatures at all heights are obtained by integration;
The age t of the ice at any height can be calculated by means of the vertical velocity relative to the surface, i.e. — A ψ (y, Y) (see Equation (11)).
The integrals in Equations (14) and (15) do not have exact solutions (see also Reference RamseierRobin, 1955). The values for T(h) and for the age t were therefore calculated on the digital computer CDC 1604 of the Rechenzentrum der E.T.H., Zürich.
V. Results
The values calculated with Equations (13) and (15) are shown in Tables I and II. In order to obtain the maximum information from Table I the product 
is given for 7 values of Y = kGH as a function of the relative heights h\H for thin (\vHH\ = 75 m2 a--1), medium (|vHH| = 375 and 750 m2
a-1) and thick (|vHH| =1 500 m2 a-1) ice sheets. ΔT is the difference between the temperatures at the bedrock and the height h/H, Go
is the standard thermal gradient taken as 1/44 deg m-1
and Ho
= 2 500 m. In order to find ΔT" for any thermal gradient G at the bedrock and any total height H of a particular point on an ice sheet, the value of Δoo, which is read from the appropriate Table, has to be multiplied by 44GH/2 500 deg-1, so that
Values of |vHH| which are between those given in Table I, a-d, are best interpolated graphically, where one has to note that, by Equation (13), for |vHH| =0, ΔT = Gh and for vHH= ∞: ΔT= o. Column (I) for Y= ∞ gives the temperature profile of Reference RamseierRobin (1955). while column (7) for Y = o represents the temperature profile obtained by Glen’s law e = const, T 3.
Table I. Temperature differences To in degrees between bedrock and the relative height h/H, for Ho. and Go .
Table 2. Age of ice in years as a function of relative heighth/H.the values refer to H/A = 2 500 a
Example 1
If we take Station Jarl Joset, Greenland, where H = 2 500 m, [uH| = 0.3 m of ice a-1, k = o.1 deg -1 and G = 1/30 deg m-1, we have to use Table Ic (|vH H| = 750), column (4;. ΔTo at 1 500 m above ground (h/H = 0.6), i.e. 1 000 m below the surface is 18.9 deg. So the real Δ T becomes
At a height of 1 500 m above the bedrock, i.e. at a depth of 1 000 in, a temperature of —30.0°C has been measured (Reference PhilberthPhilberth, 1970). Therefore the bottom temperature at Station Jarl-Joset is —30.0+27.8 = — 2.2°C, which is slightly below the pressure melting point.
Example 2
At the ice divide (Créte) in Greenland H = 3000 m, |vn| = 0.25 m of ice a-1, G = 1/44 deg/m and k =0.15 deg-1.
Table Ic gives for h/H = 1 and Y = 10.2 an interpolated value of ΔTo = 18.7 deg, so that, according to Equation (16):
which is the temperature difference between the bedrock and the surface. In comparison, Reference RamseierRobin’s (1955) paper gives ΔT= 19.4 deg.
Example 3
According to Reference HaefeliHaefeli {1961 [b]) the center of the Jungfraujoch Eiskalotte has the following values: H = 50 m, |νH| = 1-5 m a-1, G =1/44 deg m-1 and k = 0.15 deg-1, so that |vHH| = 75 m2 a-1 and Y=0.17. In order to find the temperature difference between bedrock and surface we use Table la, column (7) which yields ΔTo = 46.2 deg and ΔT=0.93 deg for the temperature difference over the whole thickness.
Example 4
In order to obtain the age at h/H = 0.1 for a station with k = O.I deg-1, G = 1/50 deg m-1, H = 400 m and A = 0.6 m of ice a-1 we have to use column (6) of Table II which yields 10 960 a. This is the value for H/A = 2 500 a. For the present values, i.e. for H/A = 400/0.6 a the age is
With Nye’s model (column (1)), an age of 1 540 years would be obtained.
Example 5
According to Hansen and Langway (1966), the temperature difference between the bedrock and the surface at Camp Century is Δ? = 11.0 deg. Let us compare this value with the one given by the present theory.
Fig. 1. Measured (by B, L. Hansen) and theoretical temperature profiles for the Camp Century bore hole,
According to Reference WeertmanWeertman (1968); H = 1400 m, k = 0.1 deg-1, G— 1/56 deg m-1 and |vH| =0.36 m a-1; thus |vHH| = 504 m2 a-1 and Y = 2.5. To obtain the ΔT 0 for these values, we have to interpolate in each Table a-d the individual ΔT 0’s between columns (5) and (6) and plot them on logarithmic paper. From the plot ΔT 0 = 24,5 deg is obtained for vHH| = 504. This yields a ΔT = 10.8 deg which is in good agreement with the measured ΔT.
Since in Camp Century there exists a measured temperature profile, we compare in Figure 1 the temperature profile obtained by the present theory with the measurement. Figure 1 also shows the temperature profiles obtained by Reference RamseierRobin (1955) and Reference WeertmanWeertman (1968). For the special case of Camp Century, Weertman used the shear stress and the temperature to calculate ∂x∂h and vx. From our Equation (5) we obtain exactly the same vx profile. But Weertman does not use this vx profile to determine vh for the calculation of the temperature and age profiles. Instead he uses—as did Robin—the simpler equation Vh = VHH/H based on Nye’s theory, which yields Vh values which are too large. Therefore, the temperatures obtained by Robin and Weertman are too low (Fig. 1) and the ages obtained by Nye’s theory are too small Fig 2
Fig. 2. Measured (by Dansgaard and Johnsen) and theoretical age profiles for the Camp Century bore hole. The distance in metres from bottom refers to a total ice value of I 370 m.
Reference Dansgaard and JohnsenDansgaard and Johnsen (1969[a]) have used the measured age profile for the determination of the vx and temperature profiles (see below). The resulting curve does not differ appreciably from our temperature profile.
Figure 2 shows the age profiles. Reference Dansgaard and JohnsenDansgaard and Johnsen (1969[b]) used for their calculation a linearized vx profile, the parameter h of which is chosen such that their age profile is identical with their measured values. For the most important range below 300 m above the bedrock, this requires νx values which differ considerably from the vz profiles given by Weertman and by the present theory. The difference between the measurement and the present theory implies that, for the conditions at Camp Century, the value for the temperature coefficient k must be smaller than 0.1 deg 1.
VI. Discussion
In previous papers the profiles for the temperature and the age have been calculated from vx-functions which have either been special cases (Dansgaard and Johnsen, loöcjfa], [b]) or rough approximations (Reference RamseierRobin, 1955)·In contrast, Equation (5) is a relatively simple function for vx which sufficiently takes into account the influences of temperature and shear stress. Possibly there is no simpler function for vx which yields satisfactory temperature and age profiles for all real values of H, A and G, because there exist ice sheets where the temperature influence is predominant (large kGH) and those where the shear-stress influence is predominant (small kGH) and those in between the two cases. The Tables and the examples show that the temperature and age profiles for these cases differ widely. In the columns (1) of the Tables the Nye flow model ∂x∂h = o) and in the columns (7) the flow according to Glen’s T 3-law are given as limiting cases (kGH = ∞ and kGH = o respectively).
Equation (9) shows that — vH
differs from —A if vxm
and ∞H/∞x are both non-zero. The temperature profile (13) contains 
, the age profile, however, 
.This difference can be viewed in the following way:
The surface receives the accumulation A dt per unit time. Therefore each ice layer in the vicinity of the surface migrates by A dt below the new surface. The increase in age di of an ice layer in the surface region is thus equal to the distance from the surface divided by A. For deeper layers the product of A with the normalized profile function 
has to be taken instead of A.
The situation is quite different in the calculation of the temperature profile. The thermal interaction of the different ice layers depends on their relative separation which is calculated from the vertical velocity v
h. It is unimportant whether the separation between the layers is shortened because of the accumulation A or because of their movement into a region with decreasing total height. Both influences together result in the quantity VH and therefore 
for the calculation of the temperature profile.
The temperature in the upper layers of ice sheets are influenced by short-time fluctuations of the climate. Furthermore, for x = o the upper layers tend to have a small negative temperature gradient in the downward direction (Reference Philberth and Federerde Quervain, 1968, p. 176; Reference WeertmanWeertman, 1968). Because of these reasons, the value of TB calculated by Equation (14)1 becomes more realistic if it is possible to fit the calculated part of the profile to one measured down to a reasonable depth. This is done for the measured 1 000 m profile at the E.G.I.G. station Jarl Joset in Greenland (to be published in Meddelelser om Gr ø nland).
MS. received 3 September 1969 and in revised form. 16 June 1970
Appendix Λ
Above a certain height the linear form of Equation (2) differs considerably from the true temperature profile. As a characteristic height where the deviation becomes important, we lake he which is the height of maximum curvature of the temperature profile. he can be calculated by use of the Equation (12a)
Differentiating with respect to h, the left-hand side becomes zero at the maximum of the curvature:
Using Equation (A1) in (A2) yields:
or
If we use the approximate expression vh = -Ah/H we arrive at the relationship described in the Introduction:
and
(The indices c denote the point of maximum curvature.)
Thus, from the bedrock to he the linear form (Equation (2)) is a very good approximation for the true temperature profile.
Above he the deviation is large. If we want to use Equation (2) in the calculations of the horizontal and vertical velocity profiles, we must prove that, above the temperature does not influence the vx profile to any appreciable extent. In order to prove this, we will show that at he - the shear velocity ∂vx/∂h has decreased by at least a factor 1/e = 0.368 with respect to its value at the bedrock (in most real situations the factor is essentially smaller than this). If this is the case the increase of Pi between he and the surface (h= H) is very small. Therefore we are not committing any important errors by using the linear form (2) which, above he , yields a lower temperature and therefore a larger rigidity of the ice than the true values. larger rigidity above hc, however, further reduces the difference of Vx between A = he and h = H..
Suppose now that the shear strain-rate (Equation (5)) has decreased by a factor less than the. Division ot Equation (5) for h = hc by the value for h = o yields:
From Equation (A5),
If condition (Λ7) is fulfilled, the linear form (2) can be used in calculating the velocity and temperature profiles. It can be seen that (Λ7) is fulfilled for every possible value of he/H if
The magnitude of the error in Equation (5) cannot very well be estimated by using the foregoing considerations. The error is equal to the difference D between the real value ∂vx/∂h and the une given by Equation (5). Tin: function D depends on the true temperature profile, which can be calculated by Equation (Λ1). In this way the following approximation is obtained:
The precise calculations show that the neglected terms of higher order (in the temperature function and in the exponential expression) nearly compensate each other, so that Equation (A9) is rather accurate for all real ire sheets.
Comparing Equation (5) with the functionf(x)(H—h)" exp (-kT), where T is the real temperature profile and n is chosen such that the values of this function are equal in the region where D has its maximum, yields the exponents n given in Table III.
Table 3. Values for n in the function f(x) (H- h)n exp (- kT) corresponding to equation (5)
Numbers in brackets are n-values for unreal ice sheets. Only Antarctica and Greenland have Y values > 3; but for these ice sheets A is smaller than 0.35 m a-1 and therefore A/(kGn) ≤ 3. Thus, for all existing ice sheets H-values between 3 and 5 are obtained.
These estimates need the following comment: Robin’s treatment is the zeroth approximation to the velocity profile, because the temperature is not taken into account and n is taken to be ∞. Equation (5) with the linear temperature profile is the first approximation if the exponent in Glen’s law is taken to be 3. On the other hand, Equation (5) is already the second approximation if it is taken to describe the true physical mechanism with the real temperature profile instead of — Gh and an exponent n taken from Table III.
Since laboratory and field experiments yield an exponent n between 3 and 4, we arrive at the following conclusion: Equation (5) gives a very accurate description of the temperature profile for all ice sheets; in many cases (especially if n ≈ 3,5) the description of the temperature profile by Equation (5) is even better than if one had substituted —Gh by the true temperature profile.
Appendix B
In this Appendix it is shown that, because of assumption (6), it is justified to putTxh
for r in Equation (1) and ∞vx
for 
In recent papers (Reference BuddBudd, 1968; Reference NyeNye, 1969) detailed studies on the stress in ice sheets have been published. Here an intuitive description of the ice flow is given. The ice sheet is considered to consist of a pile of thin horizontal layers. These are subject to the shear stress txh
and the shear strain-rate ∞vx
between the layers and to the longitudinal stress ox - oh
and the longitudinal strain-rate 
within the layers.
At the ice divide 
and τxh are zero; σ
x
–σ
h
and 
are non-zero. Therefore for the ice divide we can state that the longitudinal stress v
xh drives the longitudinal strain-rate 
and the shear stress Txh
drives the shear strain rate 
. This statement is approximately valid for the whole central region of the ice sheet. In order to prove (his fact we consider, as a first approximation, not only T but also 
to be independent of x and consequently the longitudinal stress for the maintenance of 
to be likewise independent of x This last mentioned longitudinal stress and the longitudinal stress ax— an. both being independent of. x and both being identical for x = o we can conclude that they are identical also for x≠0.
There could be an objection: for our proof we have neglected some values which have been called small. For large values of.x these could add up to an amount which could be not negligible. In principle this objection is right. It must be realized, however, that ax— ah
and 
are nearly independent of while rxh
and 
increase with x Calculations have shown that under conditions which exist on the E.G.I.G. profile in central Greenland, for. > 10 H the longitudinal force (i.e. ax— ah
= 2 ax ps
integrated over a vertical cross-section) is smaller than the total shear force (i.e. the bottom value of Txh
integrated over a horizontal cross-section). In consequence, even if an influence of 
exists, this influence remains small with respect to that of Txh-
and T in Equation (1) can therefore be understood as 
and Txh
respectively. In Reference NyeNye’s (1969) equation 
can be neglected because 
where Ps
is nearly independent of x, so that it reads
