Numerical modeling investigations of the subglacial conditions of the southern Laurentide ice sheet

Abstract

Sub- and proglacial bed conditions influence advance and retreat of an ice sheet. The existence and distribution of frozen ground is of major importance for better understanding of ice-flow dynamics and landform formation. The southern margin of the Laurentide ice sheet (LIS) was dominated by the presence of relatively thin ice lobes that seem to have been very sensitive to external and internal physical conditions. Their extent and dynamics were highly influenced by the interaction of subglacial and proglacial conditions. A three-dimensional thermomechanical ice-sheet model was coupled with a model for the thermal regime in the upper Earth crust. The model has been applied to the LIS in order to investigate the spatial distribution of thermal conditions at the bed. The evolution of the whole LIS was modeled for the last glacial cycle, with primary attention on correct reconstruction of the southern margin. Our results show extensive temporal and spatial frozen ground conditions. Only a slow degradation of permafrost under the ice was found. We conclude that there are significant interactions between the ice sheet and the underlying frozen ground and that these influence both ice dynamics and landform development.

1. Introduction

The influence of water on basal processes is a critical issue in understanding the dynamics of glaciers and ice sheets. Subglacial conditions influence not only the flow regime of ice masses, but also the formation of glacial landforms. Our knowledge of the details of physical processes at the base of the ice is rather limited. In particular, the existence and interaction of water and permanently frozen ground overridden by a glacial advance is of major importance for further understanding these processes and the past history of landscapes. The response of frozen ground to changes in surface temperature ranges from a few years to thousands of years. If there is a significant amount of unfrozen water in the soil, a much longer response time results, due to latent heat associated with change in phase.

The area of the southern margin of the Laurentide ice sheet (LIS) has been studied and mapped extensively (Colgan and others, 2003; Mickelson and Colgan, 2003). Much of the outer 500 km of the glacierized area is covered with a sediment layer. Based on our present understanding, the conditions necessary for many basal processes have a major impact on the dynamics of an ice sheet. The southern margin of the LIS was dominated by the presence of relatively thin ice lobes that seem to have been very sensitive to external and internal physical conditions, particularly whether the bed was frozen or not. Model results from Marshall and Clark (2002) suggested that 60–80% of the LIS was cold-based at the Last Glacial Maximum (LGM).

Extensive investigations with time-dependent, transient ice-sheet models have been carried out on the LIS. Three-dimensional modeling of growth and decay has been used to reconstruct the areal extent and ice volume for the LIS (Marshall and others, 2000, 2002). The impact of ice-stream dynamics in the Hudson Strait, Canada, has been investigated (Marshall and Clarke, 1997). A coupled ice-sheet model with bed thermodynamics for continental scales was introduced by Tarasov and Peltier (1999). Reconstruction of the ice lobes at the southern margin of the ice sheet has been addressed with two-dimensional flowline models (Cutler and others, 2000, 2001). The effects of permafrost on the advancing ice sheet are especially important (Winguth and others, 2004).

In this paper, we present a modeling approach to investigate the spatial distribution of subglacial thermal conditions around the southern margin of the LIS. We couple the three-dimensional thermomechanical ice-sheet model (Marshall and others, 2002) with a model for the thermal regime of the upper crust. Our numerical experiments are designed to investigate the thermal evolution at the base of the ice sheet at the LGM.

2. Model

A standard ice-sheet model approach is coupled with a model of the thermodynamics of the upper Earth crust. A detailed description of the UBC ice-sheet model used in this study is given in Marshall and others (2002). The model treats thermomechanical ice flow, a parameterization of the paleoclimatic forcing and the geodynamic calculation of isostatic rebound. The approach has been extended in this analysis by coupling with a thermodynamic treatment of a water–soil mixture based on techniques developed by Osterkamp (1987). Model simulations use a finite-difference scheme on a spherical grid over North America. The vertical coordinate z is defined to be positive upwards, with z = 0 representing the present-day sea level. The grid resolution is 1^˚ in the longitudinal and 0.5^˚ in the latitudinal direction. There are 20 evenly spaced vertical layers for the ice-sheet thickness, and 30 unevenly spaced layers for the upper 3000 m of the crust. For each time-step of the model (2–5 years) the dynamics and the mass balance are solved, while the thermodynamics in both ice and bed are updated less frequently (20 years in the presented experiments). We summarize the key factors of the model below.

Ice mechanics

The evolution of ice thickness, ice flow and ice temperature are calculated following the basic conservation equations of mass, momentum and energy together with a constitutive relationship accounting for ice rheology. Given the bed topography and a temporal climatic history, the glaciological model predicts the area, thickness and volume of the ice sheet and the distribution of ice temperature.

The three-dimensional finite-difference model follows the standard formulation for ice dynamics based on Mahaffy (1976) and the thermodynamics based on Jenssen (1977). The shallow-ice approximation (SIA) allows us to ignore longitudinal deviatoric stresses and lateral shearing (Hutter, 1983). These assumptions are valid as long as grid spacing is an order of magnitude larger than ice thickness and sufficiently small bedrock and surface slopes are present (Nye, 1969). Therefore, the effects of these limitations do not have to be considered except at the margins of the ice sheet.

Glen’s flow law, the standard ice rheology in ice-sheet models, is used. This relation for stresses and strain rates in polycrystalline ice has a temperature-dependent ice-stiffness coefficient and links the thermodynamics with the mechanics (Glen, 1955; Paterson, 1994).

Basal motion

Basal motion, sliding and bed deformation at the base of the ice sheet are an important but poorly understood aspect of ice-flow dynamics. The complex nature of subglacial processes is parameterized in standard ice-sheet models by relating basal flow rates to local driving stress (e.g. Payne, 1995). For our modeling experiments, the basal velocity u _b depends on the thermal conditions

(1)

with ρ _i being the density of ice, g acceleration due to gravity, H the local ice-sheet thickness, b _s the height of the ice surface, x _i the horizontal coordinate and B an arbitrary variable. The parameter α _T acts as a binary thermal enabling function

(2)

where T _ib is the temperature at the ice/bed interface and T _pmp is the pressure-melting temperature.

Mass balance

Mass balance is calculated from net accumulation and net ablation of ice. The ablation includes mass loss from the surface, basal and internal melt and calving of icebergs at marine margins. Snow and ice ablation at the surface are parameterized using a monthly degree-day methodology (Braithwaite, 1985). Ice accumulation is due to solid precipitation (snow) at the surface. Monthly precipitation and temperature estimates are used to calculate the fraction of total precipitation that is snow. Net annual ice accumulation is then evaluated by summing monthly snow accumulation and converted to ice equivalent. A detailed description is given in Marshall and others (2000).

While the reconstruction of the LGM ice sheet is not sensitive to the treatment of ablation terms, the calculation of snow and ice accumulation was found to be a much more sensitive matter (Marshall and others, 2002). The forcing of the model with paleoclimatic perturbations is prescribed using a glaciation index which describes the variation between interglacial (modern) and full glacial (LGM) conditions. The temperature and precipitation fields are varied based on this index. The index is based on the Greenland Ice Core Project (GRIP) δ¹⁸O record (Dansgaard and others, 1993) and atmospheric circulation model reconstructions for present-day and LGM distributions of precipitation and temperature (Marshall and others, 2000).

Thermal regime in the bed

Following Osterkamp (1987), we model the thermal evolution of a rock–water–ice–air mixture material with the one-dimensional diffusion equation

(3)

where t is time, T _b is the bed temperature, k is the thermal conductivity and C the apparent volumetric heat capacity equal to

(4)

where C _v is volumetric heat capacity, L is the volumetric latent heat of fusion for ice and θ _w is the volumetric unfrozen water content. This relation accounts for existing unfrozen water in the bed material at temperatures below freezing, which affects the heat capacity due to latent-heat effects. The second term on the righthand side of Equation (4) causes C to increase by three orders of magnitude close to the freezing point (Williams and Smith, 1989). Based on an empirical relation (Anderson and others, 1973), θ _w can be represented by

(5)

where θ _t is the porosity of the thawed material, ρ _b is the dry bulk density of the partially frozen soil/rock material, ρ _w is the density of water, T ^* is the temperature in degrees Celsius and c and d are empirically derived material-dependent constants. Cutler and others (2000) derived values for the two parameters by fitting to the empirical curves (see Table 1). Finally the bulk values for k and C _v are calculated using

(6)

(7)

where b, w, and i refer to the material components soil/rock, water and ice, respectively, and θ is the volumetric proportion. In addition, full saturation is assumed θ _b = 1 – θ _t and θ _t = θ _w + θ _i. Osterkamp (1987) summarizes the necessary relationships to calculate k and C _v.

Table 1. Model parameters

At the lower boundary the spatial distribution of the geothermal heat flux Q _G is applied at a constant rate.

3. Model Experiments

The evolution of the ice sheet was modeled through the last glacial cycle. Primary attention was focused on correct reconstruction of the southern margin of the LIS.

Our experiments were designed to address the following questions about thermal conditions:

What was the distribution of frozen and temperate basal zones at the LGM?

How thick was the frozen layer, i.e. permafrost, around the ice-sheet margin?

How rapid was the degradation of permafrost under the ice sheet?

Were there broad areas of certain thermal conditions that relate to the present distribution of landforms?

Our modeling experiments do not address the complexity of basal flow dynamics. We do not simulate any basal hydrological evolution nor sediment deformation. It has to be noted that the model does not include the physical mechanism responsible for fast flow generally termed as ice streaming.

Only a small selection of the whole set of input parameters to the model that allow us to answer the above questions have been varied. Parameters affecting basal ice movement, thermal behavior of the crust, and climate have been modified. The tuning criteria were the areal extent and timing of the LGM advance. Values of important parameters are presented in Table 1.

Two different reconstructions were modeled. The first was a reference run that did not treat thermal evolution in the crust (referred to as experiment A). Following that, the main emphasis was on the reconstructions with the improved model (B1). In order to investigate the approach by Osterkamp (1987), in a third experiment (B2) the thermal evolution was reduced to simple heat diffusion without accounting for thermal effect due to mixture-material properties of the ground (Equation (3): C and k constant).

The initial configuration at the beginning of the last glacial cycle is the present-day topography and no ice. The geothermal heat-flux distribution is kept constant during the whole cycle. In the simulation without thermal evolution in the crust, the geothermal heat flux is directly applied to the ice/bed interface. In our simulations of only the North American ice sheet, a global estimate of sea level is not achievable. Therefore, the history of sea-level changes is introduced in order to determine marine margin conditions. The climatic fields are updated every 100 years and the isostatic adjustment is considered in 1000 year intervals.

4. Results

Figure 1 shows the ice sheet resulting from the reference run A at the LGM (21 kyr BP). Ice thickness, surface flow and thermal conditions in the area of the southern LIS are shown. A continuous belt along the marginal area is warm-based. The temperate basal conditions provide enhanced flow, clearly pronounced in the Saint Lawrence River. Most marginal surface velocities are less than 200ma^–1. Some exceed this value but are less than 500ma^–1, a value observed in fast-flowing areas on the present-day ice sheets.

Fig. 1. Ice-sheet surface (contours), surface flow field (vectors) and temperate bed conditions (hatched) at the LGM (21 kyr BP) for the reference run A without thermal evolution in the crust. The dashed line represents the maximal extension based on field mapping (Dyke and Prest, 1987), and the two filled circles mark the positions of two marginal locations used for the discussion of the time-transient analysis.

Figure 2shows the same area for the extended model B1 with the coupled thermal evolution model of the crust. In contrast to Figure 1, no continuous zone of temperate conditions at the ice margin develops. The absence of a wet-bed zone along the margin area results in thicker ice and slower flow along the margin than if the bed were thawed. High surface velocities only occur in the steep outer margin. One noticeable and distinct feature is the warm-based area in the Lake Michigan basin.

Fig. 2. Ice-sheet surface (contours), surface flow field (vectors) and temperate bed conditions (hatched) at the LGM (21 kyr BP) for run B1 with thermal evolution in the crust.

At the LGM, the reduced-thermal-evolution run B2 results in continuous frozen ground, similar to B1 (shown in Fig. 2). The difference in area of the frozen base around the southern margin of the LIS is <1%. Surface topography is lower between 50 and 300 m, with the same general shape of a steep margin. This results in a similar spatial distribution of the surface velocities.

Figures 3 and 4 show the time-transient evolution of ice thickness, surface and basal temperatures, and permafrost layer thickness for two marginal locations: the Green Bay area (Fig. 3) is compared with a location farther north (Fig. 4). Basal temperatures are corrected for pressure-melting-point reduction. The LGM state (represented by maximal ice thickness) was reached between 22 and 19 kyr BP. Only minor variations of ice-margin position occurred during this period. Rapid melt of >1000m of ice within 100–300 years took place during the final retreat, when the temperature at the base reached the melting point.

Fig. 3. Transient ice-sheet evolution of B1 at 45.258 N, 87.58 W: (a) ice thickness; (b) surface (dotted curve) and bed temperature; and (c) permafrost thickness.

Fig. 4. Transient ice-sheet evolution of B1 at 50.75° N, 82.5° W: (a) ice thickness; (b) surface (dotted curve) and bed temperature; and (c) permafrost thickness.

The model’s ability to adjust to small changes in climate results in rapid development and degradation of permafrost during ice-free phases (Fig. 3c). Stable conditions occur farther north under the ice. In the central area of the ice sheet a thick permafrost layer up to 200m thick develops (Fig. 4c). The degradation of the layer during the period with thick ice cover is slow. If the base is temperate a faster degradation is clearly visible.

Figure 5 shows the comparison of the time-transient evolution of model experiments B1 and B2. The same location in the Green Bay area as in Figure 3 was chosen, because the differences are better represented than in the more centered location. While the ice sheet advances quicker in B2, B1 develops a thicker ice mass. There is no difference in the timing of the LGM. During most of the period of the last advance (–30 to –15 kyr) the base of the ice is below melting temperature. B2 is closer to thawing and starts to melt earlier. As a consequence, the frozen layer (Fig. 5c)

grows quicker and the degradation is more pronounced.

Fig. 5. Transient ice-sheet evolution at 45.25° N, 87.5˚W: comparison of ice thickness (a), bed temperature (b) and permafrost thickness (c) between the two thermal coupled model experiments B1 (black curves) and B2 (gray curves).

5. Discussion and Conclusions

As part of an ongoing investigation of physical conditions under the southern LIS, we have developed a numerical model that simulates time-dependent thermal interactions between an advancing and retreating ice sheet and the sub-and proglacial land surface. Based on our results, we suspect that the ice–permafrost interactions significantly influence ice dynamics and landform formation.

Field evidence (Mickelson and others, 1983; Pewe, 1983) that continuous permafrost was present in front of and under the ice sheet during most of the cycle except along the southernmost part of the ice sheet agrees with model results. Despite the fact that an ice-sheet model is sensitive to the poorly constrained paleoclimatic input of air temperatures and precipitation rates, the modeled permafrost distribution at the LGM is insensitive to variations of initial conditions, as well as the spatial and temporal variation of the forcing parameters. Marshall and others (2002) conclude that ice-sheet models tend to over-predict ice volumes at the LGM. The vertical downward advection of cold ice may have too great an impact. A slow degradation of permafrost even under the thick ice sheet is the result. Degradation of the frozen layer proceeds from the Earth interior toward the surface.

The thick marginal area of the ice sheet predicted by the model at the southern terminus of the ice sheet does not agree with field evidence (Mickelson and Colgan, 2003) of thin, sensitively reacting ice lobes. For the resulting conditions, the model does not allow the initiation of fast-flow instabilities that are suggested by field evidence. The thermal bed condition found in the Lake Michigan basin may be an indication of conditions needed for the initiation of such flow instabilities. Due to the scale of the model and its inability to handle flow instabilities, there is little direct indication of basal conditions that had an impact on lobe initiation. We interpret this as evidence that the lobes were caused by dynamical flow instabilities.

The modeled permafrost thickness adjusts rapidly to climate conditions in the proglacial ice-free area. A fast transition occurs during the collapse of the ice sheet at the end of the glacial cycle. The transition period with thin ice is of major importance.

We believe our results indicate some strong arguments which support a wide marginal frozen area and still allow the initiation of fast-flow instabilities. The frozen ice margin prohibits subglacial drainage and enables the storage of water. Larger subglacial water reservoirs may have formed and a sudden release may be responsible for flow instability. Cold-based ice does not generally prevent the existence of sliding, as observations in the Arctic have shown (Copland and others, 2003). The presence of water at the ice/bed interface may result from strain heating or meltwater draining from the surface. We believe a further improvement of the model should include a detailed treatment for drainage and storage of water available at the surface and at the base of the ice sheet. Water fluxes may account for a rapid energy exchange. The influence of the large amount of meltwater produced during the decay phases of an ice sheet is restricted to a relatively short period at the end of the glacial cycle and is therefore of minor interest to analysis of lobe initiation and landform formation.

One main aim of this study was to determine the broad spatial distribution of subglacial conditions. The modeled presence of continuous permafrost supports general findings of two-dimensional flowline investigations in lobe areas (Cutler and others, 2000; Winguth and others, 2004).