Ice-age simulations with a calving ice-sheet model

QUATERNARY

RESEARCH

Ice-Age

20, 30-48 (1983)

Simulations

with a Calving

Ice-Sheet

Model

DAVIDPOLLARD Climatic Research Institute. Oregon Stare University, Corvallis, Oregon 97331 Received September 9. 1982 Variations of ice-sheet volume during the Quaternary ice ages are simulated using a simple icesheet model for the Northern Hemisphere. The basic model predicts ice thickness and bedrock deformation in a north-south cross section, with a prescribed snow-budget distribution shifted uniformly in space to represent the orbital perturbations. An ice calving parameterization crudely representing proglacial lakes or marine incursions can attack the ice whenever the tip drops below sea level. The model produces a large - lOO,OOO-yrresponse in fair agreement (correlation coefficient up to 0.8) with the 8180 deep-sea core records. To increase confidence in the results, several of the more uncertain model components are extended or replaced, using an alternative treatment of bedrock deformation, a more realistic ice-shelf model of ice calving, and a generalized parameterization for such features as the North Atlantic deglacial meltwater layer. Much the same iceage simulations and agreement with the 6r*O records, as with the original model, are still obtained. The model is run with different types of forcing to identify which aspect of the orbital forcing controls the phase of the lOO,OOO-yr cycles. First, the model is shown to give a -1OO,OOO-yr response to nearly any kind of higher-frequency forcing. Although over the last 2-million yrs the model phase is mainly controlled by the precessional modulation due to eccentricity, over just the last 500,000 yr the observed phase can also be simulated with eccentricity held constant. A detinite conclusion on the phase control of the real lOO,OOO-yrcycles is prevented by uncertainty in the deep-sea core time scales before -600,000 yr B.P. The model is adapted to represent West Antarctica, and yields unforced internal oscillations with periods of about 50,000 yr.

INTRODUCTION

roughly 20,000 and 40,000 yr, and these have been well correlated with insolation variations due to the orbital perturbations of about the same periods, in accordance with the “Milankovitch” theory (e.g., Hays et al., 1976). Several simple climate and/or ice-sheet models have been able to simulate the higher-frequency oscillations as a direct response to the orbital forcing (Weertman, 1976; Schneider and Thompson, 1979; Suarez and Held, 1979; Pollard et al., 1980; Budd and Smith, 1981). Almost all of the spectral power of the orbital forcing occurs at periods around 22,000 and 41,000 yr due to variations in the precession and obliquity, respectively, and these models for the most part responded linearly with amplitudes similar to or greater than the observed minor oscillations. However, little or no sign of the dominant lOO,OOO-yr glacial cycles was produced. Recently Oerlemans (1980) and Birch-

One of the main tasks in ice-age modeling is the simulation of the fluctuations of ice sheets in the Northern Hemisphere over the last million years. Oxygen-isotope measurements in deep-sea cores are generally believed to be good proxy indicators of icesheet volume in the Northern Hemisphere, although they may also be influenced by other variables such as temperature and Antarctic ice-sheet size (Dansgaard and Tauber, 1969; Chappell, 1974; Budd, 1981). The dominant features in these records are “glacial” cycles with an average period of about 100,000 yr and with full glacial to interglacial amplitudes. Many cycles are terminated by a relatively rapid retreat from maximum ice-sheet size to conditions like those at present, as occurred, for instance, between about 18,000 and 6000 yr B.P. The records also show minor oscillations of smaller magnitude with spectral peaks at 30 0033-5894/83 $3.00 Copyright All rights

@ 1983 by the University of Washington. of reproduction in any form reserved.

ICE-AGE

31

SIMULATIONS

field er al. (1981), using simple ice-sheet models with orbital forcing, obtained results that agree with some aspects of the observed glacial cycles. This was achieved by including a realistic time lag in the bedrock deformation under the ice load, which caused some of the ice-sheet retreats to be amplified. In the model of Birchfield et al. (1981) these retreats coincided well with the observed deglaciations. Pollard (1982a) extended this mode1 by adding a crude parameterization of calving by proglacial lakes and/or marine incursions. Calving then occurred only during the rapid retreats at roughly lOO,OOO-yr intervals, further accelerating the retreats to achieve more complete deglaciations; this produced further improvements in the overall agreement with the deep-sea core records. Many of the same results were also obtained when the calving ice-sheet model was coupled to a global seasonal climate model (Pollard, 1982b). In the present paper the ice-sheet model with calving is extended and the results are analyzed further. First the model formulation is outlined briefly and the basic model results are compared with some SlsO records. Several alternate or improved treatments of the most uncertain mode1 components are described in an effort to increase confidence in the results of the model. The relationships between the model’s lOO,OOO-yr response, the orbital forcing, and the variations of eccentricity are then examined more closely. Several investigators have noted the agreement between the phases of ecc:entricity and the observed glacial cycles over the last -600,000 yr, although the magnitudes of the individual cycles are uncorrelated (e.g., Imbrie and Imbrie, 1980; Johnson, 1982). Variations in eccentricity serve mainly to modulate the -22,000-yr precessional forcing, and cause only a very small amount of lOO,OOO-yr power in the insolation variations. Consequently, sorne highly selective response would prestrmably be required to maintain a phase-lock with eccentricity. These issues

are investigated :tere by running the model with a variety 01’ imposed forcing functions to identify which components of the orbital forcing are controlling the phase of the model’s lOO,OOO-yr response. Other mechanisms that might be important in the lOO,OOO-yr glacial cycles include long-term free internal oscillations (Sergin, 1979; Ghil and Le Treut, 1981) and “tunneling” due to stochastic forcing (Benzi et al., 1982; Nicolis, 1982). Both of these are absent in the present Northern Hemisphere model. In the Appendix the model is tentatively applied to West Antarctica, where the model geometry allows long-term internal oscillations with no external forcing, along the lines suggested qualitatively by Thomas and Bentley (1978) and Thomas (1979). MODEL

FORMULATION BASIC RESULTS

AND

Only a terse formulation is given here, since the model is described more fully in Pollard (1982a). Much of the model is similar to those of Oerlemans (1980) and Birchfield et al. (1981). The overall structure of the model is shown in Figure 1. Ice thickness, h, and the elevation of the bedrock surface, h’, are predicted as functions of latitude, representing a north-south profile along a typical flow line through the Laurentide or Scandinavian ice sheets,

!$ = A$[he~a(h; h’)laa(h; h’)] + G(h + h’, x, orbit), and ah’ -=

at

v5

(1)

[h’ - hb(x) + rh].

Here t is time and x is distance to the south. Equation (1) uses a vertically integrated approximate ice-flow law with east-west flow neglected (e.g., Birchfield et al., 1981). The ice-flow coefficient A = 5.77 x 10d4 m-3 yr-I, (Y = 5, and B = 2. G is the net annual mass balance on the ice surface, representing the current distribution of snowfall

32

DAVID

POLLAKD

Latitude F’N) FIG. 1. Schematic structure of the model. Here 11 is the ice-sheet thickness. 11 ’ is the bedrock elevation above present mean sea level (full line), 116 is the equilibrium crustal topography (dotted line), and G = 0 is the equilibrium line of the prescribed snow budget (short-dashed line). A northern body of water is shaded: during deglaciations a similar water body can also form and cause calving at the southern edge of the ice sheet when the tip drops below sea level.

and ice melt; its dependence on the surface elevation h + h' follows Oerlemans (1980),

spheric channel (Walcott, 1973; Ghil and LeTreut, 1981), with the elastic response of the lithosphere neglected. However, permanent crustal topography is prescribed in - E) - b(h + h' - ET the hb (x) term. This approximately represents a North American profile running ifh + h' - ES 1500m south-southwest from Baffin Bay, pieceifh+h-E>1500m I myrm” wise-linearly joining the following (latitude (3) (ON), height (m)) points: (74,-5001, (70, 850). (66, 200), (40,400). (30,400). In Eq. (2)the where a = 0.81 x 10-3yrm’, b = 0.30 x asthenospheric flow coefficient v = 100 km? 1O-6 m-j yrl; all elevations refer to the yr-I, and Y = 0.3 is the ratio of ice density present mean sea level. E is a function of p1 to rock density pA. The most appropriate latitude denoting the equilibrium-line alti- model for the mantle response depends on tude at which the mass balance is zero (in- the vertical profiles of viscosity and density dependent of the ice-sheet profile), and is in the upper -lo3 km. which are not known taken to have a constant slope in latitude. with certainty (Peltier, 1981); this point is It is shifted uniformly in response to the discussed further below. orbital perturbations according to Ice calving is added to the basic model as described in Pollard (1982a). This crudely E = E,, + s(x - x,,) + kAQ. (4) represents accelerated flow and wastage at where AQ is the difference in the summer the southern tip of the ice sheet due to prohalf-year insolation from that of the present, glacial lakes and/or marine incursions. at the arbitrarily selected latitude of 55”N, These processes have been emphasized by calculated from the current orbital paramAndrews (1973) and Denton and Hughes eters (Berger, 1978). Unless otherwise (1981) as possibly being important during stated, k = 35.1 m/(W mm?), E,, = 550 m, the last deglaciation. Here a water body is s = IO-j, and x0 corresponds to 70”N. assumed to form outside the ice sheet Equation (2) follows from a linearized whenever the bedrock falls below sea level. flow model in a relatively thin astheno- The ice budget G in (3) is set equal to a

ICE-AGE

large negative value at ice-sheet points between the tip and the first interior ice-sheet point that cannot be floated by the hydrostatic head of the water: G(xi+l) = -20 m yr-1 if PI&~) < p&S - h’(xi)) and h’(xi+J < s

shelves in West Antarctica (e.g., Thomas and Bentley, 1978). In this one-dimensional model with no east-west resolution, the calving parameterization cannot distinguish between inland proglacial lakes above sea level and marine incursions at sea level. By choosing S = 0, Eq. (5) is envisaged as mainly representing marine incursions. However, both types of water bodies certainly existed during the last Laurentide deglaciation, and a future extension of the model to two dimensions, possibly including explicit lake water budgets, would be valuable. Equations (1) and (2) are numerically integrated forward in time using implicit Newton-Raphson estimations for all terms except those in G. Time steps of 50-100 yr are used with a latitudinal resolution of 55.5 km (= 0.5’). Latitudinal boundary conditions are taken to be that h’ = hb at the model extremities of 74”N and 30”N. Three ice-age runs are shown in Figure 2 with progressively more complex versions

(5)

where pw is the density of the water, S is the current sea level (taken to be zero), and i is the spatial grid index increasing into the ice sheet. In practice the southern tip is only attacked in this way during deglaciations and only the outermost 50 to 150 km of ice sheet is affected. The calving mechanism also attacks the northern tip and maintains it at -72”N, but does not penetrate further into the ice sheet because of the high mountain centered at 7O”N. Equation (5) is an extremely crude parameterization of processes at a floating ice-sheet tip, and a more physical parameterization is developed below based on published models of ice

700

600

.BEDROCU

500

LAG

33

SIMULATIONS

400

300

200

100

0

ADOED

AGE

I lo3

yr-

B.P.

I

FIG. 2. Ice-age simulations, i.e., total cross-sectional area of ice in the (latitude, height) plane versus time. (a) With equilibrium bedrock (V = m), and without calving. (b) With bedrock lag (v = 100 km2 yr-I), and without calving. (c) With bedrock lag (V = 100 km’ yr-I). and with calving (Eq. (5)). The dotted curve in each panel is a new composite 6180 record assembled with a new chronology by the SPECMAP group (Imbrie ef al., 1982, 1983), and used here by permission of J. Imbrie and other SPECMAP group members. Its 6180 values are shown here in units of standard deviation (0). minus its modern value.

34

DAVID

of the basic model, and are compared with an oxygen-isotope record. In Figure 3a two simplifications are made to the model as formulated above: calving is suppressed and the bedrock is always in isostatic equilibrium with the ice load [i.e., I, = 2 in (2)]. The ice sheet responds directly to the orbital forcing, in fair correspondence with the higher-frequency oscillations in the Six0 record (similarly to Weertman, 1976); however, no trace of a lOO,OOO-yr cycle is produced. In Figure 2b bedrock lag is added using Eq. (2); here as in Oerlemans (1980) and Birchfield er al. (1981), the bedrock lag in (2) acts to amplify major retreats (for instance around 10,000, 130,000,220,000, and 340,000 yr BP) by the following mechanism: at a glacial maximum the bedrock depression is relatively deep, so that any subsequent retreat is accelerated by the increased ice melt at the lower elevations near the southern tip. In Figure 2b, however, this amplifying loop is not powerful enough to

250

POLLARD

produce complete deglaciation. Some complete deglaciations can be obtained by increasing the orbital sensitivity, X. but then the amplitude of the direct -20.000- and 40,000-yr responses becomes too great (Pollard, 1982a). The calving mechanism (5) is included in Figure 2c. Calving sets in at the southern ice-sheet tip at about the midpoint of each deglaciation, enabling the ice sheet to vanish completely without resorting to an overly large value of li. The records of several other model variables for this run are shown in Figures 3 and 4. Figure 4c shows how the southern ice-sheet tip plunges below sea level during deglaciations, briefly producing model ice-sheet environments quite unlike those during the remaining -90% of the run. Visual inspection of Figure 2 suggests that the addition of calving improves the overall agreement with the 6r80 records (e.g., in the depth of the deglaciations at around 405,000 and 5000 yr B.P.). Comparisons of

150

Age

100

(103yr

B.P.)

FIG. 3. Contours of ice thickness h (m) for the last 250,000 yr of the run shown in Figure 2c using the basic model with calving. Hatching at the ice-sheet edges shows where the tip regions are attacked by calving according to Eq. (5).

ICE-AGE

4000

700

600

500

35

SIMULATIONS

400

300

200

100

a

0

B

AGE

(IO3

yr-

B.P.

/

FIG. 4. Valriation of other model parameters for the ice-age simulation with calving. The dotted curve in each panel shows the ice-sheet cross-sectional area (as in Figure 2~). (a) External forcing, shown as the elevation of the equilibrium (zero-mass balance) line E at 55”N given by Eq. (4). (b) Maximum deflection of the bedrock surface below the undisturbed topography, h; - h'. The latitude of maximum deflection generally follows the midpoint of the ice sheet. (c) Elevation of the southern ice-sheet tip above sea level.

the power spectra also bear this out (Pollard, 1982a, b). A more objective test is shown in Table 1, giving linear correlation coefficients between each model curve and various &I80 records. For each iY8O record the calving-model run yields the highest value, presumably due to the improved simulation of the lOO,OOO-yr cycles. In Table 1 the correlation coefficients between a given model curve and the various WE0 records differ considerably, and the TABLE

1. LINEAR

CORRELATION

Model Equilibrium bedrock Bedrock lag added Calving added

COEFFICIENTS THREE WO

correlation coefficients between the 6180 records themselves fall as low as 0.39. Published WO records tend to depart seriously from each other before about 400,000 yr B.P., presumably due to uncertainties in the time scales of the deep-sea cores (e.g., Broecker and Van Donk, 1970). The calvingmodel curve happens to agree best with the new SPECMAP WE0 record; as shown in Figure 2c, the model deglaciations coincide with the SPECMAP major retreats all the

BETWEEN THE MODEL ICE-AGE DEEP-SEA CORE RECORDS

SIMULATIONS

IN FIGURE

2 AND

Emiliani (1978)

Johnson (I 982)

SPECMAP

SPECMAP*

0.29 0.28 0.34

0.22 0.45 0.58

0.32 0.63 0.75

0.73 0.82

0.39

0.67

-

-

0.51

6’“O

SPECMAP

Note. The SPECMAP label refers to the new composite record generated by members of the SPECMAP group (Imbrie el al., 1982, 1983; see Fig. 2 caption). All coefficients are calculated over the last 700,000 yr, except those labeled SPECMAP*. These are calculated over the last 650.000 yr (see text).

36

DAVID

way back to -600,000 yr BP. (Although the SPECMAP time scale was fine-tuned by matching with the -2O,OOO- and -4O,OOOyr frequencies of the orbital forcing (Imbrie et al., 1982. 1983), thus guaranteeing good correlation with the model at those frequencies, the timing of the observed major retreats was not tuned at all.) If the first 50,000 yr of each model run is ignored to allow for “spin-up” from incorrect initial conditions, the correlation between the calving run and the SPECMAP record is 0.82, the highest value in Table 1. The main disagreement between the model curve and the SPECMAP record in Figure 2c is the absence of strong deglaciations in the observations before -410,000 yr B.P. Perhaps this is related to the change in character observed in other VO records at around 900,000 yr B.P., before which the lOO,OOO-yr signal is markedly reduced (Pisias and Moore, 1981). As mentioned below, this may be due to long-term processes not included in the present model. For the comparisons with the deep-sea cores above, the model ice-sheet cross-sectional area (latitude vs height) is implicitly assumed to be proportional to ice-sheet volume, which requires that the longitudinal breadth is constant. An alternative assumption is that the breadth varies linearly with the latitudinal extent; in the case of a perfectly plastic ice sheet [with vertical thickness proportional to (horizontal extent)% (Weertman, 1976)], the volume would then be proportional to (cross-sectional area)Y3. However, this transformation causes little visual difference in Figure 2 and changes the values in Table I by ~0.05. ALTERNATE MODEL

TREATMENTS COMPONENTS

OF

The basic results in Figures 2 to 4 are reasonably stable against small parameter variations (Pollard, 1982a), and this type of model sensitivity is not pursued here. It is felt that more uncertainty stems from the variety of possible representations for some of the model components. In this section

POLLARD

three alternate treatments certain model components

of the more unare described.

Bedrock Defornzution The thin-channel asthenospheric model in Eq. (2) assumes a sudden large increase in viscosity at a depth of several hundred kilometers (e.g., Walcott, 1973). However, the profiles of viscosity and density in the mantle are not known with certainty (Peltier, 1981), and the most appropriate simple model may be significantly different from the thin-channel case. For instance. another conceivable model is an infinitely deep half-space, which would yield a fundamentally different dependence of response time on horizontal wavelength (Burgers and Collette, 1958), and might inhibit the bedrock lag/ice-sheet retreat amplification. (In the thin-channel case of Eq. (2), the wide depressions under the largest ice sheets are preserved for relatively longer times during subsequent deglaciations compared to smaller features.) The deep half-space case or other more complex mantle representations have not yet been tested in the present model. However, a simpler but less-physical approach, used by Oerlemans ( 1980) and Birchfield et al. (1981), is to assume a local response with no horizontal communication in the mantle. This enables the model to be tested using at least one different type of bedrock response. Equation (2) for the bedrock elevation, h’, is then replaced by ah’ -zz at

-1

(h’ - hb(x) + rh)

(6)

7

with T = 5000 yr. With this change, some further tuning to the crustal topography hi, was necessary mainly to raise midlatitude elevations so that calving was not initiated too easily. Figure 5 shows the resulting ice-age curve. Compared to the basic model result in Figure 2c, the higher-frequency oscillations are slightly larger but the major retreats are weaker (especially at -410,000 yr B.P.). However, the main features of the lOO,OOO-yr response and the overall agree-

ICE-AGE

37

SIMULATIONS

AGE

(IO3

yr

B.P.

I

FIG. 5. Ice-age simulation using the “local” bedrock response (Eq. (6)) and 7 = 5000 yr, with hf, joining the 6atitude (“N), height (m)) points (74. - 500), (70. 600). (30. 600). and with I?,, = 500 m in (4).

ment with the observed aI80 records are preserved, indicating that the model results are not critically dependent on the form of the bedrock deformation. Ice-Shelf Model

The parameterization of calving in Eq. (5) is based only on the concept that an unconfined floating ice-sheet tongue will discharge ice much faster than a nonfloating tongue. While this appears to be generally borne out by models and observations, a more physical approach would be desirable. Fairly simple models of ice shelves coupled to marine ice sheets have been developed and applied to the Ross ice-shelf sector in West Antarctica (Weertman, 1974; Thomas and Bentley, 1978; Denton and Hughes, 1981, Ch.7), and these are adapted below to the Ipresent model. Thomas (1977) used a similar model to investigate the dynamics of a calving bay in the St. Lawrence River valley during the last retreat of the Laurentide ice sheet; here a more extensive system of ice shelves or calving tongues is envisaged around most of the southern Laurentide perimeter during deglaciations, as proposed by Andrews (1973) and Denton and Hughes ~(1981, Ch. 2). The main effect of a large water mass next to an ice sheet is to eliminate the horizontal basal shear stress near the ice-sheet tip by floating the ice, and to replace it with the horizontal component of the hydrostatic pressure of water where it contacts the ice shelf. This results in a different flow regime in the ice shelf than in the ice-sheet interior, with most of the deformation being longitudinal and vertical extensive creep, i.e., i,,

- l ZZ, assuming no transverse variation (Leertman, 1957). In the models mentioned above, the two ice-flow regimes are coupled by considerations of the continuity across a narrow transition zone at the grounding line, the point marking the farthest penetration of the water under the ice, and where the ice first floats. The creep rate at the grounding line is given by (e.g. Denton and Hughes, 1981, p. 397)

where u is the horizontal velocity southward, hGis the ice thickness at the grounding line, g is the acceleration of gravity, 4 = l/s, n = 3, and B is a flow constant given below. C represents back pressure at the grounding line due to resistance at the iceshelf sides or at grounded ice rises within the shelf; here C = 0 to represent an unconfined ice shelf. The adaptation of Eq. (7) to the present finite-difference model involved somewhat arbitrary choices for the width of the transition zone and the exact position of the grounding line within the fixed horizontal grid; Figure 6 describes the method selected. Only the outer two ice-sheet grid boxes are shown, since as the ice shelf is assumed to be freely floating (i.e., C = 0) its extent beyond the grounding line does not need to be carried explicitly. The grounding-line thickness hG is taken to be that part of the outermost ice column below sea level, and the transition zone is taken to be the width of this column, AX. (How-

38

DAVID POLLAKD

;-Ax i FIG. 6. Finite-difference representation of the iceshelf parameterization at the grounding line (Eq. (7)). Only the outermost two grid boxes of the ice sheet are shown, as are the first two grid points occupied by the water body (hatched). MaI and Niare the ice thicknesses in each box, hG is taken as the ice thickness at the grounding line. Pl’and FG are fluxes of ice between grid boxes, S is sea level, and AX is the horizontal grid size.

ever, if this ice column is floating, i.e., if reP#O < Pw ki, then it is immediately moved from the ice sheet and the grounding line advances one grid point.) The flux of ice out of the ice sheet (and implicitly into the ice shelf), FG, is given by

tween the last two ice-sheet grid boxes obtained from (I). PO’ is the value of the flux out of the ice sheet given by (I) as if there were no water body, and is used in (8) as a device to avoid vanishingly small fluxes if hc, happens to be very small. Finally the budget of the outermost ice-sheet grid box is given as JhWat = (F(l) - F,)/bu. Figure 7 shows two runs using this iceshelf model, with different values of B in (7). (These values fall within the ranges used in the Ross ice-shelf models mentioned above.) The new mechanism acts very similarly to the cruder calving parameterization (5) with ice shelves forming only during the major deglaciations which then enable the ice sheet to vanish completely. The solid curve in Figure 7, with B = 2.5 bar yr%, is almost the same as the basic model result in Figure 2c. Nearly the same results are obtained with values of B up to about 5 bar yr”’ (dashed curve), beyond which the iceshelf drawdown becomes too slow to produce complete deglaciations. These results suggest that the calving/ice-shelf parameterizations and the rates of retreat produced in these simple models are consistent with the basic knowledge of the modern Antarctic ice shelves. Generalized

Melt\rlater

ParamettTri,-ation

The calving and ice-shelf mechanisms described above both act to accelerate the major retreats so that the ice sheets can vanish completely. They are triggered only after a maximum ice-sheet stand has caused a deep bedrock depression, followed by a fairly strong negative orbital anomaly that initiates a retreat of the southern tip into

(8) where & (hG) is given by (7). and II(‘) = F”)l 0.5(h(l) + h(O)) is the velocity midway be-

AGE

1103

yr

B.P.

I

FIG. 7. Ice-age simulations using the ice-shelf parameterization (Eqs. (7) and (8)). Full curve: with B = 2.5 bar yr”’ in (7). Dashed curve: with B = 5.0 bar yr”’ in (7).

ICE-AGE

SIMULATIONS

the depression. Alternative mechanisms may exist within the real climate system which are capable of performing the same function, but any candidate mechanism would have to be influenced by the rate of change of ice-sheet size (as opposed to the static influence of the ice sheet such as albedo feedback). For instance, if the meltwater layer from ice-sheet runoff in the North Atlantic could alter the climate so as to produce more ice wastage (Adam, 1975; Ruddiman and McIntyre, 1981), this feedback loop could amplify itself during major retreats to produce complete deglaciations. (The mechanism suggested by these authors is that the runoff freshens, stabilizes, and thins the oceanic mixed layer, causing (i) lower winter temperatures and more sea ice and less continental snowfall and (ii) warmer summer temperatures with more ice-sheet melt.) This particular example is examined quantitatively using a coupled climate/ice-sheet model in Pollard (1982b). Here a generalized parameterization of meltwater is used with the present model, mainly to estimate the strength of the feedback required for good ice-age simulations. The calving mechanism is suppressed and another term is added to (4) to shift the snow-budget pattern vertically, depending on the current net rate of change of icesheet volume (h. To do this a straightforward iterative procedure is needed to solve for one year’s climate. After some experimenting, relationships like that in Figure 8 were found to be necessary to produce good ice-age simulations; this is similar to those required in the meltwater case in Pollard (1982b). Here the snow-budget pattern is unaffected above a certain “trigger” value of v, but below that value (which occurs only during major retreats) it is raised sharply. The resulting ice-age simulation shown in Figure 9 is nearly the same as the basic run with calving in Figure 2c, with the jump in Figure 8 activated only during deglaciations at about the same times as the calving mechanism was. This jump in Figure 8 assumes a raising

39

I

5; -2501 30

I

I

20 Meltwater

I

IO rate

ti (cm yi’l

0 (a

-\j)

8. Relationship between the vertical shift of the snow-budget distribution and the net rate of change of ice-sheet size (V) imposed in Eq. (4) for the generalized meltwater parameterization. ii is equivalently expressed here as the thickness of the meltwater layer added per year to the North Atlantic C&f)given by - (3 x lo-’ m-t) x d(ice-sheet cross-sectional area)/&. FIG.

of the snow-budget pattern by 500 m, or equivalently a temperature increase of -3”C, for a change of $’ corresponding to 8 cm of annual meltwater input to the North Atlantic; furthermore this sensitivity must be somehow suppressed at values of $’ above the trigger value. As noted in Pollard (1982b), this relationship seems somewhat implausible, but could be tested with more sophisticated climate models. If the jump in Figure 8 is made much less steep, the feedback loop does not self-amplify when $’ falls below the trigger value in ice-age runs, and most of the retreats are incomplete. These results suggest that mechanisms attacking the ice sheet from “below,” i.e., calving or perhaps basal surging, are more promising than those from “above,” i.e., reduced snowfall or increased surface melt, in explaining the observed rapid deglaciations. ORBITAL CONTROL OF THE PHASE OF THE lOO,OOO-yr GLACIAL CYCLES

The phase of the observed lOO,OOO-yr glacial cycles correlates well with that of eccentricity over the last 600,000 yr, despite the relatively small direct effect of eccentricity on the orbital forcing (e.g., Imbrie and Imbrie, 1980, Fig. 10). This correlation is shown in Figure 10 where the S’*O record of core V28-239 from Shackleton and Op-

40

DAVID POLLARD

FIG. 9. Ice-age simulation using the generalized meltwater parameterization

dyke (1976) is plotted using a time scale linearly interpolated between the stageboundary ages given in Pisias and Moore (1981, Table 1). The correlation is apparently lost before 600,000 yr BP., but this could easily be due to uncertainties in dating the cores (N. G. Pisias, personal communication, 1982). The magnitudes of the individual glacial cycles are not correlated with those of eccentricity, and this is perhaps partly why cross-spectral analyses between the two signals are not totally convincing (Kominz et al., 1979; Moore et of., 1982); showing coherency peaks at 80,000 to 95,000 yr rather than at the dominant lOO,OOO- to 120,000-yr peaks in the power spectra. Figure 11 shows an extended run of the

as shown in Figure 8.

basic model with calving. The --lOO,OOO-yr correlation with the phase of eccentricity is apparent, with deglaciations and periods of little or no ice coinciding with locally positive departures of eccentricity. In contrast to the V28-239 record, this correlation is maintained by the model over the last two million years. The model phase does not depend on the choice of initial conditions, since other choices converge onto the same curve within one or two glacial cycles. This correlation and other related issues are examined below by running the model with a variety of forcing functions. First, Figure 12 shows that the model’s -lOO,OOO-yr response does not depend on the exact form of the forcing. The model

FIG. 10. Full curve: oxygen-isotope record of deep-sea core V28-239 from Shackleton and Opdyke (1976), relative to its most recent Sr*O value, and with a time scale obtained from Pisias and Moore (1981) (see text). Dotted curve: variations of eccentricity (Berger, 1978).

ICE-AGE

SIMULATIONS

AGE

t IO3

yr

B.P.

I

FIG. 11. F’ull curve: ice-age simulation of the basic model. Dotted curve: variations of eccentricity (Berger. 1978).

C

100

200

400

300

500

600

700

0

TIME

c IO3

>rl

FIG. 12. (a-e) Response of the basic model to sinusoidal forcing, in which the kAQ term in (4) is replaced by H sin(2arlT) where t is time and (a) T = 10,000 yr, H = 800 m. (b) T = 20,000 yr, H = 500 m. (c) T = 30,000 yr, H = 400 m. (d) T = 40,000 yr, H 400 m. (e) T = 50.000 yr, H = 300 m. (f) As above with white noise forcing (see text) and with I& = 650 m in (4).

42

DAVID

POLLARD

ice sheet has a natural tendency to respond to nearly any kind of forcing by growing slowly to a maximum size, at which point calving and deglaciation can be triggered. (With no forcing the model does not exhibit any internal oscillations, but see Appendix.) With sinusoidal forcing as in Figures 12a-e. the period of the main long-term response is always a multiple of the forcing period falling in the 80,000- to lOO,OOO-yr range. The amplitude of the forcing at each frequency has to be adjusted so that calving can be triggered, but the amplitude of the long-term response is still mainly determined by the internal dynamics. A “whitenoise” forcing is used in Figure 12f by replacing kAQ in (4) with a random term, rectangularly distributed between t 600 m and changed every 2000 yr. The model’s main response is still a - lOO,OOO-yr cycle as above. (Although not envisioned here as real stochastic forcing, this type of white noise could possibly be due to internal climatic oscillations having a period of a few thousand years, such as those modeled by Saltzman (1982) and Ghil and Le Treut (1981) involving sea ice and/or the deep ocean.) t:::::

With the actual orbital forcing, the phase of the - lOO,OOO-yr response could conceivably be (i) quasi-random, as found in the model of Oerlemans (1980); (ii) controlled by the direct influence of eccentricity on the annual mean insolation: (iii) controlled by the modulation of the amplitude of the precessional cycle due to eccentricity: and (iv) controlled by the combined forcing of obliquity and precession independently of eccentricity (Broecker and Van Donk, 1970). According to the present model. case (i) is ruled out since the response is independent of initial conditions and small parameter variations. Case (ii) is also ruled out since in model runs with fixed precession and obliquity, and insolation averaged over the full year in (4). the very small variations in AQ (on the order of (1 - eccentricity’)? yield a negligible model response. To examine case (iii), Figure 13 shows a run with obliquity held constant. so the only forcing is the -22,000-yr precessional cycle with its amplitude modulated by eccentricity. (It should be mentioned that the relative importance of obliquity and precession in simple models is influenced by the choice of latitude, here 55”N, at which the

:::

:::::+

n, - si 1

IS r in "

4

z 0 CJ 2100

2000

1900

FIG. 13. Full curve: ice-age simulation value. Dotted curve: ice-age simulation

1900

1700

1600

1500

1400

of the basic model but with obliquity held fixed with full orbital forcing (as in Fig. 11).

at its present

ICE-AGE

insolation is computed. However, the results presented here are still valid as general examples.) Compared to the full orbital case in Figure 11, some of the deglaciations are now missing, and several others are shifted by one precessional cycle, but the overall phase of the long-term cycles is maintained. Taken at face value this implies that the precessional modulation by eccentricity is indeed responsible for the phase of the glacial (cycles and the observed correlation with eccentricity. Figure 14, however, shows a run with eccentricity held fixed so the only forcing is due to obliquity and unmodulated precession (case (iv) above). Perplexingly, over the last 500,000 yr the phase of the -1OO,OOOyr response again agrees with the full-orbital case. Before 600,000 yr B.P. this agreement is generally lost, but unfortunately the uncertainty in the deep-sea core time scales before that time is too large to determine which model result is most realistic. We are left with two possibilities: (a) The real glacial cycles are controlled by eccentricity via the precessional modulation, and with improved dating the 6180

FIG. Dotted

14. Full curve: curve: ice-age

43

SIMULATIONS

records would show a phase-lock with eccentricity over the last two million years. In this case Figure 11 is a realistic simulation over that time, and the phase agreement in Figure 14 over the last 500,000 yr is coincidental. (b) The real glacial cycles are controlled by the combined obliquity and precessional cycles, independently of eccentricity. In this case the model result in Figure 11 overestimates the influence of eccentricity, and the phase agreement between the observed glacial cycles and eccentricity over the last 600,000 yr is coincidental. These two possibilities could be distinguished by improved dating of deep-sea cores before 600,000 yr B.P. in order to see if the phase-lock with eccentricity is maintained. In the case of V28-239, this may be hindered by the change in character of the record at about 900,000 yr B.P. before which the amplitude of the lOO,OOO-yr component is markedly reduced (Pisias and Moore, 1981). This change and other long-term variations, such as the -400,000-yr signal found in some cores (Moore et al., 1982), are not simulated in the model run of Figure 11.

ice-age simulation of the basic model but with simulation with full orbital forcing (as in Figure

eccentricity 11).

held fixed

at 0.02.

44

DAVID

Suggested causes of these variations are orogeny and erosion by the ice sheets themselves in the Northern Hemisphere (Pisias and Moore, 1981), and Antarctic ice-sheet fluctuations (Moore et al., 1982). These processes are beyond the scope of the present model, and Figure 11 shows only that the response to the orbital forcing alone has been relatively constant over the last two million years. CONCLUSIONS

Some retreat mechanism triggered only during deglaciations is required to account for the major ice-sheet retreats observed in the deep-sea core records. With the bedrock lag alone only relatively weak retreats are possible, hindering satisfactory ice-age simulations. The retreat mechanism is triggered only after the ice sheets have reached maximum size and caused a large bedrock depression, followed by a fairly strong negative orbital anomaly; the retreat mechanism is turned off when the ice sheet has shrunk considerably and/or the bedrock has rebounded sufficiently. The calving or ice-shelf mechanism in the present model is triggered by the submergence of the southern ice-sheet tip as it descends into the lagged bedrock depression, whereas the generalized meltwater mechanism is triggered by the net surface wastage over the whole ice sheet. Within the limitations of the present model, the calving mechanism seems the more plausible process. The model with calving adequately simulates the sudden deglaciations and the dominant lOO,OOO-yr power observed in the 6180 deep-sea core records, and yields an overall linear correlation coefficient of up to 0.8. The form of the asthenospheric flow model and its bedrock response makes some difference to the results, especially in the triggering of the calving mechanism. Although much the same ice-age curve is still obtained with a local form of the bedrock response, further testing would be desirable

POLLARD

using alternate models for the mantle response. The model responds with a - lOO.OOO-yr glacial cycle to nearly any kind of higherfrequency forcing. The amplitude of this cycle is mainly determined by the internal dynamics, Without any external forcing the model does not exhibit internal oscillations, but this depends on the imposed crustal topography and boundary conditions (see Appendix). With actual orbital forcing, the phase of the model’s -lOO,OOO-yr response is controlled mainly by eccentricity via its modulation of the precessional cycle. However, the same phase can be produced by obliquity plus precessional forcing alone with fixed eccentricity, but only over the last 500,000 yr. A definite identification of the main phase control of the observed glacial cycles is prevented by the uncertainty in the deep-sea core time scales before about 600,000 yr B.P. As cautioned in Pollard (1982a). the simplicity of the model should be kept in mind when considering these conclusions. The reduction to one horizontal dimension of the Keewatin-Hudson Bay-Labrador-St. Lawrence complex is particularly drastic, and future model developments should include east-west structure (e.g., Budd and Smith, 1981). A more explicit treatment of the proglacial lakes and marine incursion would also be valuable. APPENDIX: SELF-OSCILLATING WEST ANTARCTICA MODEL VERSION

Thomas and Bentley (1978) and Thomas (1979) have noted the possibility of internal oscillations in a marine ice sheet such as West Antarctica, involving lagged bedrock deformation and sudden collapses due to marine incursions. Part of such a cycle may have been the recent collapse of the West Antarctic ice sheet in the Ross Sea sector from the continental shelf to its present position, which is thought to have occurred over a 5000-yr interval sometime between 16,000 and 5000 yr B.P (Thomas and Ben-

ICE-AGE

tley, 1978; Denton and Hughes, 1981, Ch. 7). Similar internal oscillations are found when the present model is applied to West Antarctica. Only a brief description is given here due to the simplicity of the model (for instance, the sea level is fixed, the Transantarctic Mountains are ignored, and the calving parameterization effectively rules out any ice-sheet advance into a body of water). Budd and McInnes (1979) also found internal oscillations in an East Antarctic icesheet model with a period of -23,000 yr, involving surges due to basal meltwater lubrication. In applying the present model to West Antarctica, the north-south grid is still Cartesian (noit circularly symmetric, since ice shelves only occupy a small fraction of the Antarctic perimeter), and extends northward from 90%. The imposed equilibrium topography hb(x) passes piecewise-linearly through the following (latitude (“S), height (m)) points: (90, 2000), (70, 0), (65, - IOO), (60, .- IOOO), (4.5, - 1000). The region between 70”s and 6.5% is meant to represent the continental shelf. Boundary conditions are no flow across 90”s [i.e., a(h + h’)l& = 0 in (1) and a(h’ - hb + rh)ldx = 0 in (2)], and h’ = hb at the northern limit of 45”s. The calving mechanism (5) is used with the sea level S tuned to yield the best results, i.e., in the runs below S = - 100 m. The snow budget G was roughly fitted

to available Antarctic observations (e.g., Thomas and Bentley, 1978, Fig. 1). G is still given by (3) except with a = 0.67 x lO-4 yr-I, b = 0.22 x lo-‘m-l yr-l, and the upper value of 0.56 replaced by 0.05 m yr-I. The equilibrium-line altitude E is still given by (4) except with E, = 1000 m, s = 0.5 x 10-3, x,, corresponding to 70”s with x positive northward, and with no orbital forcing (k = 0). Figure 15 shows results for two values of the asthenospheric flow coefficient v in (2). The cycles in this figure are essentially those envisaged by Thomas and Bentley, and four stages of one cycle are shown in Figure 16. Stage (a) takes up about half of the cycle, in which the ice sheet is steadily growing outward over dry land. The bedrock elevation is relatively high, and is in the process of rebounding from the previous cycle; however, with the smaller value of u in Figure 15b the ice advance is slowed by contact with the ocean. Stage (b) takes up nearly the rest of the cycle; it is a period of quasi equilibrium in which the ice sheet extends to the edge of the continental shelf where a “sill” blocks the ocean. However, the sill is slowly being depressed by the weight of the ice. In stage (c) the sill is broached by the ocean, and the ensuing calving causes a rapid collapse of the ice sheet back across the lagged bedrock depression. In the model this collapse takes

0 0

100

200

45

SIMULATIONS

400

300 TIME

1103

500

600

t 700

,,rl

FIG. 15. Model response (i.e., total cross-sectional area of ice versus time) for the West Antarctic version. (a) Asthenospheric flow coefficient II = 100 km2 yr-1. The period of the oscillation is 56,000 yr. (b) Asthenospheric flow coefficient v = SO km* yr-1. The period of the oscillation is 90.000 yr.

46

DAVID POI.I.AKD

4000 d

482.000 yr

3000

LATITUDE

lost

LATITUDE

1051

FIG. 16. Cross sections of latitude versus height at various times in one of the cycles in Figure ISa, corresponding to the four stages described in the text. Shown are the ice-sheet surface (dashed liner, bedrock surface (full line) and the undisturbed topography /rh (dotted line). The shaded regions show the water body that can cause ice calving according to Eq. (5).

about 4000 to 6000 yr. Finally, in stage (d) a smaller core ice sheet is left on higher ground near the pole, and bedrock to the north is starting to rebound in preparation for the next cycle. In further experiments the period of the oscillation was found to increase smoothly at first as v was reduced (i.e., values of v = 200, 100, and 50 km? yr-i produced periods of 36,000, 56,000, and 90,000 yr, respectively). With further reduction of v from 50 to 12 km” yr-I, the period values became irregular, and sometimes aperiodic fluctuations were found with time scales of 40,000 to 100,000 yr. The results shown here are only preliminary quantitative examples of the scenarios described by Thomas Bentley (1978) and Thomas (1979). As they pointed out, these Antarctic oscillations would probably be modified by the -100-m sea-level changes due to variations of the ice sheets in the Northern Hemisphere and vice versa. This would form an interactive connection between the ice sheets of the two hemispheres

that could possibly play an important in the dynamics of the ice ages.

role

ACKNOWLEDGMENTS 1 am grateful to John Imbrie and other members of the SPECMAP group for supplying a digital record of their new oxygen-isotope curve, and generously allowing its use in Figure 2 and Table 1 prior to their own publication. Robert Johnson and Niklas Pisias kindly supplied digital records of other oxygen-isotope curves used above. Also I thank Niklas Pisias for a helpful discussion. Michael Schlesinger for suggesting the white-noise experiment in Figure 12. and Naomi Zielinski for typing the manuscript. This research was supported by the Climate Dynamics Section of the NSF under Grant ATM-8019762.

REFERENCES Adam, D. P ( 1975). Ice ages and the thermal equilibrium of the Earth, II. Qu7ternar.v Resrmrcl~ 5, 161-171. Andrews, J. T. (1973). The Wisconsin Laurentide ice sheet: Dispersal centers, problems of rates of retreat, and climatic implications. Arcric crnd Alpine Research

5, 185-199.

Benzi, R.. Parisi. G. Sutera, A., and Vulpiani. A. ( 1982). Stochastic resonance in climatic change. Tellus 34, 10-16.

ICE-AGE

Berger, A. L. (1978). Long-term variations of daily insolation and Quatemary climatic changes. Journal of the Atmospheric

Sciences

35, 2362-2361.

Birchfield, G. E., Weertman, J., and Lunde, A. T. (1981). A paleoclimate model of Northern Hemispheric ice sheets. Quaternary Research 15, 126-142. Broecker, W. S., and Van Donk, J. (1970). Insolation changes, ice volumes and the OL8record in deep-sea cores. Review of Geophysics und Space Physics 8, 169-198. Budd, W. F. (1981). The importance of ice sheets in long term changes of climate and sea level. International Association of Scientific Hydrologists Publ. No. 131. 441-471. Budd, W. F., and McInnes, B. J. (1979). Periodic surging of the Antarctic ice sheet-An assessment by modelling. Hydrological Sciences Bulletin 24, 95-104. Budd. W. F.. and Smith, 1. N. (1981). The growth and retreat of ice sheets in response to orbital radiation changes. Int. Assoc. Hydrol. Sci. Publ. No. 131, 369-409. Burgers, J. M., and Collette. B. J. (1958). On the problem of the postglacial uplift of Fennoscandia I and II. Proceed’ings of the Koninklijke Nederlandse Akudemie

v’on Wetenschuppen

B 61, 221-241.

Cathles. L. M., III. (1975). “The Viscosity of the Earth’s Mantle.” Princeton Univ. Press, Princeton, N.J. Chappell, J. (1974). Relationships between sea levels, I80 variations a!nd orbital perturbations, during the past 250,000 years. Nature (London) 252, 199-202. Dansgaard. W., and Tauber, H. (1969). Glacier oxygen18 content and Pleistocene ocean temperatures. Science

166, 499-,502.

Denton. G. H., and Hughes, T J. (Eds.) (1981). “The Last Great Ice Sheets.” Wiley, New York. Emiliani, C. (1978) The cause of the ice ages. Eurth and Planetary

and Science

Letters

37, 349-352.

Ghil. M., and Le Treut. H. (1981). A climate model with cryodynalmics and geodynamics. Journal of Geophysical

Research

47

SIMULATIONS

86, 5262-5270.

Hays, J. D., Imbrie. J.. and Shackleton, N. J. (1976). Variations in the Earth’s orbit: Pacemaker of the ice ages. Science 194, 1121-l 132. Imbrie, J., and Imbrie. J. Z. (1980). Modeling the climatic response to orbital variations. Science 207, 943-953. Imbrie. J.. and other SPECMAP Project authors (1982). Chronology of ‘oceanic 6180 variations. EOS Trunsactions 63, 45, 995-996 (Abstract). Imbrie, J., and other SPECMAP Project authors (1983). Geological arguments supporting the orbital theory of Pleistocene (climates. In “Milankovitch and Climate: Understanding the Response to Orbital Forcing” (A. Eterger et al.. eds.), Reidel, Boston, in press. Johnson, R. G. (1982). Brunhes-Matuyama magnetic

reversal dated at 790,000 yr B.P. by marine-astronomical correlations. Quaternary Research 17, 135-147.

Kominz, M. A.. Heath, G. R., Ku, T. -L., and Pisias, N. G. (1979) Brunhes time scales and the interpretation of climatic change. Earth and Planetary Science Letters 45, 394-410. Moore, T. C., Jr., Pisias. N. G.. and Dunn, D. A. (1982). Carbonate time series and the Quaternary and late Miocene sediments in the Pacific ocean: A spectral comparison. Marine Geology 46, 217-233. Nicolis. C. (1982). Stochastic aspects of climatic transitions-Response to a periodic forcing. Tellus 34, l-9. Oerlemans, J. (1980). Model experiments on the lOO,OOO-yr glacial cycle. Nature (London) 287, 430-432. Peltier, W. R. (1981). Ice age geodynamics. Annual Review of Earth and Planetary Sciences 9, 199-225. Pisias, N. G.. and Moore, T. C.. Jr. (1981). The evolution of Pleistocene climate: A time series approach. Eurth and Plunetury Science Letters 52, 450-458. Pollard, D. (1982a). A simple ice sheet model yields realistic 100 kyr glacial cycles, Nature (London) 296, 334-338 Pollard, D. (1982b). A coupled climate-ice sheet model applied to the Quaternary ice ages. Climatic Research Institute. Corvallis, Oregon. Report No. 37. Also submitted to Journal of Geophysical Research. Pollard, D., Ingersoll, A. P.. and Lockwood, J. G. (1980). Response of a zonal climate-ice sheet model to the orbital perturbations during the Quaternary ice ages. Telltcs 32, 301-319. Ruddiman, W. F.. and McIntyre, A. (1981). Oceanic mechanisms for amplification of the 23.000-year icevolume cycle. Science 212, 617-627. Saltzman, B. (1982). Stochastically-driven climatic fluctuations in the sea-ice, ocean temperature, CO? feedback system. Tellus 34, 97-112. Schneider, S. H., and Thompson, S. L. (1979). Ice ages and orbital variations: Some simple theory and modelling. Quaternary Research 12, 188-203. Sergin. V. Ya.. (1979). Numerical modeling of the glaciers-ocean-atmosphere global system. Jonrnal of Geophysical Research 84, 3191-3204. Shackleton. N. J., and Opdyke, N. D. (1976). Oxygenisotope and paleomagnetic stratigraphy of Pacific core V28-239 late Pliocene to latest Pleistocene. Geological

Society

of

America,

Memoir

145,

449-464. Suarez, M. J., and Held, 1. M. (1979). The sensitivity of an energy balance climate model to variations in the orbital parameters. Journal qf Geophysical Research

84, 4825-4836.

Thomas, R H. (1977). Calving bay dynamics and ice sheet retreat up the St. Lawrence valley system. Geographie Physique et Quuternaire 31, 347-356.

48

DAVID

Thomas, R. H. (1979). The dynamics of marine ice sheets. Journal of Glac,iology 24, 167-177. Thomas, R. H.. and Bentley. C R. (1978). A model for Holocene retreat of the West Antarctic ice sheet. Quaternary Research 10, 150-170. Walcott, R. I. (1973). Structure of the Earth from glacioisostatic rebound. Annual Re~?e~~~ of Earth and Plunrtary Science 1, 15-37.

POLLARD Weertman. J. ( 1957). Deformation of floating ice shelves. Jownal ofG/ac~iology 3, 38-41. Weertman. J. (1974). Stability of the junction of an ice sheet and an ice shelf. Journal o,f‘G/rrt~io/o~y 13, 3-l I. Weertman, J. (1976). Milankovitch solar radiation variations and ice age ice sheet \ires. Nature (London) 261, 17-20.