## References

Boulton, G. S. and Clark, C. D.
1990. A highly mobile Laurentide ice sheet revealed by satallite images of glacial lineations. Nature.346(6287). 813–817.

Crank, J.
1984. Free and moving boundary problems. Oxford, Clarendon.

Dansgaard, W. and Johnsen, S. J.
1969. A flow model and a time scale for the ice core from Camp Century, Greenland. J. Glaciol., 8(53),215–223.

Firestone, J., Waddington, E. and Cunningham, J.
1990. The potential for basal melting under Summit. Greenland. J. Glaciol., 36(123) 163–168.

Greuell, W. and Oerlemans, J.
1989. The evolution of the englacial temperature distribution in the superimposed ice zone of a polar ice cap during a summer season. In Oerlemans, J., ed. Glacier fluctuations and climatic change. Dordrecht, etc., Kluwer Academic Publishers, 289–303.

Grigoryan, S.S., Krass, M. S. and Shumskiy, P. A.
1976. Mathematical model of a three-dimensional non-isothermal glacier. J. Glaciol., 17(77), 401–417.

Hooke, R.LeB
1977. Basal temperatures in polar ice sheets: a qualitative review. Quat. Res., 7(1), 1–13.

Hughes, T.
1973. Glacial permafrost and Pleistocene ice ages. In Permafrost. Second International Conference, 13 28 July 1973. Yakutsk. U.S.S.R. North American Contribution. Washington, DC. National Academy of Sciences, 213–223.

Huybrechs, P.
1994. Present evolution of the Greenland ice sheet: an assessment by modelling. Global and Planetary Change, 9(1–2), 39–51.

Jenssen, D.
1977. A three-dimensional polar ice-sheet model. J. Glaciol., 18(80), 37–389.

Jóhannesson, T., Raymond, C. and Waddington, E.
1989; Time-scale for adjustment of glaciers to changes in mass balance. J. Glaciol., 35(12), 355–369.

Kleman, J.
1994, Preservation of landforms under lce sheets and ice-caps. Geomorphology, 9. 19–32;

Kleman, J. and Borgström, I.
1994. Glacial land forms Indicative of a partly frozen bed. J. Glacol., 40(135). 255–264.

Letréguilly, A., Reeh, N. and Huybrechts, P.
1991. The Greenland ice sheet through the last glacial–interglacial cycle. Palaeogeogr.. Palaeoclimatol.. Palaeoecol..90(4). 385–394.

MacAyeal, D. R.
1993a. A low-order model of the Heinrich event cycle, Paleoceanography. 8(6). 767–773.

MacAyeal, D. R.
1993b. Binge/purge oscillations of the Lanrentide ice sheet as a cause of the North Atlantic’s Heinrich event. Paleoceanography. 8(6), 775–784.

Nye, J. F.
1951. The flow of glaciers and ice-sheets as a problem in plasticity. Proc. R. Soc. London. Ser. A, 207(1091). 554–572.

Nye, J. F.
1957. The distribution of stress and velocity in glaciers and ice-sheets. Proc. R. Soc. London. Ser. A, 239(1216), 113–133.

Nye, J. F.
1963. Correction factor for accumulation measured by the thickness of the annual layers in an ice sheet. J. Glaciol., 4(36), 785–788.

Oerlemans, J.
1991. The mass balance of the Greenland ice sheet: sensitivity to climate change as revealed by energy-balance modelling. Holocene. 1(1), 40–49.

Paterson, W. S. B.
1994. The physics of glaciers. Third edition. Oxford. Elsevier.

Ritz, C.
1987. Time dependent boundary condition for calculation of temperature fields in ice sheets. International Association of Hydrological Sciences Publication 170 (Symposium at Vancouver 1987 The Physical Basis of Ice Sheet Modelling). 207–216.

Schiesser, W. E.
1991. The numerical method of lines: integration of partial differential equations. San Diego. CA. etc., Academic Press.

Shilts, W.E.
1991. Flow patterns in the central North America ice sheet. Nature, 286(5770). 213–218.

Sugden, D. E.
1977. Reconstruction of the morphology, dynamics and thermal characteristics of the Laurentide ice sheet at its maximum. Arct. Alp. Res., 9(1). 21–47.

Trewartha, G. T. and Horn, I. H.
1980. An Introduction to climate. Fifth edition. New York. etc., McGraw-Hill.

Washburn, A. L.
1980. Geocryology. New York. etc., John Wiley and Sons.

## APPENDIX Transformation of the Governing Equation

The one-dimensional energy balance considered is:

where *z* is positive upward with its origin at the rock/ice interface, *w* is the vertical velocity and both the heat capacity, *ρc*, and the thermal conductivity, *k*, are, in general, functions of temparature. The heat flux at the base of the domain, located at some large depth *z* = –*h*
_{0}, remains fixed at *q*
_{0}:

The upper boundary, or the ice/air interface, is located at *z = h*, where h is an arbitrarily specified function of time. The surface tempaerature, *T*
_{s}, is also an arbitraarily specified function of time:

In the present case, the time-dependence of *T*
_{s} is given indirectly through its dependence on *h*(*t*).

The vertical velocity in the varies linearly through the upper three-quarters of the ice sheet and quadratically below that (Dansgaard and Johnsen, 1969):

Where *z*
_{m} is the elevation above the bed at which the upper, linear velocity profile intersects the lower, quadratic profile. (In all calculations shown in this paper, we take *z*
_{m}/*h* = 0.25) In the limit that *z*
_{m} → 0, this recovers the linear profile of Nye (1951, 1957, 1963). *W* is the magnitude of the advective velocity at the surface. Here, *W* is the volume flux of ice across the upper surface of the ice sheet, i.e. it represents the accumulation rate over and above that which goes into growth of the ice sheet. Of course, a divergent, lateral flow must balance the downward flow of ice. An assumption implicit in Equation (A4) is that the velocity field in the ice is able to adjust rapidly to increasing *h* in order to maintain *W* in the form given throughout the growth history.

The “front-fixing” method entails coordinate transformation ζ = *z*/*h*, in terms of which Equation (A1) can be written for the domain within the ice (0 ≤ ζ ≤ *h*) in the form:

where *w* (Equation (A4)) is now written in terms of ζ.

Note that, even in the absence of downward advection in the fixed reference frame (*w* = 0), the coordinate transformation has the effect of introducing into Equation (A5) an advection-like term that accounts for the moving boundary, *h(t)*. The effect of the boundary moving upward at velocity *dh/dt* and that of advection cold ice across the surface at velocity *−W* are additive. The boundary condition at the top surface (Equation (A3)) is now applied at the *fixed* point ζ = 1, which is convenient for numerical implementation. In the scheme adopted here, we solve the conduction problem in the underlying rock domain (–*h*
_{0} ≤ *z* ≤ 0) in the untransformed coordinate, *z*, while solving Equation (A5) for the ice. The temperature and heat flux are matched at the rock/ice interface, *z* = ζ = 0.

We consider two particular forms of the growth function, *h*(*t*):

Where *h*
_{x} is the final ice thickness and *t*_{r}
the rise time. In the text, we refer to Equation (A6) as an exponential growth modal and Equation (A7) as a linear modal.