Pre-lgm Northern Hemisphere Ice Sheet Topography

We here reconstruct the paleotopography of Northern Hemisphere ice sheets during the glacial maxima of marine isotope stages (MIS) 5b and 4.We employ a combined approach, blending geologically based reconstruction and numerical modeling, to arrive at probable ice sheet extents and topographies for each of these two time slices. For a physically based 3-D calculation based on geologically derived 2-D constraints, we use the University of Maine Ice Sheet Model (UMISM) to calculate ice sheet thickness and topography. The approach and ice sheet modeling strategy is designed to provide robust data sets of sufficient resolution for atmospheric circulation experiments for these previously elusive time periods. Two tunable parameters, a temperature scaling function applied to a spliced Vostok–GRIP record, and spatial adjustment of the climatic pole position, were employed iteratively to achieve a good fit to geological constraints where such were available. The model credibly reproduces the first-order pattern of size and location of geologically indicated ice sheets during marine isotope stages (MIS) 5b (86.2 kyr model age) and 4 (64 kyr model age). From the interglacial state of two north–south obstacles to atmospheric circulation (Rocky Mountains and Greenland), by MIS 5b the emergence of combined Quebec–central Arc-tic and Scandinavian–Barents-Kara ice sheets had increased the number of such highland obstacles to four. The number of major ice sheets remained constant through MIS 4, but the merging of the Cordilleran and the proto-Laurentide Ice Sheet produced a single continent-wide North American ice sheet at the LGM.


Introduction
We here address a fundamental information gap in climate science, Northern Hemisphere paleotopography during the last interglacial-to-glacial transition.At present, a good understanding of Northern Hemisphere paleotopography exists for the last glacial maximum (LGM) and the deglacial phase of the last glacial cycle (Peltier, 2004).However, there is a lack of geologically constrained data sets for defined time intervals during the 95 kyr build-up phase towards the LGM.Considerable progress has been made in understanding ice sheet dynamics in this elusive time interval (Stokes et al., 2012), but the scarcity of geological and geomorphological data (Kleman et al., 2010) that can constrain numerical models is still an impediment to our understanding of the time periods of most rapid ice sheet build-up.This situation hampers research regarding ice sheet topographical feedbacks on atmospheric circulation (Cook and Held, 1988;Roe and Lindzén, 2001;Langen and Vinther, 2009;Liakka and Nilsson, 2010), and is an obstacle to our understanding of Northern Hemisphere atmospheric circulation during periods which were very different from both interglacial and full-glacial conditions in terms of number, location and size of ice sheets.Indeed, idealized coupled atmosphere-ice-sheet models have demonstrated that the interaction between ice sheet topography and atmospheric circulation can strongly influence the spatial distribution of continental-scale ice sheets (Roe and Lindzén, 2001;Liakka et al., 2011).Thus, a number of important questions regarding ice sheet evolution and the role of climatic feedbacks can only be addressed from a foundation of improved paleotopographical reconstructions.For example, could atmospheric circulation changes have caused mass balance patterns that left Alaska largely ice-free throughout the last glacial cycle, the high-elevation Cordilleran region only ephemerally glaciated (Clague, 1989), but the almost entirely lowland Laurentide area (Fig. 1) persistently glaciated?Did downstream influence of the North American ice sheet topography on the atmospheric circulation contribute to southwestward migration of the mass center of the Eurasian Ice Sheet?Amplification of the effects of modest insolation changes appears necessary to explain ice sheet build-up in the Laurentide area (Otieno and Bromwich, 2009).In addition, the topographical and albedo changes caused by ice sheets are likely part of the non-linear response of the climate system to the insolation changes that pace glacial cycles.Considerable progress is being made in establishing the general response characteristics of the climate system by statistical analysis of the global ice volume response to insolation forcing (Lourens et al., 2010;Imbrie et al., 2011).However, to understand the impact of changing topography on atmospheric circulation and the growth of ice sheets, also the location, extent and height of individual ice sheets must be known.
Improved understanding of Northern Hemisphere paleotopography of intermediate-sized ice sheets can, if combined with the relevant array of paleoclimatic proxy data, yield a valuable calibration data set for tuning and comparison of climate models at global ice volumes smaller than the LGM case normally used (Kageyama et al., 1999;Otto-Bliesner et al., 2006).The present situation, in which research focus has been on two climatic end members, largely precludes research on non-linear responses to ice sheet build-up.Over the last decades evidence has accumulated that ice sheet extent and volume varied more drastically over shorter time scales than previously thought, and a new paradigm of highly dynamic ice sheets has emerged (Boulton and Clark, 1990).Of particular importance is the evidence for ice-free core areas in Fennoscandia (Wohlfarth, 2010;Helmens and Engels, 2010), the consequence of which is that the Fennoscandian Ice Sheet had to be fully rebuilt from zero ice volume during each of the cold phases (MIS 5d, 4 and 2).From this new paradigm, time and precipitation emerge as important limiting factors for the glacial development, and there is a need to explore the atmospheric circulation and conditions that allowed ice sheets to be fully rebuilt during the relatively short intervals of declining Northern Hemisphere highlatitude summer insolation.

Previous work
Reconstructions of the LIS at the LGM based on numerical ice sheet modeling have ranged from a thick monodome ice sheet (Denton and Hughes, 1981) to a relatively thin multidomed ice sheet (Clark et al., 1996).More recent ice sheet models (Marshall et al., 2000;Tarasov and Peltier, 2004;Álvarez-Solas et al., 2011;Gregoire et al., 2012) agree that the LIS at the LGM exhibited a multi-domed shape with a complex configuration of ice streams and inter-stream ridges in peripheral areas.Likewise there is a convergence on the existence of frozen-bed core areas under the main domes, in line with was inferred from geomorphological evidence by Kleman and Hättestrand (1999).Few numerical modeling exercises focus on, or include, the growth phase of the ice sheet (Huybrechts and T'siobbel, 1995;Charbit et al., 2007;Bintanja et al., 2002;Bonelli et al., 2009).Many of those that do so conflict with the geological evidence for ice sheet extent at the LGM, which is currently the prime available "target" period for assessing the validity and credibility of the climates used in these numerical models.Common is too much ice cover in Alaska and the Mackenzie mountains (NW sector), areas that were largely ice-free at the LGM, and also continuous and massive glaciation of British Columbia instead of the geologically indicated ephemeral ice cover (Clague, 1989).These common problems are discussed at length in Bonelli et al. (2009).Most models are also unable to generate a southward extent of the Quebec Dome during ice sheet build-up that is sufficiently large to match the geological evidence (Stea, 2004;Kleman et al., 2010).The recent ensemble-based reconstruction by Stokes et al. (2012) is focused on the build-up phase and represents a major step forward in understanding the dynamics of proto-Laurentide ice sheets.It was tuned on the basis of data for the post-LGM period, and the output for the pre-LGM stages are thoroughly discussed in relation to geological data.One of the currently leading ice sheet models (ICE-5G, Peltier, 2004) is based on isostatic recovery rather than mass-balance-driven ice physics, and is unable to reconstruct pre-LGM conditions because of the time constants involved in isostatic recovery.An early attempt at a geologically and geomorphologically constrained pre-LGM LIS model was presented by us some years ago (Kleman et al., 2002), and we here build on experience gained during that work.
Ice sheet reconstructions for the 115-21 kyr build-up phase involve additional difficulties compared to reconstructions for the LGM and deglacial time, because the amount of stratigraphic control for this time interval is much smaller than for the deglaciation period 21-8 kyr.The time interval is much longer, and can rarely be addressed with radiocarbon dates.In addition, almost all landforms and deposits created in this time interval have been subjected to overriding ice during the LGM and later stages, which in most cases has led to erosion or burial and therefore full or partial destruction and reshaping, and to severe spatial fragmentation of the record compared to the records from the LGM and deglacial periods (Kleman et al., 2010).These difficulties add up to a major methodological challenge, and require a focus on only the first-order patterns in the evolution of ice sheets in this time interval.For the pre-LGM period it is not possible to achieve the spatial and temporal precision of post-LGM reconstructions (Dyke et al., 2002;Ehlers and Gibbard, 2003).
We here present first-order paleotopographical data sets for the ice volume maxima of marine isotope stages 5b and 4, the two youngest build-up phases preceding the LGM.These reconstructions are meant to be relevant as topographic input for atmospheric circulation experiments; the data can and are intended to be used in modeling studies of the impact of these early ice sheet configurations on the large-scale atmospheric circulation.There is a considerable amount of literature concerning such studies for the LGM (e.g., Manabe and Broccoli, 1985;Cook and Held, 1988), but much less attention has been devoted to understanding the pre-LGM atmospheric response.In particular, the atmospheric stationary waves, which are topographically forced and have wavelengths of a few 1000 km (Held et al., 2002), are of key importance in this context.The bed conditions, frozen or thawed, exert some control on ice sheet thickness, but the finer-scale details of the ice sheet topography, although of importance for local climatic conditions, have a relatively weaker impact on the hemispheric-scale stationary waves.Accordingly, we have not tried to capture the details of ice sheet outline or internal ice sheet dynamics, since the large-scale atmospheric flow is not critically dependent on such information.We have used the best available geological constraints on ice sheet extent and, through numerical modeling, strived for first-order reconstructions that represent our current best understanding of location, extent and thickness in this time interval.
Our approach is not puristic.We do not regard geological data as an indisputable "truth" against which model output is evaluated.Instead, in line with Kleman et al. (2010), we argue that discrepancies between "control" data and numerical modeling results can point to flaws in interpretation of data as well as deficiencies in the numerical modeling.We employ both approaches here -geologically based reconstruction and numerical modeling -to arrive at the most probable ice sheet extent and topography for each time slice.

Methods
Our strategy is governed by the intended end use of paleoatmospheric circulation modeling.For that purpose a highprecision ice sheet outline is not crucial, but the existence or non-existence of an ice sheet in a certain location is important.It can be expected that small and thin ice sheets will not significantly affect the large-scale atmospheric circulation.Hence we are primarily concerned with capturing the main ice sheet features and ice sheets with a length of at least 1000 km.
The situation regarding data that spatially constrain pre-LGM ice sheets is very uneven for the respective glaciation centers, with a relatively good understanding in western Eurasia; some effective constraints, albeit with poor dating control, for the Laurentide Ice Sheet; and very little solid knowledge for the Cordilleran and Inuitian ice sheets in these time intervals.
We have here chosen to base our synthesis primarily on previous spatial reconstructions (Dredge and Thorleifson, 1987;Clague, 1989;Clark et al., 1993;Kleman et al., 1997Kleman et al., , 2010;;Lundqvist, 2004;Svendsen et al., 2004;Mangerud, 2004;Winsborrow et al., 2004;Lambeck et al., 2006) that cover this time interval.A full consideration of primary morphological and stratigraphic data is offered in several of the key source publications.In compiling the set of 2-D modeling targets (Fig. 2), each corresponding to a certain time interval, we have used the following approach: (i) where approximate consensus on ice sheet outline exists, we have drawn an averaging and somewhat simplified outline; (ii) where clearly dissenting views exist, we have weighed the quality and amount of geological evidence and decided on what we consider the most likely alternative; and (iii) for an outline that is considered spatially reliable, but for which alternative age assignments exist, we have weighed the age options and chosen that which we consider to be the most probable one.
The generally poor understanding of pre-LGM extent of particularly the Cordilleran and Inuitian ice sheets constitutes a "missing ice sheet" problem.There is not sufficient, nor reliable enough, data to allow us to show ice sheets' outlines in these areas on any of the panels in Fig. 2, yet embryonic or smaller ice sheets are very likely to have existed in these areas during the time intervals in question.To use the outlines shown in Fig. 2 as "hard" modeling targets (in which ice in areas marked as ice-free in Fig. 2 is regarded as indicating a poor model fit) would clearly be erroneous, and would have caused us to use inappropriate forcing of the model to keep ice out of these areas.
For both Eurasia and North America there are indications that the ice sheet extent was approximately similar during stages 5d and 5b.In view of the scarce data and dating difficulties in this time interval it is futile to separate spatial evidence from the two with any degree of certainty.Instead we have chosen to focus on stage 5b, which is closest in time LGM ice extent is indicated with stippled line.For source data and construction of outlines, see text.
to the two other major cold stages, MIS 4 and MIS 2 (the LGM).Given the similarity in global ice volume (Lisiecki and Raymo, 2005) and available sea level data (Lambeck et al., 2002) between the two stages, we consider it reasonable to regard this pattern as archetypical also for stage 5d.

North America
For the extent of North American ice sheets before the LGM, we have based our 2-D outlines primarily on the recent review by Kleman et al. (2010) and references therein, particularly Dredge and Thorleifson (1987), Clague (1989), Boulton and Clark (1990), Clark et al. (1993), andStea et al. (2004).
For the Keewatin sector of the Laurentide Ice Sheet there is good evidence for one or two, probably closely related, stages of a pre-LGM ice sheet with a dome center in the northern Keewatin or south-central Arctic region (Kleman et al., 2002(Kleman et al., , 2010;;McMartin and Henderson, 2004).However, the east-west extents, maximum southerly extents, and ages are poorly constrained.Inferences about their pre-LGM age can be made on the basis of flow pattern and relation to other flow systems, but we know of no firm dating constraints.The flow patterns do indicate that both were dynamically independent of the Quebec Dome.The lineation evidence for this ancestral northern Keewatin dome highlights the missing ice sheet problem.The lineation evidence for flow out of an ancestral northern Keewatin/central Arctic dome is solid but covers only a restricted sector or corridor which lacks any topographic funneling.Glaciological considerations therefore require this lineation system to be only a fragment of a much wider flow pattern emanating from northern Keewatin or the south-central Arctic islands.Despite its limited area it thus indicates the existence of a full ice sheet, about whose eastwest extent unfortunately not enough is known to warrant its full inclusion in Fig. 2.However, its minimum southward extent into the area of preserved traces is a viable modeling target.
In the Quebec sector, Kleman et al. (2002) assigned the "Atlantic" flow system, which indicates an ice divide closer to the Atlantic coast than any subsequent stage, to stage 5b or 5d.We consider it likely that this pattern characterized both stages and have therefore included it as a constraint for stage 5b.Regarding the Cordilleran Ice Sheet, Kleman et al. (2010) concluded that there was flow-trace evidence for "old" configurations with more easterly dispersal centers (Stumpf et al., 2000) than the ones that prevailed in near-LGM time, but constraints on age and outline are lacking.Hence, the Cordilleran Ice Sheet is one element of the missing ice sheet problem, and it appears very likely that this region was at least partly glaciated during MIS 5d, 5b, and 4, but there is simply not enough geological and geomorphological data to suggest an outline.Similarly, it appears likely that a proto-Inuitian ice sheet existed on the Queen Elizabeth Islands during these stages, but again there is not enough data to warrant the suggestion of a specific outline.

Eurasia
For northern Eurasia we mainly follow the spatial reconstructions of Svendsen et al. (2004), with consideration also of Lambeck (2006).For Fennoscandia, the main sources we have used are Mangerud (2004), Lundqvist (2004), Kleman et al. (1997), andLambeck et al. (2006).For stage 5b there is agreement that the Bothnian Bay and Gulf were glaciated, but not the Baltic.This indicates a width of the Fennoscandian Ice Sheet of less than half its LGM width.Southern Sweden was not glaciated at this time, nor were the British Isles.According to Svendsen et al. (2004), a unified ice sheet existed from southern Norway all the way to the Putorana Plateau.Lambeck et al. (2006) instead indicate a narrow gap between the Scandinavian and Barents-Kara components.We here follow the Svendsen et al. (2004) reconstruction, but regard this suggested separation as a minor issue since the overall outline of the two reconstructions are quite similar.During stage 4, there was apparently an overall shift of mass towards Fennoscandia, and a reduction of ice mass in the Kara sector.The Putorana Plateau was unglaciated or, alternatively, had an independent ice cap.The Fennoscandian Ice Sheet filled the entire Baltic depression and touched upon mainland Europe.The eastern margin of the Scandinavian Ice Sheet at this stage was approximately halfway between its stage 5b and LGM positions.The British Isles were glaciated during the LGM, but we acknowledge great uncertainty regarding the ice extent and possible confluence with Scandinavian ice during stage 4.Although volumetrically small, the southwestward lengthening of the Eurasian Ice Sheet that a British Ice Sheet may have provided was potentially important for atmospheric circulation changes.

Model
We employ the University of Maine Ice Sheet Model (UMISM), which consists of a time-dependent finite-element solution of the coupled mass, momentum, and energy conservation equations (Fastook and Chapman, 1989;Fastook, 1993) with a broad range of applications (Holmlund and Fastook, 1995;Johnson and Fastook, 2002;Kleman et al., 2002).The primary input to the model is present bedrock topography, the surface mean annual temperature, the geothermal heat flux, and the net mass balance, all defined as functions of position.The solution consists of ice thicknesses, surface elevations, column-integrated ice velocities, the temperature field within the ice sheet, the amount and distribution of water at the bed resulting from basal melting, and the amount of bed depression resulting from the ice load.
The primary solution is of the mass conservation equation.The required column-integrated ice velocities at each point in the map plane are obtained by numerically integrating the momentum equation through the depth of the ice.Stress is a function of depth and related to strain rates through the flow law of ice.The temperature, on which the flow law ice hardness depends, is obtained from a 1-D finite-element solution of the energy conservation equation, which includes internal heat generation produced by shear with depth and sliding at the bed.
Ice hardness is calculated from these modeled internal temperatures and is based on ice hardness parameterizations from Paterson (1994), modified by an appropriate flow enhancement factor designed to accommodate fabric and debris content known to differ between hemispheres (Huybrechts, 1990;Ma et al., 2010).Sliding at the bed follows Weertman (1964Weertman ( , 1969) ) and Weertman and Birchfield (1982).The amount of water at the bed is a function of the internal temperature calculation.In addition, a 2-D solution for conservation of water at the bed allows for movement of the basal water down the hydrostatic pressure gradient (Johnson and Fastook, 2002).
While this is a shallow-ice approximation model where only the basal shear stress is considered, we do include a longitudinal thinning strain rate at the grounding line that models ice drawn into a buttressing ice shelf that itself is not explicitly modeled.This thinning rate is added to the mass balance source term of the mass conservation equation only in elements that contain a grounding line, and it is proportional to grounding line thickness raised to the fourth power (Weertman, 1957).We modify this rate by a parameter we call the "Weertman" that ranges from 0 (full buttressing, no thinning into the ice shelf) to 1 (no buttressing, thinning at a rate commensurate with an unconfined floating ice slab).While not used in this exercise, and assumed to be nominally zero everywhere inside the continental shelf break, this parameter can be used to control the advance and retreat of a marine margin.

Model setup
With surface temperature and mass balance obtained from measured data, from GCM results, or from a climate parameterization, the model can be run with no more input than the specified bedrock topography, in what can be called a "free-running mode".However, it can also be tightly constrained by geological and glaciological data when such data are available.Ice margin positions can be specified and areas of basal sliding can be derived from the distribution of basal temperatures, or they can be specified by the presence of erosional features on exposed landscapes.In this exercise all of Eurasia and most of North America are free running.The exception is that in North America glaciation of the topographically high Rocky Mountains is suppressed.
We used a relatively coarse grid spacing of 95 km, suitable for the broad picture we are attempting to produce for use in atmospheric circulation modeling.A regular rectangular grid with 10 000 (100 × 100) nodal points spanning the Northern Hemisphere was used.We divided this into a Western Hemisphere grid and an Eastern Hemisphere grid so that climatic adjustment could be applied separately to the two hemispheres to obtain the best fit.The time step required at this resolution was 10 yr.
As a forcing function, we used a global temperature proxy based on the oxygen-isotope record from a spliced GRIP (Johnsen et al., 1997) and Vostok (Petit et al., 1999) ice core record (shown in Fig. 3a as a red line labeled "Composite").The latter part of the GRIP record is reliably dated, but reliability deteriorates towards the last intergacial.Also the GRIP record contains two extreme but short-lived cooling events at 70-75 kyr.The nature of these events is not well understood (Svensson et al., 2013), and the corresponding events in Antarctic records are not of a similar magnitude (EPICA community members, 2006).Because of the possibility that the large magnitude of these events in the GRIP record is a localized Greenland phenomenon, that may not well reflect regional climate in North America and Eurasia, we preferred the Vostok record for these time intervals.In previous modeling experiments (Kleman et al., 2002) these short events created brief and very large increases in ice sheet area, for which no positive geological evidence is known.In addition, one can see by comparing the SPECMAP-stacked oxygenisotope record (Imbrie et al., 1984), itself thought to represent global sea level and the basis for the definitions of the stages we are investigating, that there is better agreement with the Vostok record than with the GRIP core.This appears to also be true with regards to a terrestrial oxygenisotope record from North America, DH-11 (Landswehr et al., 1997), which may better correlate with Vostok than with GRIP.However, the Vostok core does not strongly record the Allerøod warming nor the Younger Dryas cold reversal whose impact is well documented in Northern Hemisphere ice sheet retreat records.While behavior of the retreating ice sheets is not part of this study, we did observe marginal positions at the Younger Dryas as one of our targets for goodness of fit.We therefore employed a spliced Vostok-GRIP record, Vostok prior to 20 KBP during which our target time slices occur, and GRIP during the deglaciation for comparison with Younger Dryas marginal positions.
The use of a global temperature such as our spliced GRIP-Vostok record is common in reconstruction of ice sheet glaciation cycles, and indeed provides an important timing and amplitude to whatever climate parameterization is used.Since we are fitting our modeled ice sheets to geologically constrained areal footprints, we regard it as justified to adjust this global temperature proxy with a temporally variable and regional scaling factor, shown in Fig. 4a, in order to recreate ice sheet sizes compatible with the constraints.This was done by systematically adjusting the scaling factor to yield a still stand, or maximum in the areal extent, that matched the margin positions outlined in Fig. 2 within a time range appropriate for the particular stage (we also had targets at MIS 5d and at the Younger Dryas, as well as the 5b, 4, and 2 (LGM) stages reported here).Once a target configuration has been met, the prior temporal variation is frozen, and systematic adjustment from that still stand forward proceeds to the next target configuration.
Deviations from the core proxy do not exceed 1 • C for Eurasia and 1.5 • C for North America.For both North America and Eurasia it was necessary to increase amplitudes by a factor of 1.2 (a decrease in temperature by approximately 1 • C) in the early part of the glacial cycle in order to match the target ice sheet areal footprints.The scaling amplitudes were reduced for both ice sheets during and after stage 4, in order not to overshoot documented stage 2 ice margins (1.4 • C warmer briefly for North America and 0.25 • C warmer for Eurasia almost to the LGM).
The model uses a radially symmetrical mass balance function, with mass balance being a function of elevation and the latitude relative to the climatic "pole" (Fastook and Prentice, 1992).In the mass balance scheme we have employed, the ablation is based on the positive degree day method (Braithwaite and Olsen, 1989), and the accumulation is given as a function of the saturation specific humidity and the local ice sheet slope.In the model, the temperature field has a seasonal cycle whose amplitude depends on distance from the climatic pole and a simple spatial structure: it increases linearly with distance from the climatic pole (0.5 • C degree −1 of latitude) and decreases linearly with height (with a moist adiabatic lapse rate of −5.6 • C km −1 ); i.e., a constant lapse rate is assumed.The two key control parameters of the mass balance scheme are the position of the climatic pole and the global temperature proxy that sets the general climatic state.
Moving the climatic pole to east-central Greenland is instrumental in achieving a good fit to the glacial constraints for both the Eurasian and the North American ice sheets.Figure 4b and c show the changing position of this pole throughout the modeled time interval.Green (Eurasia) and blue (North America) lines show the adjustments in the pole position that were necessary to obtain the fits.The Eurasian pole starts in the north and moves to the south and west, whereas the North American pole begins further south and west and stays stationary except for a slight shift south and east around 70 kyr.It is possible that the southward shift of the mass balance pole position for the Eurasian Ice Sheet reflects the down-stream influence of the North American ice sheets on the atmospheric circulation, an issue we plan to examine in a separate study.This fitting procedure was done in consort with the amplitude fitting procedure described above.

Modeling strategy
By tuning the "global" settings of a scaling factor for temperature and the position of the climatic pole, the latter of which can induce zonal asymmetry in mass balance (see Discussion), we have forced the model through the stages shown in Fig. 2 in runs attempting to match the "target" configurations.Having obtained an acceptable fit in areas where there is control, we accepted the model predictions for ice in areas where geological control is lacking.This is the reason why we have performed time-dependent modeling despite having only two primary target time slices.In addition to predicting the height for the ice sheets for which constraints exist, the model provides tuned and physically based extrapolation to areas for which constraints are unreliable or lacking.We modeled Eurasia and North America separately, attempting to minimize the differences in global settings of the model between the two runs.The overarching consideration was to keep differences between constraints and modeled ice sheet outlines small enough that they are unlikely to influence the results of the intended end use of the data set, atmospheric circulation modeling.
The main criterion for selecting the two runs (North America and Eurasia) was that they achieved a good fit, in ice sheet outline and location, to the available spatial constraints, with no local adjustments ("fixes") applied to achieve this fit.Fixing mismatches between model output and geological constraints through application of local adjustments would corrupt the usefulness of the model in predicting ice in areas where constraints are lacking.If over-or undershooting particular segments of a constraining ice margin was unavoidable, our main strategy was to balance over-and undershoot, i.e., give priority to correct ice sheet size, by adjusting the climatic pole position and temperature scaling.The only hard ice-cutoff function applied was the −200 m bathymetric contour.The goodness of fit in this procedure was judged through direct visual comparison of model-produced ice sheet outlines and geological constraints.An alternative approach would have been to use goodness-of-fit metrics which could then be used as the basis for an automated search algorithm.However, all meaningful metrics we know of require complete targets that can be fully numerically described.For several of our targets only incomplete outlines and first-order dome patterns are available.These incomplete outlines rep-resent real and useful constraining data, but do not allow the use of metrics for e.g.ice sheet area or center of gravity.We regard the present approach as expedient and appropriate for our ultimate goal, which is to provide reasonable first-order large-scale ice sheet reconstructions for use in climate modeling.

Fit to the constraints
Figure 5 shows ice thickness in relation to the constraints, and Fig. 6 shows ice elevations.The main deviations from a "perfect" fit to the geological constraints for stage 5b are the following.There is too little ice in Newfoundland, making the Quebec Dome shorter in the NW-SE direction than the constraints indicate, and too much ice on the western flank of the Quebec Dome.The overall size of the Quebec Dome is in line with the constraints.The model slightly undershoots the minimum size for the proto-Keewatin Dome.In Eurasia, there is too little ice in southern Norway, and the Scandinavian Ice Sheet reaches a somewhat larger easterly extent than indicated by the constraints.The Barents-Kara Dome is well centered within the constraints but slightly too small, particularly in the sector of the Putorana Plateau.In line with the situation in North America, over-and undershoots of the ice margin largely cancel out.
For stage 4, the fit is in general very good.The only significant deviations from the constraints are the somewhat too small ice extent in the Canadian Maritimes, a slightly too small extent of the proto-Keewatin Dome, and some excess ice in the easternmost part of the Eurasian Ice Sheet.The model generates a plausible Cordilleran ice sheet with restricted extent and an easterly located ice divide.

Paleotopographic evolution
In the paleotopographic plots (Fig. 6) we see that the interglacial topography is characterized by two major highlands of a size large enough to significantly affect atmospheric circulation; Rocky Mountains-Coast Range and Greenland.In contrast to the situation in North America and Eurasia, glaciation south of 60 • N is impossible in Greenland.During stage 5b (Fig. 6a), emerging ice sheets constitute four additional highlands large enough to potentially influence atmospheric circulation.These are arranged in pairs in close proximity to each other, the Barents-Kara Ice Sheet and the Scandinavian Ice Sheet in Eurasia, and the Quebec and central Arctic domes in North America.These pairs constitute two highlands, 3500 and 3000 km long, respectively, each with only a minor gap in the central parts.
In contrast to the interglacial situation with two significant and separated N-S oriented highlands located across the westerlies, by stage 5b the number has thus increased to four.In North America, a 1500-2500 km separation exists between the Rocky Mountains and the two ice sheets on the easternmost part of the continent.
By stage 4 (Fig. 6b), the number and location of ice sheets is still the same as during stage 5b, but the Quebec Dome has grown radically, constituting the largest topographic feature in the Northern Hemisphere with the exception of the Tibetan Plateau.The North American and Eurasian ice sheets are by this stage both over 3500 km long, each with a high saddle connecting two domes.The separation between the North American Ice Sheet and the Rocky Mountains has shrunk to few hundred kilometers in the north, and approximately 1500 km further south.
The time separation between stage 4 and the LGM is large, more than 40 kyr, and there is conflicting evidence (Ukkonen, 2003;Helmens and Engels, 2010;Wohlfarth, 2010) regarding the ice extent in both Eurasia and North America during stage 3. Our model indicates a significant shrinkage of ice in Eurasia, but a slow growth in North America.In stage 2 (Fig. 7c) the total ice sheet size and elevation in Eurasia is approximately similar to stage 4, but its footprint has moved 700 km to the southwest.The North American evolution leading into stage 2 is far more dramatic.Here, the two eastern and northern domes have merged with the Cordilleran Ice Sheet, to form a unified Laurentide and Cordilleran ice sheet reaching from coast to coast.There is no longer a midcontinent gap, and consequently the number of obstacles across the westerlies is down to three, although one of them has a very large east-west extent.

Evolution of mass
A surprising observation regarding the evolution of mass for Eurasian and North American ice sheets is the contrast between the behaviors of the two ice sheets.The ice volume curve (Fig. 3b) shows ice volume variation until approximately 70 kyr as being in concert on the two continents, and being of comparable magnitude.However, during stage 4 and stage 2 the volume build-up is far more dramatic in North America than in Eurasia.During stage 5b ice volume in North America is larger than in Europe by a factor of 1.3.During stage 4 this number increases to 2, and by stage 2 North America has four times the Eurasian amount of ice.

Forcing of the model
There are two main ways to force the mass balance of a stand-alone ice sheet model: either by using idealized spatial climatic fields or by using fields obtained by interpolating measured or modeled climate parameters representing endmember climates, typically the LGM and the present (e.g., Charbit et al., 2007;Langen and Vinther, 2008).However, strong zonal misfits (too much ice in the Alaska-E Siberia sector, and too little ice in Quebec, southern Scandinavia and the British Isles) is a common deficiency when simple northpole-centered mass balance schemes are employed, even when present-day zonal climate anomalies are built in (Abe-Ouchi et al., 2007).Still, in order to tune the ice sheet model to the available constraints, we have chosen the first approach because it allows the main zonal climatic asymmetries of the high latitudes to be modeled in a simple way by shifting the "climatic pole".No similar straightforward spatialtuning possibility exists for fixed climate fields, and it is un-clear what realism as representing real climates they retain if spatially manipulated.It should be noted that the glacial ice sheets, by themselves, tended to increase the zonal climatic asymmetries through effects on the surface energy balance and the large-scale atmospheric circulation (e.g., Manabe and Broccoli, 1985;Cook and Held, 1988;Charbit et al., 2007).The tunable zonal climatic anomalies in the present model are admittedly idealized and crude but yet flexible enough to yield a fair match between the geological constraints and the model results.In combination with atmospheric circulation modeling, our reconstructions can be used to examine the zonal climate anomalies that arise from interactions between the ice sheet topography and the atmospheric circulation (Roe and Lindzén, 2001;Liakka and Nilsson, 2010).Such studies should provide additional physical information, helping to refine the present first-order reconstructions of the ice sheet topography.
When comparing with the recent ensemble modeling study by Stokes et al. (2012) we notice that they agree on a firstorder pattern of ice build-up on the northeastern continental rim, but that the two models differ regarding the presence or absence of a robust Quebec Dome during stage 5b.The Stokes et al. model also indicates an early (MIS 4) full confluence between the Cordilleran and Laurentide ice sheets and complete incursion of the Canadian prairies already at that time.The volume evolution of the studies shows considerable similarity, both agreeing on small and comparable ice sheet volumes through stage 5 and a radical build-up of mass during stage 4.However, they differ in that the Stokes et al. (2012) study shows more dramatic stadial-interstadial swings in ice volume.This may be related to their inclusion of surficial and extramarginal (pro-glacial lakes) hydrology in their model, the effect of which will be that calving in such lakes is an additional fast process for ablation during interstadials.

Overall glaciation patterns during the three stadials 5b, 4 and 2
When comparing the evolution of the Eurasian and North American ice sheets through the three stages 5b, 4 and 2, some striking patterns emerge.In northeastern Eurasia, it appears that the ice sheet extent became successively smaller through the three stages, closely mimicking the diminishing amplitude of the corresponding insolation minima.In contrast, the southwestern extent of the ice sheet increased successively through these three events.Hence, the whole footprint of the ice sheet moved to the southwest.In contrast to the North American situation, a clear highland center of glaciation existed in the Scandinavian mountains throughout the glaciation cycle.Only during the LGM did the center of mass move eastwards, to eastern Sweden and the Baltic depression (Ljungner, 1949;Kleman et al., 1997).
The evolution in North America was quite different.Ice cover during stage 5b was likely to have already been established over the entire NE continental rim, from Newfoundland to the Queen Elizabeth Islands.The evolution thereafter was mainly an expansion over the interior continental areas, rapid in the east with a full expansion of the Quebec Dome to the vicinity of the LGM margin during stage 4. The expansion in the west is poorly constrained but was apparently much slower.The southern and central prairies were not infilled before stage 2.This pattern indicates that build-up of the western part of the Laurentide Ice Sheet was precipitation-limited because of poor access of moist air masses to the area.The range in glaciated area and volume through these three stages is much larger in North America than in Eurasia.We note that the volume differences in Eurasia are modest between stages 5b, 4 and 2, and that North America differs radically in response by very rapid increases towards stage 4 and stage 2, respectively.
We speculate that Fennoscandian glaciation, due to the presence of a long high-precipitation backbone which determined the basic build-up pattern, followed the same general growth and decay pattern during all glacial cycles.Such stability in glaciation pattern cannot be expected to have prevailed in North America.Due to the lack of a similar highland backbone, the growing ice sheets rapidly developed dispersal centers far removed from any topographic highlands.Hence, renewed ice growth after interstadials with incomplete deglaciation could occur on a "topography" that radically differed from the fully ice-free interglacial one.Hence, North American glaciation can be expected to have been much more prone to hysteresis effects than glaciation in Eurasia.This inference is to some extent confirmed by geological evidence for ice sheet configurations radically different from any that seem to have existed during the last glacial cycle.Magnetostratigraphy of tills indicate early Quaternary Laurentide ice sheets with a full N-S extent but with restricted ice cover on the prairies as well as in the Atlantic provinces (Barendregt and Irving, 1998).The climatology that would produce such a mid-continent glaciation pattern is to our knowledge completely unexplored.

Reconstructed ice sheet elevations
We expect reconstructions of ice elevations to be more accurate during build-up stages than for LGM and later stages, because during build-up phases the fraction of frozen bed was higher (Stokes et al., 2012), isostatic depression modest and delayed, and ice mostly terrestrial as opposed to marinebased.However, no known geological evidence allows direct evaluation of interior ice sheet elevations.An assessment of model performance in this respect therefore comes down to evaluation of the soundness of the model and the parameter settings and assumptions employed.Because of the relatively coarse grid we used, the model ice streams are rather broad and diffuse, and smaller topographically controlled streams may not be manifest at all, which may lead to an overestimate of marginal ice thickness, the effect of which will lessen inland.Figure 7 includes basal melt rates for the 3 modeled stages, as well as the resulting velocity field.Many major known ice streams are apparent here, albeit broader and less concentrated in their flow than they would be with higherresolution topography.
We note that the Eurasian Ice Sheet, despite its larger spatial extent, never attains elevations as high as the Greenland Ice Sheet, probably because the Greenland Ice Sheet over much of its perimeter has the margin on its relatively high coastal rim.The same holds true for both proto-Laurentide ice sheets in North America.In the model neither of them is as high as the Greenland Ice Sheet.Only at the LGM does our model show a Laurentide Ice Sheet that reaches above 3200 m, with one central dome.We will not here reiterate the long-standing debate over mono-dome vs. multi-dome configuration, but refer to Kleman et al. (2010) for a discussion on this topic.Our model shows an initially two-domed ice sheet that develops into a mono-dome during the LGM.The split into a deglacial two-dome configuration is well documented (Dyke et al., 2002).Hence, there is good evidence that the basic ice sheet topology changed during the glacial cycle, but the geological evidence is relatively powerless for discriminating between e.g. a high or a low saddle connecting the Quebec and Keewatin domes.Our model underestimates the LGM ice sheet extent in westernmost Canada, probably due to the long poleward distance in the model domain when the climatic pole is located in east central Greenland.In terms of area and volume this is approximately compensated by an overshoot of the ice margin on the western prairies.The LGM ice sheet reaches the geologically documented maximum southerly extent.
The total modeled Northern Hemisphere ice volume reaches approximately 100 m sea-level equivalent at LGM, thus allowing for a 20 m contribution from the Southern Hemisphere.Assuming the constraining footprints of Northern Hemisphere ice sheets to be correct, this constitutes a test on the modeled ice sheet height; any serious over-or undershoot of heights would have shown show up as implausible LGM volumes.

Conclusions
Using geologically and geomorphologically constrained numerical ice sheet modeling we have explored the evolution of Northern Hemisphere ice sheets during the build-up phase of the last glacial cycle.The model credibly reproduces the first-order pattern of size and location of geologically indicated ice sheets during marine isotope stages (MIS) 5b and 4. The differences in ice margin outline between constraints and model are considered to be small enough that they will not adversely affect the intended use; atmospheric circulation modeling based on the two paleotopographical data sets, for stages 5b and 4 respectively.
The results show that from an interglacial state of two north-south-oriented obstacles to atmospheric circulation (Rocky Mountains and Greenland), by MIS 5b the emergence of combined Quebec-central Arctic and Scandinavian-Barents-Kara ice sheets had effectively increased the number of such highland obstacles to four.This basic pattern was established already during stage 5d.This number of ice sheets remained constant through MIS 4, with the most significant change being a large increase in area and southerly extent of the Quebec Dome.By MIS 2, the con-fluence of the Laurentide and Cordilleran ice sheets had reduced the number to three, albeit one of them with a very large east-west extent.
The controls on growth pattern appear to have been distinctively dissimilar in Eurasia and North America.In Eurasia, the whole footprint of the ice sheet over the three stages moved in a southwesterly direction, towards the main precipitation source, with only modest and short-lived expansion into the continental interior at the LGM.In North America, the entire northeastern continental rim was glaciated early on, and expansion thereafter was by successive infilling of the continental interior, more rapid in the east than in the west.

Fig. 1 .
Fig. 1.The Terra Nivea and Grinnell ice caps (on the skyline) rise only about 200 m above the level 600 m-elevation plateau of Meta Incognita Peninsula, Baffin Island.The setting is typical of the inception and build-up areas of North American ice sheets in the mainly lowland Laurentide area during the early part of the last glacial cycle.

Fig. 2 .
Fig. 2. The geological and geomorphological constraints (ice margin outlines) that were used in modeling paleotopography for glacial maxima during MIS 5b (heavy solid line) and 4 (thin solid line).LGM ice extent is indicated with stippled line.For source data and construction of outlines, see text.

Fig. 3 .
Fig. 3. (a) The temperature forcing, spliced Vostok and GRIP records, used for the modeling.(b) A time-dependent scaling was used to fit ice sheet size to geological constraints.See text and Fig. 4 for full explanation.(c) Modeled ice volumes for the North American and Eurasian domains of the model, respectively.

Fig. 4 .
Fig. 4. Primary tuning parameters, (a) temperature scaling applied to the spliced Vostok-GRIP record.(b) and (c) Position of climatic pole.These two parameters were the primary controls used to achieve fit in terms of ice sheet size and location to the geologically and geomorphologically based constraints in Fig. 3. See text for a full explanation.