Introduction
The base of an ice sheet is either frozen to its bed (cold-based) or at the pressure-melting point (warm-based). A cold-based ice sheet moves by internal deformation and causes little erosion at its bed, while a warm-bused ice sheet moves by both internal deformation and basal sliding, thereby eroding the surface it over-rides. According to theory, an ice sheet becomes warm-based when the geothermal heat flux and internal heating by viscous dissipation provide enough heat to melt the base of the ice sheet, while the overlying ice mass insulates the base against the very cold air at the surface (Reference PatersonPaterson, 1994).
Determination of the basal conditions of ice sheets is of paramount importance for the reconstruction of their topography and flow charactertstics, MacAyeal (Reference MacAyeal1993a,Reference MacAyealb) suggested that a periodic change in basal conditions of the Laurentide ice sheet from cold-to warm-based could result in a dramatic thinning of the ice sheet and a massive outflow of icebergs into the North Atlantic Ocean.
The thermodynamic equations for ice divides have been analyzed by Reference Grigoryan, Krass and ShumskiyGrigoryan and others (1976). Grigoryan and others, as well as most previous ice-sheet modelers (e.g. Reference HookeHooke. 1977: Reference JenssenJenssen. 1977; Reference SugdenSugden, 1977) analyzed a steady-state temperature field, as this allows great simplification of the governing equations. A number of studies have examined the effect of climatic change on ice sheets but in all cases the ice-sheet thickness is more or less constant (e.g. Reference Greuell, Oerlemans and OerlemansGreuell and Oerlemans. 1989; Reference Firestone, Waddington and CunninghamFirestone and others, 1990; Reference Letréguilly, Reeh and HuybrechtsLetréguilly and others, 1991; Reference OerlemansOerlemans, 1991; Reference HuybrechsHuybrechts, 1994). Most of these studies examine the Greenland ice sheet, which has persisted throughout numerous glacial cycles. However, models with a steady-state temperature distribution or a constant thickness are unrealistic for the continental ice sheets of the last glaciation, which persisted only for relatively short perids of time (50 000–100 000 years). The thermal diffusion time (t D ≈ L 2/κ ≈ (107 m2)/(10–6 m2/s) = 1013s ≈ 3 × 105 years; with L the ice-sheet thickness and κ the thermal diffusivity of ice) is on the same orders as the lifespan of these ice sheets.
We use a finite-difference analysis to simulate the temperature profile at the center of an ice sheet throughout the course of a glaciation. Αt the center of the ice sheet, horizontal flow does not exist, thus simplifying our calculations in a one-dimensional model.
Model
The evolution of the thermal profile in a growing ice sheet is modeled as a one-dimensional, transient process. The governing equations are given in the Appendix. The model accounts for conduction in both the ice and the underlying rock, as well as advective transport in the ice due to the downward moition of accumulating material. The thermal conductivity and heat capacity of the ice are allowed to depend upon temperature (although this non-linearity proves to have a minimal effect in the present case), while those characteristic of the rock are assumed to be uniform and constant. To model the radial outward flow of ice in the central region of the ice sheet, the advective velocity of the ice varies from a specified value at the upper surface to zero at the base of the ice sheet. Rather than using a linear decrease in velocity (Reference NyeNye. 1951, Reference Nye1957, Reference Nye1963). we used the more realistic (personal communication from E. Waddington, 1995) model of Reference Dansgaard and JohnsenDansgaard and Johnsen (1969), in which the vertical velocity decreases linearly through the upper three-quarters of the ice sheet and then quadratically to the base. This accounts approximately for the lack of basal sliding in a cold-based ice sheet. Our model allows for a moving boundary at the top of the growing ice sheet; i.e. the time-dependent position of the boundary is specified as an input function. In our model, the moving boundary is specified independent of the downward advection, even though the two are related in real ice sheets. Our downward advection is the part of the total accumulation rate that is accomodated by divergent lateral flow in the ice sheet and therefore does not result in a thickening of the ice sheet.
At the start of a model run. the sub-surface temperature are in equilibrium with the geothermal heat flux and a given ground-surface temperature. The geothermal heal flux is applied at a depth of 4000 m. Model runs with different depths show that beyond 2000 m depth negligible changes in the temperature profiles in the ice are found. This observation is consistent with that of Reference RitzRitz (1987). Ice is deposited on the ground surface at the surface temperature. Subsequent ice is deposited at cooler temperatures, according to the adiabatic lapse rate, as the ice sheet grow s in elevation. Our model uses an adiabatic lapse rate of 8° C/1000 m, halfway between the values for dry and moisture-saturated air (Reference Trewartha and HornTrewartha and Horn, 1980). The parameters used in the model are given in Table 1. The geothermal heat flux warms bedrock and ice until basal melting occurs at –1.96°C (for a 3000 m thick ice sheet).
The model does not account for the latent heat of the melting of permafrost, This effect would be small for substrates with low porosities, and therefore a small water content, such as bare granitic bedrock as occurs in the Canadian shield. Regionally, sedimentary rocks, soil cover, etc. may result in significant ice content, which would require considerable energy for the phase change from solid to liquid. Since we model the ice sheet at its center, we do not consider the frictional heat created by ice movement, which would warm the ice. Therefore, the model results are not applicable for areas with significant horizontal ice velocities.
Numerical Analysis
The model problem involves a moving boundary, which poses a challenge for many numerical schemes typically used in condition heat transfer. In the present case, we employ a “front-fixing” method (e.g. Reference CrankCrank, 1984), in which a tranformation of the spatial coordinates fixes the moving boundary (see Appendix). Then the surface-boundary condition can be applied at a fixed nodal point in the numerical analysis. The effect of this tranformation is to introduce a new term in the govering equaiton that takes the form of an advective contribution to the transport, with the apparent velocity sealing with the rate of change of the boundary location.
The numerical analysis is carried out by the method of lines (MOL) (e.g. Reference SchiesserSchiesser. 1991). In this method, the spatial derivatives are discretized by a Central finite difference scheme, reducing the partial differential equation at each computational node to a first-order, ordinary differential equation (ODE) in time. Each ODE is coupled to the ODEs at adjacent nodes through the discretization of the spatial gradients. Ihe entire array of coupled, non-linear ODEs is passed to a highly optimized solver that performs the time integration by a backward difference scheme. The code has been tested against various exact, analytical solutions and proves to maintain excellent accuracy.
Estimate of the Time-Scale Until Basal Melting Occurs
Figure 1 shows the evolution of the temperature profile over time for a surface temperature of –12.5° C. Ice-sheet growth is exponential with fast growth at the inception and decreasing with lime. The characteristic rise time is 6000 years (i.e. at this time 63% of the final thickness has been reached) and the final ice-sheet thickness is 3000 m. Note that at 0 years no ice sheet exists, while the ground surface is at –12.5°C, which is about the current temperature in central Labrador and around the northern part of Hudson Bay (Reference WashburnWashburn. 1980). After 10 0000 years, the ice sheet has reached about 2400 m thickness and is still growing. The ice sheet between 2000 and 2500 m above the original surface actually is cooled between 10 000 years and 30 000 years by downward advection of colder ice, while the temperature signal from the basal heating has not yet propagated to the surface. The basal melting point is reached after 72 000 years. This model run will be used in all comparison as a reference.
Table 1. Ranges of parameters used to model ice-sheet temperatures
Fig. 1. Evolution of the temperature profile in an ice sheet. Initial surface temperature: –12.5°C; final ice thickness: 3000 m; downward advection: 0.025 m a−; rise time: 6000 .
Figure 2 shows the change of the basal temperature over time for ice sheets with different initial surface temperatures (and therefore initial ice and sub-surface temperatures). The warmer the ice and substrate, the more quickly the base melts. It is noteworthy that even with a relatively warm initial surface temperature of –7.5°C basal melting does not occur until 15 000 years after inception of the ice sheet. A temperature of –7.5°C may be unrealistically warm for the conditions at the inception of the Laurentide ice sheet, as this is cosiderably warmer than present ground temperatures. Reference HughesHughes (1973) pointed out that the Laurentide ice sheet developed on permafrost and probably displayed relatively, cold initial temperatures. An ice sheet with an initial surface temperature of –15°C remains cold-based throughout a 100 000 year glacial cycle.
Fig. 2. Influence of final thickness of an ice sheet on basal temperature. Note that pressure-melting point varies with ice-sheet thickness. Initial surface temperature: –12.5°C; downward advection: 0.25 m a−1; rise time: 6000 a.
The influence of the final thickness of the ice sheet on basal temperature is examined in Figure 3. A thicker ice sheet provided more insulation, resulting in a warmer base. This factor is not important for ice sheets that are thicker than 2000 m. An ice sheet thinner than 2000 m signigicantly extends the time until basal melting occurs. A 1000 m thick ice sheet under the given conditions probably would never become warm-based. Note that the pressure-melting point varies with ice-sheet thickness.
Fig. 3. Influence of final thickness of an ice sheet on basal temperature. Note that pressure-melting point varies with ice-sheet thickness. Initial surface temperature: –12.5°C; downward advection: 0.25 m a−1; rise time: 6000 a.
Figure 4 shows the influence of downward advection of cold on basal temperatures. As ice is accumlated at the surface and flows out radially from the center of the ice sheet, the downward movement of cold ice cools the ice sheet. The values commonly assumed for downward advection at the center of the Laurentide ice sheet are between 0.025 and 0.05 m a1 (Reference SugdenSugden, 1977). Downward advection of ice significantly coold the base of the ice sheet. An ice sheet with a downward advection rate of 0.075 m a11 remains frozen to its bed for 240 000 years.
Fig. 4. Influence of downward advection on basal temerature. Note that 0 m a1 means no ice flow. Final ice thickness: 3000 m, initial surface temperature: –12.5°C; rise time: 6000 a.
The growth history of the ice sheet does not significantly influence the time at which basal melting occurs (Fig. 5). For the reference calculations, we used an exponential growth curve (see Appendix) with a large growth rate at the onset, decreasing exponentially as the ice sheet approaches its final thickness. We also tested a constant growth rate of the ice sheet until the final thickness is reached and zero growth thereafter (linear growth). Calculations for linear growth ever 10 000 and 50 000 years do not show large differences in the time at which basal melting occurs, The ice sheet seems to grow faster than the propagation of the temperature signal in the ice, so that the heat sink at the ice-sheet surface only weakly affects the basal temperatures. Therefore, the specifics of the ice-sheet growth history can be neglected to a certain degree when modeling the temperature history, as long as the ice-sheet growth does not affect downward advection. In real ice sheets, the accumulation rate determines both the growth rate and the rate of downward advection (Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others, 1989). As shown above, downward advection rates exert a strong influence on basal temperatures.
Fig. 5. Influence of the growth history of an ice sheet on basal temperature. Examples of exponential and linear models (see Appendix) with various rise times shown. Note that in the model the rise time is independent of downward advection. The curces for 1 a, 6000 a (exponential) and 10 000 a (linear) are indistinguishable. Initial surface temperature: –12.5°C; final ice thickness: 3000 m; downward advection: 0.025 m a.−1
Results
For ice sheets that disappear during interglacials, the model calculations show that the basal temperature is Strongly time-dependent on the time-scale of a typical glacial cycle. Furthermore, heating due to the geothermal flux cannot account for basal melting at the center of an ice sheet during most or all of a glacial cycle if the initial ground temperature is colder than about –10°C. Basal melting at the center of an ice sheet may occur toward the end of a glacial Cycle lasting 50 000–100 000 years. Even under the most unfavorable conditions, the central part of an ice sheet takes tens of thousands of years to reach the pressure-melting point at the base. The basal temperature is determined mainly by the initial temperature of the surface (and substratum) and the advection of heat due to downward movement of ice. The thickness of the ice sheet does not influence the basal temperatures for values above 2000 m. Our results show that the growth history of the ice sheet can be chosen freely without substantially affecting the results over the time-scale of interest, as long as there is no concurrent change in downward adveciion. This uncoupling of the growth history from the downward advection is unrealistic for real ice sheets but is a useful model assumption.
Discussion
The central parts of the Pleistocene ice sheets probably were cold-based throughout a major part of a glacial cycle. Outward from the center, ice movement increases and therefore internal heating of the ice contributes to the thermal balance. Basal melting would first occur in areas of high ice velocity (ice streams, convergent flow between topographic obstacles), while other parts may never have reached the pressure-melting point. Therefore, we suggest that the basal conditions of ihe ice sheets were spatially heterogeneous and probably strongly influenced by topography. Exploration of these phenomena requires a two- or three-dimensional model which includes a transient thermal analysis.
Rcceul observations of preserved periglacial features in areas covered by the Fennoscandian ice sheet indicate that the ice sheet remained cold-based in places throughout its history (Reference KlemanKleman, 1994: Reference Kleman and BorgströmKleman and Borgström. 1994). The patchiness of the formerly cold-based areas suggests that basal melting occurred only toward the end of the glaciation and that the preserved preglacial land forms are remnants of larger areas with cold-based conditions. In other areas, despite evidence for a warm-based ice sheet, these authors observed that erosion of preglacial features is much less than expected for a glacial cycle lasting tens of thousands of years. As indicated by our results, the ice sheet may have been cold-based during most of the glaciation, while basal melting occurred in some areas toward the end of the glacial cycle, resulting only in slight erosion before the ice sheet melted completely.
Lineations, striations and other features associated with basal sliding that are used to reconstruct the topography of former ice sheets (e.g. Reference ShiltsShilts. 1980: Reference Boulton and ClarkBoulton and Clark. 1990) may have formed toward the end of the glaciation. Ice sheets reconstructed from these features may represent only the end of a glacial cycle. Alternatively, these features may have formed during earlier glaciations and have been preserved under a cold-based ice sheet.
Λ model of periodic basal melting and refreezing caused by thickening and thinning of the Laurentide ice sheet (so-called binge–purge cycles) was proposed by (Reference MacAyealMacAyeal 1993a, Reference MacAyealb) to explain periodic massive iceberg discharge into the North Atlantic Ocean (Heinrich events). Our results do not contradict this model. Once basal melting occurs, the ice sheet can oscillate between cold- and warm-based conditions. The possible occurrence of Heinrich events during the early Wisconsin could indicate that Heinrich events originate in peripheral regions of the ice sheet, while the Center remains cold-based. If the ice sheet in its entirely needs to he invoked to cause binge–purge cycles, our model results suggest that these early events were not caused by binge–purge cycles. On the other hand, the early Heinrich events could be interpreted as an indication dial surface temperatures at the inception of the Laurentide ice sheet were warmer than −7.5°C and that basal inching of the entire ice sheet occurred very early in the glacial cycle.
Acknowledgements
We thank B. Hallet and P. Clark for their comments and ideas on this project, as well as Zhilin Li and M. R. Baer for discussion of the numerical model. The comments of two anonymous reviewers. E. Waddington and D. MacAyeal were helpful to improve this paper.
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 (Reference Dansgaard and JohnsenDansgaard 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 (Reference Nye1951, Reference Nye1957, Reference Nye1963). 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 tr the rise time. In the text, we refer to Equation (A6) as an exponential growth modal and Equation (A7) as a linear modal.
