2.1. Climate and surface mass balance

The EBM climate model is a non-linear version of the original two-dimensional EBM of Reference North, Mengel and ShortNorth and others (1983) which computes a full seasonal cycle with albedo feedback. It has been validated against the predictions of coarse-resolution GGMs under LGM boundary conditions for large spatial scales in the extra-tropics (Reference Hyde, Crowley, Kim and NorthHyde and others, 1989). The EBM derives sea-level temperature from the radiative energy balance of each tropospheric and mixed-layer ocean column using a spherical harmonic representation truncated at degree and order 11 according to the following evolution equation:

The absorbed short-wave insolation is given by the product of the solar constant, Q, the solar distribution function, S/4 and the effective co-albedo a(r,l) which tracks the seasonal variation of snow and sea-ice lines. Longwave emission is represented by a linear function of sea-level temperature derived from satellite observations, [A + BT(r,l)]. Horizontal transport is modeled as a diffusion process, with diffusion coefficient D tuned to preserve the present-day observed mean meridional temperature gradient. The heat capacity, C, is a function of surface type, but does not take into account seasonal variations. The model also incorporates a parameterized time-dependent North Atlantic heat flux (NAHF) and takes into account the radiative forcing of varying atmospheric PC0₂ by assuming the Vostok chronology applies globally (Reference Barnola, Raynaud, Korotkevich and LoriusBarnola and others, 1987; Reference JouzelJouzel and others, 1987).

Given that the strength of energy-balance models is their ability to capture the impact of radiative perturbations on surface-sea-level temperature, the raw output of the EBM is not directly coupled to the ISM. Rather the EBM is employed to compute perturbations to a present-day observational climatology (based upon a 14 year mean (1982–95) of re-analyzed European Center for Medium-range Weather Forecasts monthly mean-temperature fields, ds090.0 NMG/ NGAR global re-analysis project dataset). Adjustments from sea level to surface temperature are carried out by assuming a fixed 7.5°C km⁻¹ environmental lapse rate.

The ablation model uses the now standard PDD formalism (Reference ReehReeh, 1991) with ablation parameters, unless otherwise stated for specific model runs, of 3 and 8 mm/PDD of ice equivalent for snow and ice respectively. Furthermore, 60% of melted snow is assumed to refreeze to superimposed ice. Monthly PDD are computed from monthly mean temperatures with a normal statistical model using a standard deviation of 5.5° G.

Calving represents an important component of ice-sheet mass balance, however given the small-scale controls and complexities of this process, no credible model is available, especially on "whole-ice-sheet" scales. Therefore, for the sake of simplicity, we assume that complete calving of ice occurs at the 400 m bathymetric contour.

Given its dependence on small-scale atmospheric dynamics, precipitation is much more difficult to constrain in reduced models. In the work performed to date, we have simply chosen to perturb a present-day monthly precipitation climatology (Reference Legates and WillmottLegates and Willmott, 1990) by a 3% reduction per degree cooling of gridpoint monthly mean sea-level temperature with respect to the present, in weak correspondence with the effect expected as a consequence of the temperature dependence of saturation vapour pressure. Furthermore an added desert-elevation effect (Reference Budd and SmithBudd and Smith, 1981) is also incorporated, whereby precipitation is decreased by 50% for each kilometer by which the surface elevation exceeds the height of 2 km. Rain/snow fractions are computed using a normal statistical model with a standard deviation of 4.5°C based on monthly mean temperatures, again using the logic of the PDD formalism.

2.2. Ice and bed dynamics and thermodynamics

The initial bedrock topography is derived from the ETOP05 dataset using a 2:1 weighted average of the mean and maximum gridcell elevations to capture the high-elevation nucleation zones better. Deformation of the bedrock is computed on the assumption that it occurs via a simple, local, damped return to isostatic equilibrium with an e-folding time of 4 kyr.

The 3–D thermomechanical ISM incorporates full coupling between an ISM based on Glen’s flow law and an ice-thermodynamics solver. The ice dynamics assumes the Arrhenius-type temperature dependence discussed fully in Reference PatersonPaterson (1994). The evolution equation for the ice thickness, H, is based on the vertically integrated continuity equation for ice mass as:

in which horizontal ice velocities, are determined in the shallow-ice approximation by the expression (Reference HutterHutter,1983):

Here h is the ice-surface elevation above present-day sea level, P _i is the density of ice, g is the acceleration due to gravity, T* is the ice temperature relative to the pressure-melting point and G is the mass balance. In the Glen rheology, the stress exponent n is set to 3.

According to Reference WeertmanWeertman’s (1957) original analysis, basal sliding would be expected to be controlled by a combination of two mechanisms, namely ice deformation and regelation (pressure-melting and refreezing). However, as argued by Reference PatersonPaterson (1994), it is expected that ice deformation should predominate for large-scale ice sheets (or at least for numerical models thereof). As such, the sliding law for the basal velocity of a large-scale ice sheet would be expected to be of the form:

wherever the base is at or near the pressure-melting point. This form obviously avoids the issue of hydrological controls for which no credible large-scale model presently exists. As a reference model, we have chosen to employ a simple Weertman sliding model validated for simulations of the Greenland ice sheet (Reference HuybrechtsHuybrechts, 1996), with A, = 2 × 10⁻¹³ Pa⁻³ m a"¹.

The temperature field for the ice sheet is computed on the basis of an evolution equation for the internal energy that incorporates the influence of both horizontal and vertical advection, vertical diffusion (with temperature-dependent diffusion coefficient k(T)), and heat generation due both to ice deformation, Q_d, and sliding (incorporated as a basal-boundary condition). The evolution equation for the internal temperature of the ice is therefore:

As is common in models such as this, horizontal diffusion is not included on account of the dominant influence of horizontal advection. The thermodynamic solver employs 35 equidistant vertical layers and is implicitly coupled to a local-bed-thermal model with five layers spanning the top 2 km of the bedrock. The lower boundary condition for the bed-thermal model is assumed to be provided by the digitized geothermal heat-flux map of Reference Pollack, Hurter and JohnsonPollack and others (1993), unless otherwise stated.

2.3. Model certainties and uncertainties

In considering the properties of this coupled model, it is worthwhile to clarify which components are well-constrained and which are not. The thermomechanical ISM, with only weak reliance on the basal-sliding component, has been validated in detail in the context of the European Ice Sheet Modelling Initiative (EISMINT) model-inter-comparison project (Payne and others, 1999; Ritz and others, in preparation). Analyses of the predictions made with the models demonstrated close agreement in predicted ice-sheet topographies. It should be noted that all these models were based on the same rheology and temperature dependence of the flow law. The simple relaxation model used to represent glacial-isostatic adjustment has also been validated (Reference Tarasov and PeltierTarasov and Peltier, 1997a) against the same viscoelastic gravitationally self-consistent global model of isostatic adjustment used to compute the ICE-4G reconstruction (see, e.g., Reference PeltierPeltier (1998b) for a recent review of the required theory).

The PDD formalism, employed to compute ablation, is highly reduced. However, although there are variations in the required parameter values when the model is tuned to observations from different present-day glaciers and ice sheets (see, e.g., Reference BraithwaiteBraithwaite, 1995, Reference Braithwaite1996; Reference Jóhannesson, Sigurdsson, Laumann and KennettJohannesson and others, 1995, and references therein), this sensitivity largely translates into a freedom to tune so as to provide an adequate LGM ice volume while allowing deglaciation to obtain under appropriate conditions. It is unlikely, though unproven, that a much more detailed ablation parameterization would result in a major change in ice-sheet topography that could not be obtained through simple retuning of the PDD model.

In considering the coupled model as a whole, it is therefore clear that the climate forcing and basal components of the model are in fact the least well-constrained. The EBM is a coarse-resolution energy-balance model that neglects the influence of atmospheric dynamics. Climate forcing is most important in its impact on surface mass balance, but it also contributes the upper-boundary condition required to drive the ice-bed thermodynamics. Furthermore, the precipitation parameterization is based on a highly simplified representation of the physical dynamics. Therefore, it is possible that uncertainties in the precipitation and climate components of the coupled model could of themselves account for the different flow-law enhancements required to simulate the Greenland and Laurentide ice sheets.

Basal-boundary conditions of the ISM do represent a poorly constrained and incomplete representation of actual physical processes. Though a sliding model is incorporated, there is presently little constraint on the actual form of the sliding law. Furthermore, no account is taken of basal hydrology and differing basal-surface materials. Two other processes highly sensitive to basal hydrology, namely ice streaming and deformation of subglacial till, are also not explicitly included in the model. Both sliding and till deformation reduce coupling between bed and ice and thereby enhance ice flow.

These weak constraints on modeled climate and basal processes motivated the following investigation into the possible roles that uncertainties with respect to these processes might have in explaining the aspect-ratio challenge that the IGE-4G reconstruction presents to those who would model the evolution of continental-scale ice-age ice sheets.

3.1. Impact of different sliding models

Assuming the validity of the Glen rheology, it is unlikely that the exponent in the sliding law would be below 3. However, enhancement of basal velocities through the deformation of subglacial material (Reference PatersonPaterson, 1994) is thought to have been widespread for the Laurentide ice sheet (Reference Boulton and JonesBoulton and Jones, 1979; Reference AlleyAlley, 1991) and could possibly mimic a slid-ing-law exponent as low as 1, though very large exponents approaching plastic flow have also been proposed to govern till deformation (Reference KambKamb, 1991). Lower values of the exponent in the sliding law would promote lower aspect ratios due to increased sliding away from the margins. Two simulations were carried out with linear sliding laws tuned to match the sliding velocities of the reference Weertman law for basal shear stresses of 100 kPa. As is evident from the simulated surface transects and ice-volume chronologies (Fig. 4), a sliding exponent, m, of 1 does improve both the southerly LGM extent of Tbasel0 and deglaciation, but has little impact on the elevation of the divide surface. Again, the frozen northern flank of the Laurentide ice sheet anchors the divide position and even the strongly increased sliding near the core, as compared to the exponent 3 base Weertman sliding law, associated with the 2 × enhanced linear sliding law is only able to slightly "eat away" at the divide. Stronger enhancements to the linear sliding law could have more impact on the divide. However given that the 2× enhanced linear sliding Tbasel0 model already has maximum meridional velocities of over 550 ma⁻¹ near the margin, and maximum sliding velocities of about 300 m a⁻¹, the numerical and physical acceptability of such enhanced and sustained large-scale flows would start to become doubtful. The surface-transect irregularities, associated with the linear sliding law, arise from finger instabilities (Payne and others, 1999) and from the transect not being along a flowline.

Fig. 4.

Impact of linear-sliding model (with enhancement factors 1 and 2) applied to the Tbasel0 model on (a) NA ice-volume chronology and (b) surface-elevation transect. Tbasel0 m = 1 sliding law is tuned to match the Weertman sliding law at a basal shear stress of 100 kPa.

The improvement to Tbasel0 due to exponent 1 sliding is on the optimistic side, sliding occurs everywhere that the base is near temperate, given the lack of basal hydrology and accounting for the influence of basal-surface type. The predominance of the Canadian Shield in northern Ontario would act to reduce the area of sliding and till deformation. Furthermore, the computed basal-temperature field is quite sensitive to the geothermal heat flux, which is not strongly constrained. Use of an alternative heat-flux field based on the map of Reference Blackwell and SteeleBlackwell and Steele (1992) resulted in a decrease in the percentage area of North American ice (NA) at LGM with a temperate base from 47–34% (Reference Tarasov and PeltierTarasov and Peltier, 1999). This would have a significant impact on models relying heavily on basal processes, which all require basal temperatures to be at the pressure-melting point.

When all is taken into account, these results suggest that enhanced basal flows from low-exponent sliding, possibly due to a non Glen flow ice rheology or to extensive till deformation, can only partially improve the aspect ratio of modeled Laurentide ice sheets.

3.2. Impact of mass-balance chronology uncertainties

The largest sources of uncertainty with respect to simulated mass balance are the climate fields. Given the complexity of the climate system, it is very difficult to assign clear error bars to GCM output, let alone to the much simpler EBM and precipitation parameterizations used in the coupled model. We have therefore chosen to proceed via reducto absurdum. Instead of attempting to pinpoint error bars in the time dependence of the EBM, we have chosen two extreme glacial histories to represent strongly contrasting archetypes. These contrasting scenarios are designed to produce LGM Laurentide ice sheets that are far from equilibrium, thereby offering a limited test of the suggestion of Reference Marshall, Tarasov, Clarke and PeltierMarshall and others (in press) that paleoclimate-ice-sheet fluctuations that produce far-from-equilibrium LGM Laurentide ice sheets may be able to account for much of the required flow-law discrepancy. The first scenario, referred to as "fast growth" in Figure 5, involves suppression of Laurentian glaciation from 90–40kyrBP via an ad hoc 4°C increase of surface temperatures, as computed by the EBM, followed by a fast glaciation phase aided by an applied 4°C decrease of surface temperature from 30–25 kyr BP. This modification delivers a small improvement to both the dome elevation and the southward extent. However, the final ice volume does not differ appreciably from that of Tbasel0. Furthermore, dome peak position is not altered, thereby offering little improvement in the southward extension. An even more extreme version of this scenario ("extreme growth") with an applied 5°C warming from 90–35 kyr BP followed by an applied 5°C cooling from 33–23 kyr BP does result in major "improvements" to the Tbasel0 LGM Laurentide ice sheet. Deglaciation is slightly improved to almost 50% residual (Fig. 5a). But most significantly, the surface transect of this version of Tbasel0 is a relatively close match to that of Tbase, apart from a split divide structure, slightly lower surface elevations and about a degree of latitude less southern extent. This scenario provides an example of the impact of extreme disequilibrium.

Fig. 5.

Impact of different ad hoc glacial chronologies on (a) NA ice-volume chronology and (b) surface-elevation transect. The early growth model uses the Tbasel0 base model with an imposed extra TC cooling 115–35 kyr BE Early growth 2 uses same imposed cooling in addition to a 30–25kyrBP TC imposed warming. The fast-growth model uses an imposed 4°C warming from 90–40 kyr BP and an imposed 4°C cooling for 30–25 kyr BE The extreme growth model uses a 90–35 kyr BP 5° C imposed warming followed by an applied 5°C cooling for 33–23 kyr BE All models use factor 2 sliding enhancement.

The second scenario represents an opposite extreme, early glaciation aided by an applied 2°C decrease of surface temperatures 115–35 kyr BP ("early growth" in Fig. 5). This results in a large excess of ice at LGM. Application of an ad hoc 2°C increase in surface temperature from 30–25 kyr BP ("early growth 2") brings the LGM ice volume within reasonable bounds. The same small improvement to southerly extent as for the "fast growth" scenario obtains. However, in this case peak dome elevation is not lowered. It should also be stated here that such a large pre-LGM ice volume would have left a strong isostatic signature, contrary to observations.

Both of the above non-extreme scenarios produce ice-volume chronologies that (distantly) bound the range of possible NA ice-volume chronologies that could be inferred from the spectral-mapping project (SPECMAP) δ¹⁸O record of Reference Imbrie, Berger, Imbrie, Hays, Kukla and SaltzmanImbrie and others (1984) given the continental partition of ice volumes in both Tbase (Reference Tarasov and PeltierTarasov and Peltier, 1999) and the ICE-4G reconstruction. Though SPECMAP represents a filtered record of past global ice-volume changes, it is difficult to conceive how the extremely fast growth scenario could be reconciled with it. Furthermore, the large resultant isostatic disequilibrium at LGM would imply that the ICE-4G reconstruction is severely underpre-dicting LGM NA ice volume, which would likely be difficult to reconcile with the strong sea-level constraints from the Barbados coral records (Reference FairbanksFairbanks, 1989). Therefore, possible uncertainties in the time dependence of the simulated gross temperature field or mass-balance history do not appear to offer a solution to the aspect-ratio problem or required flow-enhancement discrepancy.

3.3. Impact of the geographic pattern of mass balance

Given the strong sensitivity of ice-sheet mass balance to surface temperatures on account of ablation-rate sensitivities, a trial hypothesis we wished to test is that mass-conserving changes to the geography of accumulation within the non-ablation zones of the ice sheet, the dry-snow and percolation zones, a subarea of the accumulation zone where no runoff occurs, would not have a significant impact on the ice sheet. Model simulations using Tbase parameters were carried out with mass-conserving extremal modifications to the geography of accumulation. At one extreme, the mean accumulation for the whole non-ablation zone was used at each gridpoint. The opposite extreme involved a quadratic enhancement of the meridional accumulation gradient at each longitude. Interestingly, both of these modifications resulted in a strong reduction of the degree of glaciation. To tune the modified models for a reasonable LGM ice volume, the precipitation-reduction factor (Rp) was reduced to 1% for the enhanced-accumulation-gradient scenario (Fig. 6). However, even with the elimination of the Rp (i.e. set to 0), the average-accumulation scenario produced insufficient ice volume at the LGM, slightly less than that of Tbase10 10. Furthermore, southward extent of the ice complex of this latter simulation is even less than that of Tbasel0 (Fig. 6b). Unlike most other simulations considered in this study, mean thickness chronologies for this scenario were significantly greater than those of other scenarios with even more LGM ice volume, not shown. This is the consequence of the increase of accumulation over colder regions in this scenario which allowed increased build-up of cold ice and a subsequent increase in the modeled aspect ratio.

Fig. 6.

Impact of modifications to geographic pattern of accumulation on (a) NA ice-volume chronology and (b) surface-elevation transect. "Avg accum." applies non-ablation-region average accumulation to every point in that region. “Enhanced accum.” uses a mass-conserving quadratic enhancement of the accumulation gradient in the non-ablation zone. All models use factor 2 sliding enhancement and stated precipitation reduction factors Rp (default value = 3%).

The enhanced-gradient version of Tbasel0 (and Rp = 1%) leads to some improvements in the deglaciation phase of the simulation (Fig. 6a). It also leads to the southward displacement of the divide in better agreement with ICE-4G and Tbase. Most significantly, this scenario improves the aspect ratio of Tbasel0 both by decreasing divide elevation to within tens of meters of that of Tbase and by extending the southern margin beyond that of Tbase. This improvement relies on a major change in precipitation patterns with precipitation reduced by a third to a half or more north of 55° N compared to that computed for Tbase for transects at 86.5° W (not shown). Precipitation rates near the margin increase by a factor of two and surpass 1.5ma⁻¹. This modification appears to provide one solution to the required flow-law discrepancy. However, given that Tbase already underpredicts precipitation almost everywhere over the Laurentide ice sheet, in comparison to T32 GCM simulations using ICE-4G LGM boundary conditions (Reference Vettoretti, Peltier and McFarlaneVettoretti and others, in press), there is little support for such an extreme modification.

The above results negate our trial hypothesis. Ice-sheet evolution is sensitive to the internal geography of accumulation. Given the large uncertainties with respect to GCM precipitation fields, this suggests that definitive modeling of ice sheets is dependent on accurate precipitation fields, which are difficult to obtain even with state-of-the-art GCMs.

Surface temperature is the most important control on surface mass balance, especially given the sensitivity of ablation to temperature. The margin of LGM Tbase stays between the –2.5° and –5°C mean annual isotherms along the entire southern margin. Even the reduced ice flow in Tbasel0 does not appear to hinder its southward penetration into warmer regions significantly as its southern LGM margin stays within a degree of –4°C; a result that likely would not obtain for simulations without full two-way coupling between ice and climate. Therefore, the most obvious mechanism by which we might extend the southern margin would simply be to modify the surface-temperature field.

An ad hoc cooling in the region bounded by 70°–90° W and 40°–60°N applied to the Tbase model (Reference Tarasov and PeltierTarasov and Peltier, 1999) was found to produce a southeastern lobe corresponding with the geomorphological data (Reference Dyke and PrestDyke and Prest, 1987) used to constrain the southern margin of the ICE-4G reconstruction. This result is hardly surprising given the sensitivity of surface mass balance to surface temperature. More interesting are recent analyses that suggest this imposed regional cooling may have a physical basis. GCM simulations with LGM and ICE-4G boundary conditions (Reference Vettoretti, Peltier and McFarlaneVettoretti and others, in press) produce a cooling effect in this region that is much stronger than that predicted by the EBM, even after correcting for differing surface elevations with a 7.5˚Ckml atmospheric lapse rate. This arises as a result of the strengthening of the Rossby wave field over LGM North America which results in more northerly flows over the Great Lakes region. Applying the above mentioned ad hoc cooling to Tbasel0, however, results in a large excess of LGM ice volume and high dome center. A reduction of ice volume and dome elevation, along with maintenance of the gains in southerly extent, can be obtained with an increase of the desert-elevation gradient from a two-fold to an eight-fold reduction in precipitation for every kilometer increase in surface elevation over 2 kma.s.l. However, given that the same GCM study indicates that our reference model already underpredicts precipitation almost everywhere over the Laurentide ice sheet, other mechanisms must be invoked to explain the low surface elevation of the ICE-4G Laurentide dome.

Strongly enhanced sliding with either an 8–fold sliding enhancement of the reference sliding law or use of the linear sliding law, with sliding enhancement set to 1, coupled with an extremal value for the PDD ablation parameters (4.4 and 10.0 mm/PPD of ice equivalent), allows the model to produce southward extent better than that of Tbase and within approximately two degrees of latitude of the ICE-4G margin (Fig. 7). Furthermore, deglaciation is improved significantly with the 8× reference sliding model having even less residual ice at present than Tbase, a not unexpected consequence of the increased ablation-rate parameters. However, both divide position and surface elevation remain relatively unchanged from those of Tbasel0. Though extreme, the above PDD ablation parameters cannot be ruled out on the basis of observational data (see, e.g., Reference BraithwaiteBraithwaite, 1995, Reference Braithwaite1996; Reference Jóhannesson, Sigurdsson, Laumann and KennettJohannesson and others, 1995, and references therein) especially given the uncertainties in the surface temperatures derived from an EBM based upon the assumption of a fixed lapse rate. Slight surface irregularities, situated around 55° N in the LGM-forced, lobe-enhanced, sliding-model topographies, again arise from finger instabilities and from the fact that the transect crosses different flowlines. Stronger enhancements of the sliding law (results not shown) may extend the southern margin even further (by about a half degree of latitude for 2× linear sliding), but stronger fingering instabilities along with maximum sliding velocities >350 m a"¹ call into question the numerical, if not the physical, acceptability of the result. Furthermore, factor 2 sliding enhancements to the above do not significantly modify the position nor the surface elevation of the divide. This again supports the thesis that enhanced basal flows requiring a temperate base cannot extend the southern margin and significantly reduce surface elevation at the divide simultaneously, nor displace the divide southward when the northern flank has a frozen base.

Fig. 7.

Impact of ad hoc forced southeastern lobe development on (a) NA ice-volume chronology and (b) surface-elevation transect. Cooling proportional to the radiative perturbation from varying PCo₂ is applied to the 90°–70° W, 60°–40° N region. "d/e8" version incorporates quartic strengthening of desert-elevation effect, "m = 1 extreme" version uses linear sliding law and 4.0 and 10.0 PPD snow and ice-melt parameters respectively. "8sliding"extreme uses identical PDD parameters but 8 × enhancement of base Weertman sliding law.