Thermomechanical modelling of Northern Hemisphere ice sheets with a two-level mass-balance parameterization

Abstract

A three-dimensional time-dependent thermomechanical ice-sheet model was used together with a two-level (snow-accumulation/runoff) mass-balance model to investigate the Quaternary ice sheets of the Northern Hemisphere. The model freely generates the ice-sheet geometry in response to specified changes in surface temperature and mass balance, and includes bedrock adjustment, basal sliding and a full temperature calculation within the ice. The mass-balance parameterization makes a distinction between snowfall and melting. Yearly snowfall rates depend on the present precipitation distribution, and are varied proportionally to changes in surface temperature and the moisture content of the air. The ablation model is based on the positive-degree-day method, and distinguishes between ice and snow melting. This paper discusses steady-slate characteristics, conditions for growth and retreat, and response time-scales of ice sheets as a function of a prescribed lowering of summer temperature. Most notably, the modelled extents of the Eurasian ice sheet for a summer temperature lowering of 6–7 K and of the Laurentide ice sheet for a cooling of 9–10 K are in reasonable agreement with most reconstructions based on geological evidence, except for the presence of a large ice sheet stretching from Alaska across the Bering Strait to most of eastern Siberia. In addition, wet basal conditions turned out to be always confined to the margin, whereas central areas in these reconstructions remained always cold-based. This is of relevance for processes involving reduced basal traction.

1. Introduction

During the Quaternary, the waxing and waning of ice sheets on the continents of the Northern Hemisphere constitute a major driving force for global change and are the main source of global sea-level fluctuations. It is generally accepted that the maximum sea-level depression during the last glacial cycle was around 125 m, largely originating from ice-sheet build-up in North America and Eurasia. This corresponds to an additional global ice volume of around 50 × 10¹⁵ m³ (Fairbanks, 1989; Tushingham and Peltier, 1991). Whereas the maximum ice-sheet extent along the southern margins of the Laurentide and Fennoscandian ice sheets seems m be reasonably well documented (e.g. Demon and Hughes, 1981), the distribution of ice sheets along the Russian side of the Arctic Ocean is still largely unresolved. Here, geological reconstructions range from limited ice expansion over Scandinavia and Svalbard in the west (Boulton, 1979) to an unbroken Eurasian ice sheet extending all the way from Britain ю Scandinavia across the Barents and Kara Seas to most of eastern Siberia (Grosswald, 1980). The thickness of these ice sheets, cm the other hand, cannot be directly inferred from the geological record, but has to be estimated from models.

So far, most ice-sheet modelling studies, three-dimensional and two-dimensional (horizontal plane) have considered only a part of the Northern Hemisphere ice-sheet system (e.g. Budd and Smith. 1987; Letréguilly and Ritz, 1993) or have concentrated on the ice sheets as a lower-boundary condition for climate models (Deblonde and Peltier, 1991; Marsiat, 1992; Verbitsky and Oglesby, 1992). None of these models included thermodynamics. Nevertheless, ice temperature exerts a strong control on the ice viscosity and subsequently on the height/width ratio. Moreover, temperate ice and melting at the base are necessary preconditions for basal sediment deformation and enhanced basal sliding, which are often mentioned as mechanisms that potentially can destabilize large continental ice sheets by making them much flatter towards the end of a glacial cycle (Boulton and Jones, 1979; Oerlemans, 1982).

Here, we present results obtained with a three-dimensional thermomechanical ice-sheet model which covers all of the Northern Hemisphere where widespread continental glaciation is known to have taken place. This model is time-dependent, includes geodynamics and freely generates the ice-sheet geometry in response to changes in environmental conditions. The model has been extensively tested and validated on the present-day ice sheets of Greenland and Antarctica (Huybrechts, 1990; Huybrechts and others, 1991). The methods adopted to estimate the mass-balance components are similar to those that proved their applicability on these ice sheets. Though the prescription of climatic boundary conditions is likely to be more complex at mid-latitudes than in the polar regions, this is seen as part of an inverse modelling exercise, in which an attempt is made to investigate the climatic conditions that produced these ice sheets in the first place. This paper aims to document some basic characteristics of the Northern Hemisphere ice sheets (ice volume, ice-sheet geometry, response time-scales and basal temperature conditions) as a function of an imposed temperature lowering. Of special interest are the various thresholds required to grow and decay the individual ice sheets that composed the Quaternary glacial system. Time-dependent simulations involving glacial cycles will, however, be presented elsewhere (paper in preparation by P. Huybrechts and S. T’siobbel).

2. Model Set-Up

The ice-sheet model closely resembles the version described by Huybrechts (1993). It basically solves conservation equations for mass and heat and has a full coupling between the flow and temperature fields. Only grounded ice flow is taken into account, with a Glen-type flow law and a temperature-dependent rate factor given by an Arrhenius relationship. Longitudinal deviatoric stresses are disregarded and basal sliding (of Weertman type) is restricted to areas at the pressure-melting point. In the present formulation, the possibility of ice shelves and grounding-line dynamics was excluded and the maximum possible extent of grounded ice was constrained by specifying a coastline following the −500 m isobath. This approximately coincides with the margin of the continental platform in the Arctic basin, over which ice sheets can freely expand provided the surface mass balance is sufficiently positive. Beyond this pre-set boundary, all ice is lost as calf ice.

The thermodynamic calculations employ standard procedures. Heat transfer is considered to result from vertical diffusion, three-dimensional advection and deformational heating at layer interfaces. In the time-dependent mode, heat can also flow into an underlying rock slab of 2 km thickness. The lower-boundary condition on temperature corresponds to a geothermal heat flux of 42 m W m⁻². Bedrock adjustments are modelled as a damped return with a steady-state deflection given by local isostatic equilibrium. The e-folding time scale is set at 3000 years, with a mantle density which is 3.6 times that of ice. All calculations take place on the same grid. It comprises 193 × 193 grid cells that are laid out over a polar stereographic map projection with standard parallel at 60°N. With this standard parallel, distortions in distance range from a factor of 0.92 at 45°N to a factor of 1.07 at the pole, the effect of which has been taken into account to calculate the ice fluxes and the ice volume. The horizontal mesh size is 50 km, and there are 11 layers in the vertical, concentrated towards the base. This resolution is the minimum required to incorporate realistic geography and to capture the essential characteristics of the vertical velocity and temperature profiles.

The mass-balance model distinguishes between snow accumulation and ablation and is driven by a prescribed “background” temperature change that is uniform in space but 50% higher in January than in July. The latter feature is a characteristic common to most general-circulation model (CCM) simulations (e.g. Manabe and Broccoli, 1985) and is also supported by palaeoclimatic reconstructions based on geological and biological evidence (Frenzel and others, 1992). This feature is usually attributed to the albedo-temperature effect and a stable vertical stratification in the lower stratosphere, which prevents large-scale mixing of the air, especially during winter. The basic inputs are data sets for the mean monthly sea-level temperature and the mean monthly precipitation rate over the entire grid. These reference climatologies were taken from precipitation maps (Jaeger. 1976) and from temperature data used to initialize the Hamburg CCM model (ECHAM-1 in the T21 mode) at the 1000 hPa level. Surface temperatures were calculated using a lapse rate of 8 K km⁻¹, a value somewhat higher than the standard lapse rate for the free atmosphere. This is in line with measurements on the ice sheets of Antarctica and Greenland, and can he explained by radiative cooling at the surface, particularly d tiring the winter season. The corresponding surface elevation was taken as the maximum of cither the mean elevation of a grid cell during the calculations or the highest peak within this grid cell for the present topography. This is because the glaciation threshold is to a large extent determined by the most elevated terrain, a piece of information that is otherwise lost when values are averaged on the 50km grid.

The model considers only changes of precipitation intensity but keeps the precipitation pattern unchanged. The total precipitation is set proportional to the moisture-holding capacity of an air column, and has a sensitivity of 4% K⁻¹. This seems to be a reasonable assumption at polar latitudes, supported by analyses of Antarctic and Greenland ice cores (Yiou and others, 1985; Clausen and others, 1988), but its validity during maximum glaciation at mid-latitudes can be questioned. Here, precipitation is likely to be increasingly a function of such factors as atmospheric circulation, continentality and orography. A more sophisticated approach, however, is deemed far beyond the purpose of the present study. The fraction of precipitation falling as snow is determined by the surface temperature, and varies from one to zero for mean monthly temperatures in the range of 263–280 K. Such a relation is suggested from meteorological observations in Canada and Russia. Ice and snow ablation, on the other hand, are dealt with locally. Following a procedure previously adopted for the Greenland ice sheet (Reeh, 1989; Huybrechts and others, 1991), the melting rate is set proportional to the yearly sum of positive degree days (PDD). Degree-day factors were 3 mm a⁻¹ /PDD for snow, and 8 mm a⁻¹/PDD for ice. The mass-balance parameterization furthermore incorporates the process of super-imposed ice formation, but neglects a possible contribution from rainfall, which is assumed to run off entirely. In the model, summer temperature was taken as the relevant forcing parameter. This means that the “effective” ablation temperature roughly equals the applied temperature drop, but that the effective (i.e. mean annual) temperature for the calculation of the precipitation intensity is 25% higher. For the latter calculations, an additional temperature change resulting from elevation changes was not taken into account.

3. Results

3.1. Ice-sheet geometries

First of all, the model was run under present climatic conditions with a zero temperature perturbation starting from no ice cover at all. The result, after a total integration lime of 10⁵ years, is shown as the upper left panel in Figure hi. As demonstrated in this figure, the simulation is quite realistic. The Greenland ice sheet is well reproduced, and the smaller ice masses turn up in almost the same places as they are observed today. This is the case for the Arctic islands (northern Novaya Zemlya, Svalbard, Franz Josef Land and Severnaya Zemlya); but also the location of ice in Scandinavia (Jøstedalsbreen and Svartisen), the Alps, southeastern Iceland (Vatnajökull), Alaska and Kamchatka is in good agreement with the present situation. Apparently, the mass-balance model is able to identify the correct areas with a positive mass balance. The only real deviation from the present-day situation is found over the Canadian Arctic. Here, the model produces a unified ice cap over Ellesmere Island, Axel Heiberg Island and Devon Island of over 2000 m altitude, probably because of an inadequate representation of the coastline and smaller-scale fjords, allowing small glaciers and ice caps m merge into a larger ice cap.

The present state was then taken as an initial condition for the climate sensitivity experiments, and model calculations ran until a full steady stale was reached. Resulting geometries for summer temperature lowerings in the range of 2 to 10 K are shown in Figure 1a. It appears that major ice sheets spread from glaciation centres that are all present in proto-form today. In North America, two separate ice sheets build up over the coastal Rockies and the Canadian Arctic Islands. A cooling of more than 5 K is required for the latter ice sheet to cross Hudson Bay and cover the Ungava Peninsula and Labrador plateaus, which themselves do not act as independent inception areas. During this evolution, the summit of the Laurentide ice sheet migrates west towards the Atlantic and reaches a maximum elevation of around 4000 m. A very substantial summer cooling of 10 K (in the model this corresponds to −12.5 K in the annual mean and −15 K in winter) is required for the Laurentide ice sheet to merge with the Cordilleran ice sheet, but the respective ice domes are well conserved. The model has difficulty growing the Laurentide ice sheet far enough to the south in the Midwestern U.S.A. This may well indicate that the model underestimates precipitation rates there and that a reorganization of the atmospheric circulation produced much welter conditions during maximum glaciation in continental North America.

Fig. 1. a. Modelled steady-state ice-sheet geometries for summer temperature perturbations below present levels as indicated. Contours are for surface elevation and are spaced 500 m apart. These experiments started with a simulation of the present ice-sheet distribution (shown in the upper left panel) as initial configuration. b. Basal temperature conditions corresponding to the simulations shown in Figure 1. Dark shading indicates where the base is at the pressure-melting point, and basal sliding ran occur: intermediate shading is for frozen basal conditions; light shaded areas are ice-free.

In Eurasia, the model produces the Barents and Kara ice sheets for a moderate cooling of about 4 K. at the same time as the Fennoscandian ice sheet is still largely confined to the Norwegian mountains. This supports the view that the Siberian ice sheets originated in the Arctic basin and spread “from the sea” into the adjacent continent (Grosswald, 1980), though such an inference should be reserved because the model does not explicitly deal with marine ice-sheet dynamics. For a summer temperature lowering of 6–7 K, the simulated Eurasian ice sheet is in broad agreement with most geological reconstructions for the Last Glacial Maximum, though there are some regional deviations. Maximum elevations are about 3500 m here.

Probably the largest deviation from what is known of the Quaternary glacial cycles is the continuous ice sheet covering Alaska, the Bering Strait and northeastern Siberia, which ultimately (for a cooling below −8 K) merges with the Eurasian ice sheet and an ice cap north of the Sea of Okhotsk to form a huge circumarctic ice mass. There is no evidence to support glaciation in northern Alaska and the Bering Strait, though several reconstructions imply at least some ice in coastal eastern Siberia and the higher plateaus further south (Denton and Hughes, 1981; Tushingham and Peltier, 1991). A plausible explanation for the discrepancy would be changes in precipitation patterns, not accounted for in our mass-balance parameterization. For instance, the development of an ice sheet over southern Alaska and permanent sea ice in the Arctic Ocean may well have deprived the northern coastal areas of sufficient moisture.

3.2. Basal temperature conditions

Figure 1b shows the steady-stale thermal conditions at the base associated with these climate sensitivity experiments. The distribution of areas at pressure-melting point is important because the production of basal water is a necessary prerequisite for basal sliding and the possibility of enhanced ice flow due to the deformation of underlying sediments. Interestingly, wet-based areas in these reconstructions are confined to the outer regions, whereas central areas remain always cold-based. Despite lower accumulation rates in colder climates, implying less advection of cold ice towards the base, and central ice thicknesses above 3–4 km, contributing to a basal warming because of the “insulating effect”, the simulated Northern Hemisphere ice sheets fail to become temperate over the large areas needed for the possibility of some kind of sliding instability. Typically, mean annual surface temperatures over the higher parts of the ice sheets simulated for a cooling of 10 K are of the order of 220–230 K. with corresponding basal temperatures down to 15–20 K below the pressure-melting point. Consequently, these simulations do not support ice-dynamic theories that invoke thermal switches at the base.

3.3. Ice volume

The corresponding ice volumes are summarized in the solution diagram shown in Figure 2. Two sets of model runs are displayed, one in which the present ice-sheet distribution was used as an initial condition (“interglacial conditions”) and one in which the ice-sheet configuration for −10 K served to start the calculations up (“glacial conditions”). There are two solutions, depending on whether the ice sheets were already present or not, a behaviour well described for ice sheets where the mass-balance held is a function of surface elevation (Oerlemans and Van der Veen, 1984). Typically, the “width”‘ of the hysteresis is found to be about 2 K. An additional ice volume of about 4.5 × 10¹⁶ m³, which is a reasonable estimate for the Last Glacial Maximum contribution of the Northern Hemisphere ice sheets, is seen to be produced for a summer temperature lowering of 6–7 K. This value refers to a steady state, however. As demonstrated in Figure 3, there are long time-scales required to build ice sheets up. Typical e-folding time-scales appear to be of the order of 40 000–60 000 years, making it quite unlikely that the Xorilu-rn Hemisphere ice sheets during the LGM were in full equilibrium with the climate at that time. Consequently, the inferred cooling should be considered as a lower bound. Melting of ice sheets, on the oilier hand, is a much faster process because of the higher magnitude of ablation rates. According to the curves displayed in Figure 4, it may take as little as 5000–6000 years to entirely remove the Laurentide and Eurasian ice sheets from their continental base when the climate switches from full glacial to full interglacial conditions. In that case, only a small ice cap remains centred over northern Novaya Zemlya and the adjoining Kara Sea that would require an additional mechanism of perhaps marine origin for its decay.

Fig. 2. Solution diagram giving steady-state ice volume as a Junction of the temperature perturbation. The full squares are for model rum in which the present “interglacial” state was taken as an initial condition; the open squares started from the “glacial” ice sheet simulated far a a temperature change of −10 K. As a first estimate, 10 million km³ (i.e. 10¹⁶ m³) of ice corresponds to a worldwide sea-level lowering if about 25 m.

Fig. 3. Time-scales for ice-sheet build-up after the model was submitted to a sadden climatic cooling of the magnitude shown at the right. The present interglacial state served as an initial condition.

Fig. 4. Evolution curves for Northern Hemisphere ice volume fur various temperature rises above −10 K. These calculations started with the glacial ice-sheet configuration as an initial condition.

Another kind of (tentative) inference can be made when relating the ice-volume curves of Figure 2 to the reconstructions shown in Figure 1b. For instance, the ice-sheet reconstruction for a temperature lowering of 6 K contains roughly the right amount of additional LGM ice, but shows a Laurentide ice sheet that spreads far less to the southwest than indicated by the geological record. As the extent of the Eurasian ice sheet corresponds better with most reconstructions, this might indicate that ice sheers during the Last Glacial Maximum were significantly thinner than the model simulations suggest.

4. Conclusion

This paper discussed Northern Hemisphere ice-sheet simulations obtained with a three-dimensional thermomechanical ice-sheet model coupled with a parameterized mass-balance model. The most important findings that came out of these experiments can be summarized as follows:

• (i) The degree-day method estimating runoff based on measurements in Green hind seems to have a continent-wide validity, at least for present-day climatic conditions. This follows from a basically correct identification of ablation areas for the present state.

• (ii) The Eurasian ice sheet reaches approximately its full extent for a summer temperature lowering of 6–7 K.

• (iii) A North American ice sheet in rough agreement with LGM reconstructions is only produced for coolings of below 9–10 K.

• (iv) The Barents–Kara ice-sheet complex originated on islands in the Arctic Ocean and spread from there to the adjacent continent.

• (v) Basal temperatures in the central parts of the simulated ice sheets remain always well below the melting point, and wet basal conditions are confined to the margins.

• (vi) Several ten thousands of years are required to grow the large Northern Hemisphere ice sheets to their full size.

• (vii) Decay of full-grown ice sheets by surface melting can be accomplished in about 5000 years under an interglacial climate.

Nevertheless, this study has also highlighted some of the shortcomings of the present model. The most important ones are probably related to the influence of the ice sheets on the atmospheric circulation, in particular with respect to changes in precipitation patterns, and to the wav marine ice-sheet dynamics are treated in the model. A further logical step is then also to couple the ice-sheet model with an ice-shelf model and to incorporate an atmospheric-flow model that deals with changes in the wind field as the ice sheets grow into the troposphere. These developments are presently under way.