North American Ice Sheet reconstructions at the Last Glacial Maximum

Quaternary Science Reviews 21 (2002) 175–192

North American Ice Sheet reconstructions at the Last Glacial Maximum Shawn J. Marshalla,*, Thomas S. Jamesb, Garry K.C. Clarkec a

Department of Geography, University of Calgary, 2500 University Drive NW, Calgary, Canada AB T2N 1N4 Pacific Geoscience Centre, Geological Survey of Canada, 9860 W. Saanich Rd., Sidney, Canada B.C. V8L 4B2 c Earth and Ocean Sciences, University of British Columbia, 2219 Main Mall, Vancouver, Canada B.C. V6T 1Z4 b

Received 18 April 2001; accepted 31 August 2001

Abstract The areal extent of the last glacial maximum (LGM) ice sheets is well known in North America, but there is no direct geological proxy for ice sheet thickness or volume. Uncertainties associated with glaciological and geophysical reconstructions give widely varying estimates of North American Ice Sheet (NAIS) volume at LGM. In an effort to quantify the uncertainty associated with glaciological reconstructions, we conducted a suite of 190 numerical simulations of the last glacial cycle in North America, prescribing different climatic, mass balance, glaciologic, and isostatic treatments for the least constrained model variables. LGM ice sheet reconstructions were evaluated using the well-established geologic record of ice sheet area and southern extent at LGM (Dyke and Prest (1987)). These constraints give a subset of 33 simulations that produce reasonable LGM ice cover in North America. Ice sheet dispositions and the associated parameter settings in this subset of tests provide insight into the plausible range of NAIS thickness, form, and mass balance regime at LGM. Ice volume in this subset of tests spans a range of 28:5238:9
1015 m3 at LGM, with a predominant cluster at 32236
1015 m3 : Taking floating ice and displaced continental water into account, this corresponds to 69294 m eustatic sea level (msl). More than 75% of the accepted tests fall in the range 78288 msl: We argue that this is a plausible estimate of the volume of water locked up in the NAIS at LGM. r 2001 Elsevier Science Ltd. All rights reserved.

1. Introduction There are many outstanding questions regarding the Last Glacial Maximum (LGM) ice sheets in North America. Where were the ice sheet domes and divides? Were these stable or transitory features in the mature ice sheet? Were the deep southern forays of the ice sheet driven by basal flow or climate? What was the precipitation and mass balance regime over the ice? How thick was the ice? The areal extent of the North American ice sheets and the temporal margin history are well known over the last 21 kyr (Dyke and Prest, 1987), but a three-dimensional image of the ice sheet is difficult to construct from the geologic evidence. There is also a dearth of easily interpreted evidence to illuminate the climatic regime and subglacial environment of the LGM ice sheets. We are inevitably led to solid Earth (isostatic), glaciologic, and climate models to explore *Corresponding author. Tel.: 001-403-230-5158; fax: 001-403-2826561. E-mail address: [email protected] (S.J. Marshall).

these questions (e.g., Deblonde and Peltier, 1993; Peltier, 1994, 1996, 1998; Huybrechts and T’Siobbel, 1995; Clark et al., 1996; Jenson et al., 1996; Pollard and Thompson, 1997; Ramstein et al., 1997; Tarasov and Peltier, 1999). This paper presents a suite of simulations of North American Ice Sheet (NAIS) growth and decay through the last glacial cycle in North America. The NAIS complex in our model includes the Cordilleran, Innuitian, and Laurentide Ice Sheets in western, northern, and eastern North America. These ice sheets have separate divides and are dynamically independent through much of the glacial cycle, but they are confluent at LGM. We calculate both individual and collective LGM sea-level impacts of these ice sheets. The numerical model couples a three-dimensional representation of ice sheet thermomechanics (Marshall and Clarke, 1997) with a visco-elastic Earth isostatic model (James and Ivins, 1998). Paleoclimatic forcing is parameterized from atmospheric general circulation model (AGCM) results from Vettoretti et al. (2000). Details of the climate treatment are summarized in

0277-3791/02/$ - see front matter r 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 2 7 7 - 3 7 9 1 ( 0 1 ) 0 0 0 8 9 - 0

176

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

Section 2 and are expounded in detail in Marshall et al. (2000). The suite of simulations was designed to explore the parameter space of model uncertainties and sensitivities, with N ¼ 190 glacial cycle simulations. Our primary objective is to build an ensemble of results that pass geologic scrutiny, based on the reconstructions of NAIS area and southern margin position at LGM (Dyke and Prest, 1987). This model ensemble underlies our estimate of NAIS volume at LGM, as well as any conclusions about the associated paleoclimatic conditions and glaciologic processes. This approach also permits an estimate of uncertainty in the glaciologic reconstructions, i.e. the range of NAIS volumes and paleoclimatic conditions that produce acceptable results. Our sampling of parameter space is broad and systematic, but far from exhaustive. Furthermore, parameterizations in each model component (glaciologic, paleoclimatic, and Earth models) do not fully represent the system; there are simplifications, missing processes, and scale issues. The free parameters in the model therefore underestimate those in the natural system. Our coupled model is representative of the current state-of-the-art in Quaternary ice sheet simulation, however. We feel confident that we have explored the most sensitive free parameters at this evolutionary stage in glacial-cycle models, allowing for the best possible glaciological reconstruction of the LGM ice sheets at this time.

2. Model description A full description of the glaciological model and the paleoclimate parameterization employed in this study is given in Marshall et al. (2000). We summarize the key model ingredients below. The approach has been extended in this analysis through coupling with a full visco-elastic Earth model for isostatic rebound calculations (James and Ivins, 1998). The isostatic treatment is based on the techniques of Peltier (1974) and Wu and

Peltier (1982). Simulations use finite differences on a spherical grid, with l; y; and z denoting longitude, latitude, and vertical co-ordinates, respectively. The vertical co-ordinate z is defined to be positive upwards, with z ¼ 0 representing the modern-day global sea level datum.

2.1. Glaciological model The ice sheet model is a three-dimensional finite difference model following the formulation of Huybrechts (1990), with ice dynamics based on Mahaffy (1976) and ice thermodynamics solved after Jenssen (1977). Given the bed topography and a spatialtemporal climatic history (monthly air temperature and precipitation rates), the primary model outputs of interest here are predictions of ice sheet area, thickness, and volume evolution. Fig. 1 provides a graphical summary of model inputs and outputs. Details of the glaciological model are described in Marshall and Clarke (1997). Advection and diffusion of ice thickness and temperature fields are solved following basic conservation equations (mass, momentum, and energy), subject to a constitutive relationship for ice rheology. Numerous other processes must also be parameterized in the model, and some of the more uncertain and sensitive parameterizations are detailed below. Initial and boundary conditions for the simulations are summarized in the Inputs box in Fig. 1. Climate fields are updated every 100 years in our simulations. This climatic forcing inarguably represents the primary control on ice sheet disposition, but glaciological processes and the exact parameterization of mass balance in the model also play significant roles in dictating ice sheet evolution and form. We experiment with ice rheology and basal flow relationships in this analysis, recognizing that these are two of the most influential controls on ice sheet dynamics and thickness (Clark et al., 1996; Peltier et al., 2000).

INPUTS • Initial Bed Topography h B (λ,θ, t0) • Subgrid Topographic Attributes: h min (λ,θ), h max (λ,θ), roughness • Initial Ice Thickness H (λ,θ, t 0) • Geothermal Heat Flux QG (λ,θ) • Deformable Bed Fraction α (λ,θ) • Monthly Air Temperature and Precipitation, TA (λ,θ, t), P (λ,θ, t)

Glaciologic & Geodynamic Model Integration from t 0 to t m at 2-5 year time steps

BASIC MODEL OUTPUTS • Ice Thickness H (λ,θ, t) • Bed Topography hB(λ,θ, t) • Ice Area A (t), Volume V (t) • Ice Velocity v (λ,θ, z, t) • Ice Temperature T (λ,θ, z, t)

Fig. 1. Table of input and output fields from the coupled ice-climate-Earth model.

177

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

2.1.1. Ice rheology Glen’s flow law, the prevailing ice rheology used in ice sheet models, relates strain rates e’ to deviatoric stresses, s0 in an ice mass (Glen, 1958; Paterson, 1994), e’ ij ¼ BðTI ÞjS02 jðn1Þ=2 s0ij ;

ð1Þ

where n is the flow law exponent, TI is the ice temperature, and S02 is the second invariant of the deviatoric stress tensor, S02 ¼ s0ij s0ji =2:

ð2Þ

BðTI Þ in (1) is an ice stiffness (inverse viscosity) coefficient which follows
 Q BðTI Þ ¼ B0 exp ; ð3Þ RTI where Q is the creep activation energy, R is the ideal gas law constant, and B0 is an observationally derived rate constant. Field- and laboratory-based estimates of B0 vary by an order of magnitude (Paterson, 1994), associated with physical complexities such as crystal anisotropy and the effect of impurities (e.g., dust) on bulk ice deformation (Fisher, 1987; Paterson, 1991; Thorsteinsson et al., 1999). Disregarding details of the directional stress and strain rate components for discussion purposes, Eqs. (1) and (2) give strain rates that are proportional to the stress field to the nth power, e’pðs0 Þn : The n ¼ 1 case reduces to the constitutive relationship for a linear viscous fluid, while n ¼ 3 is the classical choice in application of Glen’s flow law for ice. Peltier (1998) has championed the idea that ice sheets subject to low deviatoric stress (thin, flat ice masses) may be governed by ice deformation mechanisms that are much close to being linear in nature, n ¼ 122 (e.g., Doake and Wolff, 1985; Peltier et al., 2000). Given the gentle topography of interior North America, a large area of the Laurentide Ice Sheet was likely in a low shear-stress state at LGM, relative to most modern analogs (e.g., the Greenland Ice Sheet, alpine glacier masses). The prospect of lower deformation-rate exponents in (1) can be entertained in a model, and we take n and B0 to be the two primary free parameters to explore in the sensitivity experiments. 2.1.2. Basal flow parameterization Basal flow F sliding and sediment deformation at the base of the ice sheet F is an interesting and poorly understood aspect of glacier dynamics. Because basal flow is governed by complex subglacial processes (i.e. basal hydrology, geological and topographic nuances), these dynamics are highly underdeveloped in glaciological models (Marshall et al., 2000). We believe this to be the largest glaciological source of uncertainty in NAIS reconstructions at LGM. Numerical experiments have demonstrated that the extent of basal motion has a large

influence on ice volume and ice sheet form (Fisher et al., 1985; Clark et al., 1996; Jenson et al., 1996). Our numerical experiments do not address the complexities of basal flow dynamics; we do not attempt to simulate basal hydrological evolution, sediment deformation mechanics, or ice stream dynamics. We are therefore far from sampling the full range of physical possibilities. Sensitivity experiments presented here make a preliminary investigation of the role of basal flow, prescribing sliding rates in proportion to local shear-stress (e.g., Payne, 1998). Basal flow is also subject to thermal and geological enabling conditions, aT ðl; y; tÞ and aG ðl; yÞ; giving the basal velocity parameterization vbj ðl; y; tÞ ¼ A0 aT ðl; y; tÞaG ðl; yÞrI gH

@hs ; @xj

ð4Þ

where Hðl; y; tÞ is the local ice sheet thickness and hs ðl; y; tÞ is the surface topography. Parameter A0 is a free variable that is used to increase the magnitude (maximum velocity) of basal flow in the experiments. We treat the thermal enabling function, aT ðl; y; tÞ; as a binary switch. Basal flow is activated when the ice sheet is at the pressure-melting point at the base, subject to the additional constraint that basal meltwater is being generated (excess energy in the local basal energy balance). The thermal switch then takes the form 8 > < 1; ice is warm-based; with local basal aT ¼ meltwater available; > : 0; ice is cold-based or no meltwater is available: ð5Þ The geologic enabling function aG ðl; yÞ is superimposed on the thermal enabling function as a time-invariant control. Basal flow is scaled to be in proportion to the fraction of deformable sediment coverage in each model grid cell (Marshall et al., 1996). Parameter aG A½0; 1 is input as an initial field, based on modern-day terrain characterization, and we assume it to remain constant for each point in the model throughout the glacial cycle. This geologic enabling condition for basal flow is turned off in some model simulations, so we treat it as a second free parameter in ð4Þ: We also experiment with the thermal control, aT ; through experiments with and without latent heat effects (thermal inertia from locally generated basal meltwater) in the basal energy balance. These experiments are ad hoc in that we do not have a realistic subglacial hydrological model to account for the transport and redistribution of basal meltwater that would certainly have occurred. Basal meltwater is assumed to pool locally, building up a free water layer at the ice–bed interface. This water layer requires a sustained energy deficit (net heat loss) at the bed to refreeze, therefore

178

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

acting to maintain warm-based conditions once they develop. Active marginal ice that is vigorously thinning refreezes to the bed more quickly without this thermal buffer. We wish to emphasize an important shortcoming of this sliding parameterization. It suggests that basal velocity increases in proportion to gravitational driving stress, linearly in this case. However, there is no clear evidence for a proportional relationship between basal sliding and driving (shear) stress with regard to largescale basal flow. In fact, the fast-flowing ice streams of West Antarctica exhibit the opposite relationship; gravitational driving stress is low in West Antarctic ice streams (e.g., Hughes, 1977; MacAyeal, 1989; Bindschadler and Vornberger, 1990). Proper exploration of basal flow effects on NAIS thickness and form should take this into account. This limits the strength of conclusions drawn from our basal flow experiments. Marshall et al. (2000) discuss a number of other glaciological process parameterizations that are either poorly understood or poorly portrayed in ice sheet models. This includes iceberg calving and grounding line dynamics at marine margins, the process of ice sheet occupation of marine channels and embayments (i.e. Hudson Bay), the numerical simulation of basal ice temperatures, representation of longitudinal stress and strain in 3D ice sheet models, the dynamics of subglacial hydrology and sediment deformation processes, and the representation of discrete fast-flow elements (ice streams, surge lobes) in continent-scale models. We will not return to these important processes here, but note that they represent the current frontiers in ice sheet modelling, with consequent uncertainties in our conclusions. 2.2. Mass balance model

ð7Þ

where Tc is a threshold temperature, below which precipitation falls as snow. We take this threshold to ! be 11C after Johannesson et al. (1995). Temperature T is assumed to be normally distributed about the monthly mean, Tm ; with standard deviation s: Net annual snowfall is then aðl; y; tÞ ¼

12 X m¼1

am ¼

12 X m¼1

fm Pm ðl; y; tÞ

rI ; rw

ð8Þ

where the final factor, the ratio of densities of ice and water, converts from m of water to m of ice equivalent. Snow and ice ablations are parameterized as a function of positive degree days,

 Z N tm ðT  Tm Þ2 Dm ¼ pffiffiffiffiffiffi T exp dT; ð9Þ 2s2 s 2p 0 where tm is the length of the month in days. We sum the monthly P positive degree days to give the annual total, D ¼ 12 m¼1 Dm : The heat energy associated with the net annual D goes first towards melting of new snowfall, aðl; y; tÞ; ms ðl; y; tÞ ¼ max½ds Dðl; y; tÞð1  rf Þ; aðl; y; tÞ ;

ð10Þ

where ds is the degree-day melt factor for snow and rf is the refreezing factor, which accounts for superimposed ice formation (refrozen meltwater). Surplus D; after melting off this year’s snowfall, is used for ablation of superimposed and glacial ice: ( ms dI ½Dðl; y; tÞ  ds ð1r ð1  rf Þ; ms > a; fÞ mI ðl; y; tÞ ¼ 0; ms pa: ð11Þ

We parameterize ice sheet mass balance using monthly degree-day methodology (Braithwaite, 1984, 1995). Mass balance is calculated from the net accumulation minus ablation of ice, bðl; y; tÞ ¼ aðl; y; tÞ  mðl; y; tÞ;

to fall as snow from

 Z Tc 1 ðT  Tm Þ2 dT; exp fm ¼ pffiffiffiffiffiffi 2s2 s 2p N

ð6Þ

where a and m are expressed in meters ice-equivalent. The ablation term m includes mass loss from surface, basal, and internal melt and from iceberg calving at marine margins. Calculation of the last three ablation terms is discussed in Marshall et al. (2000). We do not find LGM ice sheet reconstructions to be sensitive to the treatment of these terms, at least within the range of reasonable calving parameterizations, so we do not present sensitivities to these processes here. Calculation of snow and ice accumulation and surface melt is a much more sensitive matter. Given monthly precipitation and temperature information, Pm ðl; y; tÞ and Tm ðl; y; tÞ; we calculate the fraction of precipitation

Note that degree-day factors are distinct for snow and ice due to differences in radiative absorption (albedo). We identify five free parameters within this degreeday model, Tc ; rf ; s; ds ; and dI : These parameters have been examined at various field sites, and they do vary with latitude and/or region, e.g., North vs. South Greenland. It is not at all clear that the melt rate parameters, d; should extrapolate to lower latitudes, where summer insolation is much more intense. We therefore explore this parameter space in the tests described below. We do not investigate any other family of mass balance models, however, so our results on mass balance model sensitivity cannot be seen as exhaustive. Degree-day methodology is appealing for paleo-ice sheet reconstruction because the only climatic fields required are monthly (or annual) mean temperatures and precipitation rates. In reality, the actual processes and conditions that control snowfall and surface melt are much more meteorologically and energetically

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

complex. A full energy balance model, with knowledge of winds, relative humidity, cloud cover, etc., is needed to rigorously simulate the melt process. One can calculate the full energy balance in an atmospheric general circulation model (AGCM), but it is not yet computationally feasible to apply AGCMs to millennial-scale integrations. Simplified energy-moisture balance models or combined temperature-radiation melt models offer more physically based alternatives to degree-day models (e.g., Hock, 1999). While we do not experiment with alternative treatments here, we emphasize the importance of this avenue for increased confidence in mass balance modelling. 2.3. Isostatic model Ice sheet simulations reported here calculate the isostatic response to surface loading using an Earth model with an elastic lithosphere and Maxwell viscoelastic mantle, following standard methods discussed by Wu and Peltier (1982) and Peltier (1985). The equations of mass and momentum conservation and Poisson’s equation for gravitational potential are Laplace transformed, yielding an equivalent elastic problem. The elastic problem is solved by numerical integration of six coupled first-order ordinary differential equations obtained from scalar decomposition of the Laplacetransformed conservation equations. Specifically, the response to an impulsive point load is determined, and then inverse Laplace transformed to give the timedomain response. For application to the glaciological modelling described here, the surface response to a finite-radius spherical cap (disk) load with a Heaviside history is derived from the point load impulsive response using methods described by James and Ivins (1998). The area of the disk is chosen to agree with the nominal grid spacing of the glaciological model (55 km by 55 km), although response functions are scaled to match the exact load area Aðl; yÞ over the variable-area grid. The Heaviside disk response depends on the distance from the disk center and on the time elapsed after the load is imposed. We calculate response functions at 500- or 1000-yr intervals, ranging from the instantaneous elastic response to a time lag of 150 kyr; giving the viscoelastic Green’s functions that are convolved at each point in the model grid for each isostatic model update. Response functions were derived for compressible Maxwell viscoelastic Earth models with a specified density and elastic structure. The thickness of the lithosphere and viscosity layering of the mantle vary from simulation to simulation. For sensitivity experiments, we simplify this to three free variables in the Earth model: lithosphere thickness, upper mantle viscosity, and lower mantle viscosity. Some simulations in the sensitivity experiments also make use of a simple isostatic treatment that assumes

179

a damped, local viscous relaxation to surface loads (Deblonde and Peltier, 1993; Tarasov and Peltier, 1999),

 @hb hb  hb0 rI H þ rw dHw ¼ þ ; ð12Þ @t t rb t where t is the exponential response time of the underlying mantle, hb is the bed topography, hb0 is the initial (isostatically equilibrated) bed topography, and rI ; rw ; and rb are the densities of ice, water, and shallow Earth crust, respectively. Note that the surface forcing, both in (12) and for the Maxwell-Earth load history, includes both the ice load, Hðl; y; tÞ and any changes in waterlayer thickness associated with sea level fluctuations at marine points, dHw ðl; y; tÞ: 2.4. Paleoclimate forcing Paleoclimate perturbations are prescribed using the Greenland Ice Core Project (GRIP) d18 O record (Dansgaard et al., 1993) in conjunction with AGCM results from the Canadian Centre for Climate Modelling and Analysis (CCCma) (McFarlane et al., 1992). AGCM temperature and precipitation reconstructions for present day and LGM (Vettoretti et al., 2000) are linearly combined in a GRIP-based interpolation of modern and glacial climatologies. This parameterization is based on a ‘‘glacial index’’ IðtÞ which varies between 0 and 1 for interglacial (modern) and full glacial (LGM) conditions. Assumptions and limitations in this simple paleoclimatic forcing are discussed in Marshall et al. (2000). Paleoclimate fields Pðl; y; tÞ and Tðl; y; tÞ are based on perturbations to modern-day climatology, Pðl; y; tÞ ¼ fP ðl; y; tÞPobs ðl; yÞ;

ð13Þ

Tðl; y; tÞ ¼ Tobs ðl; yÞ þ DTðl; y; tÞ  bDhs ðl; y; tÞ;

ð14Þ

where fP is the fractional precipitation scaling and DT is the climatic air temperature perturbation. The second term in (14) incorporates atmospheric lapse-rate effects associated with the change in surface elevation, Dhs ; where b is the lapse rate. The paleo-precipitation scaling fP is prescribed based on the harmonic combination of LGM and present-day CCC fields,



IðtÞ fP ðl; y; tÞ ¼ ðP  EÞCCC ðl; y; 0Þ ðP  EÞCCC ðl; y; 21Þ 1 1  IðtÞ ; ð15Þ þ ðP  EÞCCC ðl; y; 0Þ where ðP  EÞCCC denotes the AGCM precipitation less evaporation (available moisture) at each point, t ¼ 0 refers to the modern control run, and t ¼ 21 denotes the LGM snapshot, 21 kyr BP.

180

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

This default model treatment of precipitation generates spatial variations as predicted by the climate model, with glacial conditions that are wetter in some locations (e.g., the southern margins) and drier in others. We examine model sensitivity to this treatment through two means: (a) truncation of the maximum possible fractional increase in wetness at LGM (I ¼ 1), through fP ¼ max½Eq:ð15Þ; fPmax ;

fPmax A½1:0; 2:0 ;

ð16Þ

and (b) elevation-drying effects when the surface elevation increases above a threshold altitude, ht ; following Pðl; y; tÞ ¼ P0 ðl; y; tÞ exp½bP max½hs ðl; y; tÞ  ht ; 0 : ð17Þ The variable P0 in this expression refers to the precipitation calculated from Eq. (13), while bP is an elevation lapse rate for precipitation. This parameterization reduces paleo-precipitation exponentially above the elevation-desert threshold ht : For completeness here, note that we have also experimented with spatial precipitation patterns that are based on present-day distributions. Under this treatment, local precipitation rates diminish exponentially with local atmospheric cooling, reflecting the increased aridity that can be expected under glacial conditions (Tarasov and Peltier, 1999). Paleo-precipitation under this parameterization has the form Pðl; y; tÞ ¼ Pobs ðl; yÞð1 þ dP ÞDTðl;y;tÞ
exp½bP max½hs ðl; y; tÞ  ht ; 0 :

ð18Þ

The parameter dP in this equation represents the percentage of drying per 1C; Tarasov and Peltier ð1999Þ choose a value of 3% per 1C; dP ¼ 0:03: All told then, our paleo-precipitation experiments involved four free parameters: fPmax ; bP ; ht ; and dP : Paleo-temperatures over the ice sheet are calculated from (14) with spatially variable perturbations DTðl; y; tÞ based on a linear weighting of glacial and interglacial temperatures from the CCC model, DTðl; y; tÞ ¼ fT IðtÞ½TCCC ðl; y; 21Þ  TCCC ðl; y; 0Þ :

ð19Þ

The coefficient fT has been introduced to add a free parameter to this temperature forcing; it allows us to experiment with fractional proportions of the CCCbased temperature perturbation (e.g., 80% of the predicted LGM cooling). 2.5. Numerical model A conventional finite difference model is used for the ice sheet dynamics, with a grid resolution of 11 longitude (15–75 km) and 0:51 latitude (55 km) over North America. The model grid does not include Greenland. There are 20 vertical layers in the model, with a stretched vertical grid giving vertical resolutions from

ca. 10–300 m over the range of ice sheet thicknesses. The glacial cycle integration extends from 120 kyr BP to present, with model time steps of 2–5 years for the ice dynamics and mass balance solution, 5–20 years for the ice thermodynamics update, 100 years for the climate updates, and 500–1000 years for the isostatic model solution. Bed topography is taken from TerrainBase DEM (Row and Hastings, 1994). Modern mean monthly air temperatures are interpolated from 14-year averages (1982–1995) of the NOAA NCEP-NCAR CDAS-1 2-m temperature reanalysis dataset (Kalnay et al., 1996), modern precipitation rates are from Legates and Willmott (1990), and the geothermal heat map for North America is taken from Blackwell and Steele ð1992Þ: 2.6. Model experiments We carried out a series of 190 simulations of the last glacial cycle in North America, essentially sieving for combinations of parameters that give ice sheet reconstructions in accord with geological reconstructions of LGM ice sheet area. Simulations run from 120 kyr BP to the present, beginning with no initial ice, and we define LGM as 21 kyr BP. Fig. 2 depicts a typical glacial cycle simulation in North America, from a case that comes close to the mean predicted LGM ice volume for the suite of ‘‘reasonable’’ simulations (see below). Fig. 2(a) illustrates the mean air temperature perturbation over the ice sheet, % ¼ DTðtÞ

5

1 SSðtÞ Aðl; yÞ½Tðl; y; tÞ  Tðl; y; 0Þ ; SðtÞ

ð20Þ

(a)

Mean Air Temperature over the Ice Sheet ( C)

0 -5 -10 -15 120 35 30 25 20 15 10 5 0 120

100

80

60

40

(b)

Ice Volume (1015 m3)

100

80

0

20

60

40

20

0

Time (kyr BP)

Fig. 2. Time series of the last glacial cycle integration: (a) Mean air temperature perturbation over the ice sheet, 3 C; and (b) NAIS ice volume, 1015 m3 : Results from the median case in the suite of acceptable results (cf. Fig. 4).

181

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

where Aðl; yÞ is the area of a 0:51 by 11 grid cell in the model. The sum over SðtÞ refers to an integration/ summation over all model grid cells that contribute to the total area of the ice sheet at time t: For latitude yj ; longitude li ; and Earth radius RE ; cell area is calculated from Z yjþ1=2 Z liþ1=2 R2E sin y dl dy Aðli ; yj Þ ¼ yj1=2

li1=2

¼ R2E Dlðcos yjþ1=2

 cos yj1=2 Þ ð21Þ P P and total ice sheet area SðtÞ ¼ i j Aðli ; yj Þ: The air temperature time series in Fig. 2(a) should be thought of as a model forcing rather than a climate model prediction. It is the result of the CCC/GRIP paleoclimatology (a pure forcing) and the ice sheet surface elevation evolution (a model feedback) for a given model realization. Fig. 2(b) plots the total NAIS volume VðtÞ ¼ SSðtÞ Hðl; y; tÞAðl; yÞ:

ð22Þ

Units are m3 of ice; conversion of ice volume to eustatic sea level is not straightforward and is discussed below.

The range of parameter space explored in the full suite of numerical experiments is summarized in Table 1. The 190 simulations focus on sensitivity to model treatment of ice sheet rheology, bed isostatic treatment, basal flow, mass balance parameterization, and climatic forcing. Parameters have been varied one-at-a-time in approximately two-thirds of the experiments, holding all other parameters at the default value. In the remaining experiments, parameters were varied simultaneously in a preliminary exploration of the multivariate parameter space. Simulated LGM ice extent in North America is very sensitive to elevation-desert effects (aridity in interior portions of the ice sheet), and the joint sensitivity experiments were centered on two focal treatments of the elevation-drying threshold: ht ¼ 1200 m and 2000 m:

3. Results Fig. 3 plots the distributions of predicted ice sheet and paleoclimate diagnostics for LGM snapshots from the full simulation suite (N ¼ 190). A summary of the

Table 1 Parameter range for glacial cycle simulations Parameter

Default

Range

Ice rheology experiments Flow law exponent B0 ð103 CÞ; n ¼ 3 B0 ð103 CÞ; n ¼ 2:5 B0 ð103 CÞ; n ¼ 1:7

3 4:2
1017 Pa3 yr1 9:5
1013 Pa2:5 yr1 7:8
1011 Pa1:7 yr1

1.7–3 1.4–14
1017 Pa3 yr1 1.9–47
1013 Pa2:5 yr1 3.9–78
1011 Pa1:7 yr1

Isostatic model experiments Local model, relaxation time t Upper mantle viscosity, mu Lower mantle viscosity, ml Lithosphere thickness Maxwell Earth model updates

4000 yr 1021 Pa s 4
1021 Pa s 120 km 500 yr

100–20; 000 yr 5–10
1020 Pa s 2–8
1021 Pa s 70–170 km 100–1000 yr

0:008 m yr1 Pa1

0–0:01 m yr1 Pa1 Pure aT regulation aG ; aT controls Basal melt/latent heat

Mass balance experiments Maximum refreezing, fr ðmaxÞ Degree day melt factor for snow, ds Degree day melt factor for ice, dI Degree day model, s Elevation-aridity threshold, ht Accumulation lapse rate, bP

0.6 3 mm yr13 C1 8 mm yr13 C1 53 C 1500 1 m yr1 km1

0.3–0.8 2–4 mm yr13 C1 7–9 mm yr13 C1 3–83 C 1000–2000 m 0.5–1:2 m yr1 km1

Climate experiments Atmospheric lapse rate, b Maximum LGM precipitation ratio, fPmax Maximum LGM DT fraction, fT Temperature aridity, dP

73 C km1 2 0.7 0.02

5–103 C km1 1–2 0.6–0.9 0–0.03

Sliding experiments Sliding law coefficient, A0 Basal regulation (switch)

182

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

Fig. 3. Last Glacial Maximum (21 kyr BP) diagnostics from the full set of glacial cycle simulations: (a) NAIS volume, 1015 m3 ; (b) NAIS area, 1012 m2 ; (c) Basal melt percent, (d) Mean mass balance over the ice sheet, m yr1 ice-equivalent, (e) Mean precipitation rate over the ice sheet, m yr1 water, and (f) Mean LGM air temperature depression over the ice sheet, 3 C: All histograms plot the number of occurrences of each class, for a total sample size of N ¼ 190:

Table 2 LGM (21 ka BP) summary statistics, North American ice sheets. Full simulation suite (N ¼ 190)

15

3

Ice volume (10 m ) Ice area (1012 m2 ) Basal melt fraction Mass balance (m yr1 ice) Mean precipitation (m yr1 ) Mean DT ð3 CÞ Southern limit (3 N)

Mean

Std. Deviation

Minimum

Maximum

37.7 14.3 0.37 0.26 0.93 11:6 40.2

9.5 2.0 0.12 0.07 0.09 1.7 2.1

15.2 9.4 0.02 0.15 0.62 16:6 38.0

56.7 18.2 0.75 0.60 1.09 7:3 46.25

means and ranges of these diagnostics is given in Table 2. There are two main modes to the predicted LGM NAIS volume, corresponding to the two focal elevation-drying thresholds. Climate and mass balance models have been perturbed from standard or previous treatments (e.g., Marshall et al., 2000) so have already been ‘‘pre-tuned’’ to a large degree; the default parameters give a complete 120-kyr ice age cycle in North

America, with the area and southern extent of the LGM ice sheets predicted reasonably well. The scatter and uncertainty in model experiments is very evident from Fig. 3a, however. Dramatically different glacial cycles are realized over the range of sensitivity experiments. Fig. 3b plots the areal extent of the LGM ice sheets in North America for the full set of numerical experiments and Fig. 3c illustrates the range of predictions for the

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

area of the ice sheet that is at the pressure-melting point at the base. The basal temperature field is a complicated function of the advective regime in the ice (flow law, precipitation rates, extent of basal flow), the long-term ð104 -year) air temperature forcing, the thickness of the ice, and the residence time of ice at a site (new, thin ice cover is typically cold-based). The fraction of the base that is at the basal melting point is of interest for meltwater production and basal flow regulation, and Fig. 3c indicates the general uncertainty in this important variable. This percentage increases throughout the glacial period, peaking during late deglaciation stages; at LGM, however, results indicate a mode of only 30–40% of the base able to support basal flow. Figs. 3d–f plot the distributions of mean mass balance and climate model diagnostics for the experiments. All three variables are averaged over all North American grid cells containing ice at LGM (hence, a different area for each simulation). Mean mass balance in Fig. 3d depicts the mean accumulation minus mean ablation (including surface melt, basal melt, and calving) for the entire ice sheet, % ¼ 1 SSðtÞ Aðl; yÞbðy; l; tÞ: bðtÞ SðtÞ

ð23Þ

Figs. 3e and f illustrate the range in paleoclimatic forcing for the experiments. Mean air temperature perturbations are calculated over the ice sheet area and are equivalent to the temperature forcing time series shown in Fig. 1a. This field has a strong contribution from ice sheet elevation; the coldest cases correspond to those with thick ice. Mean precipitation rate in Fig. 3e is expressed in m/yr water-equivalent and is the mean annual value over the ice sheet. This varies between model experiments under the control of the CCC/GRIP forcing (LGM/present maximum precipitation ratio) and the strength of the elevation-precipitation lapse rate parameterization (threshold elevations and the lapse rate itself). The CCC spatial patterns of air temperature and precipitation fields at LGM are presented in Vettoretti et al. (2000) and Marshall et al. (2000). This AGCM-based spatial pattern of glacial precipitation rates is implicit in all of the model runs that use the CCC/GRIP forcing. Other runs use the modern spatial distribution of precipitation, reducing net rain/snowfall on a local basis in proportion to local cooling (Eq. (18)). The driest cases in Fig. 3e are associated with these models. 3.1. Model assessment There is dismayingly little paleoclimatologic information available to evaluate and constrain the model. A mean cooling of 10–151C over the ice sheet area in North America is consistent with the ca. 201C cooling in

183

central Greenland at LGM (Cuffey and Clow, 1997; Dahl-Jensen et al., 1998) and with development of periglacial features in the Laurentide Ice Sheet (LIS) forefield (e.g., Johnson, 1990). Paleohydrological evidence supports the likelihood that regions south of the LIS were more moist at LGM, probably a result of both (a) reduced evaporation and (b) increased precipitation, associated with jet stream splitting under the orographic influence of the NAIS. These paleoclimatological inferences can be used for consistency checks but do not have the quantitative footing to give them discriminating power in the model experiments. Geological and geophysical controls therefore represent the primary means of evaluating model results. We apply the geologic record of NAIS area and southern extent (Dyke and Prest, 1987) as a control on the glaciological reconstructions. LGM ice sheet area in North America, including Greenland, is estimated to be 15:98
106 km2 (A.S. Dyke, personal communication, 2001). Separate simulations of the Greenland Ice Sheet through the last glacial cycle give a Greenland Ice Sheet area of 2:61
106 km2 at 21 kyr BP (Marshall and Cuffey, 2000). Greenland is excluded from the NAIS simulations in this study, so we subtract this to give an estimate of 13:37
106 km2 for the remaining North American ice at LGM. The map sheet reconstructions of Dyke and Prest (1987) also give exceptional control on the LGM southern margin position as a function of longitude, e.g., 41:51N at 761W; etc. We include this control in a simple way by examining the maximum LGM southern extent (minimum LGM latitude) and the longitudinal position of this southern foray. To apply these constraints, we sieve the full suite of model experiments for those simulations with (a) LGM NAIS areas within 1
106 km2 of the expected geologic value, 13:4
106 km2 ; and (b) maximum southeastern advances beyond 42:51N: Fig. 4 illustrates the LGM diagnostics from the subset of 33 simulations that result from this filtering. Results embodied in Fig. 4 will be discussed extensively in Section 3.3. Fig. 5 plots the LGM surface topography from four different NAIS realizations. The first of these, Fig. 5a, is a result that was rejected based on the area and southern extent criteria. This is a case with a relatively dry glacial climate. It is similar in appearance to many cases that ‘‘came up short,’’ due to low basal flow, not enough moisture or temperature depression in the climate parameterization, or insufficient ice flux in the thin ice sheets that result from low-exponent ice rheology in our tests. The LGM ice sheets in Figs. 5b–d were all accepted. These illustrate a general sample of ice divide structures and ice sheet thicknesses in our experiments. We typically find multiple Laurentide Ice Sheet (LIS) domes, with divides evident to the east, south, and west of Hudson Bay.

184

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

Fig. 4. Simulated Last Glacial Maximum (21 kyr BP) diagnostics from the set of acceptable glacial cycle simulations (acceptability based on LGM area and southern ice extent): (a) NAIS volume, 1015 m3 ; (b) NAIS area, 1012 m2 ; (c) Basal melt percent, (d) Mean mass balance over the ice sheet, m yr1 ice-equivalent, (e) Mean precipitation rate over the ice sheet, m yr1 water, and (f) Mean LGM air temperature depression over the ice sheet, 3 C: All histograms plot the number of occurrences of each class, for a total sample size of N ¼ 33:

Cordilleran and Innuitian ice complexes are confluent with the LIS at LGM, but they maintain distinct divide structures. 3.2. Additional model constraints While the ice sheets in Figs. 5b–d were accepted based on the LGM area and southern extent, additional constraints could probably be applied to these reconstructions. There is almost certainly too much ice west of Hudson Bay in Figs. 5b and d, for instance, where the model suggests a dominant Keewatin ice divide and thick ice at the confluence of Laurentide and Cordilleran ice. Ice advance into the western interior plains is thought to have been late, ca. 30–20 kyr BP (Clark et al., 1993; Jackson et al., 1997); thick ice in the southwestern plains predicted by many of the model simulations is in conflict with this. The late-glacial flow indicators described by Shilts (1980) may offer support for a major Keewatin ice dome, but the results of Veillette et al. (1999) indicate that the main ice dispersal center was east of Hudson Bay through much of the glacial period. Examination of earlier snapshots in the glacial cycle will

be useful to further the model evaluation, wherever firm geological insights are available. The LGM constraints that we do apply could also be tightened; NAIS area at LGM is geologically known to better than 1
106 km2 ; for instance, and the detailed NAIS margin structure could be used for model discrimination. Unfortunately, the glaciological model is not sophisticated enough for this as yet. Southern lobe structure is poorly captured in the model, as seen in Fig. 5; margin positions are diffuse and uniform relative to the Dyke and Prest (1987) reconstruction. The complicated southern margin structure as geologically inferred was almost certainly a result of glaciological processes (subglacial sediment deformation and subglacial hydrological regulation) that are neglected in the ice sheet model. Improved physical treatments of basal flow controls and a physical description of the longitudinal stress/strain regimes appropriate to thin, fastflowing lobes (Clark, 1992, 1994) persist as high priorities for improved glaciological models (Payne, 1998; Marshall et al., 2000). In addition to areal and geographical constraints, some model discrimination is possible based on LGM

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

185

Fig. 5. LGM (21 kyr BP) ice sheet reconstructions in North America for four sample simulations. (a) Cold, dry case with enhanced sliding. This case failed to meet the acceptance criteria (area, southern extent too limited at LGM). (b–d) Example LGM reconstructions under different climate model configurations (see text for a detailed description of parameter settings for these cases). Cases (b–d) were all accepted based on LGM ice area and southern extent, although detailed considerations of divide structure can be invoked to argue against cases (b) and (d), while no model simulations adequately represent LGM southern lobe structure.

NAIS volume. Natural brackets can be imposed from the global sea-level record; grounded NAIS volume must account for less than the total eustatic sea-level lowering of 115–135 m at LGM (Bard et al., 1990; Fleming et al., 1998; Yokoyama et al., 2000). Beyond this we have not introduced a volume constraint, as we wish to avoid a predisposition and reserve ice volume as an objective result of the analysis. Existing estimates of LGM ice volume in North America vary from ca. 33% to 90% of this total eustatic signal (Peltier, 1994; Clark et al., 1996; Marshall et al., 2000; discussion in Lambeck et al., 2000). While it is likely that North America harbored a large percentage of the global LGM ice, based on known areas of LGM global ice cover, NAIS volume must be considered a poorly constrained quantity.

Peltier (1994, 1996) has used the postglacial rebound record to develop convincing geophysical controls on late-glacial ice sheet thickness. It is possible, in principle, to use relative sea-level observations in conjunction with simulated isostatic rebound to help evaluate LGM ice volume predictions from the glaciological model. Our initial examination of this possibility, using relative sealevel data from northeastern North America, was inconclusive. This is because the detailed history of deglaciation is the primary control on modern and Holocene isostatic rebound rates; the model is relatively insensitive to the system state at 21 kyr BP (Lambeck et al., 2000). Inferences about LGM ice thickness are therefore weak in our simulations. Relative sea-level data certainly offer a potential window to the LGM ice

186

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

sheets, but insights require careful consideration and constraint on the ice sheet deglaciation. While this is well-known geologically (i.e., the areal history of ice retreat is well known), we do not have the glaciological/ paleoclimatic modelling skill to address this as yet.

3.3. Summary of acceptable results Table 3 summarizes the mean and range of model results from the ‘‘acceptable suite’’ of simulations shown in Fig. 4. There is a considerable spread in LGM ice volumes within the limits of reasonable LGM ice area and southern ice extent, 28.5–38.9
1015 m3 ice. The majority of acceptable tests cluster around 32–36
1015 m3 : LGM basal melt percent for this subset of results has a mean of 32%; this indicates that only one-third of the NAIS is able to support basal flow at this time, a surprisingly low figure. At this point in the glacial cycle, cold climatic conditions have prevailed for several tens of millennia and ice has not been in position for long enough to reach warm-based conditions in many regions. Subglacial melt fraction increases steadily and rapidly from 20 kyr BP until the final stages of deglaciation (rising to ca. 80% by 11 kyr BP). This basal temperature evolution may play a key role in facilitating late-glacial fast flow in the ice sheet, therefore operating as an important feedback in deglacial ice sheet collapse (Fisher et al., 1985; Clark et al., 1999). This is beyond the scope of the LGM focus in this paper. An exclusively positive LGM mass balance is associated with this suite of results, with a mode of b% ¼ 0:2–0:22 m=yr ice-equivalent. The mean climatologic state associated with this range of experiments is an air temperature depression DT ¼ 111C relative to modern conditions and a mean precipitation rate of 0:89 m yr1 water-equivalent, slightly drier than modern. The interior regions of the ice sheet are very dry but precipitation in the vicinity of the southern margin is

Table 3 LGM (21 ka BP) summary statistics, North American ice sheets. Range of acceptable results based on LGM area and southern limit (N ¼ 33) Mean Ice volume (1015 m3 ) 34.3 Ice area (1012 m2 ) 13.75 Basal melt fraction 0.32 Mass balance (m yr1 ice) 0.22 Mean precipitation (m yr1 ) 0.89 Mean DT ð3 CÞ 11.0 42.0 Southern limit (3 N)

Std. Minimum Maximum Deviation 2.3 0.44 0.06 0.05 0.08 0.5 0.2

28.5 12.9 0.08 0.16 0.70 12.3 41.25

38.9 14.4 0.42 0.39 1.02 10.1 42.5

enhanced, compensating for the interior dryness to a large degree. 3.4. Regional ice volumes and conversion to sea-level equivalent The raw ice volumes presented above are relevant to ‘‘isotopic’’ ice locked up on the continent; all of this ice plays a role in LGM d18 O increases in the ocean. This ice volume cannot be directly translated to sea-level equivalent, however. Floating ice does not contribute to LGM sea-level change; nor does that which displaces modern-day marine bodies such as the Hudson Bay. To convert from m3 of ice to sea-level equivalent, we postprocess the archived ice sheet thickness and bed topography to remove floating and marine-based contributions. The latter is accounted through consideration of the modern-day water column in marine channels and embayments that were ice-occupied at the LGM. Applying this correction, the sea-level contributing ice volume Vadj is reduced by about 1.5%; that is, roughly 0.5
1015 m3 of ice is floating or displacing modern marine waters at LGM (cases with well-developed LGM ice sheets). As an explicit example, a case that falls in the middle of the acceptable results (Fig. 4) has a total LGM ice volume of 34:32
1015 m3 ice, with 0:48
1015 m3 of floating/marine ice, giving Vadj ¼ 33:83
1015 m3 : This sea-level impacting ice is converted to global sealevel equivalent, hsl ; through

 r Vadj hsl ¼ w ; ð24Þ rI 0:7ð4pR2E Þ where the denominator is an estimate of the global ocean area. The result for this case is a eustatic sealevel drawdown of 82:7 m (msl) associated with the NAIS at LGM. Revisiting the ice volume summary in Table 3, we calculate a range of global sea-level predictions of 68.7–93:7 msl; with a mean of 82:7 msl and a dominant mode of 77.8–87:6 msl (76% of acceptable tests). It is interesting for comparison with other reconstructions to isolate the geographic contributions to this total. For the representative case above, we calculate individual contributions of 30:0
1015 m3 ice from the Laurentide Ice Sheet, 2:9
1015 m3 ice in the western Cordilleran, and 1:4
1015 m3 ice in the Innuitian Ice Sheet. This corresponds to 87% of the ice in the Laurentide Ice Sheet, 8.5% in the Cordillera, and 4.5% in the Innuitian complex. The floating/marine ice that does not contribute to global sea-level draw-down is concentrated in the Laurentide and Innuitian Ice Sheets, 0:38
1015 and 0:1
1015 m3 ice, respectively. With these corrections applied, the sea-level contributions from the Laurentide, Cordilleran, and Innuitian complexes are 72:4; 7:1; and 3:2 msl; respectively.

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

4. Discussion and interpretation 4.1. Reconciliation with geophysical inferences of low NAIS volume Results suggest a global sea-level impact of 78–88 m associated with the NAIS at the LGM. This estimate goes against the grain of recent geophysical and glaciological reconstructions that suggest ca. 30% less ice than this in North America. Clark et al. (1996) modelled a LIS that was subject to widespread basal sediment deformation at LGM, creating a WAIS-like ice sheet containing 19:7
1015 m3 of ice, or 49 msl; in contrast with our estimate of 72 msl: The ICE-4G reconstruction of Peltier (1994, 1996) suggests a LIS that contains 50–60 msl equivalent: again, significantly less ice than that suggested from this study. Previous studies related to this discussion (e.g., Marshall and Clarke, 1999; Marshall et al., 2000), have taken the ICE-4G reconstruction as an objective constraint or ‘‘target’’ with respect to LGM ice volume in North America. This leads to the conclusion that the glaciological model predicts excessive ice, presumably due to deficiencies in basal flow modelling and/or inappropriate ice rheology. We have retreated from this position here, and argue that a 78–88 msl impact is in fact possible from the NAIS. This can be reconciled with the relative sea-level evidence under the condition that the LGM ice sheets thin substantially in the several millennia following the LGM, prior to large-scale areal retreat (Lambeck et al., 2000). This point, discussed earlier, results from the fact that relative sea-level data speak primarily to the deglaciation period subsequent to LGM. This is particularly true for the LIS, as the late evacuation of ice from northeastern Canada does not allow for uplift data extending beyond the Holocene (e.g., 6- to 7-kyr records from Hudson Bay and Ungava Bay). The ice thickness at 21 kyr BP is not critical to these uplift records, relative to the ice thickness history in the period 21–8 kyr BP (Lambeck et al., 2000; preliminary analyses in our simulations, not presented here). Hence, the relative sea-level data that are available are not definitive with respect to LGM ice thickness in North America. Furthermore, data relevant to the LIS are confined to coastal regions in eastern North America and Labrador/Hudson Bay, providing little constraint on ice thickness in interior regions. We therefore argue that geophysical reconstructions permit this ‘‘restoration’’ of ice in North America at LGM. This creates a new problem of course; the isostatic record demands substantial ice thinning subsequent to LGM, at a time (14–20 kyr BP) when there is no strong geologic signal of ice sheet retreat (Dyke and Prest, 1987). Model results suggest that this is possible via the increasing role of fast basal flow in this period, as more of the ice sheet

187

becomes warm-based (Fig. 3c) and basal meltwater accumulates at the bed. This essentially argues for a transition from thicker, largely cold-based ice sheets at LGM to a thin and mobile, more WAIS-like ice sheet through the deglaciation (Clark et al., 1999). 4.2. Sensitivity experiments Discrimination of results based on LGM ice area and southern extent gives some suggestion as to the range of ice volume and climatic states associated with glaciological reconstruction of the North American ice sheets. The subset of acceptable model results includes a wide spectrum of parameter settings. Following is a brief synopsis of the phase space of acceptable results, as summarized in Fig. 4 and Table 3. 4.2.1. Climate forcing Cases in Fig. 4 cluster around a CCC/GRIP climate forcing with air temperature perturbations of 67–70%, with the range from 62% to 78% (fT ¼ 0:62–0.78 in Eq. (19)). Predicted LGM ice sheet disposition is not sensitive to the precipitation thresholds within the CCC/ GRIP climatology. None of the cases with modern spatial patterns of precipitation were accepted, however. Ice sheets typically fail to advance far enough south without the high-precipitation belt that results from the polar jet stream diversion and the orographic effects of the advancing southern margin. 4.2.2. Mass balance model We found the main climate/mass-balance sensitivity to be the strength and elevation threshold of the precipitation lapse rate. With exponential drying (Eq. (17)) beginning 2000 m above current topography, LGM ice volume was consistently 10–12
1015 m3 (ca. 30%) higher than in cases with ht ¼ 1000 m; with a concomitant dynamical expansion of ice area. Most cases that have been classified as acceptable have ht ¼ 1000–1200 m: The nonlinearity in this result is interesting; this parameter emerged as the single most sensitive parameter in our simulations. This sensitivity is large because lowering of the drying threshold acts to suppress ice thickness in interior regions, without discouraging new ice growth in low-lying or marginal areas. Cases with increased interior ice thickness lead to excessive ice flow to the west and south, giving areal expansion of ice beyond the cutoff used for model evaluation in this paper. The contrasting ice sheet evolution in these two families of simulations, with ht ¼ 1200 m and 2000 m; is graphically represented in Figs. 6 and 7, respectively. Fig. 6 plots ice surface topography at 40 kyr BP in two cases, showing the early penetration and growth of thick ice cover in the western interior plains for ht ¼ 2000 m (Fig. 6b). This evolution is at odds with geologic

188

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

Fig. 6. Modelled ice sheet surface topography at 40 kyr BP: (a) A case with relatively low glacial precipitation, suppressing ice expansion into the western plains of North America until late in the glacial cycle. (b) A case with higher glacial precipitation (a less severe elevation-desert effect), resulting in early ice occupation in the western plains and excessive western ice at LGM (cf. Fig. 5d).

Fig. 7. Duration of ice occupation in each model grid cell for the 120-kyr glacial cycle simulation, kyr. (a) A case with relatively low glacial precipitation, suppressing ice expansion into the western plains of North America until late in the glacial cycle. (b) A case with higher glacial precipitation (a less severe elevation-desert effect), resulting in early ice occupation in the western plains and excessive western ice at LGM (cf. Fig. 6d).

expectation, as discussed in Section 3.2, and it gives rise to LIS cover at LGM that is thick and too far south in the continental interior. The higher precipitation rates in the ice sheet interior in the ht ¼ 2000 m case therefore lead to too much LGM ice in southern and western sectors of the LIS, excess Cordilleran ice (Clague, 1989; Clark et al., 1993), and a premature confluence of Laurentide and Cordilleran ice (Jackson et al., 1997). This contrast and the sensitivity to precipitation– elevation threshold is also seen in Fig. 7. Cumulative occupation time of ice during the last glacial cycle is plotted at each point in the model grid, again for the two

families ht ¼ 1200 m and 2000 m: Fig. 7 illustrates the ice sheet inception and spreading centers that are typical in model simulations (blue areas, which experience ice occupation for more than 100 kyr). These areas are similar in each case, but ice sheet history diverges in the southwestern interior plains and in the northwestern interior, where the dry-climate case shown in Fig. 7a does not see western LIS expansion until late in the glaciation. We do not have a good feel for the elevation-drying that actually took place over the mature ice sheets in North America, if any. However, these results suggest

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

that this is an extremely important factor in ice sheet mass balance and thickness. Synoptic circulation and precipitation patterns during the glacial period play a pivotal role in dictating LGM ice extent. This is particularly true in western North America, where LIS expansion and Cordilleran ice sheet growth appear to be largely precipitation-controlled. In fact, all of our model simulations overpredict pre-LGM Cordilleran ice (Clague, 1989; Clark et al., 1993), suggesting that the GCM climatology is missing an important paleoclimatic (paleoprecipitation) shift in western North America. This should be a priority for future paleoclimatic reconstructions concerning North American glaciations. 4.2.3. Ice rheology The suite of acceptable results included 28 simulations with a Glen-like n ¼ 3 ice rheology and five simulations with n ¼ 2:5: For the n ¼ 3 cases, interestingly, both high and low values of B0 (very fluid and very stiff ice sheets) failed to produce acceptable results, due to insufficient ice area and southern ice extent. This was also true of the n ¼ 1:7 tests. Ice sheets that are too thin and fluid have broader ablation zones and are not necessarily associated with greater ice flux; with low gravitational driving stress, shear-stress driven flow is much reduced. We have not yet found a positive result for the high-enhancement or low-exponent rheologies, although there is likely to be some combination of mass balance and climate model parameters that would yield this. Within the acceptable n ¼ 2:5 and 3 cases, ice sheets with the low flow-law exponent were thinner, with less LGM ice volume, as one would expect. This is also true for increasing enhancements. Beyond a certain threshold in each rheologic treatment, the ice sheet thins to the point where the Laurentide Ice Sheet does not extend far enough south and west at LGM, resulting in LGM areas and geometries that are not representative of the geologic record. It is unclear whether this decreased ice mobility can be climatically overcome. Our most serious concern for the gaunt ice sheets that result from large enhancement factors and low stress exponents is associated with the modelled basal temperature fields. Thin ice flowing in a low-stress regime is extremely cold-based, due to (a) high diffusive heat loss from the bed in thin ice (less insulation), (b) strong advective cooling (both vertical and horizontal), and (c) low strain heating. The latter is especially pronounced in low-exponent flow laws, as internal deformational heating varies in proportion to snþ1 : All of these thermodynamic influences promote cold-based ice, such that less than 10% of the LGM ice sheets in North America are warm-based in the experiments with n ¼ 1:7 (with most of this warm ice in the Cordillera). This is a serious concern with respect to the ample geological

189

record of warm-based environments in the last glaciation. It also prohibits basal meltwater production and basal flow, which further limits the ability of these models to expand southward. 4.2.4. Sliding law parameterization Fisher et al. (1985) and Boulton et al. (1985) have pointed out the potential importance of basal flow in dictating NAIS form and thickness. Model studies reported in Clark et al. (1996) and Jenson et al. (1996) illustrate and emphasize this, with LIS volume reductions of 30–40% under the influence of widespread basal flow. The importance of this control on ice sheet thickness is also manifest in the contemporary West Antarctic Ice Sheet, where fast-flowing outlets have created a thin, low-sloping ice mass (e.g., MacAyeal, 1989; Bindschadler and Vornberger, 1990). Some contribution from basal flow appears to be necessary to bring the ice far enough south in our simulations. No tests with A0 o0:008 m yr1 Pa1 were accepted. The thermal vs. geologic controls on basal flow enabling impacted ice sheet divide structure and transient evolution to some extent, but had no large effect on the total ice volume or area in these tests. As a qualifier, these experiments do not adequately explore subglacial flow dynamics. The basal sliding that we specify, (Eqs. (4) and (5)) is a simple formulation that is not likely to represent the large-scale basal flow that occurs in nature. In a separate study, for instance, we have not been successful in simulating contemporary West Antarctic ice sheet thickness. This implies that the modelled LGM ice volume is likely to be an upper limit for North America; increased basal flow and longitudinal ice deformation would produce thinner ice sheets (e.g., Clark, 1994; Clark et al., 1996). We suggest that this effect would not be fully developed at LGM, as much of the NAIS complex is still cold-based at this time. However, a more physically based treatment of basal flow (WAIS physics) would still be expected to lower the LGM ice volume in North America. 4.2.5. Isostatic model We did not find any strong or systematic effect of isostatic treatment on LGM NAIS reconstructions in these experiments. Ice volumes and divide structures do change with variations in the Earth model, but the isostatic model does not act as a primary determinant of LGM ice sheet disposition, as evaluated by first-order constraints from ice sheet area and LGM margins. This is consistent with the findings of Tarasov and Peltier ð1999Þ: Earth model impacts can be important at ice sheet margins and during deglaciation. Furthermore, ice-sheet Earth models in conjunction with relative sea-level curves can provide compelling additional constraints

190

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

on ice sheet history, as Peltier (1994, 1996) has shown. We have focused on LGM reconstructions here, for which sea-level curves offer only weak model discrimination in our model; Holocene and modern-day uplift rates are primarily controlled by the deglacial history, with much less memory of LGM conditions. Because the deglaciation is not well represented in the model, particularly with respect to complexities and asynchronies in the spatial retreat pattern, we do not yet feel confident discriminating between LGM ice sheet reconstructions based on relative sea-level records. We suggest that this avenue should be carefully investigated in followup studies, however, as the relative sea-level record offers unique and independent insights into LGM ice sheet form and thickness.

5. Conclusions The suite of 190 glacial-cycle simulations produced a wide spread of scenarios for the LGM ice sheets in North America. A subset of 33 cases came close to the geologically inferred ice sheet area and southern margin position at LGM. This subset of tests, summarized in Table 3, had a mean LGM ice volume of 34:3
1015 m3 ; corresponding to 82:7 msl after the correction for floating ice and displaced bodies of water (ice grounded below present-day sea-level). This includes contributions from the Cordilleran, Laurentide, and Innuitian Ice Sheets, 7:1; 72:4; and 3:2 msl; respectively. There is a considerable spread in ice volume in the accepted tests, but more than 75% of the accepted simulations vote for a eustatic sea-level change in the range 78–88 msl: The accepted results represent a fairly consistent climatic range, with a positive net mass balance of 0.2–0:22 m yr1 over the ice sheet at LGM and a mean temperature depression of ca. 111C over the ice sheet, relative to modern conditions. This temperature depression includes both lapse-rate and climatic effects. NAIS area, volume, divide structure, and transient evolution are very dependent on glacial-period precipitation rates and spatial patterns in the model. We failed to achieve a successful simulation using modern-day spatial precipitation patterns. This emphasizes the important role of enhanced precipitation on the ice sheet margins, particularly in the south, as a result of orographic effects and synoptic circulation shifts caused by the ice sheets. The single most sensitive model parameter in our simulations was the threshold elevation ht in the elevation-drying parameterization, Eq. (18). This parameter appears critical to depress precipitation rates in interior portions of the ice sheet, amplifying the elevation-desert or interior plateau effect that is only weakly present in the GCM simulations for LGM. NAIS reconstructions are substantially improved, parti-

cularly in western Canada and the southern interior plains, when a strong elevation-desert effect is prescribed. We conclude that synoptic circulation and precipitation patterns during the glacial period play a critical role in dictating LGM ice extent in North America. Ice rheology and basal flow represent the glaciological processes with the largest impact on NAIS reconstructions. Both processes influence the ice sheet aspect ratio (horizontal/vertical extent), with thinner ice sheets resulting from the low-exponent flow laws. We had difficulty achieving sufficient southern ice sheet extent with the low-exponent flow laws (n ¼ 1:7) and without fairly vigorous basal flow. One of the greatest problems with the former class of runs is the propensity for pervasive cold-based conditions in the model, with LGM ice sheets having more than 85% of the ice sheet frozen to the bed in the n ¼ 1:7 simulations. This is a predictable result of thin ice sheets with diminished strain heating, more diffusive heat loss at the base, and high flow rates. Enhanced advection increases the downward transfer of cold surface ice, as well as the horizontal (downglacier) transfer of cold ice from interior portions of the ice sheet. Collectively, these thermodynamic effects promote cold-based ice, which prohibits basal flow and inhibits marginal ice advance. While these findings do not support the likelihood of low-exponent flow laws, we cannot rule them out. Determined exploration of parameter space will doubtlessly reveal a combination that produces plausible LGM ice sheets for any given ice rheology. Modelled basal temperature fields are also uncertain (e.g., Payne and Huybrechts, 2000). The basic physical controls of basal temperature conditions (strain heating, advection, and diffusion) nevertheless present an important consideration. The geologic record testifies to a much more active subglacial environment than one would expect from widespread frozen conditions (e.g. Clark, 1997). The model ensemble examined here is far from exhaustive. It provides some insight into the glaciological possibilities, however. We feel confident that we have adequately sampled the parameter range of glaciological, paleoclimatological, and mass balance variables within the current level of model sophistication. Climatic patterns and glaciological processes beyond the current model capabilities are unexplored of course; improved treatment of basal flow processes, for instance, will be likely to lower the range and mean of model simulations. We suggest that 78–88 msl be taken as a maximum estimate of LGM ice volume in North America, with future work to improve the representation of basal flow (WAIS physics) expected to bring this down. Because much of the LIS is still cold-based at LGM, however, we believe that 80 msl may prove to be a reasonable estimate.

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

Acknowledgements This paper is a contribution to the Climate System History and Dynamics Program (CSHD), Phase II, sponsored by the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Meteorological Service of Canada. SJM also acknowledges support from the U.S. National Science Foundation Earth System History program ESH NSF 9905535. Martin Siegert provided a perceptive and helpful review. We thank A.S. Dyke at the Geological Survey of Canada, Ottawa, for sharing his ongoing insights into North American Ice Sheet history. He generously provided LGM ice sheet area estimates and relative sea level data from northeastern Canada. We are grateful to P.U. Clark, E. Bard, and A.C. Mix for the chance to participate in the EPILOG meeting at Mt. Hood; this was an extremely stimulating and productive workshop. Geological Survey of Canada contribution 2001087.

References Bard, E., Hamelin, B., Fairbanks, R.G., 1990. U-Th ages obtained by mass spectrometry in corals from Barbados: sea level during the last 130,000 years.. Nature 346, 456–458. Bindschadler, R.A., Vornberger, P.L., 1990. AVHRR imagery reveals Antarctic ice dynamics. EOS 71 (23), 741–742. Blackwell, D.B., Steele, J.L., 1992. Geothermal Map of North America. Decade of North American Geology, Map 006, Gelogical Society of America, Boulder, Colorado. Boulton, G.S., Smith, G.D., Jones, A.S., Newsome, J., 1985. Glacial geology and glaciology of the last mid-latitude ice sheets. Geological Society of London Journal 142, 447–474. Braithwaite, R., 1984. Calculation of degree-days for glacier-climate research. Zeitschrift fur . Gletscherkunde und Glazialgeologie 20, 1–8. Braithwaite, R., 1995. Positive degree-day factors for ablation on the Greenland ice sheet studied by energy-balance modelling. Journal of Glaciology 41, 153–160. Clague, J.J., 1989. Cordilleran ice sheet. In: Fulton, R.J. (Ed.), Quaternary Geology of Canada and Greenland. Geological Survey of Canada, Geology of Canada, Vol. 1, pp. 40–42. Clark, C.D., 1997. Reconstructing the evolutionary dynamics of former ice sheets using multi-temporal evidence, remote sensing, and GIS. Quaternary Science Reviews 16, 1067–1092. Clark, P.U., 1992. Surface form of the southern Laurentide Ice Sheet and its implications to ice-sheet dynamics. Geological Society of America Bulletin 104, 595–605. Clark, P.U., 1994. Unstable behaviour of the Laurentide Ice Sheet over deforming sediment and its implications for climate change. Quaternary Research 41, 19–25. Clark, P.U., Clague, J.J., Curry, B.B., et al., 1993. Initiation and development of the Laurentide and Cordilleran ice sheets following the last interglaciation. Quaternary Science Reviews 12, 79–114. Clark, P.U., Licciardi, J.M., MacAyeal, D.R., Jenson, J.W., 1996. Numerical reconstruction of a soft-bedded Laurentide ice sheet during the last glacial maximum. Geology 24 (8), 679–682.

191

Clark, P.U., Alley, R.B., Pollard, D., 1999. Northern Hemisphere ice-sheet influences on global climate change. Science 286, 1104–1111. Cuffey, K.M., Clow, G.D., 1997. Temperature, accumulation, and ice sheet elevation in central Greenland through the last deglacial transition. Journal of Geophysical Research 102 (C12), 26, 383–26,396. Dahl-Jensen, D., Mosegaard, K., Gundestrup, N., Clow, G.D., Johnsen, S.J., Hansen, A.W., Balling, N., 1998. Past temperatures directly from the Greenland ice sheet. Science 282, 268–271. Dansgaard, W., Johnsen, S.J., Clausen, H.B., Dahl-Jensen, D., Gundestrup, N.S., Hammer, C.U., Hvidberg, C.S., Steffensen, . J.P., Sveinbjornsd! ottir, A.E., Jouzel, J., Bond, G.C., 1993. Evidence for general instability of past climate from a 250 kyr ice-core record. Nature 364, 218–220. Deblonde, G., Peltier, W.R., 1993. Pleistocene ice age scenarios based upon observational evidence. Journal of Climate 6, 709–727. Doake, C.S.M., Wolff, E.W., 1985. Matters arising. Flow law for ice in polar ice sheets. Nature 318, 83. Dyke, A.S., Prest, V.K., 1987. Late Wisconsinan and Holocene history of the Laurentide ice sheet. G!eographie Physique et Quaternaire 41, 237–264. Fisher, D.A., 1987. Enhanced flow of Wisconsin ice related to solid conductivity through strain history and recrystallization. International Association of Hydrological Sciences 170, 45–51. Fisher, D.A., Reeh, N., Langley, K., 1985. Objective reconstructions of the late Wisconsinan Laurentide Ice Sheet and the significance of deformable beds. G!eographie Physique et Quaternaire 39, 229–238. Fleming, K., Johnston, P., Zwartz, D., Yokoyama, Y., Lambeck, K., Chappell, J., 1998. Refining the eustatic sea-level curve since the Last Glacial Maximum using far- and intermediate-field sites. Earth and Planetary Science Letters 163, 327–342. Glen, J.W., 1958. The flow law of ice: a discussion of the assumptions made in glacier theory, their experimental foundations and consequences. International Association of Hydrological Sciences 47, 171–183. Hock, R., 1999. A distributed temperature-index ice- and snowmelt model including potential direct solar radiation. Journal of Glaciology 45 (149), 101–111. Hughes, T.J., 1977. West Antarctic ice streams. Reviews of Geophysics and Space Physics 15 (1), 1–46. Huybrechts, P., 1990. A 3-D model for the Antarctic ice sheet: a sensitivity study on the glacial–interglacial contrast. Climate Dynamics 5, 79–92. Huybrechts, P., T’Siobbel, S., 1995. Thermomechanical modelling of northern hemisphere ice sheets with a two-level mass-balance parameterisation. Annals of Glaciology 21, 111–117. Jackson Jr., L.E., Phillips, F.M., Shimamura, K., Little, E.C., 1997. Cosmogenic 36 Cl dating of the Foothills erratics train, Alberta, Canada. Geology 25, 195–198. James, T.S., Ivins, E.R., 1998. Predictions of Antarctic crustal motions driven by present-day ice sheet evolution and by isostatic memory of the Last Glacial Maximum. Journal of Geophysical Research 103, 4993–5017. Jenssen, D., 1977. A three-dimensional polar ice sheet model. Journal of Glaciology 18, 373–389. Jenson, J.W., MacAyeal, D.R., Clark, P.U., Ho, C.L., Vela, J.C., 1996. Numerical modeling of subglacial sediment deformation: implications for the behaviour of the Lake Michigan Lobe, Laurentide Ice Sheet. Journal of Geophysical Research 101 (B4), 8717–8728. ! Johannesson, T., Sigurdsson, O., Laumann, T., Kennett, M., 1995. Degree-day glacier mass-balance modelling with application to glaciers in Iceland, Norway, and Greenland. Journal of Glaciology 41, 345–358.

192

S.J. Marshall et al. / Quaternary Science Reviews 21 (2002) 175–192

Johnson, W.H., 1990. Ice-wedge casts and relict patterned ground in central Illinois and their environmental significance. Quaternary Research 33, 51–72. Kalnay, E., et al., 1996. The NCEP/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society 77, 437–471. Lambeck, K., Yokoyama, Y., Johnston, P., Purcell, A., 2000. Global ice volumes at the Last Glacial Maximum and early Lateglacial. Earth and Planetary Science Letters 181, 513–527. Legates, D.R., Willmott, C.J., 1990. Mean seasonal and spatial variability in gauge-corrected global precipitation. International Journal of Climate 10, 111–127. MacAyeal, D.R., 1989. Large-scale ice flow over a viscous basal sediment: theory and application to Ice Stream B, Antarctica. Journal of Geophysical Research 94 (B4), 4071–4087. Mahaffy, M.W., 1976. A three-dimensional numerical model of ice sheets: tests on the Barnes Ice Cap, Northwest Territories. Journal of Geophysical Research 81 (B6), 1059–1066. Marshall, S.J., Clarke, G.K.C., 1997. A continuum mixture model of ice stream thermomechanics in the Laurentide Ice Sheet 1. Theory. Journal of Geophysical Research 102, 20,599–20,613. Marshall, S.J., Clarke, G.K.C., 1999. Modeling North American freshwater runoff through the Last Glacial cycle Quaternary Research, 52, 300–315. Marshall, S.J., Cuffey, K.M., 2000. Peregrinations of the Greenland Ice Sheet divide in the last glacial cycle: implications for central Greenland ice cores. Earth and Planetary Sciences Letters 179, 73–90. Marshall, S.J., Clarke, G.K.C., Dyke, A.S., Fisher, D.A., 1996. Topographic and geologic controls on fast flow in the Laurentide and Cordilleran ice sheets. Journal of Geophysical Research 101 (B8), 17,827–17,840. Marshall, S.J., Tarasov, L., Clarke, G.K.C., Peltier, W.R., 2000. Glaciological reconstruction of the Laurentide ice sheet: physical processes and modelling challenges. Canadian Journal of Earth Sciences 37, 769–793. McFarlane, N.A., Boer, G.J., Blanchet, J.-P., Lazare, M., 1992. The Canadian Climate Centre second generation general circulation model and its equilibrium climate. Journal of Climate 5, 1013–1044. Paterson, W.S.B., 1991. Why ice-age ice is sometimes ‘‘soft’’. Cold Regions Science and Technology 20, 75–98. Paterson, W.S.B., 1994. The Physics of Glaciers, 3rd Edition. Elsevier Science Ltd, New York. Payne, A.J., 1998. Dynamics of the Siple Coast ice streams, West Antarctica: results from a thermomechanical ice sheet model. Geophysical Research Letters 25 (16), 3173–3176. Payne, A.J., Huybrechts, P., 2000. Results from the EISMINT model intercomparison: the effects of thermomechanical coupling. Journal of Glaciology 46 (153), 227–238.

Peltier, W.R., 1974. The impulsive response of a Maxwell Earth. Reviews of Geophysics and Space Physics 12, 649–668. Peltier, W.R., 1985. The LAGEOS constraint on deep mantle viscosity: results from a new normal mode method for the inversion of viscoelastic relaxation spectra. Journal of Geophysical Research 90, 9411–9421. Peltier, W.R., 1994. Ice age paleotopography. Science 265, 195–201. Peltier, W.R., 1996. Mantle viscosity and ice-age ice sheet topography. Science 273, 1359–1364. Peltier, W.R., 1998. Postglacial variations in levels of the sea: implications for climate dynamics and solid earth geophysics. Reviews of Geophysics 36 (4), 603–689. Peltier, W.R., Goldsby, D.L., Kohlstedt, D.L., Tarasov, L., 2000. Ice age ice sheet rheology: constraints from the Last Glacial Maximum form of the Laurentide ice sheet. Annals of Glaciology 30, 163–176. Pollard, D., Thompson, S.L., 1997. Climate and ice-sheet mass balance at the last glacial maximum from the GENESIS version 2 global climate model. Quaternary Science Reviews 16 (8), 841–863. Ramstein, G., Fabre, A., Pinot, S., Ritz, C., Joussaume, S., 1997. Icesheet mass balance during the last glacial maximum. Annals of Glaciology 25, 145–152. Row, L.W. III, Hastings, D.A., 1994. TerrainBase: Worldwide Digital Terrain Data. CD-ROM. National Oceanic and Atmospheric Administration, National Geophysical Data Center, Boulder, Colorado. Shilts, W.W., 1980. Flow patterns in the central North American ice sheet. Nature 286, 213–218. Tarasov, L., Peltier, W.R., 1999. Impact of thermomechanical ice sheet coupling on a model of the 100 kyr ice age cycle. Journal of Geophysical Research 104 (D8), 9517–9545. Thorsteinsson, T., Waddington, E.D., Taylor, K.C., Alley, R.B., Blankenship, D.D., 1999. Strain rate enhancement at Dye 3, Greenland. Journal of Glaciology 45, 338–345. Veillette, J.J., Dyke, A.S., Roy, M., 1999. Ice-flow evolution of the Labrador sector of the Laurentide Ice Sheet: a review, with new evidence from northern Quebec. Quaternary Science Reviews 18, 993–1019. Vettoretti, G., Peltier, W.R., McFarlane, N.A., 2000. Global water balance and atmospheric water vapour transport at Last Glacial Maximum: climate simulations with the CCCma atmospheric general circulation model. Canadian Journal of Earth Sciences 37, 695–723. Wu, P., Peltier, W.R., 1982. Viscous gravitational relaxation. Geophysical Journal of the Royal Astronomical Society 70, 435–485. Yokoyama, Y., Lambeck, K., De Deckker, P., Johnston, P., Fifield, L.K., 2000. Timing of the Last Glacial Maximum from observed sea-level minima. Nature 406, 713–716.