A numerical investigation of ice-lobe–permafrost interaction around the southern Laurentide ice sheet

Abstract

Permafrost existed around and under marginal parts of the southern Laurentide ice sheet during the Last Glacial Maximum. The presence of permafrost was important in determining the extent, form and dynamics of ice lobes and the landforms they produced because of influences on resistance to basal motion and subglacial hydrology. We develop a two-dimensional time-dependent model of permafrost and glacier-ice dynamics along a flowline to examine: (i) the extent to which permafrost survives under an advancing ice lobe and how it influences landform development and hydrology, and (ii) the influence of permafrost on ice motion and surface profile. The model is applied to the Green Bay lobe, which terminated near Madison, Wisconsin, during the Last Glacial Maximum. Simulations of ice advance over permafrost indicate that the bed upstream of the ice-sheet margin was frozen for 60–200 km at the glacial maximum. Permafrost remained for centuries to a few thousand years under advancing ice, and penetrated sufficiently deep (tens of meters) into the underlying aquifer that drainage of basal meltwater became inefficient, likely resulting in water storage beneath the glacier. Our results highlight the influence of permafrost on subglacial conditions, even though uncertainties in boundary conditions such as climate exist.

1. Introduction

This paper reports on our investigation of the influence of permafrost on subglacial conditions under the southern margin of the Laurentide ice sheet (LIS) during the period leading up to the Last Glacial Maximum (LGM). We introduce a coupled numerical ice-flow–permafrost model to assess the time-dependent feedback between an ice sheet and any underlying permafrost, and perform a suite of sensitivity tests to explore the range of possible scenarios for interactions between the ice sheet and permafrost. We focus on subglacial conditions during ice advance, as it is under such conditions that subglacial permafrost is likely to be most extensive (Mooers, 1990), most dynamic and most influential on landform development and ice motion.

Permafrost existed in Wisconsin, Illinois, Minnesota and North Dakota around the time of the LGM (Moran and others, 1980; Attig and others, 1989; Johnson, 1990; Mooers, 1990; Colgan, 1996). Differences in the nature of glacial deposits and former ice dynamics around the southern margin of the LIS have been attributed to spatial variations in the degree of permafrost development (Mickelson and others, 1983). In the southern Great Lakes region, dates on spruce wood are numerous. In the northern Great Lakes region, the paucity of radiocarbon-dateable material (especially wood) with ages from 30 to 15.5 kyr before present (BP) suggests that tundra conditions prevailed, and supports the idea that ice lobes overrode permafrost as they advanced to their maximum extent at 22–18 kyr BP (Attig and others, 1989). Frozen-bed conditions may have also been influential in the genesis of features including tunnel valleys/channels (Wright, 1973; Attig and others, 1989), drumlins (Mickelson and others, 1983; Patterson and Hooke, 1995), ice-walled lake plains (Ham and Attig, 1996) and ice-thrust features (Moran and others, 1980; Mooers, 1990). These landforms are noticeably less common, or absent, in the southern Great Lakes region: perhaps the result of limited influence of permafrost on ice dynamics in that area.

Time-dependent numerical ice-sheet models can be used to investigate influences on ice dynamics and place constraints on paleoglaciological conditions during the formation of glacial landforms. Any model aimed at assessing conditions in the northern Great Lakes region of the LIS should include treatment of permafrost, as the temperature distribution in the ground must have influenced ice-lobe dynamics through (i) prevention of basal sliding at the termini of some ice lobes due to frozen-bed conditions, (ii) alteration of the basal heat budget by release or uptake of latent heat of fusion during phase changes, and (iii) alteration of subglacial hydraulic conditions. Steady-state analyses of Mooers (1990) and Colgan (1996) suggest that the frozen-bed zone upstream from an ice-sheet margin may have been restricted in horizontal extent to the first few kilometers at most. However, the extent of this zone is likely to vary under transient conditions. Wright (1973) hypothesized the existence of a broad marginal frozen zone under the Superior lobe of the LIS, and Boulton and others (1995) showed that a frozen zone approximately 100 km wide (throughout this paper “wide” refers to the horizontal extent along the flowline) probably developed upstream of the southern margin of the Scandinavian ice sheet at the LGM. The presence of a broad marginal frozen zone is potentially important for landform development (Mooers, 1990; Shoemaker, 1994; Boulton and Caban, 1995) through its impact on subglacial water flow.

Most models of ice-sheet dynamics utilize the assumption of drained-bed conditions, thus removing the need to treat latent heat (e.g. Payne, 1995). A more likely bed condition is complete saturation of the groundwater aquifer (Pio-trowski, 1997). The latent heat associated with phase changes in such an aquifer influences the subglacial heat balance and may be important to ice dynamics on sub-millennial time-scales. Greve and MacAyeal (1996) note that the latent heat released by freezing of a 25 m thick layer of water-saturated substrate was sufficient to stop thermally regulated cycles in their model of ice dynamics through Hudson Strait. Such cycles have been proposed as a likely mechanism for generation of Heinrich layers in the North Atlantic (MacAyeal, 1993), with possible feedbacks to ice-lobe fluctuations on the southern margin (Mooers and Lehr, 1997). Clearly, the time-dependent interaction between permafrost and ice dynamics needs to be further investigated.

2. Geologic setting

We are interested in conditions under the southern LIS in the Great Lakes region (Fig. 1). In this paper we explore possible ice–permafrost interactions under the Green Bay lobe. Ice in this lobe probably originated from a divide near the western shore of James Bay (Dyke and Prest, 1987). From here flowlines were directed south into the Great Lakes region. For the initial ∼400 km, flow was over Paleozoic carbonate rock units (Karrow and Geddes, 1987) that are considered “soft” in terms of resistance to basal motion (Boulton and others, 1985; Fisher and others, 1985). The term “soft” implies an ability to generate tills that may facilitate rapid basal motion. The transition to “hard” Precambrian crystalline Shield lithologies at 400 km (henceforth termed S–H) is smoothed by a train of carbonate till that overlies crystalline rocks for > 100 km south of the carbonate bedrock (Hicock and others, 1989). This till is often sufficiently deep to mask the underlying bedrock topography. For example, it is up to 52 m thick near Geraldton, Ontario (Kristjansson and Thorleifson, 1987). Whilst the existing carbonate till has been linked with post-LGM readvances into the area (Hicock and others, 1989), it is probable that earlier advances also deposited trains of carbonate till on the Shield. Thus, we follow the recommendation of Hicock and others (1989) that this area should be considered soft-bedded. South of the till-covered and exposed Shield lithologies, a second transition in the bed surface, from hard, crystalline rock to soft, sedimentary rock (henceforth called H–S) occurs near the southern margin of Lake Superior (Farrand and others, 1984; Licciardi and others, 1998). These lithologic transitions have hydrologic, thermal and ice-dynamic significance, as they delimit permeable from impermeable substrates with differing thermal conductivity and differing potential for basal motion (Boulton and others, 1985; Fisher and others, 1985; Alley, 1991; Clark and others, 1996; Jenson and others, 1996; Marshall and others, 1996).

Fig. 1. Location of the Green Bay lobe in the southern LIS. D marks the local ice dome; the lines passing through this point are ice divides. Dashed line indicates path of the modeled flowline. C, Chippewa lobe; W, Wisconsin Valley lobe; L, Langlade lobe. The solid grey area represents land beyond the late-Wisconsinan ice limit. Areas with diagonal grey stripes are soft-bedded. Following Hicock and others (1989), a belt of carbonate till on the Shield is considered soft (unbounded diagonal grey stripes immediately north of Lake Superior. Sources: Dyke and Prest (1987) Attig and others (1989); Hicock and others (1989); Johnson and others (1992); Jenson and others (1996); Licciardi and others (1998).

The vertical distribution of lithologies must also be considered when incorporating permafrost into a numerical simulation. Approximations of the depth of the crystalline basement are drawn from Johnson and others (1992) for the Moose River Basin in Ontario and from the USGS (1996) in Wisconsin. For our model, depth to basement in the Moose River Basin is assumed to be 400 m compared to 200 m along the flowline axis in Wisconsin. At this point we lump all Paleozoic sedimentary lithologies into a single category. We consider unconsolidated Quaternary deposits as a separate unit with a thickness of 10 m. This allows distinctions to be made in the porosity of lithified and unlithified material. It is assumed that the present-day distribution of Quaternary material is representative of pre-LGM conditions.

3. Outline of model

General features

Here we summarize the primary features of the model. The reader is referred to MacAyeal (1997) and Parizek (2000) for further details. The model is a two-dimensional (vertical and horizontal coordinates) flowline model with a domain that extends to a specifiable depth into bedrock. An elastic lithosphere and relaxation asthenosphere (Le Meur and Huybrechts, 1996) are used to determine isostatic adjustment. The model is time-dependent, and is constructed using high-order Hermite-polynomial basis functions within a finite-element framework. Node spacing in the horizontal direction, Δx, is 17 km, and it is assumed that the flowline is fixed in space with respect to time. Each vertical column within the ice is composed of 50 evenly spaced nodes.

Ice flow in the model is fully thermomechanically coupled, and includes a first-order treatment of polythermal ice conditions where ice temperature is allowed to reach, but not exceed, the local pressure-melting point. A constant geothermal flux is prescribed at the lower boundary of the model, and thermal conduction through the bedrock is calculated through time. Climate (precipitation and temperature) acts on the upper boundary. Temperatures within the ice are computed with vertical conduction, strain heating (applied at source rather than solely at the bed), vertical and horizontal advection and sliding friction. The thermal profile in a column of ice that has advanced into a previously barren nodal location 17 km downstream assumes the temperature profile of its upstream neighbor. In the unlikely case where ice advances over several nodes in one time-step (requiring an ice-front velocity of at least 1360 m a⁻¹), the temperature profiles in newly ice-covered nodes that were not located at the previous terminus are set to the local surface temperature.

Model performance has been tested against the European Ice Sheet Modelling Initiative (EISMINT) level 1, 2 and 3 benchmarks (Huybrechts and others, 1996; unpublished information from C. Ritz). Of importance to the findings presented here, the flowline simulations produced both dynamic and thermodynamic results that are consistent with ice-sheet model predictions for trends in the Greenland ice-sheet evolution over the past glacial cycle (Dahl-Jensen, 1993; Johnsen and others, 1995; Ritz and others, 1997). The most obvious deviation from three-dimensional simulations, as expected, occurred in the dynamic profile. The flowline produced a thicker-than-observed present-day ice sheet along the Greenland Ice Core Project (GRIP) transect, as the model does not treat horizontal divergence. This overestimation decreases significantly when simulations are conducted along true ice-sheet “flowlines”.

Basal motion

The term “basal motion” includes both sliding and bed deformation. For model runs in which basal motion, u, is nonzero we use

(1)

(Payne, 1995; Greve and MacAyeal, 1996), where B is an adjustable parameter, ρ _i is the density of ice, H is the ice thickness, z _b is the bedrock elevation, T is temperature and T _pmp is the pressure-melting temperature (Table 1 lists definitions of all notation used herein). We utilize the simple viscous relationship expressed in Equation (1) because the proper partitioning of motion between basal sliding (Iverson and others, 1995; Engelhardt and Kamb, 1998) and till deformation (Jenson and others, 1996; Iverson and others, 1997) is not understood. Hindmarsh (1997) has argued that a viscous relation such as Equation (1) is applicable on length scales rather greater than a few meters even though till possesses a plastic rheology at smaller length scales (Kamb, 1991; Iverson and others, 1998). Considering the uncertainties surrounding components of basal motion, and the simple nature of experiments herein, we feel that Equation (1) is adequate for our purposes.

Table 1. Notation, parameter values and units

Thermal enabling of basal motion is internally solved by the model (cf. Payne, 1995; Greve and MacAyeal, 1996; Marshall and Clarke, 1997a): once the bed becomes thawed at two or more neighboring nodes, basal motion may occur in that region if B ≠ 0. The value of B is governed by the nature of the substrate. In our initial sensitivity experiments, B is assigned a value of zero everywhere, to maintain simplicity. Later, we retain B = 0 for hard-bedded areas (Jenson and others, 1996), and assume values of 2 × 10⁻³ to 5 × 10⁻³ m a⁻¹ Pa⁻¹ in soft-bedded areas (Payne, 1995; Greve and MacAyeal, 1996). The latter value allows maximum sliding speeds of 200–300 m a⁻¹. By retaining B = 0 in hard-bedded areas, we assume that sliding speeds are unimportant for overall ice–permafrost interaction compared to speeds over soft-bedded areas. For example, making B finite, but small, in hard-bedded areas (one-tenth the value of B in soft-bedded areas) does not alter conclusions on ice–permafrost interaction: computed frozen-bed widths, maximum permafrost depths and ice thickness are changed by <10%.

Spatial variation of climate and mass balance

Mean annual temperature, T _a, and precipitation, P, are the variables in the climate parameterization. Our desire to explore the fundamental interactions between ice dynamics and permafrost (and a lack of detailed proxy evidence) leads us to assume constant values for P at a given latitude, φ, and elevation, H + z _b, for all experiments described herein. Modern P declines northwards (Walter and Leith, 1967; Korzun, 1977), though trends in the past are unknown (Jenson and others, 1996; Marshall and Clarke, 1997b). Thus, we adopt the following relation:

(2)

where P _sl is the annual precipitation at sea level in m a⁻¹, Π _z is an elevation correction, z _s is surface elevation and Π _φ is a latitude correction (Table 1). Evidence from permafrost features and fossil ostracodes in Illinois (Johnson, 1990; Curry, 1998) and speleothem growth rates in Missouri (personal communication from J. Dorale, 1999) suggests that P was less than today throughout marine oxygen-isotope stages 3 and 2. These data support results of paleoclimate simulations by Kutzbach and others (1998) showing that P in the Midwestern sector of the LIS was 20–30% lower in simulations for 21 and 16 kyr BP than modern P. We therefore assume a value of 0.6 m a⁻¹, or 75% of modern P at Madison, Wisconsin, as our baseline value around which sensitivity tests are performed.

Mean annual air temperature, T _a, varies in space as:

(3)

where T _b is mean annual temperature at a fixed location, Γ _z is the adiabatic lapse rate and Γ _φ is the latitudinal lapse rate (Table 1). We performed runs with either fixed or time-dependent values of T _a. In the time-dependent runs the value of T _b is prescribed through time, as discussed in the next section.

Rather than assume a fixed net accumulation pattern that operates independently of the specified temperature regime (Hooke, 1977; Mooers, 1990; Payne, 1995; Jenson and others, 1996), we use the degree-day method to determine net annual accumulation (Huybrechts and others, 1991; MacAyeal, 1997). An example of the simulated mass-budget curve for the Green Bay lobe flowline at the end of a run with steady air temperature is provided in Figure 2. In the absence of modern analogs with which to compare our simulated mass balance (cf. Jenson and others, 1996), the curve in Figure 2 can be compared with a measured profile from Antarctica and a simulated profile from Greenland (Boulton and others, 1984). Our profile mimics that for Greenland with regard to the steepening balance gradient near the equilibrium-line altitude. This arises primarily because we link mean annual accumulation to T _a. The negative mass balance at lower elevations also seems intuitively reasonable (cf. Hughes, 1997), as ice was advancing into increasingly warm surroundings in the Great Lakes region.

Fig. 2. Distribution of mass balance with elevation in Greenland and Antarctica (Boulton and others, 1984) compared with simulated distributions for the Lake Michigan lobe (Jenson and others, 1996) and the Green Bay lobe (this paper). The Green Bay lobe mass-balance curve is a snapshot of conditions after 30 kyr in a run with steady air temperature. Adapted from Jenson and others (1996).

Time-dependent temperature boundary conditions and initial conditions

In choosing temperature boundary conditions and initial conditions we strike a compromise between the complexity of the global climate system over the last glacial cycle (Bond and others, 1993; Dorale and others, 1998) and our goal to understand the fundamental interactions between ice and permafrost under simple climatic scenarios. Our investigation focuses on ice advance culminating in LGM conditions during marine oxygen-isotope stage 2 because that is the period for which we have convincing evidence of ice–per-mafrost interactions. Thus, the following constraints are considered: (i) favorable conditions for ice advance within stage 3 began at 55 kyr BP (Dorale and others, 1998) and culminated at ∼21 kyr BP (Pollard and Thompson, 1997); (ii) within the last glacial period (Clark and others, 1993), ice did not advance into the upper Mississippi valley until stage 3 (Johnson and Follmer, 1989; Curry and Follmer, 1992; Curry and Pavich, 1996; Dorale and others, 1998); (iii) the likely maximum southern extent of the south-central LIS at 55 kyr BP was north of Lake Superior (Clark and others, 1993), whilst the likely minimum southern extent at 55 kyr BP was at the north end of our flowline in the James Bay lowland; (iv) the Green Bay lobe reached Madison, Wisconsin, at the LGM, depositing the Johnstown moraine (Attig and others, 1989); and (v) temperatures must have been sufficiently warm at 55 kyr BP to prevent glaciation south of the Hudson Bay lowland, and sufficiently cool at the LGM for permafrost to evolve at latitudes where field evidence for its former presence exists.

We simplify the temperature boundary condition for our simulations in the following manner. Two sets of runs are performed. In the first, ice advance over permafrost is examined under steady temperature. This permits investigation of basic information such as time-scales for thawing of subglacial permafrost and possible magnitudes of frozen-zone width in response to changes in variables such as precipitation rate, geothermal flux and substrate porosity. The second temperature scenario introduces an idealized version of possible temperature trends through stages 3 and 2, based on trends in reported climate proxies. Using this temperature scenario as input to the model, as well as the above-mentioned constraints on ice extent through time, we draw conclusions on the evolution of ice–permafrost interactions in the Great Lakes region.

The idealized time-dependent temperature scenario for our second group of runs is summarized in Figure 3. This scenario is constructed from a range of sources. The principal source is a U/Th-dated record of δ ¹⁸O in speleothems from Missouri (Dorale and others, 1998) that is assumed to indicate relative air-temperature variations between 55 and 27 kyr BP. This record originates from only 200 km south of the LGM ice margin and is accurately dated throughout. Temperature dropped by ∼4°C in Missouri between 55 and 41 kyr BP (Dorale and others, 1998), and subsequently rose ∼2°C by ∼35 kyr BP before roughly stabilizing until the end of the record.

Fig. 3. Idealized temperature variation at Madison, Wisconsin, 55–21 kyr BP (see text for details).

It is possible that changes in source δ ¹⁸O related to ice-sheet expansion after ∼31 kyr BP masked the effect of declining air temperature on local δ ¹⁸O values in Missouri (personal communication from J. Dorale, 1999). This question of timing of the final decline in temperatures to the LGM is one of three that must be addressed when incorporating the Missouri record into a complete reconstruction of the temperature history for the period 55–21 kyr BP. The other two issues are: (i) the magnitude of temperature change in this final decline, and (ii) the initial temperature at 55 kyr BP. We address these three questions in order. At the recent end of the time range it is clear that temperatures dropped for at least 5 kyr prior to the LGM (Pollard and Thompson, 1997). It is also evident that sea level was beginning to drop distinctly by ∼31 kyr BP (Chappell and others, 1996), and temperatures in the North Atlantic and Antarctica were decreasing by 31 kyr BP (Ruddiman and others, 1986; Jouzel and others, 1987). We therefore test the impact on results of assuming that the final decline in temperatures to LGM conditions endured for 7–10 kyr. Related to the magnitude of the final temperature drop, evidence of ice-wedge polygons around the margin of the LIS in Wisconsin suggests LGM temperatures of −5°C or lower (Washburn, 1980). Total cooling in the Vostok core between 31 and 21 kyr BP is slightly greater than between 55 and 41 kyr BP (Jouzel and others, 1987); thus, by analogy, cooling in the Great Lakes region may have been approximately 5°C (when compared to the ∼4°C drop in Missouri between 55 and 41 kyr BP). Using this information to work backwards in time from the LGM yields a mean annual temperature of ∼2°C at Madison at 55 kyr BP. The overall reconstruction yields a drop of 7°C between 55 and 21 kyr BP. The total glacial–interglacial change is then at least 13°C, using the modern-day temperature at Madison of 8°C. This concurs with differences of 12–14°C in the same area from global circulation model simulations (Kutzbach and others, 1998).

Treatment of permafrost

Permafrost growth and decay is modeled by numerically solving the heat-transport equation for mixtures of parent material, water and ice. We use relationships outlined by Osterkamp (1987) for the time-dependent simulation of permafrost dynamics. In one dimension, the heat-transport equation is

(4)

where z is the vertical coordinate, K is the thermal conductivity, t is time and C is the apparent volumetric heat capacity equal to

(5)

where C _v is the volumetric heat capacity, L is the volumetric latent heat of fusion for ice and θ _u is the volumetric unfrozen water content. The second term on the righthand side of Equation (5) causes C to increase by ∼3 orders of magnitude close to the freezing point (Williams and Smith, 1989). Anderson and others (1973) expressed θ _u in the readily differentiable form

(6)

where ρ _b is the dry bulk density of the parent material, ρ _u is the density of water, T _* is the temperature in degrees Celsius below freezing, θ _t is the thawed porosity of the material, a and b are empirically derived material-dependent constants, and ΔT is a small offset that ensures a smooth transition between the two regimes. Values of these parameters for unlithified material were derived from data of Williams and Smith (1989) using a least-squares fit. We obtained a = 0.011 m³ H₂O m⁻³ dry sediment °C ^b , and b = −0.660, respectively, for “sandy loam”. For lithified sedimentary rocks we initially use a = 0.015 m³ H₂O m⁻³ °C ^b and b = −0.415, with a porosity of 0.2. Model sensitivity to lower porosity is tested later on. Strictly, Equation (6) applies with no overburden pressure (Anderson and others, 1973). However, we assume that the relationship holds in the subsurface if T _* is substituted by the homologous temperature, or temperature below the pressure-melting point (also expressed in degrees Celsius below freezing), T _H*, in Equation (6).

There are 75 subsurface nodes in each column of the flowline model. It is assumed that each subsurface element has homogeneous values of K, C, ρ _b, ρ _u and θ _u for a given T, though these can vary between geologic units. Nodes are closely spaced near the top of the subsurface. Node spacing, Δz, is 1 m for the upper 10 m, 5 m for the next 200 m, and 20 m for the lowest 500 m. Three subdivisions of the subsurface are delineated: (i) an unlithified upper sediment, (ii) an intermediate lithified sedimentary rock, and (iii) a lower impermeable crystalline rock. The horizontal boundaries between geological materials are not necessarily coincident with changes in Δz, as the thickness of each unit may vary in space.

In order to obtain bulk estimates of K in Equation (4) and C _v in Equation (5) it is necessary to evaluate the contributions of the individual thermal properties of rock, water, air and ice. By assuming a saturated system at all times, we exclude the contribution from air. The bulk value of K for a mixture is represented using a weighted geometric mean (Lunardini, 1995):

(7)

where subscripts “b”, “u” and “i” refer to parent material, water and ice, respectively, and θ is the volumetric proportion of each component. It is assumed that θ _b = 1 − θ _t, where θ _t = θ _u + θ _i. Osterkamp (1987) summarizes the necessary relationships to calculate K _i and K _u. K _b is assumed to be constant with T, though it can vary with parent material. Values of C _v are computed using

(8)

where the contributions of individual volume heat capacities of ice, water and parent material are weighted by their volumetric contribution (Osterkamp, 1987).

Thawing of permafrost due to heat advection in groundwater

Anomalously high basal heat flow has been detected under some valley glaciers, probably as a result of potential energy released from groundwater (Clarke and others, 1984; Echelmeyer, 1987). Though the magnitude of this energy source is likely reduced under ice sheets (Echelmeyer, 1987), it may become important under thinner ice near the terminus (where permafrost may persist). The model therefore includes a simple treatment of potential energy, following Clarke and others (1984) and Paterson (1994). Energy is released at the rate β _a = Q(−dϕ/ds) per unit length in the direction of water flow, s. Here Q is discharge per unit width and φ is water potential. The gradient in water potential can be expressed as

(9)

where z _s is the elevation of the ice surface (Paterson, 1994). An important assumption behind Equation (9) is that subglacial water pressure is equal to the ice overburden pressure (Paterson, 1994). The presence of a permafrost cap introduces an additional influence on the potential gradient, but this is insufficient to undermine the conclusion that surface slope dominates the determination of the gradient. Calculation of the gradient in the proglacial area accounts for the onset of discontinuous permafrost, as pressurized flow is not sustained beyond this point. Following Boulton and others (1995), we assume a threshold permafrost thickness of 50 m for continuous permafrost in the proglacial area. When it is shallower on average, discontinuous permafrost exists. This threshold reflects the effect of heterogeneity in the Earth’s surface that would have caused accelerated deterioration of permafrost in some areas.

In order to determine β _a , Q(s, t) is computed from the cumulative volume of basal meltwater generated upstream from the node in question (Q _b(j)). When downstream nodes are frozen, the value of (Q _b(j)) may not exceed Q _max(j), the maximum possible discharge through the subglacial aquifer. This is determined from Q _max = (K _h D/ρ _u g)(−dϕ/ds), where K _h is the hydraulic conductivity of the substrate, and D is the depth over which water flows in the aquifer. The available flow depth between the base of the permafrost and the base of the aquifer varies as permafrost evolves. We track this thickness over time and space. If Q _b(j) > Q _max(j), then β _a is determined using Q _max. The hydraulic consequences of this condition will be examined in future work; in this paper we consider only the contribution of groundwater flow to the basal thermal regime.

We implement the heat flux, β _a, as a point source at the base of the ice sheet or under the permafrost, if present. It is assumed that all energy is released at the ice–bed, or permafrost–bed interface. Not all of β _a is available to melt permafrost, however. An amount −Qkρ _i g∂(z _s − z _b)/∂s is used to maintain water at the pressure-melting point as it moves down the potential gradient. We use k = 0.410 (Hooke, 1998), a value for air-saturated water. Typically, β _a peaks at <0.001 W m⁻² in our experiments, because only a small fraction of basal meltwater can escape through the aquifer.

4. Outline of Model Runs

Table 1 contains values used in the model runs. We performed runs to examine the following questions.

1. 1. How deep was permafrost around the margin of the Green Bay lobe?

2. 2. How wide does the frozen zone become during an ice advance?

3. 3. How quickly does permafrost thaw under an advancing ice front?

4. 4. What are the potential interactions between basal meltwater, permafrost and ice dynamics?

The model runs are divided into two groups. In the first we explore basic interactions between ice and permafrost under constant-temperature boundary conditions for a total run time of 30 kyr. The sensitivity of model results to the following variables is tested: (i) mean annual precipitation, P; (ii) mean annual temperature, T _a; (iii) the standard deviation of the daily temperature range, σ; (iv) geothermal flux; and (v) bedrock porosity. Table 2 summarizes details of each run. For this first set of runs the prescribed initial ice extent assumes ice with a parabolic profile associated with a basal shear stress of 20 kPa (Paterson, 1994) that terminates at 550 km along the flowline. In some cases the prescribed climate favors instant snow accumulation beyond this initial ice extent, but initial ice volume is identical in all cases. Initial temperatures in the bedrock and ice are in equilibrium with the geothermal flux, and all runs are conducted with a time-step, Δt, of 25 years.

Table 2. Parameter values in steady-climate runs

The second group of experiments (Table 3) explores the evolution of ice–permafrost interactions under variable temperature boundary conditions (Fig. 3). In these experiments we consider the constraints on variations in ice extent between 55 and 21 kyr BP, as outlined earlier.

Table 3. Parameter values in variable-climate runs

5. Results

I. Experiments with steady air temperature

The first set of results is derived from runs with a constant-air-temperature boundary condition and zero sliding (Table 2) over a period of 30 kyr. Note that time is in model years, not years BP. The latter will be used later, in the second set of runs with time-variable air temperature. From our initial tests we develop a feeling for the sensitivity of permafrost depth and frozen-zone width to individual climatic and lithologic influences. Figure 4 orients the reader to a typical situation in which ice is advancing, by internal deformation, over proglacial permafrost. T _a is steady at 10°C below the modern value (T _a at the latitude and elevation of Madison is −2°C), and precipitation is 0.6 m a⁻¹, or 75% of the modern value. Under these conditions, the terminus has migrated to 1180 km from its initial position at 550 km over a period of 11.5 kyr. Permafrost has formed in the unglaciated forefield, is deepest immediately beyond the ice margin (∼120 m thick), and thins under the advancing ice lobe because of insulation by the overlying ice from the colder atmospheric conditions. The width of the marginal frozen zone is about 100 km in Figure 4, though it varies through time. For example, as ice advances further south it encounters shallower permafrost that takes less time to thaw. Accordingly, the width of the frozen zone decreases with further advance under a steady climate.

Fig. 4. Cross-section of an ice lobe advancing over permafrost. North is to the left. Isotherms indicate ice temperature. The glacier bed is shaded grey. Permafrost within the bed is shaded black. Small black squares indicate nodes that are at the pressure-melting point. The extent of the frozen zone under the ice is measured from the terminus to the point where the bed is first at the melting point. Small circles indicate nodes where sliding could take place if B>0 in Equation (1). Madison (“Mad”) is located at 1260 km along the flowline. H–S is the hard-to-soft bed transition, and S–H is the soft-to-hard bed transition

Results from runs 1–3 (Table 2) are shown in Figure 5. The figure illustrates the time variation over a period of 30 kyr of three parameters: (i) ice-marginal position (Fig. 5a), (ii) maximum depth of permafrost (Fig. 5b), and (iii) width of the frozen zone upstream from the terminus (Fig. 5c). Note that (i) the position of deepest permafrost may migrate during an ice advance, and (ii) Figure 5b shows the time variation of maximum permafrost depth only in areas where sedimentary bedrock is present (the freezing front may penetrate more deeply into igneous bedrock, but this is assumed to be of no significance to groundwater flow because of the relative impermeability of igneous bedrock). The stepping character of lines in Figure 5 is caused by the 17 km horizontal spacing of nodes in the model.

Fig. 5. Variation of (a) ice extent, (b) maximum permafrost depth in the southern aquifer (Fig. 4) and (c) frozen-zone width for runs 1–3 (Table 2). Note that maximum permafrost depth (b) may vary in space as ice advances over frozen ground. Note also that frozen-zone width (c) is measured from the terminus to the subglacial location where unfrozen ground first appears. In (a) the thin horizontal line marked “S–H” is the position of the soft-to-hard bed transition, the line marked “H–S” is the position of the hard-to-soft bed transition, and the line marked “Mad” is the position of Madison. Run 1 has a geothermal heat flux of 0.04 W m ⁻², compared to 0.035 W m⁻² in run 2 and 0.045 W m⁻² in run 3.

The three runs in Figure 5 demonstrate model sensitivity to choice of geothermal flux. The geothermal flux was adjusted by ±0.005 W m⁻² (runs 2 and 3) around the “standard” value of 0.040 W m⁻² (run 1). Such a variation is within the uncertainty of field measurements (Touloukian and others, 1981). Clearly, the general trends in Figure 5b and c are unaffected by choice of geothermal flux. The maximum permafrost depth is 130–145 m, and the maximum width of the frozen zone is 100–140 km. Not surprisingly, permafrost reaches greater depths with a weaker geothermal flux (run 2). When ice reaches the latitude of Madison after ∼16 kyr of the model run, the width of the frozen zone is 60–90 km, and maximum permafrost depth is 100–120 m (approximately half of the depth of the aquifer). As suspected by Mooers (1990), and demonstrated by Boulton and others (1995), the width of the frozen zone can be much greater than a few kilometers in the transient situation of an ice advance. Only as the runs approach steady state, after ∼25 kyr does the frozen-zone width permanently drop to less than the spatial resolution of the model. Not shown in Figure 5 is the time taken for permafrost to thaw under the advancing ice. Thawing takes 5 kyr at 1040 km along the flowline, and 4 kyr at 1220 km.

Figure 6 is the first of two figures demonstrating model sensitivity to aspects of the climate parameterization. This figure contains results from runs 1 and 4, in which the standard deviation of daily temperature variation, σ, is 5° and 6°C respectively. A value of ∼5°C has been used in mass-balance calculations for Greenland (Huybrechts and others, 1991; Reeh, 1991), and we suggest that a slightly higher value could be applicable in a mid-continental climate. As seen in Figure 6a, σ impacts the ice-marginal position because warmer peak daily temperatures associated with higher σ cause more melting. Notice that although both runs have the same prescribed initial ice extent, the climate in run 1 favors initial snow accumulation further south than run 4 (Fig. 6a). T _a is the same in both runs; therefore differences in maximum thickness of permafrost (Fig. 6b) are controlled solely by the timing of ice advance. For example, the less hostile climate (with respect to ice growth: σ = 5°C, run 1) favors early ice cover at H–S, which in turn prevents permafrost from penetrating so deeply into the sedimentary aquifer (140 m, compared with 210 m in run 4 in which ice takes an additional 12 kyr to reach H–S). Longer thaw times for deeper permafrost influence the trend in frozen-zone width (Fig. 6c). In contrast to run 1, in which a single peak occurs after 10 kyr of the model run, a broad initial peak at 10 kyr is followed by a more distinct peak at 22 kyr in run 4. In this case, maximum frozen-zone width is 260 km. The explanation for the double peak is as follows: Permafrost thaw rates on the crystalline bedrock are higher than in the sedimentary aquifer because of the additional latent heat that must be released in the latter case before temperature can rise. Consequently, upon ice advance onto the sedimentary rocks, the trend for decreasing frozen-zone width is reversed until the rate of ice advance once again becomes smaller than the now-reduced rate of migration of the unfrozen zone.

Fig. 6. Variation of (a) ice extent, (b) maximum permafrost depth in the southern aquifer, and (c) frozen-zone width, for runs 1 and 4 (Table 2). In run 1, the standard deviation of daily temperature variation, σ, is 5°C, compared to σ = 6°C in run 4.

Figure 7 illustrates the relationship between mean annual precipitation and permafrost dynamics (runs 1, 5 and 6 in Table 2). The ice margin reaches the latitude of Madison after 9 and 16 kyr in runs 6 and 1, respectively. It never reaches that point in run 5, due to insufficient supply of mass. Not shown in Figure 7 is run 7 in which snow immediately accumulates at Madison due to favorably high precipitation and cool temperature. As with runs in Figure 6, the maximum depth of permafrost in Figure 7b is regulated by how quickly ice covers the cooler northern portion of the flowline. The range is 100 m in run 6, with P = 0.7 m a⁻¹, to 200 m in run 5, with P = 0.5 m a⁻¹. In Figure 7c, peak frozen-zone width increases with decreasing P. As with results in Figure 6, the longer the permafrost formation time, the wider the frozen zone.

Fig. 7 Variation of (a) ice extent, (b) maximum permafrost depth in the southern aquifer, and (c) frozen-zone width for runs 1, 5 and 6 (Table 2). The influence of precipitation is tested: P = 0.5 m a⁻¹ in run 5, 0.6 m a⁻¹ in run 1 and 0.7 m a⁻³ in run 6.

The impact on permafrost dynamics of raising or lowering steady-state mean annual temperature by 2°C (runs 8 and 9) is predictable and will be discussed only briefly. Warmer conditions result in slower rates of ice advance, shallower permafrost (maximum 130 m) and, hence, a narrower frozen zone that peaks at 80 km in width. Cooler temperatures promote earlier ice advance, but not before permafrost has penetrated to 190 m. The frozen zone achieves a maximum width of 260 km. One contrast between results from runs 8 and 9 and those relating to precipitation (runs 5–7) is the association of shallower permafrost with slower-moving ice advance in the former, and of deeper permafrost with slower-moving ice in the latter.

Results in Figure 8 show model sensitivity to the porosity of sedimentary rocks that comprise the aquifer (runs 1, 10 and 11). Maximum permafrost depth increases as porosity decreases, due to the reduced volume of water making the phase change to ice (Fig. 8a). The range is 140 m for a porosity of 0.2, to 190 m for zero porosity. For the geologic setting of the Green Bay lobe, the latter almost equals the mean depth of the aquifer along the flowline (USGS, 1996). Despite reaching greater depths with zero porosity, permafrost also thaws more rapidly upon coverage by advancing ice because no latent-heat exchange occurs during vertical migration of the permafrost front. This results in a narrower frozen-zone width than in cases with non-zero porosity (Fig. 8b). The maximum width in this case is 50 km, compared to 100 and 120 km for the porosities of 0.1 and 0.2, respectively.

Fig. 8. Variation of (a) maximum permafrost depth in the southern aquifer, and (b) frozen-zone width, for runs 1, 10 and 11 (Table 2). Sensitivity to aquifer porosity is tested: porosity is 0.2 in run 1, 0.1 in run 10 and 0.0 in run 11.

So far we have evaluated individual influences on permafrost dynamics and seen that model sensitivity to climatic factors is significant. It is possible that effects of some parameters may compound or counteract each other. Nevertheless, order-of-magnitude generalizations about model results for the Green Bay lobe can be made:

1. 1. Permafrost lasts for centuries to a few thousand years under advancing ice.

2. 2. The maximum width of the frozen zone is probably about 100 km.

3. 3. Permafrost almost certainly prevents basal meltwater drainage through unconsolidated sediments that only reach tens of meters in depth, and permafrost may reach the base of the lithified sedimentary aquifer under some circumstances. Maximum permafrost depths are achieved under conditions of slowly advancing ice.

II. Experiments with time-dependent air temperature

The second set of runs is driven by the air-temperature scenario outlined in Figure 3, and therefore runs from 55 to 21 kyr BP. Details of runs are listed in Table 3. Figure 9 contrasts the evolution of lobe profiles in runs 12 (Fig. 9a) and 13 (Fig. 9b). Run 13 includes basal motion north of 550 km and south of 860 km if the bed is thawed. The effect of basal motion on the profile is two-fold. First, the profile flattens in areas where basal motion is possible. The profile is dissimilar from those of Jenson and others (1996) for the Lake Michigan lobe, however, as the frozen margin prevents development of a low profile in that zone. Second, lower surface elevations in areas with basal motion, relative to those in the non-sliding case (run 12), result in increased mass supply from precipitation (Equation (2)). This, in combination with increasing P with decreasing latitude, allows the ice to penetrate further south in the allotted time. At the very least, this result illustrates the limitations of a simple precipitation parameterization (Equation (2) ) in which P varies as a function of elevation. Future refinements of the parameterization could account for orographic lifting, perhaps as a function of surface slope.

Fig. 9. Evolution of the cross-sectional profile (including bed elevation ) of an ice lobe from 55 to 21 kyr BP with (a) zero sliding (run 12) and (b) sliding over unfrozen soft-bedded areas (run 13). In (b), the value of B (Equation (1)) is 2 × 10⁻³ m a⁻¹ Pa⁻¹. Times of the snapshots, working from left to right, are 48, 41, 38, 35, 31, 27, 23 and 21 kyr BP. In both panels, profiles are initially ice-free at 55 kyr BP. In (a) a small snowpatch forms at the terminus between 23 and 21 kyr BP: the result of southern progression of the equilibrium line outpacing ice advance.

An apron of snow at the terminus of the most extensive profile in run 12 (Fig. 9a) is caused by the rate of southerly progression of the equilibrium line temporarily exceeding that of the ice-lobe terminus between 21.5 and 21 kyr BP. This is evident in Figure 10a, which shows ice-marginal position through time for runs 12–14. Also evident are two earlier occurrences of the same behavior, at 53 and 47 kyr BP. In general, the air-temperature scenario in Figure 3 causes ice to reach the latitude of Madison close to the LGM: 27 kyr BP in run 14, 25 kyr BP in run 13 and 21.5 kyr BP in run 12 (Fig. 10a). Noteworthy is the monotonic advance of ice from zero initial ice cover in all three runs, despite a warming trend in air temperature between 41 and 35 kyr BP (Fig. 3). The accelerated advance due to sliding in runs 13 and 14 begins once ice is south of the hard–soft transition, H–S (Fig. 10a). Similar behavior was not possible when ice was north of the soft–hard transition, S–H (Fig. 10a), because the bed was frozen everywhere. This frozen condition can be deduced from an additional dataset in Figure 10a (the thin dotted line), which depicts the total area of the bed that is frozen (the trends in total frozen area from runs 12–14 are very similar, so only that from run 13 is shown). Until 43.5 kyr BP, the entire bed is frozen, whilst most of the bed is thawed after ∼35 kyr BP.

Fig. 10. Variation of (a) ice extent, (b) maximum permafrost depth in the southern aquifer, and (c) frozen-zone width, for runs 12–14 (Table 3). Also shown in (a) is the area of bed that is frozen in run 13 (thin dotted line that descends between 43 and 30 kyr BP). Sensitivity to the sliding constant, B (Equation (1)), is tested. The value of B over unfrozen soft-bedded areas is zero in run 12, 2 × 10⁻³ m a⁻¹ Pa⁻¹ in run 13, and 5 × 10⁻³ m a⁻¹ Pa⁻¹ in run 14.

The variation of maximum permafrost depth south of H–S is shown in Figure 10b. Deepest permafrost (210 m) occurs after 39 kyr BP; ∼2 kyr after the air-temperature minimum in Figure 3. Because of the greater extent of ice, maximum permafrost depth is shallower under cooler conditions at the LGM (60 m deep in run 14 and 160 m in run 12). Until 35 kyr BP, trends in maximum depth are identical between runs because ice has not advanced beyond H–S. After 35 kyr BP, the most rapid advance (run 14) covers the deeper permafrost earlier, thus accelerating demise of the permafrost.

As noted in the discussion of Figure 10a, the lobe is completely frozen at its bed until 43.5 kyr BP; thus, the width of the frozen zone (Fig. 10c) equals the total extent of ice until this time. Thereafter, a thawed patch develops 150 km upstream from the terminus, and until 31 kyr BP the slow rate of ice marginal progression under warmer atmospheric conditions (Fig. 3) leads to a declining frozen-zone width. Once cooling air temperatures allow accelerating ice advance after 31 kyr BP, the frozen-zone width begins to increase again. By the time ice reaches Madison, the frozen-zone width is 80 km in run 14, 100 km in run 13 and 240 km in run 12. The large value of the latter is partly due to proglacial snow-patch formation at 21.5 kyr BP, however. Not shown in Figure 10 is run 15, in which final cooling to LGM conditions is assumed to occur over 7 kyr instead of 10 kyr. As a result of the additional 3 kyr of warmer conditions, permafrost thins to 50 m, and subsequent ice advance results in a narrower frozen-zone width of 60 km instead of 100 km (run 13).

A similarity between runs 12–15 (Fig. 10) is the assumption of zero initial ice. This follows the “minimum” case of Clark and others (1993, p. 105). Figure 11 (displaying runs 16–18) demonstrates the impact on results of assuming an initial ice extent of 275 km (one-sixth of the model domain). Runs 16 (Fig. 11) and 13 (Fig. 10) differ only in their initial conditions. Clearly, choice of initial condition strongly affects final ice-marginal position and, as a result, permafrost dynamics. One “improvement” in the trend of ice-marginal position with time is the lack of periods when the equilibrium-line advance outruns the rate of ice advance (cf. Fig. 10a). However, ice reaches Madison at 27 kyr BP, only 4 kyr after the beginning of cooling temperatures (Fig. 3), and advances to ∼400 km south of Madison by the end of the run. Such a situation is unreasonable, though we cannot discount the possibility of non-zero initial ice extent at 55 kyr BP, or the possibility that lack of treatment of flow divergence towards the margin will allow ice to penetrate further south than if it were considered.

Fig. 11. Variation of (a) ice extent, (b) maximum permafrost depth in the southern aquifer, and (c) frozen-zone width, for runs 16–18 (Table 3). Initial ice extent in all three runs is 275 km, compared with zero for runs in Figure 10. σ = 6°C in run 17, compared to 5°C in run 16, and P = 0.55 m a⁻¹ in run 18, compared to 0.6 m a−¹ in run 16.

Guided by results from the first set of runs (1–11), we tested model response to lowering of P to 0.55 m a⁻¹ (run 17), increasing σ to 6°C (run 18) and raising T _a by 1°C throughout the 34 kyr range (run 19). These adjustments are reasonable within the context of our limited knowledge of climate from 55 to 21 kyr BP. Trends in Figure 11a show that the ice-marginal position in runs 17 and 18 is closer to Madison at the LGM than in run 16. Similarly, the final ice marginal position in run 19 (not shown) equaled that in run 18. Slower southerly migration of the ice in runs 17 and 18 than in run 16 allows permafrost to remain at greater depths from 40 kyr BP onwards (Fig. 11b). In LGM conditions, maximum permafrost depth is 160 m in run 18 and 180 m in run 17. The width of the frozen zone never exceeds 270 km in runs 16–18 (Fig. 11c). On reaching or approaching Madison at the LGM, frozen-zone widths are 180 km in runs 17 and 18.

In the final runs (20 and 21 in Table 3) we tested model sensitivity to the latitudinal gradient in P. Predictably, the steeper gradient resulted in slower ice advance, and conversely. As with previous runs, slower-moving ice allowed deeper permafrost. Frozen-zone widths at the LGM position were 120 km in run 20 and 200 km in run 21.

Rather than toil further with what is at this point a poorly constrained problem with respect to climatic inputs and initial conditions, we now summarize some general concepts revealed by our experiments. First, permafrost probably existed in portions of the lithified sedimentary aquifer along the Green Bay flowline throughout the period 55–21 kyr BP. Second, after initially approaching depths equal to the aquifer thickness, permafrost was always a maximum of a few tens of meters thick up to the LGM. Such depths were well into the lithified sedimentary rocks. Third, the subglacial permafrost, which equates to the marginal frozen zone, varied in extent through time, and was probably 60–200 km wide at the LGM. The robustness of these results in spite of large uncertainties in climate and initial conditions leads us to suspect that untreated effects such as orographic influences on precipitation and divergence along the flowline would not alter our broad conclusions on ice–permafrost interaction.

6. Discussion

Based on the above results, we suspect that ice–permafrost interactions significantly influenced landform development and ice dynamics around the southern margin of the LIS. We now discuss implications of our results for various aspects of geologic interest.

Hydrology

In the event that permafrost developed to depths below the base of unlithified sediments, basal meltwater would have been unable to drain through any form of channelized network (e.g. canals of Walder and Fowler, 1994, or R channels of Röthlisberger, 1972). Instead, water would have been forced to drain through underlying lithified sedimentary rocks. Like Piotrowski (1997) and Fountain and Walder (1998), we find that reasonable aquifer hydraulic conductivities (10⁻⁶ m s⁻¹ in our case) and the limited available vertical thickness of unfrozen aquifer (∼100 m) are too small to allow significant drainage of basal meltwater through the groundwater system. Even if aquifer thickness were an order of magnitude greater under the southern Green Bay lobe, as it was under portions of the southern Scandinavian ice sheet (Boulton and others, 1995), groundwater discharge would remain grossly inadequate. Only by increasing aquifer hydraulic conductivity to unreasonable values of 10⁻³ m s⁻¹, close to that of gravel (Matthess and Ubell, 1983), is it possible to evacuate all basal meltwater through the groundwater system. One result of this, also noted by Piotrowski (1997), is that water could have been temporarily stored subglacially, and was released upon later thawing of permafrost. It is likely that conditions conducive to drainage probably became more prevalent as the ice margin stabilized after the glacial maximum and as the width of the frozen zone declined. Boulton and Caban (1995) have noted the likelihood for hydrofracturing in the vicinity of the ice margin when permafrost is present. Outburst events may have originated through such a process. We see numerous tunnel channels around the LGM margin of the Green Bay lobe that may have resulted from such drainage events (Attig and others, 1989; Clayton and others, 1999).

Drumlin formation

A striking feature of drumlins behind the Johnstown (glacial maximum) moraine of the Green Bay lobe is the systematic decrease in length-to-width ratio from the upstream end of the field from 25:1 to 1.25: 1 (Borowiecka and Erickson, 1985; Attig and others, 1989; Colgan and Mickelson, 1997). This has been interpreted as indicating decreasing basal sliding velocity toward the ice margin (Chorley, 1959). A further idea arises from results herein: following ice advance, the frozen zone began to thin from its upstream end such that sculpting of erosional–depositional landforms became possible in that region first. Because sculpting was possible at the upstream end for the longest time, the landforms in that area became most streamlined. Based on results from our experiments with steady air temperature, time-scales for wastage of the frozen zone are on the order of thousands of years. Radiocarbon dates suggest that the Green Bay lobe may have been close to its maximum extent for a few thousand years (Maher and Mickelson, 1996), thus allowing sufficient time for gradual thawing of the frozen zone. Mickelson and others (1983) and Stanford and Mickelson (1985) suggested that drumlins form in a patchily thawing zone. Certainly sub-kilometre-scale heterogeneities in the landscape would promote uneven thawing of the frozen zone, and this effect could be superimposed on that of a regionally thinning frozen margin. Sparsity of drumlins in the last few kilometers of the flowline (Attig and others, 1989; Colgan and Mickelson, 1997) may be due to retreat of the ice margin before thawed-bed conditions developed in that zone. Smaller drumlin fields created during minor ice readvances after the LGM (Colgan and Mickelson, 1997) may result from readvance over temporarily re-forming permafrost.

Proglacial lakes

One detail explored in a subsequent paper (P. M. Cutler, unpublished information) is the influence of proglacial lakes on ice–permafrost interactions. In large basins, such as that containing Lake Michigan, the presence of a lake would prevent development of permafrost (Mackay and Mathews, 1964; Boulton and others, 1995). A regional-scale warm thermal anomaly would exist under the lake, with lake-bottom temperatures of up to 4°C. Ice advancing into the lake would thus not encounter deep permafrost that would favor development of a steep profile. Instead, rapid sliding of grounded ice might ensue immediately, depending on the hypsometry of the lake bed. Ice-surface profiles in such situations would be gentle, as proposed by Beget (1986), Clark (1992) and Jenson and others (1996) for the Lake Michigan lobe. For reasons relating to basal thermal conditions, the distribution of large proglacial lake basins in the Midwestern U. S.A. may thus have had a significant impact on the heterogeneity of ice dynamics and landforms in the region. For example, landforms possibly relating to thawing frozen zones and inefficient subglacial drainage (drumlins and tunnel channels) would not be expected to form under a lobe advancing into a large lake basin. Indeed, such features are strikingly absent in the area covered by the Lake Michigan lobe (Mickelson and others, 1983).

Ice-surface profile reconstructions

A number of authors have reconstructed ice-surface profiles for glacial maximum extent around the southern margin of the LIS using trends in moraine elevations and flow-direction indications (Wright, 1972; Mathews, 1974; Clark, 1992; Colgan, 1999). Clarks (1992) reconstruction of the Green Bay lobe stood out in comparison to those of surrounding lobes as one with higher basal shear stress (17–22 kPa vs <5 kPa in all other cases). Colgan (1999) confirmed these higher values, but also noted that basal shear stress associated with shallow ice-surface gradients at 50–150 km from the terminus was inconsistent with higher reconstructed values within the last 15 km of the lobe in the vicinity of the Baraboo Hills. This pattern is, however, consistent with a profile impacted by frozen-bed conditions near the terminus (see, e.g., Fig. 9b). An important question for future work relates to the evolution of such a profile as the subglacial permafrost thaws after the LGM. Compared to the initial stable advance over frozen ground, the post-LGM retreating lobe may have had a much lower surface profile and higher sliding speeds because its bed was thawed and the profile was no longer influenced by a broad frozen zone near the margin. Such a sequence might explain the general synchronicity of ice advance to the LGM around the southern LIS, followed by asynchronicity of lobes during deglaciation (personal communication from H. Mooers, 1999).

7. Conclusions

As part of an ongoing investigation of physical conditions under the southern LIS, we have developed a numerical flowline model that simulates time-dependent interactions between an advancing ice lobe and proglacial permafrost. It is applied to the Green Bay lobe which terminated near Madison, Wisconsin. Sensitivity tests demonstrate that our current knowledge of initial conditions and the spatial and temporal variation of air temperature and precipitation limits our ability to draw detailed quantitative conclusions on lobe–permafrost interactions. Nonetheless, the following broad conclusions are considered robust: (i) permafrost was at least tens of meters deep in the forefield of the Green Bay lobe at the LGM; (ii) the frozen-zone width was 60–200 km at the LGM; and (iii) time-scales for thawing of subglacial permafrost are centuries to a few thousand years. Permafrost was sufficiently deep and widespread at the LGM to remove the possibility for steady channelized drainage through unconsolidated sediments under the margin of the Green Bay lobe. Consequently, water may have been stored subglacially until it could be released at a later time, possibly exiting through tunnel channels. Progressive thawing of the frozen zone following the LGM may have allowed erosional–depositional processes to act longest on the upstream end of drumlin fields, perhaps partly explaining the decreasing length-to-width ratio of individual drumlins toward the ice margin. The absence of permafrost in the Lake Michigan basin due to the presence of a lake may explain differences in landforms and ice-surface profiles with the adjacent Green Bay lobe.