Introduction

One of the striking features of glacial–interglacial cycles is the formation and disintegration of some ice sheets and the apparent stability of others which survive the warm periods. The question of which climatic conditions can start ice sheets growing has been addressed in several ways by using mathematical modelling. The first type of approach uses ice-sheet models with prescribed bed topography and climatic input (Reference Andrews and MahaffyAndrews and Mahaffy, 1976; Reference Budd and SmithBudd and Smith, 1981; Reference Kuhle, Herterich and CalovKuhle and others, 1989; Reference Payne and SugdenPayne and Sugden, 1990; Reference Letréguilly, Huybrechts and ReehLetréguilly and others, 1991). The second type of approach combines global climatic energy balance models with ice-sheet models (Reference Birchfield and WeertmanBirchfield and others, 1982; Reference Deblonde and PeltierDeblonde and Peltier, 1991). The third group uses atmospheric general-circulation models (GCMs) for calculating the ice mass balance at given geographical locations (Reference Rind., Peteet and KuklaRind and others, 1989; Reference OglesbyVerbitsky and Oglesby, 1992). These studies try to simulate the initiation of past and present ice sheets realistically, but the resulting ice sheets also reflect the chosen climatic scenarios, especially the snow accumulation rate, which are subject to uncertainties.

In order to understand the recurrence of glaciations it is necessary to answer the following questions concerning the initiation processes and the growth of ice sheets: given an area of positive mass balance, (1) what is its necessary size in order to act as a seed for the growth of a large ice sheet; (2) how does the final ice-sheet size depend on the snow accumulation rate; and (3) how does the growing speed depend on the accumulation rate?

In this paper we address the problem of ice-sheet initiation by applying a thermodynamically coupled ice-sheet model to a prescribed bed topography and mass balance at its surface, similar to the first group of approaches. In contrast to their attempts to simulate realistic initiation situations, we perform a series of systematic sensitivity studies. A wide range of mass-balance forcing was applied in order to find the critical climatic conditions for ice-sheet initiation.

The Ice-Sheet Model

We used a time-dependent mathematical ice-sheet model for the plane flow (2-D) approximation (Abe-Ouchi and others, in press), mainly following the methods of Reference HuybrechtsHuybrechts (1990) and Reference Esch and HerterichEsch and Herterich (1990). The calculation of the evolution of the free surface takes the surface mass balance into account, and a local isostatic adjustment with a prescribed time lag of 3–5 × 10³ years is used for calculating the current bed topography. The velocity field is calculated by solving the stress balance equation. The shallow-ice approximation is applied (Reference HutterHutter, 1983) but the velocity field is coupled to the evolving temperature field through a temperature-dependent flow-law parameter. No basal sliding is taken into account. The time-dependent temperature field in the ice is calculated, including vertical conduction, horizontal and vertical advection and strain heating. For transient situations, the temperature in the underlying rock is also calculated to ensure that the thermal conditions at the lower boundary remain stationary for the time span in consideration (Reference RitzRitz, 1987; Reference HuybrechtsHuybrechts, 1992).

In order to achieve numerical stability a time step of 0.5–1 year was used for calculating the evolution of the free surface and the bed topography, but five-year steps were sufficient for the temperature calculation. Not only the speed of growth of the ice sheet, but also its final size, can critically depend on the choice of the horizontal grid resolution. No horizontal stretching was applied, therefore advance or retreat occurred in steps of the given horizontal grid width. It was found that a grid size of less than 10 km is safe. Terrain-following coordinates were used (Reference HaltinerHaitiner, 1971) with the number of vertical grid points between 50 and 100, depending on the model run.

At the free surface the mass balance and the ice temperature are prescribed. The idealized surface boundary conditions are parameterized in a simple way with parameters as few as suitable for the sensitivity studies. The surface ice temperature is chosen as a linear function of height with a prescribed lapse rate and a given temperature at sea level. The mass balance is also taken as a linear function of height with a prescribed lapse rate and is zero at the prescribed equilibrium-line altitude (ELA). Above the height where this function exceeds a prescribed positive accumulation, A_c, the mass balance is taken to be constant and equal to A_c.

For the study of the initiation of ice sheets a range of values for ELA and A_c was chosen for varying climatic conditions. In nature, the initiation and growth of an ice sheet also changes the local climate through changes in the albedo (Reference KoernerKoerner, 1980; Reference Birchfield, Weertman and LundeBirchfield and Weertmann, 1983) and in the atmospheric boundary layer and the katabatic-wind regime, but also through changes in the larger-scale atmospheric circulation. However, in the calculations no feedback due to the changing climate above the ice sheet is taken into account. The only resulting climatic feedback works indirectly through the changing topographic elevation. Only the accumulation rate, A_c, and the ELA were varied for the sensitivity studies, and the lapse rate of the mass balance was fixed at 0.005 a⁻¹ in all calculations. The results of this study do not depend critically on the choice of this lapse rate.

In all cases the lower boundary for the temperature calculation was fixed at 3000 m below sea level, where a geothermal heat flux of 0.042 W m⁻² was prescribed. The ice temperature at the free surface was always prescribed with a lapse rate of 0.008 Κ m⁻¹ and −5°C at sea level. Two different topographical situations were used. The first corresponds to an east-west cross section of the Greenland ice sheet near the profile of the Expédition Glaciologique Internationale au Groenland (EGIG) between 70° and 72° Ν. The initial bed topography without an ice sheet was calculated as the isostatically rebound topography of the present bed under the present ice load. The second topography is an idealized single-mountain island of sinusoidal shape ( Fig. 1).

Fig.1 . Initial topography chosen in the model runs: a, EGIG profile; b, idealized plateau-type island. Dashed lines indicate sea surface.

The performance of the model was tested for the EGIG profile, which will be presented by Abe-Ouchi and others (in press). The model reproduced the general features of the present geometry quite well, despite its limitation in reproducing details for the three-dimensionality of the problem. On the other hand, the reduction to two dimensions allows for systematic sensitivity studies over large ranges of the parameters in question and is very flexible in handling and moderate in CPU time consumption.

The Growth of an Ice Sheet: Greenland

An example of the evolution of ice-sheet geometry in 4000 a steps starting with the initiation on the mountain top in the east of the EGIG profile is shown in Figure 2. The chosen surface boundary conditions, with ELA = 1600 m a.s.l. and A_c = 300 mm a⁻¹, correspond to a warmer situation than the present climatic conditions. This shows that the present Greenland ice sheet is not only a remnant of the ice-age climate but it would even start growing in a warmer climate than the present, which is in good agreement with the three-dimensional calculations of Letréguilly and others (1991). The initial growth to almost full size is achieved in about 40 000 a, but for reaching equilibrium in a stationary climate more than 10⁵ a are necessary. This is a consequence of the slow response of the ice temperature.

Fig.2 . Evolution of the ice sheet for the EGIG profile in 4000 a steps. The ice mass starts growing on the mountain top in the east and almost reaches the west coast after 40 000 a. The dashed lines show the bed topography corresponding to the transient surface topography. The ELA is kept at 1600 m a.s.l. and the accumulation rate is 300 mm a⁻¹.

Figure 3 shows that even a small area of about 50 km diameter with a positive mass balance on a mountain top can be sufficient to initiate the growth of an ice sheet. The graph also illustrates the dependence of the growth speed and the final size on the chosen horizontal grid size. In order to be safe, a grid width as small as 5 km should be applied. For small accumulation rates this problem is even more critical. If the grid size chosen is too large, the model not only tends to underestimate the growth speed, but also fails to reach the proper equilibrium state.

Fig.3 . Width of the EGIG cross-section versus time after initiation. This result was achieved with horizontal grid resolutions of 2 km (solid line), 5, 10 and 20 km (dashed lines from top to bottom). With 20 km resolution the model fails to grow the full-size ice sheet.

Multiple Steady-State Solutions

In this section we discuss some aspects of steady-state sizes of ice sheets. Various papers (Reference WeertmanWeertman, 1976; Reference OerlemansOerlemans, 1981; Reference HindmarshHindmarsh, 1990) dealing with the initiation, growth and decay, and stability of ice sheets have pointed out the dependence of steady-state solutions on the initial conditions. Such intransitive behaviour is attributed to the non-linear coupling between climatic boundary conditions and ice-sheet dynamics with largely different characteristic time scales. In order to study the dependence of the ice-sheet growth and its final steady state to different climatic conditions, we performed a series of model runs with a range of values for ELA and

A_c. We used again the EGIG geometry and varied the ELA between 1100 and 2000 m a.s.l. for two values of A_c, 200 and 500 mm a⁻¹, representing a continental and a maritime climate situation, respectively. The steady-state solutions under the prescribed climatic boundary conditions are obtained by calculating the transient ice-sheet evolution starting from some initial state, either from no ice or from a steady-state solution, until the ice sheet reaches the final steady state.

The pattern of the multiple solution branches displays a clear dependence on the accumulation rate. Figure 4 shows the size of the steady-state ice sheets for the corresponding climatic conditions as a function of the ELA, but for different initial conditions. In the first case of the larger A_c at least two solutions with finite ice mass exist for a common range of ELAs: the “small ice-cap branch” and the “large ice-sheet branch”. The situation is different for the second case with smaller A_c. The two branches have no overlapping range for the ELA, but near ELA = 1550 m a.s.l. a steep transition between the two solutions can be observed.

Fig.4 . Area of the cross-section of the steady-state ice sheets versus ELA for two accumulation rates (A_c). Branches a, b and c correspond to large ice-sheet solutions, small ice-cap solutions and no-ice solutions, respectively. The arrows indicate the possible directions of transition from an initial steady-state solution on one branch to a final steady-state solution on another branch.

A third “no-ice branch” can be defined for situations where the ELA lies above the ice-free mountain top. When the ELA touches the mountain top (ELA = 1660 m a.s.l.) we find a sharp transition from the “no-ice solution” to a solution with a finite ice mass. Again we see a different pattern for the different accumulation rates. In the case of large accumulation the solution jumps from the “no-ice branch” directly up to the “large ice-sheet branch”. On the other hand, with the small accumulation rate this transition occurs from the “no-ice branch” only to the “small ice-cap branch”, because the transition between small ice caps and large ice sheets occurs at an ELA which is almost 200 m below the top of the mountain. This clearly illustrates the importance of the amount of annual snowfall for the initiation of ice sheets.

Time Scales of Ice-Sheet Evolution

Glacial cycles are assumed to be triggered by astronomical forcings of 20 000 and 40 000 a periods. Therefore, not only climatic conditions themselves, but also speed of growth, determines whether an ice sheet can grow to full size in the time span available. We performed a series of numerical experiments to investigate the duration of growth as a function of ELA and the annual amount of snowfall. In these experiments we used a steady-state climate and the idealized plateau-type topography as described above (Fig. 1b). The assumed island has a width of 1100 km and an elevation of 800 m a.s.l. at the top of the central mountain. This shape produces large initial areas with snow accumulation except for cases where the ELA is chosen near the top.

Figure 5 shows the growth curves of the ice sheet for different accumulation rates with an ELA of 600 m a.s.l. The result suggests a nearly linear relation between growth rate and accumulation rate. Only for a large A_c, can the ice sheet reach full size (> 1000 km) in 5000–10 000 a: short enough to grow completely during the period of favourable conditions within Milankovitch cycles. Small values of accumulation, such as occur in parts of Arctic Canada or ice-age Greenland, make growing times of more than 50 000 a necessary, longer than Milankovitch cycles.

Fig.5 . Area of the cross-section of the ice sheet on the idealized topography vs time after initiation. The crosses indicate the time when the ice reaches the coast of the island.

The dependence of growth time on ELA is illustrated in Figure 6, where growth curves for three ELAs are plotted for two A_cs. The initial stage of growth with small accumulation area (high ELA) seems to delay final growth, but final size does not depend critically on ELA. The influence of englacial temperature on growing speed was also tested by increasing the surface ice temperature by about 10 K, however, growing speed was accelerated only by about 10%. This illustrates how strongly accumulation rate controls growing speed, and therefore, the possibility for ice-sheet initiation.

Fig.6 . Area of cross-section of the ice sheet on the idealized topography vs time after initiation, calculated with two accumulation rates for each ELA.

Discussion and Conclusions

The model calculations illustrate the fact that even a small area of less than 50 km diameter with a positive mass balance on a mountain top can be sufficient to initiate the growth of an ice sheet. Such conditions can easily be missed by GCMs trying to find climatic conditions for initiating ice sheets at the proper geographic locations; for example, Reference Verbitskiy and OglesbyVerbitsky and Ogles-by (1992) find Siberia and Tibet more favourable locations for ice-sheet growth than the locations of the ice-age Fennoscandian and Laurentide ice sheets. We suspect that their model produces too much snowfall in Siberia and Tibet (which seems to be a general trend in GCMs); on the other hand, the coarse grid of 200–500 km mesh size and the smoothing of the topography lead to failure in producing positive mass budgets on the mountains in Scandinavia and Canada, which are believed to have functioned as “seeds” for the formation of these former ice sheets.

Additionally, the smoothing of topography in GCMs may lower mountain heights below the critical threshold, and may again miss the initiation conditions for an ice sheet. Positive feedback, such as albedo feedback and climatic feedback (ice-sheet climate), may further enhance ice-sheet growth and stability. Therefore, for GCMs the ice sheet initiation should be considered as a sub-grid phenomenon which needs an appropriate treatment.

For the initiation of a large ice sheet, it is not enough to find areas with occurrence of positive mass balance alone, e.g. Reference Rind., Peteet and KuklaKRind and others (1989) and Reference OglesbyOglesby (1990). Although a glacier will always form from an area of persistent positive mass balance, our steady-state solutions show that a large ice sheet does not always develop: sometimes it ends up as a small ice cap. Using an idealized ice-sheet model, Reference WeertmanWeertman (1976) pointed out the possibility of multiplicity of such steady-state ice-sheet solutions for a range of climatic conditions. This was confirmed for Greenland (Reference Letréguilly, Huybrechts and ReehLetréguilly and others, 1991) and discussed more systematically by Reference HindmarshHindmarsh (1990), but in addition to Weertman’s solutions, multiple stable solutions with finite-size ice masses may occur. In our model runs this was found again, and seemingly, it is the amount of positive mass balance that critically determines the pattern of solution branches: small ice caps or large ice sheets. This result may have important implications for the different growth and evolution of ice sheets in maritime (Reference KerrKerr, 1990; Reference Payne and SugdenPayne and Sugden, 1990) and continental climates.

The amount of snow accumulation is important for yet another reason: ice-sheet growth and decay occur in a changing climate with typical time scales of 10⁴ a. This may be too short to build up a full-size ice sheet if small annual snowfall leads to slow growth of the ice mass. This aspect needs further consideration for ice-sheet initiation and growth in an unsteady climate, such as glacial–interglacial cycles. The transient response of an ice sheet even to very simple periodic climatic forcing displays a wide variety of regular to almost chaotic patterns, strongly depending on the time lag for isostatic adjustment, which will be one aspect of our future work.

Acknowledgements

The work on this topic was initiated by discussions with A. Ohmura. We wish to thank him for continuous encouragement, P. Huybrechts for sharing his experience with ice-sheet modelling, K. Hutter for stimulating discussions, and G. M. Sigut for helpful advice on programming. All model calculations were done on the Cray YMP of ΕΤΗ in Zürich. The first author was supported by the Swiss National Science Foundation (grant no. 20-25271.88).