Numerical Simulation of the Evolution of Ice Covers Using the Scandinavian and Laurentide Ice Sheets as Examples

Abstract

A non-linear parabolic equation describing the evolution of an isothermal linearly viscous ice sheet is numerically solved in non-dimensional coordinates obtained by normalization over the horizontal size of a glacier. The horizontal size of the ice sheet is defined from the solution of an ordinary differential equation, the integral mass balance. For simple climate models, approximate relations describing the evolution of glaciers are proposed. These relations and palaeogeographical data are used to estimate changes in the mass balance on the surface of the Scandinavian and Laurentide ice sheets during retreat of the last glaciation.

Simulation numérique de l'évolution den couvertures glaciaires en utilisant les calottes Scandinaves et laurentides comme exemples. Une équation parabolique non linéaire décrivant l'évolution d’une calotte glaciaire isotherme à viscosité évoluant linéairement est résolue numériquement en coordonnées adimensionnelles obtenues en les rapportant à l’etendue en plan du glacier. L’étendue en plan d’une calotte glaciaire est définie à partir de la solution d’une équation différentielle ordinaire, le bilan de masse intégral. Pour des modèles sous climat simple, on propose des relations approchées décrivant l’évolution des glaciers. Ces relations et des données paléogéographiques sont utilisées pour estimer les changements intervenus dans le bilan de masse de la surface des calottes glaciaires scandinaves et laurentides durant la période de retrait de la dernière glaciation.

Zusammenfassung

Numerische Simulation der Entwicklung von Eisdecken am Beispiel des skandinavischen und laurentidischen Eisschildes. Eine nichtlineare, parabolische Gleichung für die Entwicklung eines isothermen, linear viskosen Eisschildes wird numerisch in dimensionslosen Koordinaten, die sich durch Normierung über den Grundriss eines Gletschers ergeben, gelöst. Der Grundriss des Eisschildes folgt aus der Lösung einer gewöhnlichen Differentialgleichung für die integrale Massenbilanz. Für einfache Klimamodelle lassen sich Näherungsbeziehungen für die Entwicklung von Gletschern angeben. Diese Beziehungen und paläogeographische Daten werden zur Abschätzung der Massenbilanzänderungen an der Oberfläche des skandinavischen und des laurentidischen Eisschildes während des Rückzugs aus der letzten Vereisung herangezogen.

Mathematical reconstruction of glaciations in the geological past (Kvasov and Verbitsky, 1981) is based on calculations using a thermohydrodynamical model of a stationary ice sheet (Verbitsky and Chalikov, 1980[a], [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 (Weertman, 1964, 1976; Birchfield and Weertman, 1978) and the so-called “transport” (Sergin, 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 (Budd, 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

(1)

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

(2)

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 (Oerlemans, 1981) makes Equation (1) somewhat “poorer”.)

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

(3)

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

(4)

At the edge of the glacier, we assume:

(5)

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

(6)

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

(7)

(8)

(9)

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

(10)

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

(11)

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

(12)

Equation (6) with conditions (7)–(9) is solved numerically using the grid method jointly with Equation (12) (Verbitsky, 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:

(13)

(14)

(15)

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:

(16)

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

(17)

and, finally, Equation (13) yields

(18)

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:

(19)

or

(20)

Hence, setting L = L ₀ at t = 0

(21)

The simple expressions obtained, similar, to a certain extent, to those in Weertman (1964, 1976), 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 (Gerasimov, 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 (Andrews 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 (Kvasov, 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:

(22)

the solution of which is as follows:

(23)

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