Summary of further model features

Model equations based on the shallow-ice approximation (Hutter, 1983; Morland, 1984).

Two-dimensional flowline model incorporating flow-band width variation, based on finite differences.

Polythermal calculation (as described above).

Lithosphere; evolution of temperature calculated in a layer of thickness H _r (this requires knowledge of lithosphere density ρ_r, specific heat c _x and heat conductivity K _r); reaction on varying ice load accounted for by assuming local isostasy (based on ice density ρ and asthenosphere density ρ_a) with time lag τ_V.

Horizontal resolution: 17.5 km.

Time step: 1a for the calculation of velocity and ice topography, 10a for the calculation of ice temperature and water content.

In the vertical; σ coordinates for each of the three regions (lithosphere: 11 equidistant grid points, temperate ice: 11 equidistant grid points, cold ice: 51 grid points with densification towards the bottom), i.e. vertical columns are mapped to [0, 1] intervals.

Temperature and water-content equation: implicit discretization of z derivatives, explicit discretization of x derivatives (upwind scheme for the horizontal advection terms).

Height-evolution equation: implicit discretization of x derivatives.

Glen’s flow law with n = 3 and an enhancement factor E = 3 (accounting for the apparent stiffness reduction of ice accumulated during glacial periods); rate factor A depending on homologous temperature in the case of cold ice (commonly used exponential law A(Tʹ) = A ₀exp(−Q/R(T ₀ + Tʹ)) with two different activation energies Q for Tʹ < − 10°C and Tʹ > − 10°C, respectively and the absolute temperature T ₀ = 273.15 K, see Paterson (1994)) and on water content in the case of temperate ice (A _t(ω) = A(Tʹ = 0°C) × (1 + 184ω), according to Lliboutry and Duval (1985)).

The values for the physical quantities used in the model are listed in Table 1.

3.1. Flowline topography

Our main interest is the behavior of the Laurentide ice sheet in Hudson Bay and Hudson Strait. We therefore consider an approximate flowline, running from the position (ϕ, λ) = (40° N. 99° W), south of the ice margin determined for the Last Glacial Maximum 21 000 calendar years ago, through Hudson Bay and the center of Hudson Strait. Its eastern margin is situated at (ϕ, λ) = (60.9° N, 60.8° W), beyond the sill separating Hudson Strait from the open sea. Figure 1 shows the course of this flowline together with the regions of crystalline bedrock and unlithified sediment, determining the basal sliding regime. Figure 2 depicts the present bedrock topography z= b ₀, assumed to be approximately in isostatic equilibrium with no ice load.

To correct for flow-band width denoted by s) variation, we introduce a term proportional to ds/dx in the mass-balance equation governing the evolution of ice thickness H, so that

(5)

(q _x , horizontal mass flux; surface accumulation-ablation function; , basal melting rate), and estimate a flow-band width s(x) around the constructed flowline. We assume s(x) to be approximately constant in the southern part and use the present margin contours of Hudson Bay and Hudson Strait as flow-band limits in the northern part (Fig. 1). This affects mainly the transition zone between Hudson Bay and Hudson Strait, where the correction accounts for the channeling effect of Hudson Strait on the ice flow toward the North Atlantic.

Table 1.

Values for the distinct physical quantities used in the model calculations

Fig. 2.

Present bedrock topography along the modelled flowline (see Fig. 1), assumed to represent the isostatic equilibrium with no ice load, z = b₀, x = 0 represents the southern flowline end, x = 3500 km represents the north-eastern end where icebergs are assumed to discharge into the North Atlantic.

3.2. Standard glacial boundary conditions

All model experiments simulate the Laurentide ice sheet under time-independent, glacial climate conditions. For the mean annual air temperature, T _ma, above the ice we apply a function of latitude, ϕ, and height above sea level, h, derived by Reeh (1991) for the Greenland ice sheet under present climate conditions:

(6)

with θ₁ = 38.38°C, θ₂ = −0.007924°C m⁻¹ and θ₃ = −0.7512°C × (deg N. lat.)⁻¹. The value of θ₁ given here is 10°C below Reeh’s value to account for glacial conditions.

The accumulation—ablation function consists of two parts, namely surface snowfall, S ₀, and surface melting, M, from which it can readily be obtained via . As for the snowfall, ice-core analyses indicate that glacial accumulation rates were about one-hall to one-fourth of present values. We thus choose S ₀ = 0.3 ma⁻¹.

The ablation rate (surface melting) is parameterized according to Calov (1994) following the degree-day model (Braithwaite and Olesen, 1989; Reeh, 1991). This method couples the melting rate linearly to the air-temperature excess above 0°C; it is distinguished between (i) transformation of snow via water to superimposed ice (i.e. ice above firm and (ii) melting of ice.

The ice-surface temperature T _s is assumed equal to the mean annual air temperature, T _ma, unless the formation rate of superimposed ice, M*, exceeds the ice-melting rate, M. In this circumstance, an empirical firn-warming correction due to the latent-heat release of the remaining superimposed ice is applied (Reeh, 1991);

(7)

with the firn-warming correction coefficient μ_fwe = 24.206°C m⁻¹ a. Naturally, T _s must not exceed the melting point of ice. Hence, we constrain T _s to a maximum of −0.001°C; the slightly negative value is artificially introduced to avoid a temperate ice-sheet surface. This circumvents problems associated with the upper-boundary condition for water content and has little effect on the outcome of our study.

The remaining boundary condition affects the lower boundary of the model domain, the lithosphere base. Here, the geothermal heat flux is required by the mode and we use the standard value for Precambrian shields, namely (Lee, 1970).

3.3. Model experiments

In the following, we wish to discuss ten model experiments, all of them aimed at building up the Laurentide ice sheet under idealized, time-independent glacial climate conditions, starting from an ice-free initial condition. Run rlp01 is our control experiment and is carried out subject to the standard conditions defined in sections 2 and 3.2: in runs rlp02–rlplO several properties are varied, namely

Run rlp02: subglacial meltwater storage in the unlithified sediment regions included, acting as reservoir of latent energy. Hereby, the evolution of the thickness of the meltwater column, H _w, is calculated by

(8)

H _w corresponds to an amount of latent energy per unit area

(9)

at the ice-sediment interface. It is furthermore assumed that the restriction H _w ≤ 25 m holds and any surplus is regarded as draining away. Because of temperature steadiness, the effect of the presence of subglacial meltwater is that at any position the ice base remains temperate as long as H _w > 0 is fulfilled; in other words, the entire amount of meltwater must refreeze before the basal temperature can switch from temperate to cold.

Run rlp03: θ₁ = 37.38°C (mean annual air temperature decreased by 1°C).

Run rlp04: c _sl = 50 a⁻¹ (basal-sliding coefficient cut down in half).

Run rlp05: θ₂ = −0.007 C m⁻¹ (lapse rate for mean annual air temperature increased by approximately 0.001°C m⁻¹).

Run rlp06: S ₀ = 0.2 ma⁻¹ (snowfall decreased by one-third).

Run rlp07: E = 2.5 (enhancement factor in the flow law decreased by 0.5).

Run rlpOB: τ_v = 10 000 a (time lag for bedrock sinking increased by 7000 a).

Rim rlp09: polythermal mode switched off. Instead solution of the temperature equation for cold ice in the entire ice-sheet domain, afterwards resetting of temperatures that exceed the pressure-melting point on the pressure-melting point.

Run rlp10: calculation without horizontal advection in the temperature and water-content equations.

The runs cover an iteration time of 60 000 years, representing the approximate time the Laurentide ice sheet had for growth during the Wisconsin glacial period.

3.4. Discussion of model results

In Figure 3, the total ice volume of the modelled transect V _tot, the ratio of temperate to total ice volume P _temp the ice-covered area A _i,b, the area covered by temperate ice A _t,b, and the iceberg discharge from Hudson Strait qHS are depicted for all runs as functions of time. Here, qHS assumed equal to the mass flux q _x through the last grid cell x ϵ [3482.5 km, 3500 km] of the model domain immediately adjoining the ocean (see Fig. 2). This is indeed a good measurement for the iceberg discharge, because it represents the amount of ice flow towards the ocean that is not compensated by melting, i.e. the calving rate. If the second-to-last grid point at x = 3482.5 km is ice-free and therefore the mass flux through the last grid cell is equal to zero, the seaward ice flow from the inner ice-sheet region is annihilated by surface and basal melting before reaching the coast, and thus calving does not appear (no iceberg discharge).*

In all cases except rlp02. P _temp and A _t,b, undergo an oscillatory behaviour. In rlp02, the A _t,b oscillations are much smaller and vanish after 30 ka but p _temp, remains slightly oscillating. Of particular interest in our effort to understand Heinrich events is the iceberg discharge through Hudson Strait, qHS. This variable is seen to be oscillatory in all runs except rlp02 and rlp06, demonstrating that the presence of discharge events is not particular to a special set of parameters. The onset of iceberg-discharge oscillations is in all cases delayed 15–25 ka from the inception of the ice sheet. This is because the ice sheet needs a finite amount of time to reach the North Atlantic by growing down the Hudson Strait channel. Iceberg-discharge in run rlp02 is predominantly smooth (we are unable to judge the significance of two tiny negative peaks during the last 4 ka of the model time). Also of note is the iceberg-calving flux of rlp06 which does not become oscillatory until the last 1000 years of the model time. This is explained by the fact that lower snowfall in this run means that the ice sheet does not grow sufficiently fast to reach its iceberg-calving terminus until late in the model run.

The oscillations of A _t,b, p _temp and qHS coincide; high .A _t,b and P _temp values entail high qHs values and vice versa. This is in agreement with the idealization of the binge/purge theory: extensive temperate ice at the base provides extensive basal sliding and, as a consequence, extensive iceberg discharge. Due to the high ice velocities, advection of cold surface ice towards the bottom is enhanced during such an event, resulting in a decrease of A _t,b, until basal sliding can no longer sustain a significantly elevated iceberg-calving flux. During the quiet phase, when little temperate ice exists at the bed, the ice sheet grows. As it grows, the base becomes warmer and eventually temperate basal conditions build up sufficiently to initiate the next discharge event.

Why do the oscillations of qHS not appear in rlp02? This can be understood by the influence of the meltwater layer on the basal temperature conditions. In section 3.3, we stated that with a meltwater layer included, a temperate ice base must remain temperate until the underlying meltwater has refrozen completely. This condition causes a thermal “inertia effect” on the meltwater layer and its stored latent energy, hampering the switch from a predominantly temperate base to a predominantly cold one after a discharge event. Therefore, as opposed to the other runs with no meltwater layer taken into account, temperate basal conditions prevail, leading to a continuous iceberg discharge.

Fig. 3.

Fig. 3. Time evolution of the total ice volume V_tot (10³ km²), the ratio of temperate to total ice volume p_temp (per mille), the ice coverage A_b, (10³ km) along the modelled flowline (dashed line; total ice-covered basal area A_i,b, solid line: basal area covered by temperate ice A_t,b and the iceberg discharge into the North Atlantic qHS (km² a⁻¹ for the ten runs rlp01–rlp10 (rlp01 control, rlp02 subgtacial meltwater reservoir included. rlp03 mean annual air temperature decreased. rlp04 basal sliding decreased, rlp05 lapse rate for mean annual air temperature increased, rlp06 snowfall decreased, rlp07 enhancement factor decreased, rlp08 time lag for bedrock sinking increased, rlp09 polythermal model switched off. rlp10 no horizontal advection).

In rlp09, the polythermal calculation has been switched off. This has a major effect on the ratio of temperate ice volume to total ice volume, P _temp, whose absolute values and oscillations are much bigger than for the polythermal simulation rlp01 (note the different axis scale for p _temp of run rlp09 in Figure 3). This is due to the fact that solving the temperature equation for cold ice in the whole domain provides unrealistic temperatures above the pressure-melting point in pans of the ice sheet. Because of the implicit solution technique, this leads to a temperature field slightly too warm in the entire ice sheet, so that the amount of temperate-ice regions is over-predicted. However, the effect on A _t,b, and qHS is only marginal; the main characteristics of the oscillations are the same in both cases. Evident differences occur in the time interval between 30 and 40 ka, where the oscillations provided by rlp09 appear to be more regular than the ones of rlp01. This shows that the oscillations are mainly due to the switch between basal sliding and basal adhesion; the presence of a soft (due to the strong dependence of the rate factor A _t(ω) on water content temperate-ice layer above the base is of minor importance.

Run rlp10 was carried out without accounting for horizontal advection in the evolution equations for temperature (cold regions; Equation (1)) and water content (temperate regions; Equation (2)). The impact on the dynamic behaviour of the modelled ice-sheet flowline is remarkable. The ratio P _temp is more than 10 times bigger than in the control run rlp01, and the oscillations of A _t,b and qHS are of higher frequency and lower amplitude. In particular, after the onset of iceberg discharge at t = 29.2 ka it does not break down anymore, i.e. this run does not provide distinct discharge events interrupted by quiet periods. This is due to the fact that horizontal advection is responsible for the transport of cold surface ice from the inner ice-sheet regions towards the margin, where it enforces the switch from temperate basal conditions (basal sliding) to cold basal conditions (basal adhesion). However, since vertical advection also transports cold surface ice to the base, the binge/purge mechanism as explained above can still be sustained, so that some oscillatoric behaviour remains.

The range of model integration time steps that provide stable runs is very small. In order to test the possible influence on results, we carried out run rlp01 two more times with (i) △t = 0.5 a/5 a and (ii) △t = 2 a/20 a (instead of the standard set △t = 1 a/10 a (see “summary of farther model features”). Both of these runs crash between t = 50 ka and t = 55 ka; since then, they provide almost the same dynamic behaviour of the modelled ice sheet as standard run rlp01. Furthermore, the simulated periodicity (> l ka) is much bigger than the integration time steps, so that a purely numerical reason for the simulated oscillations can in all probability be excluded. However, we have not conducted three-dimensional simulations to check the consequences of the flowline restriction; with a horizontal resolution fine enough to include adequately the topography of Hudson Strait, this would require a large amount of CPU lime. Payne (personal communication (1995)) carried out three-dimensional simulations on the Laurentide ice sheet with a somewhat coarser resolution (△x = △y = 50 km) and observed subsequent on- and off-switch events of draining tee streams, indicating that oscillatoric behaviour still holds in the three-dimensional case.

Fig. 4.

Run rlp01; ice-sheet topography (z = h free surface, z = b ice base) along the modelled flowline for t = 58.7 ka (solid line, corresponds to high A_t,b and maximum qJHS) and t = 59.1 ka (dashed line, corresponds to low A _t,b and qHS = 0).

Figure 4 shows the ice-sheet geometries for rlp01 at two stages in the cycle of iceberg discharge. At an active point in the cycle, specifically at t = 58.7 ka, corresponding to the second-to-last qHS peak in Figure 3, the seaward margin has advanced to the coast and the ice surface near the dome has lowered. For this time, qHS amounts to a total value of 2.71 km² a⁻¹, from which 0.99 km² a⁻¹ result from basal sliding according to Equation (4) and 1.72 km² a⁻¹ from internal deformation. At a quiet stage in the cycle, specifically at t = 59.1 ka, where qHS has come to zero afterward, the seaward margin has retreated and the surface elevation has grown again.