Introduction

It is the purpose of this paper to review work on ice-sheet modeling of the Pleistocene ice sheets that has been done since the mid 1970s. The volume of the Pleistocene ice sheets is, of course, a major “output” sought in the astronomical theory of the ice ages. The “input” of the theory is the insolation variations that occur over time because of small changes in the Earth’s orbital parameters.

Perfect Plasticity Ice-Sheet Model

All essential information about the behavior of ice sheets under any time-variable forcing can be obtained from the perfect plasticity ice-sheet model. That is, the dimensions of an ice sheet are found under the assumption that the constitutive equation for plastic flow of ice is that of perfect plasticity. The reason that this assumption works is that flow of ice is governed by the power-1 aw creep equation (Reference GlenGlen 1955)

(1)

where de/dt is the strain-rate, S is the stress, S* is a constant, C is a temperature-dependent constant, and n is the power exponent. If n were set equal to infinity, Equation (1) describes the flow of a perfectly plastic solid. For ice, n is approximately equal to 3. This value is sufficiently large to enable ice to approximate the constitutive flow law of perfect plasticity. For the perfect plasticity law given by Equation (1) no deformation occurs until S = S*; unlimited flow occurs at this value. The stress cannot exceed the value S*. A two dimensional ice sheet of fixed half-width L has a profile given by

(2)

where h is the elevation of the upper surface at x, H is the elevation at the ice-sheet center, and x is measured from the ice-sheet center in the horizontal direction (see Fig.1). The elevation H is given by

(3)

where T is the basal shear stress and has a constant value, ρ is the density, and g the gravity acceleration, and B is a constant equal to 2 (if the ice-sheet base is horizontal) or approximately equal to 4/3 (if isostatic sinking has taken place). (The land surface would be perfectly horizontal If the ice sheet were removed and full rebound took place.)

If the perfectly plastic law is replaced with Glen’s creep law for ice (Equation (1)), it is easily shown that H is proportional to a^z, where a is the average accumulation rate over an ice sheet with no ablation zone and which does not slide over its bed, and z = 1/2(n + 1) ~ 1/8. The ice sheet is “chopped off” at I × I = L. If sliding occurs the 1/8th power is raised to a somewhat higher power, depending upon the sliding law used. In any case the dependence on accumulation rate is a very weak one. If both accumulation and ablation zones exist then the value of H is a weak function of the average ablation rate in the accumulation zone, a, and average ablation rate in the ablation zone, a' (both a and a’ are positive quantities). The perfectly plastic ice sheet can take this weak dependence on a and a' into account by considering the basal shear stress T to be a weak function of these quantities. Account also can be made, through alterations in the value of T, for changes in ice temperature or a change from no sliding to sliding or vice versa. By varying the value of B the amount of isostatic sinking or rebound can also be accounted for with the perfect plasticity model.

Fig. 1

Cross-section of a two-dimensional ice sheet whose northern and southern borders terminate at the water’s edge.

Perfect Plasticity Model Coupled with a Snow Line

Suppose the snow line rises linearly with distance from the northern edge of an ice sheet on a continent in the northern hemisphere. Let the snow-line equation be given by

(4)

where h_s is the elevation of the snow line, s is its slope, and h_So is its value at x = −L.

The firn Tine in the southern half of the two dimensional ice sheet is given by the condition that h(x) of Equation (2) is equal to h_s(x) of Equation (4). Let x = R be the value of x when h − h_s(x). For a steady-state ice sheet aR = a'(L−R) in its southern half. Inserting this relationship into Equation (2) results in

(5)

where h_R is the elevation of the ice surface at the steady-state firn line. Figure 2 shows schematic plots of h_R versus R and h_s versus R. The h_R curve in this figure is shown for the case where a and a' are held constant. Hence the curve is parabolic because More likely, both a and a' vary with R. The actual h_R curve thus is modified. However, this modification can only alter the intersection points but not the qualitative situation of Figure 2. Steady-state solutions exist at the intersections of the two curves. If h_s0 is negative there is one solution and it is a stable one. (By stable we mean that, if it is not initially in a steady state, in time the ice sheet will either grow or shrink towards the steady-state size. If it. is unstable, then, if perturbed slightly from the steady-state size, it will not approach the steady configuration but in time will either shrink until it disappears or grow until it reaches the stable steady-state size.) If h_so is positive there is always one solution in which no ice sheet exists. If an ice sheet does exist, in general there are two solutions. The smaller steady-state ice sheet is one of unstable equilibrium. The larger is a stable steady-state ice sheet.

Fig. 2

Schematic plot of ice-surface elevation at steady-state equilibrium firn line and snow-line elevation versus x = L+R.

A more complete discussion of the material in these last two sections can be found in several of our papers (Reference WeertmanWeertman 1961, Reference Weertman1976, Reference Birchfield, Weertman and LundeBirchfield and others 1981).

Ice Sheet of Fixed Half-Width and Ice Sheet of Variable Half-Width

One fact that should be appreciated with regard to Pleistocene ice sheets is that the southern edges were never limited by encountering the edge of a continent. Because the position of their southern edges was free of any physical constraint, paradoxically these ice sheets were inherently more unstable to changes in their mass budget than are the existing ice sheets in Greenland and Antarctica. The outlet glaciers and ice streams that drain these ice bodies terminate at the waters’s edge. The major effect from large changes in the accumulation on present-day ice sheets is to change the rate at which icebergs are discharged from them. The ice-sheet thickness, as already pointed out, is a very insensitive function of accumulation rate when the horizontal dimensions of an ice sheet are held fixed. The horizontal dimensions are indeed fixed because the ice sheet cannot extend beyond the land mass upon which it is resting. Only if the areas of their ablation zones (virtually non-existent in Antarctica) were to be greatly enlarged would these ice sheets behave somewhat like a Pleistocene ice sheet.

A Pleistocene ice sheet would behave like the Antarctic ice sheet if its continent were rifted apart, east to west, at a high latitude. The major ice loss would be from iceberg discharge at the continental border; ice loss in the ablation zone would be a minor fraction of the total loss. A very large increase of the elevation of the snow Tine would be required to increase the area of the ablation zone to the extent that ice loss in this region becomes important in the mass-balance inventory. This hypothetical stable Pleistocene ice sheet on a rifted continent actually exists; it is the Greenland ice sheet. (Of course, the primary ice-flow directions of this ice sheet are in east-west directions rather than north-south ones.)

In brief, if the horizontal extent of an ice sheet is limited by iceberg calving, then the ice thickness of that ice sheet is also approximately constant when the ice sheet is at its maximum possible width. If the snow line is such that when the ice sheet is at its maximum width the ablation zone becomes of marginal importance, then the ice sheet will be very stable because very large increases in the snow-line elevation are required in order to make ablation processes more important than iceberg production. When the horizontal dimensions of an ice sheet are not restricted, as is the case for the Pleistocene ice sheets, what happens? Figure 2, as already discussed, presents the possibilitfes.

Time Dependence of Build-Up and Destruction of Nonsteady-State Ice Sheets

The change with time of the volume of an ice sheet whose profile is determined with the perfect plasticity constitutive equation can be calculated in a straightforward manner (Reference WeertmanWeertman 1964). The growth time to steady-state equilibrium size from a small nucleus is generally much larger than the destruction time. The qualitative reason for this is that when the ice sheet is small there is only a small area upon which to deposit ice. The rate of change of volume dV/dt is proportional to and approximately equal to āV^2/3(ρg/T)^1/3. Here ā is the average accumulation rate over both the accumulation and ablation areas of the southern half of the ice sheet and is considered to be a positive quantity.

If, on the other hand, ā is negative, the ice sheet shrinks. In this case, particularly if the area of the ablation zone is large, the ice in the ablation zone becomes stagnant because the basal shear stress is smaller than T. Then dV/dt is of the order of āV/H. The destruction time can be very much faster than the build-up time.

We suggest that under orbital forcing ice sheets rarely reach their steady-state condition. That is, during the growth stage they never become as large as they could if the growth rate were very large and the orbital forcing frequencies were very small. During the waning stage the ice sheet at some instant may have a half-width corresponding to a steady-state half-width, but the profile of the ice sheet will not necessarily be that of steady-state ice sheet. During shrinkage, part of the ice sheet will not be flowing; for the steady-state ice sheet, however, the basal stress is large enough that ice flows through every cross-section.

Astronomical Orbital Perturbation Forcing of Pleistocene Ice Sheets

If the snow line, the ablation rate, or the accumulation rate varies over time periods comparable to or larger than the characteristic response times of ice sheets, the ice-sheet dimensions too will vary with time. A relatively small change in the accumulation rate or elevation of the snow line or in the ablation rate can lead to large changes in the steady-state dimensions of an ice sheet. Whatever ultimately drives snow line elevation changes etc., drives the ice sheet. If the time periods of the mechanism that ultimately drives the Pleistocene ice sheets are not very much larger than the characteristic time periods of the ice sheet itself, the ice sheet never will attain a steady-state condition. It will either be growing or shrinking. Its only stationary period occurs when it switches between these two states. The Pleistocene ice sheets appear never to have attained a true steady-state condition.

The simplest assumption that can be made in a model calculation of a Pleistocene ice sheet driven by insolation variations produced by the Milankovitch (astronomical) orbital perturbations, is that the snow line goes up and down by amounts that are proportional to the variation. This simple-minded picture was used by us (Reference WeertmanWeertman 1976, Reference BirchfieldBirchfield 1977, Reference Birchfield and WeertmanBirchfield and Weertman 1978) in an attempt to see if the insolation swings, which actually are not all that small (about 10% of the average value during the high radiation half of the year at 50°N), are potent enough to create and destroy the huge Pleistocene ice sheets. The answer from this calculation is that they are. The snow line in this calculation was made to go up and down by an amount proportional to the insolation change. In other words, the elevation term h_so of Equation (4) varies by an amount proportional to these insolation variations. The proportionality constant was set by noting that the insolation and the snow-line elevation both increase towards the equator. Thus a change in snow-line elevation can be related to a change in insolation. It was assumed that this spatial change gives a reasonable order of magnitude estimate of the temporal changes.

Figure 3 is one plot of Pleistocene ice-sheet half-width as a function of time before present, obtained in this simple manner. Also shown in this plot are two projections of future Pleistocene ice-sheet half-widths made for two very slightly different accumulation rates. In one case no ice sheets form for over the next 50 ka. In the other case very large ice sheets form in the immediate future. This result is an indication that it may not be reasonable to expect any Pleistocene ice-sheet model calculation, no matter how sophisticated, to match perfectly some geologic record, such as the deep-sea core isotope record. At some critical time bifurcation may be a possibility when noise can push the system into either a glacial or an interglacial path that runs for a significantly long period of time.

Fig. 3

Half-width L of a Pleistocene ice sheet as a function of time and its prediction into the future for two slightly different accumulation rates (Reference WeertmanWeertman 1976).

The 100 ka Signal

It is well known that if the response of ice sheets to the orbital forcing were linear, the time variation of ice-sheet volume should show virtually no 100 ka signal. But a strong 100 ka signal exists in the deep sea core record of ice volume. The obvious explanation of how this signal arises is that the response of the ice sheets is not linear. Reference BirchfieldBirchfield (1977) pointed out that half-wave rectification of two forcing frequencies close to each other produces an output signal containing a strong signal at the difference in the two frequencies. This, in effect, is what an ice sheet can do.

Consider Figure 4 which shows schematically a modulated 20 ka signal with a modulation at 100 ka. Suppose the mean signal position, as indicated by the horizontal dashed line in Figure 4(a), that separates growing from shrinking of an existing ice sheet is not the mean signal position of zero. The ice sheet might grow and then shrink as indicated in Figure 4(b) with rapid destruction, rather than a relatively minor retreat, of the ice sheet occurring only where the modulated amplitude becomes unusually large in the negative sense. Clearly this “half-wave” rectification effect of Birchfield can produce an ice-volume history with a strong signal at the modulation frequency. To us, this still appears as a likely mechanism of producing power at 100 ka. In the early Pleistocene the conditions for the rectification may not have been optimal. (Reference Watts and HayderWatts and Hayder (1983) have discussed a possible switch at 700 ka BP which may be the equiva-lent of going from a non-rectification to one of rectification.) Natural internal mechanisms of this frequency are hard to come by. It is true (Reference PeltierPeltier 1982) that decay modes within the Earth’s mantle that arise at layers where phase transformations have produced a density contrast, can have long-time constants of this magnitude. As shown in the appendix, the trouble with such modes is that the total rebound or isostatic sinking they produce is only a very minor fraction of the total (of order of 5% or less, or 50 m or less). It is difficult to see how such a small amplitude effect could ever be significant in the waxing and waning of ice sheets.

Fig. 4

(a) A 20 ka signal with a modulated amplitude of 100 ka period used as forcing for an ice sheet, (b) A schematic history of ice-sheet volume change if ice-sheet growth is favored above dashed curve and destruction below it in (a) and a rapid (crash) destruction occurs under unusually large amplitudes.

Adding Complexities to the Ice-Sheet Dating

The description that has been summarized up to this point, we believe, presents the essential physical features which are involved in the waxing and waning of ice sheets under a time-dependent forcing and which have been considered in numerous papers by ourselves and others (see, for example, Reference Budd and SmithBudd and Smith 1981, Reference Bindschadler and GoreBindschadler and Gore 1982, Reference Birchfield and WeertmanBirchfield and Weertman’ (1982), Reference OerlemansOerlemans 1982[a], Reference Oerlemans[b]). It does not consider the forcing itself (other than for the most simple case of forcing in which the snow-line elevation moves up and down by amounts proportional to the insolation changes at a fixed latitude). Much more sophisticated forcing can be placed into the analysis. The degree of sophistication is in a continuum ranging from the simple-minded forcing already described to that that cannot be handled in a reasonable time by an existing computer.

For the ice sheet, the essential response needed from the atmosphere is the mean annual accumulation rate. (The ablation rate too is needed.) This accumulation rate is one part of the entire hydrological cycle, and hence quite complex. Several models of the ice sheets have incorporated models of the atmosphere and ocean to some degree of complexity, for the purpose of estimating the accumulation on the ice sheet. For example Reference PollardPollard and others (1980) and Pollard (Reference Pollard1978, Reference Pollard1980, Reference Pollard1982) include a one-parameter energy-balance model (EBM) which is used to calculate the position of the snow line as a function of time; Reference Källën, Crafoord and GhilKällén and others (1979) and Reference Ghil and Le TreutGhil (1981[a]) have made a similar calculation. Reference Birchfield and WeertmanBirchfield and others (1982) briefly introduce a more complex EBM which, although it does not include a closed hydrological cycle, does explicitly calculate the mean annual accumulation rate at the ice-sheet surface. This is a result of fairly detailed calculation of energy and mass fluxes at’the surface of the ice sheet.

It is apparent that, for the accurate prediction of the net growth or decay of the ice sheet, there must be a careful treatment of the hydrological cycle. There exists a fundamental problem, in that, although the mass and energy flows involved in the growth or decay of the ice sheet are large, the time scale is so very large that the fluxes are extremely small when averaged globally and hence may be very difficult to predict accurately. The existence of this inherent limitation on modeling has been proposed and explored by Reference SaltzmanSaltzman (1983) and Reference Saltzman and SuteraSaltzman and Sutera (in press). The possibility of internal oscillations within the forcing system itself (Reference Ghil and Le TreutGhil 1981[b], Reference Ghil and Le TreutGhil and Le Treut 1981, Reference Le Treu and GhilLe Treut and Ghil 1983) also may bring in additional complications.