General features

Here we summarize the primary features of the model. The reader is referred to Reference MacAyealMacAyeal (1997) and Reference ParizekParizek (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 (Reference Le Meur and HuybrechtsLe 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 (Reference Huybrechts and PayneHuybrechts 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 (Reference Dahl-Jensen and MosegaardDahl-Jensen, 1993; Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others, 1995; Reference Ritz, Fabre and LetréguillyRitz 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

(Reference PaynePayne, 1995; Reference Greve and MacAyealGreve 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 (Reference Iverson, Hanson, Hooke and JanssonIverson and others, 1995; Reference Engelhardt and KambEngelhardt and Kamb, 1998) and till deformation (Reference Jenson, MacAyeal, Clark, Ho and VelaJenson and others, 1996; Reference Iverson, Baker and HooyerIverson and others, 1997) is not understood. Reference HindmarshHindmarsh (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 (Reference KambKamb, 1991; Reference Iverson, Hooyer and BakerIverson 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.

Thermal enabling of basal motion is internally solved by the model (cf. Reference PaynePayne, 1995; Reference Greve and MacAyealGreve and MacAyeal, 1996; Reference Marshall and ClarkeMarshall 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 (Reference Jenson, MacAyeal, Clark, Ho and VelaJenson and others, 1996), and assume values of 2 × 10⁻³ to 5 × 10⁻³ m a⁻¹ Pa⁻¹ in soft-bedded areas (Reference PaynePayne, 1995; Reference Greve and MacAyealGreve 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 (Reference Walter and LeithWalter and Leith, 1967; Reference KorzunKorzun, 1977), though trends in the past are unknown (Reference Jenson, MacAyeal, Clark, Ho and VelaJenson and others, 1996; Reference Marshall and ClarkeMarshall and Clarke, 1997b). Thus, we adopt the following relation:

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 (Reference JohnsonJohnson, 1990; Reference CurryCurry, 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 Reference Kutzbach, Gallimore, Harrison, Behling, Selin and LaarifKutzbach 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:

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 (Reference HookeHooke, 1977; Reference MooersMooers, 1990; Reference PaynePayne, 1995; Reference Jenson, MacAyeal, Clark, Ho and VelaJenson and others, 1996), we use the degree-day method to determine net annual accumulation (Reference Huybrechts, Letréguilly and ReehHuybrechts and others, 1991; Reference MacAyealMacAyeal, 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. Reference Jenson, MacAyeal, Clark, Ho and VelaJenson and others, 1996), the curve in Figure 2 can be compared with a measured profile from Antarctica and a simulated profile from Greenland (Reference Boulton, Smith and MorlandBoulton 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. Reference HughesHughes, 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 (Reference Boulton, Smith and MorlandBoulton and others, 1984) compared with simulated distributions for the Lake Michigan lobe (Reference Jenson, MacAyeal, Clark, Ho and VelaJenson 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 Reference Jenson, MacAyeal, Clark, Ho and VelaJenson 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 (Reference BondBond and others, 1993; Reference Dorale, Edwards, Gonzalez and ItoDorale 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 (Reference Dorale, Edwards, Gonzalez and ItoDorale and others, 1998) and culminated at ∼21 kyr BP (Reference Pollard and ThompsonPollard and Thompson, 1997); (ii) within the last glacial period (Reference ClarkClark and others, 1993), ice did not advance into the upper Mississippi valley until stage 3 (Reference Johnson and FollmerJohnson and Follmer, 1989; Reference Curry, Follmer, Clark and LeaCurry and Follmer, 1992; Reference Curry and PavichCurry and Pavich, 1996; Reference Dorale, Edwards, Gonzalez and ItoDorale and others, 1998); (iii) the likely maximum southern extent of the south-central LIS at 55 kyr BP was north of Lake Superior (Reference ClarkClark 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 (Reference Attig, Mickelson and ClaytonAttig 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 (Reference Dorale, Edwards, Gonzalez and ItoDorale 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 (Reference Dorale, Edwards, Gonzalez and ItoDorale 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 (Reference Pollard and ThompsonPollard and Thompson, 1997). It is also evident that sea level was beginning to drop distinctly by ∼31 kyr BP (Reference ChappellChappell and others, 1996), and temperatures in the North Atlantic and Antarctica were decreasing by 31 kyr BP (Reference Ruddiman, Shackleton, McIntyre, Summerhayes and ShackletonRuddiman and others, 1986; Reference JouzelJouzel 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 (Reference WashburnWashburn, 1980). Total cooling in the Vostok core between 31 and 21 kyr BP is slightly greater than between 55 and 41 kyr BP (Reference JouzelJouzel 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 (Reference Kutzbach, Gallimore, Harrison, Behling, Selin and LaarifKutzbach 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 Reference OsterkampOsterkamp (1987) for the time-dependent simulation of permafrost dynamics. In one dimension, the heat-transport equation is

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

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 (Reference Williams and SmithWilliams and Smith, 1989). Reference Anderson, Tice and McKimAnderson and others (1973) expressed θ _u in the readily differentiable form

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 Reference Williams and SmithWilliams 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 (Reference Anderson, Tice and McKimAnderson 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 (Reference LunardiniLunardini, 1995):

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. Reference OsterkampOsterkamp (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

where the contributions of individual volume heat capacities of ice, water and parent material are weighted by their volumetric contribution (Reference OsterkampOsterkamp, 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 (Reference Clarke, Collins and ThompsonClarke and others, 1984; Reference EchelmeyerEchelmeyer, 1987). Though the magnitude of this energy source is likely reduced under ice sheets (Reference EchelmeyerEchelmeyer, 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 Reference Clarke, Collins and ThompsonClarke and others (1984) and Reference PatersonPaterson (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

where z _s is the elevation of the ice surface (Reference PatersonPaterson, 1994). An important assumption behind Equation (9) is that subglacial water pressure is equal to the ice overburden pressure (Reference PatersonPaterson, 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 Reference Boulton, Caban and van GijsselBoulton 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 (Reference HookeHooke, 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.

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 (Reference Touloukian, Judd and RoyTouloukian 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 Reference MooersMooers (1990), and demonstrated by Reference Boulton, Caban and van GijsselBoulton 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 (Reference Huybrechts, Letréguilly and ReehHuybrechts and others, 1991; Reference ReehReeh, 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 Reference Jenson, MacAyeal, Clark, Ho and VelaJenson 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 Reference ClarkClark 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.

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 Reference Walder and FowlerWalder and Fowler, 1994, or R channels of Reference RöthlisbergerRöthlisberger, 1972). Instead, water would have been forced to drain through underlying lithified sedimentary rocks. Like Reference PiotrowskiPiotrowski (1997) and Reference Fountain and WalderFountain 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 (Reference Boulton, Caban and van GijsselBoulton 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 (Reference Matthess and UbellMatthess and Ubell, 1983), is it possible to evacuate all basal meltwater through the groundwater system. One result of this, also noted by Reference PiotrowskiPiotrowski (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. Reference Boulton and CabanBoulton 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 (Reference Attig, Mickelson and ClaytonAttig and others, 1989; Reference Clayton, Attig, Mickelson, Mickelson and AttigClayton 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 (Reference Borowiecka and EricksonBorowiecka and Erickson, 1985; Reference Attig, Mickelson and ClaytonAttig and others, 1989; Reference Colgan and MickelsonColgan and Mickelson, 1997). This has been interpreted as indicating decreasing basal sliding velocity toward the ice margin (Reference ChorleyChorley, 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 (Reference Maher and MickelsonMaher and Mickelson, 1996), thus allowing sufficient time for gradual thawing of the frozen zone. Reference Mickelson, Clayton, Fullerton, Borns and WrightMickelson and others (1983) and Reference Stanford and MickelsonStanford 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 (Reference Attig, Mickelson and ClaytonAttig and others, 1989; Reference Colgan and MickelsonColgan 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 (Reference Colgan and MickelsonColgan 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 (Reference Mackay and MathewsMackay and Mathews, 1964; Reference Boulton, Caban and van GijsselBoulton 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 Reference BegetBeget (1986), Reference ClarkClark (1992) and Reference Jenson, MacAyeal, Clark, Ho and VelaJenson 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 (Reference Mickelson, Clayton, Fullerton, Borns and WrightMickelson 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 (Reference Wright, Sims and MoreyWright, 1972; Reference MathewsMathews, 1974; Reference ClarkClark, 1992; Reference Colgan, Mickelson and AttigColgan, 1999). Reference ClarkClarks (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). Reference Colgan, Mickelson and AttigColgan (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).