Introduction

With the publication of the paper by Hays and others (1976), most of the scientists interested in ice-age modelling were convinced that the past changes in the global ice volume were related to variations of the Earth’s orbital parameters. Encouraged by this result, a hierarchy of ice-sheet climate models was developed, intended to explain the link between radiation and ice volume. Only a few ice models were coupled to explicit atmospheric models, from which the snow budget driving the ice sheet could be derived. The most convincing results were obtained by Pollard (1983). He used an ice model following Birchfield and others (1981), and coupled it to an atmospheric model derived from a zonally-averaged and seasonal energy-balance equation (but with land-ocean contrast). Most of the ice-sheet climate models did not include a rigorous calculation of ice temperature, except the model of Morland-Hutter (Morland (1984), and following papers), which considers temperature-dependent planar flow in ice sheets but without application to specific ice-age conditions.

The aim of this work is not to make a further attempt to match data and model results, but rather to inspect the dynamics of the model without tuning the parameters to ice-age conditions. Especially, we are interested in the effect of ice temperature on the evolution of the ice-age ice sheets. Therefore, our ice-sheet model contains a fully coupled flow-temperature calculation and the parameteri-zations used in the atmospheric part of the model are fitted to modern conditions. In the next section we present a brief description of the coupled ice-atmosphere-continent model to simulate the build-up of the Laurentide ice sheet and to study stationary-state dynamics. In the final section we show model results, together with our conclusions.

Model Formulation

The ice-sheet model predicts the ice thickness h along a meridian. The northern boundary of the model ice sheet lies at 70°N (the Arctic coast). The rate of change of ice-sheet thickness h follows from the mass-balance equation:

where b is the annual snow balance and q is the vertically-integrated horizontal ice velocity u (see Fig. 1). Delayed bedrock sinking (or rising) is described in the manner proposed by Weertman (1976), assuming a response time of 10 000 years. For a more detailed review of the basic equations used for ice-sheet modelling see Paterson (1981) and Oerlemans and Van der Veen (1984).

The rate of change of ice temperature T_I, within the ice sheet is controlled by the conservation of energy:

where u, v are the horizontal and vertical ice-velocity components, respectively, k = 36 m a^–1 the thermal diffusivity of ice, and d the production of deformational heat (see below).

Fig. 1. Coordinate system of the coupled model and definition of model variables used in the text. The resolution of the ice-sheet model is 100 km in the horizontal direction and 200 m in the vertical direction.

The surface temperature of the ice sheet is determined by the atmospheric model. At the ice-sheet bottom, a geothermal heat flux of G = 5.02 – 10^–2Wm^–J enters the ice. If the bottom temperature reaches the melting temperature T_M, T = T_M replaces the flux boundary condition. The melting temperature is corrected for pressure: T_M = −a_M h where h is the ice thickness and a_M = 0.87 K km^–1. The values of the physical constants used in the ice model were all taken from Paterson (1981).

The prognostic Equations (1) and (2) determine the time evolution of the ice sheet, provided the ice velocity is known. It can be derived from a set of diagnostic equations, the force balance and the empirical flow law for ice. In the shallow-ice approximation (Reference HutterHutter, 1983), an analytical expression for the horizontal ice-velocity component u can be derived (Reference HerterichHerterich, 1988):

where u_b is the sliding velocity at the ice-sheet bottom h _b(being zero in the present model version) and h_s is the height of the ice surface. In the model, the density of ice (p = 910 kg m^–3) is constant, g = 9.8ms^–2 is the acceleration by Earth’s gravity, and n = 3. The temperature-dependent coefficient A(T′), where T′ is the ice temperature measured above pressure-melting point, is based on measurements compiled by Paterson (1981).

The vertical velocity component v follows from the in-compressibility condition. Finally, the production of defor-mational heat d, defined by

, where and σ _ik are the components of the strain-rate tensor and stress tensor, respectively, can be expressed in terms of the shape of the ice sheet:

with c = 2 – 10³ m² s^–1 JT¹, the heat capacity of ice.

The above equations were formulated on a finite-difference grid. To prevent numerical instabilities, an up-wind scheme was used. The movement of the ice-sheet margin, where the shallow-ice approximation breaks down, was determined by mass conservation (Equation (1)).

The atmospheric model, forcing the model ice sheet, is based on a zonally-averaged energy-balance equation (cf. Reference NorthNorth, 1975), including a seasonal cycle, which is solved analytically for the atmospheric sea-level temperature T_A as a function of latitude. At the ice-sheet surface, the temperature is reduced with respect to sea-level temperature (cf. Reference BowmanBowman, 1982), using a lapse rate of γ = −6deg/km. The mean annual surface temperature T_s of the ice sheet is obtained by integrating over the year.

The model precipitation as a function of latitude φ in mid-latitudes is approximated by a Gaussian distribution:

where P₀ is a proportionality factor, determining the precipitation maximum at latitude φ₀, and c a constant to give the best fit between model and present data. The position φ₀ is related to a critical temperature gradient, derived from the theory of baroclinic instability (Reference HoltonHolton, 1979). In the model, φ₀ is shifted by 10° in latitude during the year, in accordance with observations.

Some additional effects, also influencing precipitation, are contained in the proportionality factor P₀:

By Equation (6), precipitation increases with the surface slope

of the ice sheet thus sirnulatin 8 orographic rain. P₁, is the average precipitation in mid-latitudes and P₂ was chosen such that the model precipitation fits the observations in mountain areas of North America. The form (6) includes the “elevation-desert effect” (Budd and Smith, 1979) with h_s1 , the sea-level height, and a factor τ(T₀), which crudely models the (linear) reduction in precipitation with decreasing global temperature T₀.

The snowfall s in the model depends on precipitation and surface temperature (s = P for T_s < 0, and s = 0 otherwise). Following Pollard (1980), snow melt is parameterized linearly in terms of surface temperature and solar insolation. The annual balance b of snowfall and snow melt follows by integrating over the year.

Results and Conclusions

The aim of the first numerical experiment was to build up a typical ice-age ice sheet within a realistic time. Using the parameterization (6) for precipitation fitted to present-day observations, the Laurentide ice sheet reaches an extent of 3300 km and a maximum height of 3770 m within a build-up time of 40 000 years.

Fig. 2. Ice-sheet shape and isolines of temperature (°C). Dotted areas indicate temperate ice (a) after 35 000 model years integration (with temperature-flow coupling), (b) after 35 000 model years integration (without temperature-flow coupling), (c) after 200 000 model years integration (with temperature-flow coupling).

Figure 2a shows the temperature field within the model ice sheet after 35 000 model years. The temperature increases from the ice surface towards the base. Near the base, however, a temperature inversion is visible. This inversion is due to temperate regions (dotted areas), where the ice is at the pressure melting point everywhere. These temperate regions started to develop after 15 000 model years near the southern edge of the ice sheet. The existence of temperate regions within an ice sheet requires an internal energy source, which is the heat released by deformation. The occurrence of temperate regions in ice-age ice sheets has not yet been reported by other modelers. Since we did not expect to have such temperate regions occurring in our model, we were also not prepared to treat the ice deformation in these regions properly. In the model, the temperature was simply not allowed to rise above the pressure melting point and the deformation was calculated as if it were ice at the pressure melting point, without considering the complexities introduced by melting.

Figure 2b shows another model run, also over 35 000 model years, but without temperature-flow coupling. The flow was calculated assuming a mean ice temperature of −10°C. This build-up experiment produces a higher ice sheet with a shorter southern extent compared to the build-up experiment with temperature-flow coupling. The thickness-to-length ratio is now unrealistically high (increased by a factor of about 1.5). This difference in ice-sheet shape reflects the strong temperature dependence of ice flow. The flow increases by a factor 10 from −10° to 0°C. Since the ice-sheet shape influences its own mass balance, we conclude that ice temperature is an indispensable model variable which has to be included in a realistic ice-sheet model.

An almost-stationary state of the ice sheet is achieved after about 200 000 model years integration (see Fig. 2c). For a stationary state, the amount of accumulated snow over the year has to be equal to the amount of ice melting away. In our model, the ice sheet needs more than 100 000 years build-up time to arrive at a nearly stationary shape. Its length is then approximately 5000 km. The ice sheet needs another 100 000 years to achieve a stationary temperature distribution. We infer from this result that the Laurentide ice sheet was probably never in a stationary state. The temperature distribution of this stationary state differs from the temperature distribution during the initial phase of ice-sheet build-up. In the build-up phase, heat diffusion produces an almost constant vertical temperature gradient. In the stationary state advection of cold ice has flattened the gradient in the upper half of the ice sheet, with a steepening below. Temperate regions are still present, but their vertical extent is now reduced.

In a final experiment, the almost-stationary state of the ice sheet was forced stochastically. Our intention was to study the effect of short time-scale weather fluctuations on the evolution of the ice sheet. Weather fluctuations occur on time scales of hours to days, which is small compared to the response time of ice sheets (on the order of 10 000 years and longer). The weather fluctuations will therefore be treated as a white-noise process, formally introduced in the model by splitting the proportionality factor P₀ in Equation (5) into a mean and a fluctuating part, with a variance equal to the square of the mean. For interpretation of the resulting model time series (and spectrum) we used standard methods of spectral analysis. The resulting time constant is near 18 000 years. The induced stochastic changes of the ice volume had an amplitude of about 10% of the total volume. This leaves 90% of the change between minimum and maximum ice volume to be explained by deterministic processes. In the case of a stronger positive feedback between the ice sheet and the atmosphere than incorporated in our model, this amplitude may be much larger (Reference OerlemansOerlemans, 1979). The stochastic component, however, cannot be too large. Otherwise, it would be hard to explain the observed high coherence between the changes of solar insolation and the ice-volume record (Reference ImbrieImbrie and others, 1984).