Mathematical reconstruction of glaciations in the geological past (Reference Kvasov and VerbitskyKvasov and Verbitsky, 1981) is based on calculations using a thermohydrodynamical model of a stationary ice sheet (Reference Verbitsky and ChalikovVerbitsky and Chalikov, 1980[a], Reference Verbitsky and Chalikov[b]) for prescribed regions of its existence. This model can be used for ice sheets where an important component of mass balance is ice calving, i.e. such as the glaciers of Antarctica and Greenland. For analysis of the history of the Quaternary glaciation it is essential to study the behaviour of glaciers extending on land and changing their sizes. This problem was considered on the basis of simple balance (Reference WeertmanWeertman, 1964, Reference Weertman1976; Reference Birchfield and WeertmanBirchfield and Weertman, 1978) and the so-called “transport” (Reference SerginSergin, 1979) models. This paper is concerned with analysis of the evolution of ice sheets based on solution of complete non-linear hydrodynamical equations.

Here we study the dynamical properties of a large ice sheet using as an example a two-dimensional glacier composed of isothermal linearly viscous ice. If the assumption about the linear ice flow law is reasonably justified for large plain ice covers with characteristic values of stresses not exceeding 10⁴–10⁵ Pa (Reference BuddBudd, 1969), then the isothermal condition is, in effect, determined by the “strong” assumption about the two-dimensional character of a glacier, after which it is hardly appropriate to allow for the “fine” temperature structure.

The equation for calculating the dynamics of the free surface of an ice sheet is the kinematic condition on its upper boundary, which, allowing for the vertically integrated continuity equation, takes the form

where h is a free surface ordinate, U is the horizontal velocity, a is a function giving the mass influx distribution, and x, z are rectangular Cartesian coordinates with the origin at the centre of the glacier base.

Substitution into Equation (1) of the values of U derived from solution of the motion equations obtained using the assumption about relative thinness of the glacier, reduces Equation (1) to the form

where α = ⅔Kρg, K, ρ being respectively the fluidity and the density of ice and g the gravity acceleration. (The assumption that U = Bτ ^m where τ is the bottom stress (Reference OerlemansOerlemans, 1981) makes Equation (1) somewhat “poorer”.)

The initial condition for the non-linear parabolic Equation (2) is:

Requiring the continuity derivative h _x at the origin of the coordinates gives

At the edge of the glacier, we assume:

where L is the glacier horizontal size. As the ice sheet approaches the edge of the continent, a situation may occur when h ≠ 0. But the values of the edge thicknesses of glaciers are small compared with their characteristic thicknesses, which also enables condition (5) to be used with a high degree of accuracy.

The fact that L is dependent on time produces certain difficulties in solving of Equation (2). To overcome them, we introduce non-dimensional coordinates x = Lξ and t = T ₀ τ where T ₀ = H/a ₀ and H and a ₀ are the scales of quantities h and a. In the new coordinates, Equation (2) will be written in the form

(here h and a are nondimensional quantities and β = αH ⁴/a ₀ L ²); conditions (3)–(5) will take the form:

To calculate L, we integrate Equation (1) horizontally from 0 to L:

The expression obtained is the integral mass balance and can be transformed as follows:

where

We divide both parts of Equation (11) by h^⋆L and solve the differential equation obtained with the condition h ^⋆ L = h ₀ ^⋆ L ₀ for t = 0, and obtain

Equation (6) with conditions (7)–(9) is solved numerically using the grid method jointly with Equation (12) (Reference VerbitskyVerbitsky, 1981).

At time t = 0 let a = a ₁ = const. for ξ∈ [0, 1 — ε], ε ⪡ 1, a ^⋆ = 0. When t > 0 a = a ₁ + a _m sin (2πτ/T) for ξ ∈ [0, 1 – ε], a ^⋆ = a _m sin (2πτ/T). Results of the calculations are given in Figure 1. The amplitudes of oscillations of h ^⋆ and L appear to be strongly damped with T decreasing from 10 (T ₀ T = 10¹³ s) to 3 (T ₀ T = 3 × 10¹² s). Changes in h ^⋆ and L at T = 1 (T ₀ T = 10¹² s) do not exceed 200 m and 50 km respectively, for a _m = 2 (a ₀ a _m = 2 × 10⁻⁹m s⁻¹) and are not shown in the figure. As a _m increases to 6 (a ₀ a _m = 6 × 10⁻⁹ m s⁻¹), the oscillations of h ^⋆ and L (curves 5, 6) differ from sinusoidal ones. Here L ₀ = 10⁵ m, K = 3 × 10⁻⁴ Pa⁻¹ s⁻¹, ρ = 917 kg m⁻³, g = 9.8 m s⁻² , H = 10³ m, a ₀ = 10⁻⁹ m s⁻¹, and a ₁ = 3. h ₀ is calculated using the expression

which is the solution of Equation (6) without non-stationary terms where

Fig. 1.

Changes in the mean thickness of the ice cover h* (1,3,5) and horizontal size L (2,4,6) for periodic changes in precipitation. 1,2-a_m = 2; T = 10, 3, 4-a_m = 2; T = 3, 5, 6-a_m = 6; T = 3.

It is interesting that in simple climate models it is possible to avoid a quite laborious solution of Equation (6) using the methods of similarity theory. Consider the glacier motion equations and the continuity equation:

where p is the pressure and μ = ½K is the ice viscosity. We choose the value of L as the horizontal and of h ^★ as the vertical scale, then, from Equation (14) the pressure scale 〈p〉 will be defined:

Let 〈W〉 = a ₀ then, allowing for Equation (15),

and, finally, Equation (13) yields

where k is an empirical constant of order 1. Substituting Equation (18) into Equation (11) gives the approximate equation describing the evolution of the ice-sheet edge:

or

Hence, setting L = L ₀ at t = 0

The simple expressions obtained, similar, to a certain extent, to those in Reference WeertmanWeertman (1964, Reference Weertman1976), allow one to perform approximate calculations of the regime (mass balance) of the Scandinavian and Laurentide ice sheets. Their dynamics has been thoroughly studied for the time of the retreat of the last glaciation which reached a maximum 18000 years b.p.

Consider the retreat of the Scandinavian ice sheet from edge formations during its maximum stage at Vyshniy Volochek to the Salpausselkä ranges north of Vyborg (Reference GerasimovGerasimov, 1973). From 18000 to 11000 years b.p. the sheet edge became 500 km closer to its dynamical centre in the northern part of the Gulf of Bothnia and the radius decreased from 1000 to 500 km. The calculations thus far will not allow for the non-uniform rate of retreat of the glacier which was higher during the inter-stage warmings and slower during the stages. In a range of 18000-12500 years b.p. these warmings were small and brief; but 12500-11000 years b.p. the inter-stages Bölling and Alleröd caused a rapid retreat of the glacier edge. Between 11000 and 10000 years b.p. the ice sheet was in a state close to stationary and then disintegrated rapidly.

The retreat of the Laurentide sheet occurred at a higher rate. Beginning from 18000 to 8000 years b.p. its edge retreated from the Ohio River almost to the southern part of Hudson Bay and became 1200 km nearer to the dynamical centre (in Hudson Bay) and the radius decreased from 2400 to 1200 km. In North America there were no inter-stages similar to Bölling and Alleröd; the retreat of the glaciers continued to edge formations of the Cochrane stage, which were formed 8000 years b.p. Subsequent disintegration of the Laurentide sheet cannot be described using the model proposed here. The decisive factor in this case was the formation of a “calving bay” in Davis Strait and penetration of sea-water into the central part of Hudson Bay, which, under the effect of isostasy, was lowered by several hundred metres. As a result, the Laurentide ice sheet formed a great number of icebergs and disintegrated. This phenomenon is referred to as the collapse of an ice sheet (Reference Andrews and PeltierAndrews and Peltier, 1976).

Calculation of the value of a ^⋆ according to Equation (20) for a linear presentation of L ₀ – Ψt is summarized in Table 1.

Table I.

Calculated ablation rates of ice sheets

Thus, for the period of retreat of the Scandinavian ice sheet (18000-11000 B.P.) we obtained a value for the excess of ablation over accumulation within the range 125–170 mm/year (this value is averaged over the sheet area). This value appears to be sufficient to provide a retreat of the glacier edge of some 70 m per annum, on average. The retreat, apparently, was mainly caused by a reduction of atmospheric precipitation (Reference KvasovKvasov, 1978). The Laurentide ice sheet 18000-8000 years b.p. retreated as a result of an excess of ablation over accumulation amounting to 130-200 mm/year at a mean rate of 120 m per year. Relatively small changes in the mass balance on the surface of the ice sheets led to a rapid change in their sizes. As a result, during the time of the order of ten thousand years the Earth passed from glaciation maximum to the beginning of interglaciation.

In conclusion, we note that in a more general case when γ = γ(t) Equation (19) reduces to an equation of the Bernoulli type:

the solution of which is as follows:

When γ = const. Equation (23) is transformed into Equation (21).