The Mathematical Model of Ice Sheets and the Calculation of the Evolution of the Greenland Ice Sheet

Abstract

An evolutionary mathematical model of ice sheets is presented. The model takes into account the basic climatic and geophysical parameters, with temperature parameterization. Some numerical data derived from experiments on the Greenland ice sheet are received. At present the Greenland ice sheet is found to be in a state essentially different from a stationary one corresponding to modern climatic conditions.

Résumé

On expose un modèle mathématique de l’évolution des calottes polaires prenant en compte des éléments climatiques et géophysiques avec paramétrisation des températures. Les applications numériques à l’inlandsis groenlandais montrent que ce dernier est dans un état très différent de l’état stationnaire qui correspondrait aux conditions climatiques actuelles.

Zusammenfassung

Das mathematische Entwicklungsmodell für Eisdecken wird dargelegt. Es berücksichtigt die grundlegenden klimatischen und geophysikalischen Parameter und deren Temperaturabhängigkeit. Ergebnisse einiger Berechnungsversuche für das grönländische Inlandeis werden vorgelegt. Derzeit befindet sich das grönländische Inlandeis in einem Zustand, der sich vom stationären Zustand wesentlich unterscheidet, entsprechend den herrschenden klimatischen Bedingungen.

Introduction

It has now become fairly urgent to achieve a reliable and scientifically valid forecast of the mechanics of large ice masses, both on account of Man’s anticipated active use of currently glacierized regions and in order to add to our knowledge. A scientifically valid forecast of the changes in size and temperature of continental ice sheets will enable us to predict global variations in the overall system of interactions between glaciers, climate, and lithosphere, and to carry out palaeoreconstructions of old glacierizations.

Recent progress in computing technology has made it possible to carry out calculations for modelling complex natural processes. Thus, there has arisen a desire for numerical modelling of large ice sheets. The mathematical theory permits us to specify and analyse the effect of external factors, such as climate and thermal flux from the Earth’s interior, on the character of variations of glacierization, as well as to throw light on the relationship with regional tectonics.

The simplified case of an axisymmetric ice sheet on a flat base with ice viscosity dependent on temperature according to the Arrhenius law was discussed by Strukov (1976). Later there have been proposed either isothermal (Verbitskiy, 1981) or stationary (Verbitskiy and Govorukha, 1979; Verbitskiy and Chalikov, 1980; Verbitskiy, 1981) models of ice-sheet glaciations. In all these models the geothermal flow was assumed to be constant over the whole area.

The response of an ice sheet to seasonal and climatic changes has been discussed by Weertman (1961), Nye (1960), Morland and Johnson (1980), Hutter ([^c1983]), and others. As a rule only the response of an isothermal ice sheet to changes of the ice balance at the surface was studied.

A model of a large ice sheet, based on the mathematical model of a three-dimensional non-stationary glacier (Grigoryan and Shumskiy, 1975) taking into account the basic climatic and geophysical factors was presented in a simplified form by Krass (1981) and extended by allowing for an advective penetration of heat through the depth, by Krass (1983).

A forecast of thermo-hydro-mechanics of a large ice sheet has for its purpose:

• (1) specification of the relationship between the development and existence of glaciation and regional tectonic and climatic conditions;

• (2) determination of the principal stages in the evolution of glaciation by selecting determining parameters with respect to its present state;

• (3) description of the modern tectonic pattern concealed beneath the ice and of its relationship with the tectonics of the surrounding areas.

Hence the principal requirements for the model are the following. First, it must be evolutionary; in particular, the various stages with a stationary state should be derived from it (rather than postulated) as being those established in a time of stable conditions. Secondly, the calculated characteristics must explain the main observed facts and their relationships with external effects. Thirdly, with respect to the overall confidence level of initial information and to forecast scale, the model must as far as possible be simple and easy to analyse. These considerations being taken into account, we shall proceed below to elaborate such a model. Furthermore, the mathematical model of a large ice sheet has served as a basis for an attempt to model the dynamics of the Greenland ice sheet, since sufficiently full data are available for it, i.e. surface and bed relief, accumulation and ablation at the surface, and temperature at the base of the active layer.

Basic Equations

A mathematical model of the evolution of large ice sheets is constructed according to the latest papers on the mathematical theory of glaciers (Grigoryan and Shumskiy, 1975; Krass, 1981, 1983; Buyanov, 1983). The rheological law of ice flow is ordinarily assumed to be a power law

(1)

where K and n are constants, T is dimensionless temperature normalized with respect to the ice-melting temperature T_m ,

are deformation-rate tensor components, σ _ik are stress-tensor deviator components, σ is the shear-stress intensity, E is the activation heat of deformation and mechanical relaxation, and R is the universal gas constant. Arrhenius-type relations are known to fit data poorly between 263 and 273 K. Smith and Morland (1981) and others proposed the relations with time factors, but experimental evidence is insufficient (Hutter, [^c1983]). In order to keep the model simple, we here assume the validity of Equation (1) for all temperatures up to the ice-melting point, 273 K. We shall further make use of a Cartesian coordinate system with the Oz axis directed vertically upwards, the xOy plane corresponding to sea-level. Since ice sheets obey the relation

(2)

where H is a characteristic thickness and L a characteristic horizontal dimension of the ice sheet, the expressions for horizontal velocity vector components have the form (Grigoryan and Shumskiy, 1975)

(3)

Here, h1 (x, y) is the z-coordinate of the bed, Z(x, y, t)is the z-coordinate of the upper glacier surface, v _i (h ₁) is the sliding velocity of the ice over the bed; x, y, z are dimensionless coordinates normalized with respect to H, and

are functions of the inclination angles of the bed and the surface. To the first approximation

We shall further everywhere assume that n = 1 (Newtonian medium), which is acceptable as an approximation to Equation (1) when considering ice-sheet dynamics. The sliding velocity of the ice over the bed is one of the most difficult unknown parameters to define. Special theories of sliding have been developed to set a connection between the sliding velocity of ice and the shear basal stress (Kamb, 1970; Nye, 1970; Morland, 1976; Lliboutry, 1979; Weertman, 1979; Hutter, [^c1983]; Shumskiy and Krass, 1983). For simplicity the sliding velocities are supposed to be zero.

To describe the evolution of a dimensionless thickness of an ice sheet

we shall make use of the equation

(4)

which can be derived from the continuity equation. Here v _h is the rate of change of the mass balance at the glacier surface and v ₀ is the rate of sinking due to bottom melting. The quantity v _h (x, y, t) is characteristic of the climate.

To determine the temperate field we have the equation of thermal conductivity or heat transfer (Grigoryan and Shumskiy, 1975)

(5)

where ∇, ∆ are, respectively, the Hamiltonian and Laplacian operators; N is the heat-generation parameter; τ = t/t ₀ is dimensionless time; t ₀ = μ/2ρgH; μ and ρ are the viscosity and density of ice; g is the acceleration due to gravity, and k is the thermal diffusivity.

Certain conditions have to be added to Equation (5), as follows.

Near the glacier surface as a result of the effect of climate beneath the layer in which seasonal temperature variations attenuate (its thickness being about 10 m) there is formed a temperature corresponding to the long-period climate variations. We shall assign this temperature to be a boundary condition at the glacier surface as a function of areal coordinates and of time, i.e.

(6)

The condition (6) is one of the basic relations interconnecting the ice sheet with climate.

The second boundary condition determines a geophysical relationship between the ice sheet and the thermal regime of the Earth’s interior, the geothermal flux affecting the ice sheet at its lower bounding surface. This condition under low inclination angles of the bed takes the form (the surface of the bed being assumed not to vary with time)

where all the quantities with the sign ~ are dimensional, q _g is geothermal flux, and λ is the thermal conductivity of the ice. In dimensionless form, this condition can be written down as follows:

(7)

If we allow for the possibility of bottom melting once the ice temperature has attained the melting point at the bed (barring heat generation in sliding), the condition (7) can be written down in the form

(8)

The second term on the right-hand side of Equation (8) determines the heat dissipated in the phase transition.

Equations (4) and (5) being evolutionary, it is necessary to assign conditions describing the initial state:

(9)

The problem of solving Equations (1) and (3) to (9) is a closed one. Its effective solution even by numerical methods is possible only for particular cases, such as for instance for axisymmetric dome-like glaciers.

Formulation of the Mathematical Model of the Ice Sheet

In the case of ice sheets, the ratio of characteristic scales of the phenomenon δ has the order δ ≈ 0.001 i.e. their relative thickness is extremely low. In studying the evolution of ice sheets particularly, great interest is devoted to their reaction to variations of the main geophysical conditions. In other words, we are above all interested in the pattern of variation of the ice sheet in plan, whilst the structure of the temperature field in a relatively thin ice layer is an effect of a higher order of smallness. The main thing is to take account in the first approximation of the dependence of temperature distribution in the ice thickness upon outside effects (geophysical relationships): climate, Equation (6), the assignment of V _h, and geothermics Equation (7). Accordingly, it is proposed that we avoid solving the intricate temperature problem of Equations (5)–(9) associated with the kinematic Equation (4), but rather to achieve the parameterization of the model by a priori assigning the temperature of the ice mass and by assigning its profile along depth, account being taken of the boundary conditions (6) and (7). Such a parameterization was discussed in the monographs by Krass (1981, 1983).

Let us introduce the three-step dependence of temperature upon depth by

(10)

Here

is the penetration depth of the surface temperature T _s due to the vertical ice mass-transfer. The value of h_* is determined from an “empirical” formula (empirics being numerical modelling of numerous variants),

(11)

where

c being the specific thermal capacity of the ice, and λ its thermal conductivity. The value of C ₁ = 0.06 has been determined from the statistical processing of the results of a non-linear unidimensional thermal problem allowing for non-linear mass-transfer and heat generation, as well as for the phase transition ice–water (Krass, 1983). a ^* _h is to be computed according to the rule:

(12)

Here a _h is a dimensionless rate of mass exchange or balance (ablation plus accumulation) at the ice-sheet surface.

In the parameterization (10), the temperature T is a function of the coordinates x, y, and time τ, since the parameters entering in the equation are dependent upon these variables. By way of normalizing, the value of T = 1 corresponds to the attainment of the ice melting point. We shall now explain the other quantities entering into Equation (10). In a monograph by Krass (1983), a parameterization of the type of Equation (10) comprises only the parameter of geothermal flux Q _g. Generally speaking, the ice sheet stretches with a sufficiently large ice thickness and is characterized by a significant dissipative heat release, so that the temperature gradient in the bottom ice layer may be much more than the geothermal one (Budd, 1969; Krass, 1983). Therefore, we shall consider the dimensionless gradient of temperature in its linear dependence upon depth in Equation (10) as a sum

(13)

where Q is the temperature gradient, Q _g being the contribution from the geothermal flux at the glacier bottom and Q _d being a temperature gradient averaged over thickness due to dissipative heat generation (Budd, 1969), which in a dimensionless form is represented by the expressions

(14)

where α is the inclination of the surface

to the horizontal, and J is the mechanical equivalent of heat.

Equations (10) involve h _m – the z-coordinate of the upper boundary of the bottom ice layer at the melting point (bulk melting zone).

The value of h _m will be determined from the parametric temperature assignment in Equations (10) and the condition T(h _m) = 1:

(15)

Depending on the values of h_* , T_s, and Q, the value of h _m may be either superior or inferior to h ₁. This corresponds either to the presence of absence of a bulk melting zone, e.g. in the case of h _m < h ₁ the temperature T < 1.

We shall now introduce one more simplification of the model. The dimensionless temperature T in the mass of the ice being different from unity by not more than 0.25 deg, the following relationship is valid

so that

(16)

An error due to the substitution indicated in Equation (16) does not amount to more than 50%, which is quite acceptable for the purpose of mathematical modelling of such a large-scale phenomenon as the evolution of an ice sheet.

Substituting Equation (16) in Equation (3), with the condition v _i (h _j ) = 0 (basal sliding of ice is supposed to be negligible), we shall then perform the integration and afterwards substitute the expressions obtained for velocities into Equation (4). As a result we obtain for the thickness h of the ice sheet

(17)

where v ₀ is the dimensionless rate of bottom melting. Should melting in the ice volume be taking place in the bottom layer of thickness h _e = h _m − h ₁ > 0, then according to Equation (8) and the third of Equations (10), for z = h ₁, ∂t/∂z = 0, and v ₀ will be determined from

(18)

For the coefficient a in the equation of evolution (17) we have (according to Equation (10))

(19)

(20)

It can be shown that for h _m = h ₁ the coefficient a, as determined by Equations (19) and (20), retains continuity and smoothness.

Dimensionless parameters of external geophysical relations of ice sheets enter into the coefficients of the evolution Equation (17) in the following manner:

• (a) Climate – through two functions: the distribution of a _h(x, y, τ) of the mass balance at the surface and the distribution of temperature T _s (x, y, τ) of the ice surface. T _s represents the dimensionless temperature at the bottom of the attenuation layer of seasonal fluctuations, which is formed under the impact of long-period air-temperature variations.

• (b) Geothermal flux – through a single function: the distribution of Q _g (x, y, τ). Since geotectonic cycles have lengths some 2–3 orders of magnitude higher than those of climatic variations, it is assumed in the model that the geothermal flux, like the function describing the profile of the bed, is independent of time, and so is considered to be a function of the coordinates x and y. Thus, the vertical gradient Q in the temperature parameterization (10) consists of two components: the geothermal flux parameters Q _g which is independent of time and the dissipative heat-generation parameter Q _d which is dependent on time and is determined by ice dynamics.

At the places of transition of land ice into ice shelves we have modelled the ice discharge (formation of icebergs) with the ice-shelf thicknesses at the sea–land boundary attaining the “maximally permissible” value to be empirically chosen from the observed average iceberg thickness h ₁:

(21)

where Г_* is the sea boundary.

We have to add to Equation (17) the initial condition (9) assigning the initial form of the ice sheet:

(22)

Along with Equation (17) this condition determines the Cauchy problem for a non-linear equation of parabolic type for the unknown function of ice-sheet thickness h(x, y, τ); in this case the position of the land ice-sheet margin Г will be determined by the condition

(23)

The system of Equations (10)–(23) completely determines the mathematical model of ice-sheet evolution (involving a parametric interrelationship between the mechanics and thermal state of the ice) from which we can compute all basic characteristics, including the movement of the glaciation boundary. Certain conclusions can be drawn from a general analysis of the structure of this problem as follows.

A non-stationary solution for the entire xOy plane is available given a total positive balance of solid precipitation at the surface.

In the presence of ice discharge across the coastline or of ablation there is the possibility of a stablized regime with the ice-sheet boundary stable in time.

Scheme of Numerical Solution

To achieve a numerical solution of the Cauchy problem (10)–(23) describing the evolution of the ice sheet, Buyanov (1983) proposed a finite-difference “hopscotch” scheme (Roach, 1976; Gourlay, 1970) for which calculations were made for both test and model to establish its functional reliability. The ice-sheet evolution calculation program “EVGLAC” was realized in the language FORTRAN. Calculations were made on the BESM-6 computer. As a result of preliminary model-estimate calculations, the dimensionless time step Δτ was chosen equal to 500.

Calculation from the initial state was performed initially with a smaller time step, to secure stability of counting and to avoid an initial “driving” of the solution, and is then increased exponentially to the above value. Error estimation was achieved by undertaking a control calculation with a smaller step Δτ = 250; the difference in the results amounted to less than 1%.

Initial Modelling Object. Principal Sources

We chose the Greenland ice sheet as the object of modelling. The initial calculation data – the glacier surface and bed relief, the distribution of the average annual temperature of the glacier surface, and the mass balance at the surface were taken from the books by Bader (1961) and Weidick (1975).

Bader’s work is a summary of research results on the Greenland ice sheet obtained by a number of expeditions prior to 1961. Weidick’s (1975) work contains, apart from a description of the present state of the Greenland ice sheet, a review of researches on the history of the Quaternary glaciation in Greenland. It is to be noted that the degree of trustworthiness of the initial information for the calculations is not great (as is noted, in particular, by Bader (1961)). This particularly applies to the ice-bed centerography. The magnitudes of the average annual accumulation are likewise inaccurately known. The values of the mass budget calculated by different authors differ greatly from one another not only in absolute values but even in sign (Weidick, 1975). Thus, on Bauer’s evidence the annual glacier budget is 84 km³, on Benson’s evidence this quantity is +13 km³, whilst Bader’s estimate lies between +120 and +270 km³; Loewe estimated the budget to be −20 km³ (in 1936) and +158 km³ (in 1964). Such different estimates have for their cause, first, inaccurate information on the magnitude of accumulation at the glacier surface and, secondly, the different evaluations of the area of the ablation region and the rate of ablation. Thus, Bader estimated the area of the ablation region from the position of the firn line, while accepting the ablation intensity as constant throughout the ablation area.

The magnitude of the geothermal flux at the bed beneath the ice sheet is at present not measured at all, its distribution over the area covered by ice is unknown.

Initial Formulation

Numerical calculations were made on a grid of 77 × 37 nodes with a 32 km step for the spatial variables, upon which all of the data presented above was placed. To achieve this, the glacier surface and bed profiles given by Bader (1961) were taken and isolines of the mean surface temperatures and accumulation were used to assign initially a discrete set of parameters. Then the values of the initial distribution of parameters were interpolated into the calculation grid nodes (triangulation process). For ablation we accepted an average value of 110 cm/year throughout the ablation area (Bader, 1961). The function a _h (the mass balance at the surface) was calculated as an algebraic sum of the accumulation and ablation values at the calculation points.

The system of coordinates in plan was chosen as follows: the x-axis is directed from north to south, the y-axis from west to east, the xOy plane corresponds to sea-level. The origin is placed in the left-hand upper corner of the calculation template.

Calculation Variants

To achieve the modelling and forecast of ice-sheet evolution under different climatic conditions we made calculations for several variants in which we assigned different distribution of the mass-balance function at the surface and of the mean temperature at the glacier surface. The distribution of geothermal flux was assumed to be the same for all the calculation variants. The geothermal flux density was assumed to be 13 µJ cm⁻² s⁻¹ beneath outlet glaciers, traced according to the data of Bader (1961) and equal to the normal value of 4 µJ cm⁻² s⁻¹ at the rest of the calculation nodes. This inhomogeneity is introduced into the distribution of geothermal flux because, as shown by Krass (1981, 1983), outlet glaciers which are, in effect, gigantic ice “rivers” in the body of the ice sheet with high ice velocities, coincident with the zones of neotectonic fractures in the Earth’s crust with increased flux of heat from the Earth’s interior. In other words, glaciers and ice sheets happen to be indicators of inhomogeneities in the distribution of the geothermal flux. This was, moreover, confirmed by test calculations by Buyanov (1983).

The various climatic conditions which were adopted in carrying out numerical model calculations for the Greenland ice sheet*; are given below:

1. Forecast of ice-sheet evolution on the assumption of present climatic conditions being stationary.

   At present there is no established relationship between the precipitation changes (Δa _h ) and temperature at the glacier surface (ΔT _s ). Therefore, to be able to throw light on the possible reaction of the ice sheet to a variation in the climatic conditions we calculated some other variants, as follows.

2. Warming-up with a higher mass balance at the glacier surface, as compared to the present balance field a _h .

   The balance rate at the surface was raised by 30%, i.e.

   

   (24)

   The mean annual temperature at the surface T ’ _s was given a linear increment increasing from north to south

   

   (25)
   i.e. in the north the rise in temperature of the ice surface is 5 deg, in the south 10 deg; T _s is the present distribution of temperature.

3. Warming-up, as in the previous variant, but with a decline in the mass-balance by 30%, i.e.

   

   (26)

4. An overall fall in temperature with a decline in the mass balance at the surface by 30%, similarly as in Equation (26) and with the mean annual surface temperature given in the form

   

   (27)
   where ΔT _s is the same as in Equation (25); in other words, cold is increasing from north to south: in the north by 5 deg, in the south by 10 deg.

5. With an unchanged surface temperature we assumed an increase in accumulation in the northern half of the ice sheet by 30% and a 30% decrease in accumulation in the southern half of the island, i.e. we introduced a decrease in the observed difference between the quantities of solid precipitation between the two regions.

   By retaining the present distribution of mass exchange at the glacier surface we calculated two more variants with different stepped distributions of mean surface temperature, as follows.

6. The surface temperature in the northern half of the island is raised by 10° C; in the southern half of the island the surface temperature is decreased by 10° C.

7. An opposite change in temperature conditions is assigned, i.e. a 10° C increase in temperature in the southern half of the island and its 10° C decrease in the northern half.

For all the above variants a calculation of the ice-sheet evolution was carried out through 13 000 years of dimensional time.

The results of the calculations are maps showing the distribution of ice-sheet thickness at the final moment of its calculated evolution, maps showing the distribution of ice-bed temperature calculated using Equation (10) of the model, maps showing the distribution and thickness of melting zones at the glacier bed, calculated using Equation (15); we also determined the variability in time of integral characteristics of the ice sheet: its area, volume, average ice thickness, and temperature with respect to the ice-sheet volume. The principal parameters and their values as used in the calculations are indicated in Table I.

Table I. Basic Parameters

Model Testing

Model testing was carried out for the following calculated variants from the basic hypothesis that the present climatic conditions are invariable. We calculated the ice-sheet evolution:

• (a) without special zones of increased geothermal flux;

• (b) disregarding dissipative heat generation in the mass of the glaciers; and

• (c) with dissipative heat generation averaged over the 64 × 64 km² area. In this case, the value of the additional temperature gradient Q _d in the calculation grid node was computed as an average from its values in this and in eight neighbouring calculation nodes. The results of this calculation were compared with a calculation using the accepted model with present climatic conditions (variant 1).

As revealed by calculations, ice outflow in the zones 4ith increased geothermal flux for the conditions in Greenland is a relatively small part of the total balance; for overall ice-sheet characteristics the results differed insignificantly from the variant with identified inhomogeneities of q _g, the differences amounting to about 2–3%. Such a conclusion is not unexpected: compared to the Antarctic, where 80% of the ice transport from the centre of the sheet to its periphery is accounted for by outlet glaciers, the Greenland ice sheet is “warmer”, the region itself being tectonically more placid (with relatively fewer tectonic fractures).

Results of calculations with an averaged dissipative heat generation and of those without averaging hardly differ from each other. The results of the calculations allowing for dissipative heat generation and those disregarding it differ by not more than 7%. Allowing for (or disregarding) dissipative heat generation does have a significant effect on the temperature at the glacier bed and on the distribution and thickness of melting zones at the bed, this being due to linear parameterization of the temperature profile in the mass of the ice. In the areas with a large ice thickness a relatively insignificant increase of the temperature gradient due to dissipation leads to a perceptible increase of temperature at the bed. Nevertheless, it is to be noted that the temperature at the bed in the area of Camp Century, calculated from the model using present conditions (h ₀, T _s , a _h ) and the adopted distribution of q_g, was found to be −16° C, whilst measurements in the bore hole yielded a value of −13° C (Weertman, 1968). The calculated values of ice-travel velocities at the surface correspond to the observed ones, including those for outlet glaciers.

Results and their Interpretation

The results of calculations of response of the Greenland ice sheet to climatic changes corresponding to the variants mentioned above are shown in the set of illustrations.

Figure 1 shows the variability of integral parameters of the Greenland ice sheet, area S and volume V, as functions of time. The curves in the plots are labelled with numbers which correspond to the numbering of the variants in the text (see above). One general qualitative property of the response of the ice sheet to all supposed climatic changes is observed: the area of glaciation initially decreases as the ice volume increases.

Fig. 1. Time variation of integral characteristics of the Greenland ice sheet: V – volume, S – ice-sheet area.

An intensive decrease of the Greenland ice sheet would take place in the case of a cold spell; with a decrease of mass balance at the surface by 30%; the area covered by ice diminishes by 12% (curve 4). The ice sheet becomes less mobile and the volume of ice increases by 35% approximately. The increase of ice volume along with the growth of the average thickness and temperature of the ice sheet tends to produce a fig. 1 gradual increase of the deformability of ice, as a result of which the ice-sheet area increases again after 1200 years. A stationary phase of the integral parameters V and S and the average parameters H and T (Fig. 2) is due to conformity of accumulation of the ice balance at the surface to the ice discharge across the coastline. During the first 1000 years the area covered by ice would decrease by approximately 11% (from 1.92 × 10⁶ to 1.76 × 10⁶ km²), and then the area of ice sheet would increase up to its present value because of accumulation of ice. The period of approach to the steady state lasts some ten thousand years. In this case the ice volume would increase by 44% (from 2.7 × 10⁶ to 3.95 × 10⁶ km³); average thickness of ice would increase by 42% and average ice temperature by 2 deg. The calculation results of the first variant (curves 1) differ from values of the “cold” variant (curves 4): response of the ice sheet is more rapid because the ice is more deformable in this case. The period of approach to minimum area of glaciation is less than for the “cold” one as is the period of approach to the steady state.

Fig. 2. Time variation of the mean temperature T and mean thickness H of the Greenland ice sheet.

A warm spell of climate influences in particular the ice-sheet area (dashed curves 2 and 3 in Figure 1). The ice sheet expands again over its present boundaries (in some 600–800 years) because the mobility of ice increases. The volume of ice and the average parameters of ice sheet are essentially dependent on the balance of ice on the surface: with an increase of mass-balance rate by 30%, the volume and average thickness of ice increase by 45 and 42%, respectively. With a decrease of mass-balance rate by 30% V and H increase only by 17% and 9%, respectively. The same values occur for variant 1, as for the “warm” variant 2 in 13 000 years, but for a warm spell of climate the steady state is approached more rapidly.

Calculations of ice-sheet evolution under hypothetical changes of climatic conditions have brought to light some general laws in the reaction of the Greenland ice sheet to these variations. As mentioned above, a characteristic feature is the initial reduction of glacier area from 900–1500 years with a simultaneous growth of ice volume (Fig. 1). After this the growth of ice-sheet area takes place until it attains a stationary state. In this case the quantitative variations of the integral parameters and the distribution of ice-sheet thickness are dependent upon specific combinations of climatic conditions. Relatively smaller (curves 2 and 3 in Figure 1) variations of the ice-sheet area, correspond to warm spells, as compared to cold spells (curve 4 in Figure 1) and to a shorter time of the degradation phase. This is accounted for by the higher ice deformability, which is due to the higher temperature and, correspondingly, to a more rapid build-up of ice inflow from the central ice-sheet regions. With a similar temperature regime at the surface (variants 2 and 3), a larger inflow of ice mass to the surface leads to a more rapid growth of ice volume (curves 2 and 3 in Figure 1): the ice sheet is characterized by a larger thickness. A similar effect on ice-sheet dynamics is produced by a lower surface temperature under a similar mass-exchange regime at the glacier surface (variants 3 and 4). A cold spell contributes to a decrease in the outflow of ice from central ice-sheet regions and to a rapid growth of its volume (Fig. 1, curves 3 and 4).

Now let us look at some regional effects of the evolution of the Greenland ice sheet. Maps of the present state are shown in Figures 3 and 4: with the distribution of ice-sheet thickness (Fig. 3) according to Bader (1961), bottom melting stretches are shown in Figure 4a and temperature at the bed in Figure 4b as calculated from the model. Two ice domes appear: the central ice dome with a thickness of about 3000 m and the southern one with a thickness of 2000 m. The bottom zones of temperate ice (Fig. 4a) correspond to these ice domes; the zone of bottom-melting ice which stretches from the central region of the ice sheet to the western coast corresponds to an outlet and is a consequence of the assumed higher geothermal flux (see point 6). The thickness of the bottom temperate ice layer is calculated by Equations (15) and (11). These calculations have a rather qualitative character because the formulae are based on the temperature parameterization of the model; therefore Figure 4a shows zones of ice melting at the bed where thickness of temperate ice exceed 200 m. The ice is frozen to the bed over an area which consists of approximately 60% of the total area of Greenland ice sheet (Fig. 4b), the lowest bottom temperatures occur in the northern and north-eastern parts of the island. As mentioned above, the temperature at the bed in the area of the Camp Century calculated from the model was found to be −16° C, whereas measurements in the bore hole yielded the value −13° C (Weertman, 1968). This fact proves the validity of our parameterization of the model.

Fig. 3. Map of the present thickness of the Greenland ice sheet: 1 – island boundary, 2 – ice-sheet boundary, 3 – ice-thickness isolines marked in metres, 4 – land area.

Fig. 4. Maps showing the distribution of melting zones (a), bottom ice temperature (b), and the mass balance at the surface in the present state (c); (a) melting zones, 1 – thickness of melting ice layer in metres; (b) ice temperature at the bed, 2 – temperature isolines in °C; (c) the mass balance at the ice-sheet surface, 3 – mass-exchange isolines in m/year.

The map of present mass-balance rate v _h on the surface of the Greenland ice sheet is shown in Figure 4c (Bader, 1961, Weidick, 1975). With the preservation of the present climatic conditions (variant 1) the initial decrease of the ice-sheet area (in northern and south-western parts of the island) is due to the presence at the island periphery of zones of sufficiently intensive fig. 3 ablation which are not compensated by the inflow of ice from the central ice-sheet region. However, as ice accumulates at the centre and its temperature rises (Fig. 2) along with the growth of ice-sheet thickness, the deformability of ice also increases, as a result of which the rate of ice inflow into the periphery rises too.

The prognostic map of the ice thickness in 13 000 years is shown in Figure 5. Part of the island surface, which had been free from ice, has again been covered by ice sheet. This “return” of ice cover is very distinct in the south-west of the island, where higher surface temperatures of ice define the greater activity of this part of the ice sheet. The average ice-sheet thickness increases from 1380 to 2000 m, the average ice-sheet temperature rising from −18.5° C to −17° C (Fig. 2, curves 1). The area of maximum ice-sheet thickness (above 3000 m) has been displaced to the north-west. The southern ice dome has been fused with the northern one. Degradation of the ice sheet is quite distinct in the north and north-east of the island and on the south-western coast. An expansion of the ice sheet is observable on the eastern coast, where ice covers all the land area.

Fig. 5. Map of ice-sheet thickness as calculated using the first variant. 1 – ice-thickness isolines in metres, t = 13 000 years after present.

The results of calculations of bed conditions for variant 1 are shown in Figure 6. There is a significant increase in the area of regions in which the temperature of the bottom layer of ice attains the melting point (Fig. 6a). The temperature in regions of the bed with frozen ice would rise (Fig. 6b). This is due to an overall increase of the ice-sheet thickness. The rate of increase of the ice-sheet volume at the initial moment of calculation (present time) amounted to 200 km³/year, which is in satisfactory agreement with the data given by Weidick (1975).

Fig. 6. Map of the distribution of melting zones (a) and bottom ice temperature (b). Calculation using the first variant: (a) melting zones, 1 – thickness of melting ice layer; (b) ice temperature at the bed, 2 – temperature isolines in °C. t = 13 000 years after present.

A map showing the distribution of ice-sheet thickness corresponding to the minimum area of ice is given in Figure 7 (in 900 years approximately, variant 1). The central ice dome would considerably decrease, the southern one would actually disappear. The northern and south-western parts of the island are free from ice, at the same time the ice sheet is expanding into the fig. 4 south-eastern part of the island. The ice-sheet area diminishes by 10% as compared with the initial value. Nevertheless the movement of the isoline for 2000 m of ice thickness is insignificant (see Fig. 3). Therefore the central part of the Greenland ice sheet is more passive than its periphery. The ice-sheet boundary displacement during the first 1000 years is caused mainly by the reaction of the periphery of the Greenland ice sheet of the present distribution of mass balance on the surface (see Fig. 4c).

Fig. 7. Map of ice-sheet thickness corresponding to minimal glaciation area: 1 – thickness isolines in metres. t = 900 years after present.

Let us regard the results of calculations of the evolution of the Greenland ice sheet as a reaction to an assumed hypothetical climatic change. The map of ice thickness in 13 000 years is given in Figure 8 for variant 2: warm spell of climate with increase of mass-balance rate on ice surface (Equations (24) and (25)). The ice sheet becomes more mobile when the average temperature of ice increases due to growth of ice thickness. The northern part of land is free from ice, but ice covers almost all other areas of the island with the exception of a small part in the south-west. An area of more than 70% of the total for the ice sheet is covered by an ice layer with thickness exceeding 2000 m. The central ice-dome area expands three times as much, the southern ice dome is absorbed by the growing central part of the ice sheet and disappears. The distribution of ice thickness is similar to that of variant 1, and this fact proves the general steady state for the Greenland ice sheet. The mass balance on the surface in 8000 years is almost compensated by ice discharge across the coastline.

Fig. 8. Map of ice-sheet thickness as calculated using the second variant: 1 – thickness isolines in metres. t = 13 000 years after present.

Figure 9 shows the forecast map of the ice sheet in 13 000 years time under climatic changes according to Equations (25) and (26): a warm spell of climate with a decline of mass-balance rate on the surface of the ice (variant 3). In this case the temperature and deformability of ice increase due to the effects of the warm spell on the surface only. The ice-sheet area is the same as in previous variants, but the average thickness of ice is essentially less due to this deficit of fig. 6 mass balance on the surface. Spreading of the ice sheet and smoothing take place; the height of the central ice dome decreases by 500 m.

Fig. 9. Map of ice-sheet thickness as calculated using the third variant. 1 – thickness isolines in metres. t = 13 000 years after present.

Reaction of the ice sheet to the climatic change given by Equations (26) and (27) is presented in Figure 10. This is variant 4: cold spell of climate with decrease of mass-balance rate on the surface. In 12 000 years the ice-sheet dynamics would become stationary. The northern part of the island is free from ice too; a western land area of the southern half of the island, which initially was free from ice, is subject to slight spreading. The ice sheet expands into eastern and south-eastern regions which initially were free from ice. The total area of glaciation diminshes slightly. The relatively low temperature of ice provides only a small mobility of the central part of the ice sheet. The central ice-dome area decreases, and the southern ice dome disappears because of the increase of the region of the ice sheet with thickness exceeding 2000 m (up to 60% of total area).

Fig. 10. Map of ice-sheet thickness as calculated using the fourth variant. 1 – thickness isolines in metres. t = 13 000 years after present.

Figure 11 shows the prognostic map of the Greenland ice sheet under the conditions of variant 5: temperature of the ice surface is supposed to be invariable (the present distribution is suggested), but the observed difference between the qualities of solid precipitation between the two regions of the island (northern and southern) decreases. Variations with time of integral ice-sheet characteristics for variants 1 and 5 proved to be sufficiently close-lying. However, the distribution of ice thicknesses for these variants are significantly different. A larger inflow of ice mass to the northern, cold half of the island, and lower inflow if ice mass to the warm southern half, accounts for the fact that the glaciation with > 3000 m thickness increases and further advances north-westerly. The area free from ice in the north decreases. On the contrary, the island area free from ice in the south-west becomes fig. 10. significantly larger. In other words, the ice sheet shows a tendency for an accumulation of the main ice mass in the cold part of its accumulation area.

Fig. 11. Map of ice-sheet thickness as calculated using the fifth variant. 1 – thickness isolines in metres. t = 13 000 years after present.

Figure 12 gives the results of calculations for variants 6 and 7. It shows sections across the Greenland ice sheet approximately along the meridian long. 43 °W: 1 – present state; 2 – the calculation result for variant 1; 3 – the calculation result for variant 6; 4 – the calculation result for variant 7. The figure clearly shows the effect of a variation in the temperature regime at the surface on the ice-sheet dynamics. Higher temperature leads to rising ice deformability, to an intensified ice outflow from central regions towards the “periphery”, to “smoothing out” of the surface profile. The ice dome perceptibly shifts over into the cold zone. It is clearly seen that the calculation result for unperturbed climatic conditions occupies an intermediate position (curves 2–4 in Figure 12).

Fig. 12. Greenland ice-sheet profile at long. 43 °W. 1 – present state, 2 – calculation using the first variant, 3 – calculation using the sixth variant, 4 – calculation using the seventh variant. t = 13 000 years after present.

Conclusions

• 1. At present the Greenland ice sheet is in a state essentially different from stationary corresponding to modern climatic conditions. Since no researches have been undertaken into the preceding history, the reason for the appearance of non-stationarity can still only be pointed out hypothetically. An abrupt change in external conditions had apparently taken place and what we now observe is a process of reconstruction of the ice sheet, and its reaction to this change. The larger part of the ice-sheet margin lies on land and, therefore, the reason for the reduction of the ice-sheet area cannot be a eustatic rise of the ocean level or an epeirogenetic sinking of the land, but must principally be accounted for by an intensified ablation in the marginal glaciated zone.

• 2. At present there have not been identified any connections between temperature variations at the ice-sheet surface and precipitation variations; nor has any reciprocal connection been established between climate and glacierization. Therefore, at the present stage we can only investigate the possibility of predicting the type of research into the reaction of ice sheets to a change in the parameters of external impacts. As such relationships are brought to light, they can be included into the model.

• 3. The model can be used as a unit in the system “glaciers-ocean-atmosphere-land”.

• 4. The model is convenient in making calculations to throw light on the basic laws of the Quaternary glaciation. In this sense it may serve as a filter of geographical and glaciological hypotheses, an instrument for evaluating the trustworthiness of their relationships.