Numerical simulations of the Cordilleran ice sheet through the last glacial cycle

After more than a century of geological research, the Cordilleran ice sheet of North America remains among the least understood in terms of its former extent, volume, and dynamics. Because of the mountainous topography on which the ice sheet formed, geological studies have often had only local or regional relevance and shown such a complexity that ice-sheet-wide spatial reconstructions of advance and retreat patterns are lacking. Here we use a numerical ice sheet model calibrated against field-based evidence to attempt a quantitative reconstruction of the Cordilleran ice sheet history through the last glacial cycle. A series of simulations is driven by time-dependent temperature offsets from six proxy records located around the globe. Although this approach reveals large variations in model response to evolving climate forcing, all simulations produce two major glaciations during marine oxygen isotope stages 4 (62.2–56.9 ka) and 2 (23.2– 16.9 ka). The timing of glaciation is better reproduced using temperature reconstructions from Greenland and Antarctic ice cores than from regional oceanic sediment cores. During most of the last glacial cycle, the modelled ice cover is discontinuous and restricted to high mountain areas. However, widespread precipitation over the Skeena Mountains favours the persistence of a central ice dome throughout the glacial cycle. It acts as a nucleation centre before the Last Glacial Maximum and hosts the last remains of Cordilleran ice until the middle Holocene (6.7 ka).


Introduction
During the last glacial cycle, glaciers and ice caps of the North American Cordillera have been more extensive than today.At the Last Glacial Maximum (LGM), a continuous blanket of ice, the Cordilleran ice sheet (Dawson, 1888), stretched from the Alaska Range in the north to the North Cascades in the south (Fig. 1).In addition, it extended offshore, where it calved into the Pacific Ocean, and merged with the western margin of its much larger neighbour, the Laurentide ice sheet, east of the Rocky Mountains.
More than a century of exploration and geological investigation of the Cordilleran mountains have led to many observations in support of the former ice sheet (Jackson and Clague, 1991).Despite the lack of documented end moraines offshore, in the zone of confluence with the Laurentide ice sheet and in areas swept by the Missoula floods (Carrara et al., 1996), moraines that demarcate the northern and southwestern margins provide key constraints that allow reasonable reconstructions of maximum ice sheet extents (Prest et al., 1968;Clague, 1989, Fig. 1.12;Duk-Rodkin, 1999;Booth et al., 2003;Dyke, 2004).As indicated by field evidence from radiocarbon dating (Clague et al., 1980;Clague, 1985Clague, , 1986;;Porter and Swanson, 1998;Menounos et al., 2008), cosmogenic exposure dating (Stroeven et al., 2010(Stroeven et al., , 2014;;Margold et al., 2014), bedrock deformation in response to former ice loads (Clague and James, 2002;Clague et al., 2005), and offshore sedimentary records (Cosma et al., 2008;Davies et al., 2011), the LGM Cordilleran ice sheet extent was short-lived.However, former ice thicknesses and, there-  (Dyke, 2004, red line) and the modelling domain used in this study (black rectangle).The background map consists of ETOPO1 (Amante and Eakins, 2009) and Natural Earth Data (Patterson and Kelso, 2015).fore, the ice sheet's contribution to the LGM sea-level lowstand (Carlson and Clark, 2012;Clark and Mix, 2002) remain uncertain.
Our understanding of the Cordilleran glaciation history prior to the LGM is even more fragmentary (Barendregt and Irving, 1998;Kleman et al., 2010;Rutter et al., 2012), although it is clear that the Pleistocene maximum extent of the Cordilleran ice sheet predates the last glacial cycle (Hidy et al., 2013).In parts of Yukon and Alaska, and in the Puget Lowland, the distribution of tills (Turner et al., 2013;Troost, 2014) and dated glacial erratics (Ward et al., 2007(Ward et al., , 2008;;Briner and Kaufman, 2008;Stroeven et al., 2010Stroeven et al., , 2014) ) indicate an extensive marine oxygen isotope (MIS) stage 4 glaciation.Landforms in the interior regions include flow sets that are likely older than the LGM (Kleman et al., 2010, Fig. 2), but their absolute age remains uncertain.
In contrast, evidence for the deglaciation history of the Cordilleran ice sheet since the LGM is considerable, albeit mostly at a regional scale.Geomorphological evidence from south-central British Columbia indicates a rapid deglaciation, including an early emergence of elevated areas while thin, stagnant ice still covered the surrounding lowlands (Fulton, 1967(Fulton, , 1991;;Margold et al., 2011Margold et al., , 2013b)).This model, although credible, may not apply in all areas of the Cordilleran ice sheet (Margold et al., 2013a).Although solid evidence for late-glacial glacier readvances has been found in the Coast, Columbia, and Rocky mountains (Reasoner et al., 1994;Osborn and Gerloff, 1997;Clague et al., 1997;Friele and Clague, 2002a, b;Kovanen, 2002;Kovanen and Easterbrook, 2002;Lakeman et al., 2008;Menounos et al., 2008), it appears to be sparser than for formerly glaciated regions surrounding the North Atlantic (e.g.Sissons, 1979;Lundqvist, 1987;Ivy-Ochs et al., 1999;Stea et al., 2011).Nevertheless, recent oxygen isotope measurements from Gulf of Alaska sediments reveal a climatic evolution highly correlated to that of Greenland during this period, including a distinct Late Glacial cold reversal between 14.1 and 11.7 ka (Praetorius and Mix, 2014).
In general, the topographic complexity of the North American Cordillera and its effect on glacial history have inhibited the reconstruction of ice-sheet-wide glacial advance and retreat patterns such as those available for the Fennoscandian and Laurentide ice sheets (Dyke and Prest, 1987;Boulton et al., 2001;Dyke et al., 2003;Kleman et al., 1997Kleman et al., , 2010;;Stroeven et al., 2015).Here, we use a numerical ice sheet model (the PISM authors, 2015), calibrated against fieldbased evidence, to perform a quantitative reconstruction of the Cordilleran ice sheet evolution through the last glacial cycle and analyse some of the long-standing questions related to its evolution: Asthenosphere  2) with parameters from Table 3.
Although numerical ice sheet modelling has been established as a useful tool to improve our understanding of the Cordilleran ice sheet (Jackson and Clague, 1991, p. 227;Robert, 1991;Marshall et al., 2000), the ubiquitously mountainous topography of the region has presented two major challenges to its application.First, only recent developments in numerical ice sheet models and underlying scientific computing tools (Bueler and Brown, 2009;Balay et al., 2015) have allowed for high-resolution numerical modelling of glaciers and ice sheets on mountainous terrain over millennial timescales (e.g.Golledge et al., 2012).Second, the complex topography of the North American Cordillera also induces strong geographic variations in temperature and precipitation (Jarosch et al., 2012), thus requiring the use of high-resolution climate forcing fields as an input to an ice sheet model (Seguinot et al., 2014).However, evolving climate conditions over the last glacial cycle are subject to considerable uncertainty and still lie beyond the computational reach of atmosphere circulation models.
Our palaeoclimate forcing therefore includes spatial temperature and precipitation grids derived from a present-day atmospheric reanalysis (Mesinger et al., 2006).These were previously tested against observational data and shown to best reproduce the steep precipitation gradients previously identified as necessary to model the LGM extent of the Cordilleran ice sheet in agreement with its geological imprint among four atmospheric reanalyses available over the study area (Seguinot et al., 2014).To mimic climate evo-lution through the last glacial cycle, these grids are simply supplemented by lapse-rate corrections and temperature offset time series.The latter are obtained by scaling six different palaeotemperature reconstructions from proxy records around the globe, including two oxygen isotope records from Greenland ice cores (Dansgaard et al., 1993;Andersen et al., 2004), two oxygen isotope records from Antarctic ice cores (Petit et al., 1999;Jouzel et al., 2007), and two alkenone unsaturation index records from northwestern Pacific ocean sediment cores (Herbert et al., 2001).
Although these proxy records were all obtained outside the model domain, more regional palaeotemperature reconstructions spanning over the last glacial cycle are lacking.For instance, the Mount Logan ice core oxygen isotope record covers only the last 30 000 years (30 ka) and has been interpreted as a proxy for source region rather than for palaeotemperature (Fisher et al., 2004(Fisher et al., , 2008)).Seasurface temperatures have been reconstructed offshore Vancouver Island from alkenone unsaturation indices over the last 16 14 C cal ka (Kienast and McKay, 2001), and from the Mg / Ca ratio in planktonic foraminifera over the period from 10 to ca.50 14 C cal ka (Taylor et al., 2014(Taylor et al., , 2015)), but these records cover only parts of the last glacial cycle.
After testing the model sensitivity to these climate forcing time series (Sect.3) and to some of the most influential ice flow parameters (Sect.4), we then proceed to compare the model output to geological evidence and discuss the timing and extent of glaciation and the patterns of deglaciation (Sect.5).

Overview
The simulations presented here were run using the Parallel Ice Sheet Model (PISM, version 0.7.2), an open source, finite difference, shallow ice sheet model (the PISM authors, 2015).The model requires input on basal topography, sea level, geothermal heat flux, and climate forcing.It computes the evolution of ice extent and thickness over time, the thermal and dynamic states of the ice sheet, and the associated lithospheric response.
Basal topography is bilinearly interpolated from the ETOPO1 combined topography and bathymetry data set with a resolution of 1 arcmin (Amante and Eakins, 2009) to the model grids.It responds to ice load following a bedrock deformation model that includes local isostasy, elastic lithosphere flexure, and viscous asthenosphere deformation in an infinite half-space (Lingle and Clark, 1985;Bueler et al., 2007).A relatively low viscosity value of ν m = 1 × 10 19 Pa s is used for the asthenosphere (Table 1) in accordance with the results from regional glacial isostatic adjustment modelling at the northern Cascadia subduction zone (James et al., 2009).Sea level is lowered as a function of time based on the Table 2. Palaeotemperature proxy records and scaling parameters yielding temperature offset time series used to force the ice sheet model through the last glacial cycle (Fig. 5).f corresponds to the scaling factor adopted to yield last glacial maximum ice limits in the vicinity of mapped end moraines, and [ T TS ] 22 32 refers to the resulting mean temperature anomaly during the period 32 to 22 ka after scaling.A  1. Basal sliding follows a pseudo-plastic law where the yield stress accounts for till deconsolidation under high water pressure (Sect.2.3).Ice shelf calving is simulated using an ice thickness threshold and principle components of the strain rate tensor (Sect.2.4).Surface mass balance is computed using a positive degree-day (PDD) model (Sect.2.5).Climate forcing is provided by a monthly climatology averaged from 1979 to 2000 from the North American Regional Reanalysis (NARR, Mesinger et al., 2006), perturbed by time-dependent offsets and lapse-rate temperature corrections (Sect.2.6, Table 2).A sensitivity study to some of the most influential ice rheology (Sect 2.2) and basal sliding (Sect.2.3) parameters is performed (Table 3).
Each simulation starts from assumed ice-free conditions at 120 ka and runs to the present.Our modelling domain of 1500 by 3000 km encompasses the entire area covered by the Cordilleran ice sheet at the LGM (Fig. 1).The simulations were run on two distinct grids, using a lower horizontal resolution of 10 km and a higher horizontal resolution of 5 km.

Ice rheology
Ice sheet dynamics are typically modelled using a combination of internal deformation and basal sliding.PISM is a shallow ice sheet model, which implies that the balance of stresses is approximated based on their predominant components.The shallow shelf approximation (SSA) is combined with the shallow ice approximation (SIA) by adding velocity solutions of the two approximations (Winkelmann et al., 2011, Eqs. 7-9 and 15).Although this heuristic approach implies errors in the transition zone where gravitational stresses intervene both in the SIA and SSA velocity computation, this hybrid scheme is computationally much more efficient than a fully three-dimensional model with which the simulations presented here would not be feasible.
Ice deformation is governed by the constitutive law for ice (Glen, 1952;Nye, 1953): where ˙ is the strain-rate tensor, τ the deviatoric stress tensor, and τ e the effective stress defined in our case by τ e 2 = 1 2 tr(τ 2 ).
where T pa is the pressure-adjusted ice temperature calculated using the Clapeyron relation, T pa = T − βp.R = 8.31441 J mol −1 K −1 is the ideal gas constant, and A c , A w , Q c , and Q w are constant parameters corresponding to values measured below and above a critical temperature threshold T c = −10 • C (Paterson and Budd, 1982;Cuffey and Paterson, 2010, p. 72).The water fraction, ω, is capped at a maximum value of 0.01, above which no measurements are available (Lliboutry and Duval, 1985;Greve, 1997, Eq. 5.7).Finally, E is a non-dimensional enhancement factor which can take different values: E SIA in the SIA and E SSA in the SSA.In all our simulations, we set constant the power-law exponent, n = 3, according to Cuffey and Paterson (2010, p. 55-57), the Clapeyron constant, β = 7.9 × 10 −8 K Pa −1 , according to Lüthi et al. (2002), the water fraction coefficient, f = 181.25,according to Lliboutry and Duval (1985), and the SSA enhancement factor, E SSA = 1, according to Cuffey and Paterson (2010, p. 77).These fixed parameter values are summarised in Table 1.
Contrarily, we test the model sensitivity (Sect.4) to different values for the two creep parameters, A c and A w , the two activation energies, Q c and Q w , and the SIA enhancement factor, E SIA , as follows.
-Our default configuration used in the control run of the sensitivity study and all other simulations in the paper includes rheological parameters, A c , A w , Q c , and Q w , derived from Paterson and Budd (1982) and given in Bueler and Brown (2009, Eq. 5), and E SIA = 1.
-Our hard ice configuration includes rheological parameters, A c , A w , Q c , and Q w , derived from Cuffey and Paterson (2010, p. 72 and 76), and E SIA = 1, which correspond to a stiffer rheology than that used in the control run.
-Our soft ice configuration includes rheological parameters from Cuffey and Paterson ( 2010) and E SIA = 5, the recommended value for ice age polar ice (Cuffey and Paterson, 2010, p. 77).
An additional simulation using the ice rheology from Cuffey and Paterson (2010) and E SIA = 2, the recommended value for Holocene polar ice (Cuffey and Paterson, 2010, p. 77), has been performed but is not presented here, since its results were very similar to that of our default run.Actual parameter values for A c , A w , Q c , Q w , and E SIA used in our simulations are given in Table 3, while the effect of the three different parametrisations on temperaturedependent ice softness, A, is illustrated in Fig. 2. Surface air temperature derived from the climate forcing (Sect.2.6) provides the upper boundary condition to the ice enthalpy model.Temperature is computed in the ice and in the bedrock to a depth of 3 km below the ice-bedrock interface, where it is conditioned by a lower boundary geothermal heat flux of q G = 70 mW m −2 .Although this uniform value does not account for the high spatial geothermal variability in the region (Blackwell and Richards, 2004), it is, on average, representative of available heat flow measurements.In the low-resolution simulations, the vertical grid consists of 31 temperature layers in the bedrock and up to 51 enthalpy layers in the ice sheet, corresponding to a vertical resolution of 100 m.The high-resolution simulations use 61 bedrock layers and up to 101 ice layers with a vertical resolution of 50 m.

Basal sliding
A pseudo-plastic sliding law, relates the bed-parallel shear stresses, τ b , to the sliding velocity, v b .The yield stress, τ c , is modelled using the Mohr-Coulomb criterion: where cohesion, c 0 , is assumed to be 0. Effective pressure, N , is related to the ice overburden stress, P 0 = ρgh, and the modelled amount of subglacial water, using a formula derived from laboratory experiments with till extracted from the base of Ice Stream B in West Antarctica (Tulaczyk et al., 2000;Bueler and van Pelt, 2015): where δ sets the minimum ratio between the effective and overburden pressures, e 0 is a reference void ratio, and C c is the till compressibility coefficient (Tulaczyk et al., 2000).
The amount of water at the base, W , varies from 0 to W max , a threshold above which additional meltwater is assumed to drain off instantaneously.In all our simulations, we set constant the pseudo-plastic sliding exponent, q = 0.25, and the threshold velocity, v th = 100 m a −1 , according to values used by Aschwanden et al. (2013), the till cohesion, c 0 = 0, whose measured values are consistently negligible (Tulaczyk et al., 2000;Cuffey and Paterson, 2010, p. 268), the till reference void ratio, e 0 = 0.69, and the till compressibility coefficient, C c = 0.12, according to the only measurements available to our knowledge, published by Tulaczyk et al. (2000).These fixed parameter values are summarised in Table 1.
We also use a constant spatial distribution of the till friction angle, φ, whose values vary from 15 to 45 • as a piecewise-linear function of modern bed elevation, with The Cryosphere, 10, 639-664, 2016 the lowest value occurring below the modern sea level (0 m above sea level, m a.s.l.) and the highest value occurring above the generalised elevation of the highest shorelines (200 m a.s.l.; Clague, 1981, Fig. 5).This range of values span over the range of measured values for glacial till of 18 to 40 • (Cuffey and Paterson, 2010, p. 268).It accounts for frictional basal conditions associated with discontinuous till cover at high elevations and for a weakening of till associated with the presence of marine sediments (cf.Martin et al., 2011;Aschwanden et al., 2013, Supplement;the PISM authors, 2015).
An additional simulation with a constant till friction angle, φ = 30 • , corresponding to the average value in Cuffey and Paterson (2010, p. 268), has been performed but is not presented here, since the induced variability was small.
Additionally, we test the model sensitivity (Sect.4) to different values for the minimum ratio between the effective and overburden pressures, δ, and the maximum water height in the till, W max , as follows.
-Our default configuration used in the control run of the sensitivity study and all other simulations in the paper includes δ = 0.02 and W max = 2 m as in Bueler and van Pelt (2015).
-Our soft bed configuration uses δ = 0.01 and W max = 1 m.
-Our hard bed configuration uses δ = 0.05 and W max = 5 m.
The effect of the three different parametrisations on the effective pressure on the till, N, in response to water content, W , is illustrated in Fig. 3.All parameter choices are listed in Table 3.

Ice shelf calving
Ice shelf calving is computed using a double criterion.First, a physically realistic calving flux is computed based on eigenvalues of the horizontal strain rate tensor (Winkelmann et al., 2011;Levermann et al., 2012).This allows floating ice to advance in confined embayments but prevents the formation of extensive ice shelves in the open ocean.Second, floating ice thinner than 50 m is systematically calved off.A subgrid scheme by Albrecht et al. (2011) allows for a continuous migration of the calving front.This formulation of calving has been applied to the Antarctic ice sheet and has shown to produce a realistic calving front position for many of the present-day ice shelves (Martin et al., 2011).

Surface mass balance
Ice surface accumulation and ablation are computed from monthly mean near-surface air temperature, T m , monthly standard deviation of near-surface air temperature, σ , and monthly precipitation, P m , using a temperature-index model (e.g.Hock, 2003).Accumulation is equal to precipitation when air temperatures are below 0 • C and decreases to 0 linearly with temperatures between 0 and 2 • C. Ablation is computed from PDD, defined as an integral of temperatures above 0 • C in 1 year.
The PDD computation accounts for stochastic temperature variations by assuming a normal temperature distribution of standard deviation σ around the expected value T m .It is expressed by an error-function formulation (Calov and Greve, 2005), which is numerically approximated using week-long subintervals.In order to account for the effects of spatial and seasonal variations of temperature variability (Seguinot, 2013), σ is computed from NARR daily temperature values from 1979 to 2000 (Mesinger et al., 2006), including variability associated with the seasonal cycle, and bilinearly interpolated to the model grids (Fig. 4).Degree-day factors for snow and ice melt are derived from mass-balance measurements on contemporary glaciers from the Coast Mountains and Rocky Mountains in British Columbia (Table 1; Shea et al., 2009).
The present-day monthly climatology was bilinearly interpolated from near-surface air temperature and precipitation rate fields from the NARR, averaged from 1979 to 2000.Modern climate of the North American Cordillera is characterised by strong geographic variations in temperature seasonality, timing of the maximum annual precipitation, and daily temperature variability (Fig. 4).Although the ability of the NARR to reproduce the steep climatic gradients is limited by its spatial resolution of 32 km (Jarosch et al., 2012), it has been tested against observational data in our previous sensitivity study and identified as yielding a closer fit between the modelled LGM extent of the Cordilleran ice sheet and the geological evidence than other atmospheric reanalyses (Seguinot et al., 2014).Temperature offset time series, T TS , are derived from palaeotemperature proxy records from the Greenland Ice Core Project (GRIP; Dansgaard et al., 1993), the North Greenland Ice Core Project (NGRIP; Andersen et al., 2004), the European Project for Ice Coring in Antarctica (EPICA; Jouzel et al., 2007), the Vostok ice core (Petit et al., 1999), and Ocean Drilling Program (ODP) sites 1012 and 1020, both located off the coast of California (Herbert et al., 2001).Palaeotemperature anomalies from the GRIP and NGRIP records were calculated from oxygen isotope (δ 18 O) measurements using a quadratic equation (Johnsen et al., 1995), while temperature reconstructions from Antarctic and oceanic cores were provided as such.For each proxy record    2) used as palaeoclimate forcing for the ice sheet model (top panel), and modelled sea-level relevant ice volume (bottom panel) through the last 120 ka, expressed in metres of sea-level equivalent (m s.l.e.).Grey fields indicate marine oxygen isotope stage (MIS) boundaries for MIS 2 and MIS 4 according to a global compilation of benthic δ 18 O records (Lisiecki and Raymo, 2005).Hatched rectangles highlight the time-volume span for ice volume extremes corresponding to MIS 4 (61.9-56.5 ka), MIS 3 (53.0-41.3ka), and MIS 2 (LGM, 23.2-16.8ka).Dotted lines correspond to GRIP-and EPICA-driven 5 km resolution runs.used and each of the parameter set-up used in the sensitivity tests, palaeotemperature anomalies were scaled linearly (Tables 2 and 3) in order to simulate comparable ice extents at the LGM (Table 4) and realistic outlines (Fig. 6).
Finally, lapse-rate corrections, T LR , are computed as a function of ice surface elevation, s, using the NARR surface geopotential height invariant field as a reference topography, b ref : thus accounting for the evolution of ice thickness, h = s − b, on the one hand, and for differences between the basal topography of the ice flow model, b, and the NARR reference topography, b ref , on the other hand.All simulations use an annual temperature lapse rate of γ = 6 K km −1 .In the rest of this paper, we refer to different model runs by the name of the proxy record used for the palaeotemperature forcing.
3 Sensitivity to climate forcing time series

Evolution of ice volume
Despite large differences in the input climate forcing (Fig. 5, upper panel), model output presents consistent features that can be observed across the range of forcing data used.In all simulations, modelled sea-level relevant ice volumes remain relatively low during most of the glacial cycle, except during two major glacial events which occur between 62.2 and 56.9 ka during MIS 4 and between 23.2 and 16.9 ka during MIS 2 (Fig. 5, lower panel).An ice volume minimum is consistently reached between 55.9 and 42.9 ka during MIS 3.However, the magnitude and precise timing of these three events depend significantly on the choice of proxy record used to derive a time-dependent climate forcing (Table 4).Simulations forced by the Greenland ice core palaeotemperature records (GRIP, NGRIP) produce the highest variability in modelled ice volume throughout the last glacial cycle.In contrast, simulations driven by oceanic (ODP 1012, ODP 1020) and Antarctic (EPICA, Vostok) palaeotemperature records generally result in lower ice volume variability throughout the simulation length, resulting in lower modelled ice volumes during MIS 4 and larger ice volumes during MIS 3. The NGRIP climate forcing is the only one that results in a larger ice volume during MIS 4 than during MIS 2.
While simulations driven by the GRIP and the two Antarctic palaeotemperature records attain a last ice volume maximum between 19.1 and 16.9 ka, those informed by the NGRIP and the two oceanic palaeotemperature records attain their maximum ice volumes thousands of years earlier.Moreover, the ODP 1012 run yields a rapid deglaciation of the modelled grounded ice extent prior to 17 ka.The ODP 1020 simulation predicts an early maximum in ice volume at 20.4 ka, followed by slower deglaciation than modelled using the other palaeotemperature records.Finally, Table 5. Extremes in Cordilleran ice sheet grounded ice extent and sea-level relevant ice volume corresponding to MIS 4, 3, and 2 using the GRIP palaeoclimate forcing with each parameter configuration (Fig. 3).Relative differences (R. diff.)give rough error estimates related to varying selected ice rheology and basal sliding parameters ( whereas model runs forced by the Antarctic palaeotemperature records result in a rapid and uninterrupted deglaciation after the LGM, the simulation driven by the GRIP palaeotemperature record also results in a rapid deglaciation but in three steps, separated by two periods of ice sheet regrowth (Fig. 5).

Extreme configurations
Despite large differences in the timing of attained volume extrema (   (Lisiecki and Raymo, 2005).attained during MIS 2 (Table 4; Fig. 6, upper panels).Corresponding maximum sea-level potentials also differ significantly between model runs and vary between 3.66 and 8.71 m sea-level equivalents (m s.l.e.; Table 4).
During MIS 3 ice volume minimum reconstructions, a central ice cap persists over the Skeena Mountains (Fig. 6, middle panels).Although this ice cap is present in all simulations, its dimensions depend sensitively on the choice of the applied palaeotemperature record.Modelled sea-level relevant ice volume minima spread over a wide range between 1.54 and 2.81 m s.l.e.(Table 4).
Modelled ice sheet geometries during the LGM (MIS 2; Fig. 6, lower panels) invariably include a ca.1500 km long central divide above 3000 m a.s.l.located along the spine of the Rocky Mountains.As an indirect result of the choice of scaling factors applied to different palaeotemperature proxy records (Table 2), modelled maximum sea-level potentials also fall within a tight range of 8.56 to 8.75 m s.l.e.(Table 4).

Sensitivity to ice flow parameters
Using the GRIP ice core palaeotemperature record as climate forcing time series, the model shows a significant sensitivity to selected ice rheological and basal sliding parameters (Fig. 7; Table 5) in terms of modelled grounded ice extent and even more so in terms of modelled sea-level relevant ice volumes.
As a direct result of the different scaling factors applied (Table 3), the resulting simulations show very little difference in the modelled glaciated areas corresponding to the maximum stage during MIS 2, as well as during MIS 4 (Table 5).However, this cannot be said of the modelled glaciated area corresponding to the minimum stage during MIS 3. In fact, the extent of the remnant ice cap which persists over the Skeena Mountains during this stage shows a significant sensitivity to ice rheology of 25 % and an even higher sensitivity to basal sliding parameters of 57 % (Table 5).
The modelled sea-level potentials show a stronger variability than the modelled glaciated areas (Table 5, Fig. 7).As one could expect, softer ice and weaker till both result in a thinner ice sheet, while harder ice and stronger till result in a thicker ice sheet.For instance, the modelled peak sea-level relevant ice volume during MIS 2 (LGM) varies by 25 % between the two parametrisations of ice rheology used and by 14 % between the two parametrisations of basal sliding used.The differences in ice volume are greatest during MIS 3 (Table 5, Fig. 7) where both the areal and thickness contributions add up.

Comparison with the geologic record
Large variations in the model responses to evolving climate forcing reveal its sensitivity to the choice of palaeotemperature proxy record.To distinguish between different records, geological evidence of former glaciations provides a basis for validation of our runs, while the results from numerical modelling can perhaps help to analyse some of the complexity of this evidence.In this section, we compare model outputs to the geologic record, in terms of timing and configuration of the maximum stages, location and lifetime of major nucleation centres, and patterns of ice retreat during the last deglaciation.

Timing of glaciation
Independently of the palaeotemperature records used to force the ice sheet model, our simulations consistently produce two glacial maxima during the last glacial cycle.The first maximum configuration is obtained during MIS 4 (62.2-56.9ka) and the second during MIS 2 (23.2-16.9ka; Figs. 5, 6; Table 4).These events broadly correspond in timing to the Gladstone (MIS 4) and McConnell (MIS 2) glaciations documented by geological evidence for the northern sector of the Cordilleran ice sheet (Duk-Rodkin et al., 1996;Ward et al., 2007;Stroeven et al., 2010Stroeven et al., , 2014) ) and to the Fraser Glaciation (MIS 2) documented for its southern sector (Porter and Swanson, 1998;Margold et al., 2014).There is stratigraphical evidence for an MIS 4 glaciation in British Columbia www.the-cryosphere.net/10/639/2016/The Cryosphere, 10, 639-664, 2016 (Clague and Ward, 2011) and in the Puget Lowland (Troost, 2014), but their extent and timing are still highly conjectural (perhaps MIS 4 or early MIS 3; e.g.Cosma et al., 2008).The exact timing of modelled MIS 2 maximum ice volume depends strongly on the choice of applied palaeotemperature record, which allows for a more in-depth comparison with geological evidence for the timing of the maximum Cordilleran ice sheet extent.In the Puget Lowland (Fig. 1), the LGM advance of the southern Cordilleran ice sheet margin has been constrained by radiocarbon dating on wood between 17.4 and 16.4 14 C cal ka (Porter and Swanson, 1998).These dates are consistent with radiocarbon dates from the offshore sedimentary record, which reveals an increase of glaciomarine sedimentation between 19.5 and 16.2 14 C cal ka (Cosma et al., 2008;Taylor et al., 2014).Radiocarbon dating of the northern Cordilleran ice sheet margin is much less constrained but straddles presented constraints from the southern margin.However, cosmogenic exposure dating places the timing of the maximum northern ice sheet margin extent during the McConnell glaciation close to 17 10 Be ka (Stroeven et al., 2010(Stroeven et al., , 2014)).A sharp transition in the sediment record of the Gulf of Alaska indicates a retreat of regional outlet glaciers onto land at 14.8 14 C cal ka (Davies et al., 2011).
Among the simulations presented here, only those forced with the GRIP, EPICA, and Vostok palaeotemperature records yield Cordilleran ice sheet maximum extents that may be compatible with these field constraints (Fig. 5, lower panel; Table 4).Simulations driven by the NGRIP, ODP 1012, and ODP 1020 palaeotemperature records, on the contrary, yield MIS 2 maximum Cordilleran ice sheet volumes that pre-date field-based constraints by several thousands of years (about 6, 6, and 4 ka respectively).Concerning the simulations driven by oceanic records, this early deglaciation is caused by an early warming present in the alkenone palaeotemperature reconstructions (Fig. 5, upper panel; Herbert et al., 2001, Fig. 3).However, this early warming is a local effect, corresponding to a weakening of the California current (Herbert et al., 2001).The California current, driving cold waters southwards along the south-western coast of North America, has been shown to have weakened during each peak of global glaciation (in SPECMAP) during the past 550 ka, including the LGM, resulting in paradoxically warmer sea-surface temperatures at the locations of the ODP 1012 and ODP 1020 sites (Herbert et al., 2001).
Because most of the marine margin of the Cordilleran ice sheet terminated in a sector of the Pacific Ocean unaffected by variations in the California current, it probably remained insensitive to this local phenomenon.However, the above paradox illustrates the complexity of ice-sheet feedbacks on regional climate and demonstrates that, although located in the neighbourhood of the modelling domain, the ODP 1012 and ODP 1020 palaeotemperature records cannot be used as a realistic forcing to model the Cordilleran ice sheet through the last glacial cycle.Similarly, the simulation using the NGRIP palaeotemperature record depicts an early onset of deglaciation (Fig. 5) following its last glacial volume maximum (22.9 ka, Table 4) attained about 6 ka earlier than dated evidence of the LGM advance.There is a fair agreement between the EPICA and Vostok palaeotemperature records, resulting in only small differences between the simulations driven by those records.These differences are not subject to further analysis here; instead we focus on simulations forced by palaeotemperature records from the GRIP and EPICA ice cores that appear to produce the most realistic reconstructions of regional glaciation history, yet bearing significant disparities in model output.To allow for a more detailed comparison against the geological record, these two simulations were re-run using a higher-resolution grid (Sect.2; Fig. 5, lower panel, dotted lines).A supplementary animation presenting their results is available online.

Ice configuration during MIS 2
During maximum glaciation, both simulations position the main meridional ice divide over the western flank of the Rocky Mountains (Figs. 6,lower panels and 8).This result appears to contrast with palaeoglaciological reconstructions for central and southern British Columbia with ice divides in a more westerly position, over the western margin of the Interior Plateau (Ryder et al., 1991;Stumpf et al., 2000;Kleman et al., 2010;Clague and Ward, 2011;Margold et al., 2013b).These indicate that a latitudinal saddle connected ice-dispersal centres in the Columbia Mountains with the main ice divide (Ryder et al., 1991;Kleman et al., 2010;Clague and Ward, 2011;Margold et al., 2013b).A latitudinal saddle does indeed feature in our modelling results, however, in an inverse configuration between the main ice divide over the Columbia Mountains and a secondary divide over the southern Coast Mountains (Fig. 8).
Such deviation from the geological inferences could reflect the fact that the NARR has difficulties with simulating orographic processes in some areas of steep topography (Jarosch et al., 2012).In our previous study (Seguinot et al., 2014), we have evaluated an applicability of the climate forcing bilinearly interpolated from the NARR to constant-climate simulations of the Cordilleran ice sheet during the LGM against that of an observation-based data set (WorldClim, Hijmans et al., 2005).Indeed, the use of NARR in these simulations produced slightly different patterns of glaciation relative to WorldClim, including a more extensive ice cover on the Columbia and Rocky mountains (Seguinot et al., 2014, Figs. 6-7).Our simulations have shown that these differences are mainly caused by disparities between the precipitation fields of the two data sets (Seguinot et al., 2014, Figs. 13-14).It has been shown that over the southern part of our model domain that temperature and precipitation downscaling can potentially address these limitations (e.g.Jarosch et al., 2012).However, extending this downscaling method to the entire model domain used in our study is challenging, because the northern part of the model domain is characterized by stronger precipitation gradients (Fig. 4) and fewer weather stations than the southern part where the previous analysis has been performed (Hijmans et al., 2005).
Since the model does not include feedback mechanisms between the ice sheet topography and the regional climate, the modelled easterly positions of the ice divide and eastern ice sheet margin may also be sensitive to the assumption of fixed modern-day spatial patterns of air temperature and precipitation.In fact, it is reasonable to think that the cooling during the last glacial cycle was greater inland than near the coast, prohibiting melt at the eastern margin.However, our simulations already produce an excess of ice inland.Including such a temperature continentality gradient in the model while keeping the precipitation pattern constant would thus cause an even greater mismatch between the model results and the geologically reconstructed ice margins during the LGM.
Consequently, the mismatch between the modelled and reconstructed LGM ice margins is likely due to the assumption of the fixed modern-day precipitation patterns rather than the assumption of the fixed modern-day temperature patterns.Firstly, during the build-up phase preceding the LGM, rapid accumulation over the Coast Mountains enhanced the topographic barrier formed by these mountain ranges, which likely resulted in a decrease of precipitation and, therefore, a decrease of accumulation in the interior.Secondly, latent warming of the moisture-depleted air parcels flowing over this enhanced topography could have resulted in an inflow of potentially warmer air over the eastern flank of the ice sheet, thereby counterbalancing the potential continentality gradient discussed above through increasing melt along the advancing margin (cf.Langen et al., 2012).Because these two processes, both with a tendency to limit ice-sheet growth, are absent from our model, the eastern margin of the ice sheet and the position of the main meridional ice divide are certainly biased towards the east in our simulations (Seguinot et al., 2014).
However, field-based palaeoglaciological reconstructions have struggled to reconcile the more westerly centred ice divide in south-central British Columbia with evidence in the Rocky Mountains and beyond that indicates the Cordilleran ice sheet invaded the western Interior Plains, where it merged with the south-western margin of the Laurentide ice sheet and was deflected to the south (Jackson et al., 1997;Bednarski and Smith, 2007;Kleman et al., 2010;Margold et al., 2013a, b).Ice geometries from our model runs do not have this problem, because the position and elevation of the ice divide ensure significant ice drainage across the Rocky Mountains at the LGM (Fig. 8).
During MIS 2, the modelled sea-level potential peaks at 8.62 m s.l.e.(19.1 ka) in the GRIP simulation and at 8.56 m s.l.e.(17.4 ka) in the EPICA simulation.However, these numbers are subject to significant uncertainties in the ice flow parameters embedded in the reference model set-up.The range of parameter values tested in our sensitivity study (Sect.4) yielded relative errors of 25 % with regard to ice rheological parameters and 14 % regarding basal sliding parameters (Fig. 7; Table 5).

Ice configuration during MIS 4
The modelled ice sheet configurations corresponding to ice volume maxima during MIS 4 are more sensitive to the choice of atmospheric forcing than those corresponding to ice volume maxima during MIS 2 (Figs. 6, upper panels, and 9).The GRIP simulation (Fig. 9, left panel) results in a modelled maximum ice sheet extent that closely resembles that obtained during MIS 2, with the only major difference of being slightly less extensive across northern and eastern sectors.In contrast, the EPICA simulation produces a lower ice volume maximum (Fig. 5), which translates in the modelled ice sheet geometry into a significantly reduced southern sector, more restricted ice cover in northern and eastern sectors, and generally lower ice surface elevations in the interior (Fig. 9, right panel).Thus, only the GRIP simulation can explain the presence of MIS 4 glacial deposits in the Puget Lowland (Troost, 2014) and that of ice-rafted debris in the marine sediment record offshore Vancouver Island at ca. 47 14 C cal ka (Cosma et al., 2008).
During MIS 4, the modelled sea-level potential peaks at 7.43 m s.l.e.(57.6 ka) in the GRIP simulation and at 4.84 m s.l.e.(61.9 ka) in the EPICA simulation, corresponding to respectively 86 and 57 % of modelled MIS 2 ice volumes.These estimates show little sensitivity to the range of parameter values tested in our sensitivity study (Sect.4), which yielded relative errors of 12 % to ice rheological parameters and 19 % to basal sliding parameters (Fig. 7; Table 5).

Transient ice sheet states
Palaeoglaciological reconstructions are generally more robust for maximum ice sheet extents and late ice sheet configurations than for intermediate or minimum ice sheet extents and older ice sheet configurations (Kleman et al., 2010).However, these maximum stages are, by nature, extreme configurations, which do not necessarily represent the dominant patterns of glaciation throughout the period of ice cover (Porter, 1989;Kleman and Stroeven, 1997;Kleman et al., 2008Kleman et al., , 2010)).
The resulting maps show that, during most of the glacial cycle, modelled ice cover is restricted to disjoint ice caps centred on major mountain ranges of the North American Cordillera (Fig. 10, blue areas).A 2 500 km long continuous expanse of ice, extending from the Alaska Range in the north-east to the Rocky Mountains in the south-west, is only in operation for at most 34 ka, which is about a third of the timespan of the last glacial cycle (Fig. 10, hatched areas).However, except for its margins in the Pacific Ocean and in the northern foothills of the Alaska Range, the maximum extent of the ice sheet is attained for a much shorter period of time of only few thousand years (Fig. 10, red areas).This result illustrates that the maximum extents of the modelled ice sheet during MIS 4 and MIS 2 were both short-lived and therefore out of balance with contemporary climate.
A notable exception to the transient character of the maximum extent of the Cordilleran ice sheet is the northern slope of the Alaska Range, where modelled glaciers are confined to its foothills during the entire simulation period (Fig. 10, AR).This apparent insensitivity of modelled glacial extent to temperature fluctuations results from a combination of low precipitation, high summer temperature, and large temperature standard deviation in the plains of the Alaska Interior (Fig. 4) which confines glaciation to the foothills of the mountains.This result could potentially explain the local distribution of glacial deposits, which indicates that glaciers flowing on the northern slope of the Alaska Range have remained small throughout the Pleistocene (Kaufman and Manley, 2004).

Major ice-dispersal centres
It is generally believed that the Cordilleran ice sheet formed by the coalescence of several mountain-centred ice caps (Davis and Mathews, 1944).In our simulations, major icedispersal centres, visible on the modelled ice cover duration maps (Fig. 10), are located over the Coast Mountains, the Columbia and Rocky mountains, the Skeena Mountains, and the Selwyn and Mackenzie mountains.
The location of the modelled ice-dispersal centres is potentially biased by the present-day ice volumes contained in the ETOPO1 basal topography data.The most problematic part of the model domain in this respect is that of the Wrangell and Saint Elias mountains, where ice thicknesses of up to 1200 m have been measured by a low-frequency radar (Rignot et al., 2013).In this area, located over the USA-Canada border just north of 60 • N, temperate ice, glacier surging dynamics, and deep subglacial depressions in the ice-field interior pose a fundamental challenge to the reconstruction basal topography for the entire ice cap (Rignot et al., 2013).Recent ice thickness reconstructions available for all glaciers around the globe (Huss and Farinotti, 2012), for Canadian glaciers south of 60 • N (Clarke et al., 2013), and for Alaska tidewater glaciers (McNabb et al., 2015) have not been implemented in our simulations.However, this drawback seem to have little effect on modelled Cordilleran ice sheet dynamics.In fact, the Wrangell and Saint Elias mountains, heavily glacierized at present, host an ice cap for the entire length of both simulations, but that ice cap does not appear to be a major feed to the Cordilleran ice sheet (Fig. 10, WSEM).With this exception of the Wrangell and Saint Elias mountains ice field, present-day ice volumes are small relative to the ice volumes concerned in our study.
Although the Coast, Skeena, Columbia, and Rocky mountains are covered by mountain glaciers for most of the last glacial cycle, providing durable nucleation centres for an ice sheet, this is not the case for the Selwyn and Mackenzie mountains, where ice cover on the highest peaks is limited to a small fraction of the last glacial cycle.In other words, the Selwyn and Mackenzie mountains only appear as a secondary ice-dispersal centre during the coldest periods of the last glacial cycle.The Northern Rocky Mountains (Fig. 10, NRM) do not act as a nucleation centre, but rather as a pinning point for the Cordilleran ice sheet margin coming from the west.
Perhaps the most striking feature displayed by the distributions of modelled ice cover is the persistence of the Skeena Mountains ice cap throughout the entire last glacial period (ca.100-10 ka) and its predominance over the other icedispersal centres (Figs. 6 and 10,SM).Regardless of the applied forcing, this ice cap appears to survive MIS 3 (Fig. 6, middle panels) and serves as a nucleation centre at the onset of the glacial readvance towards the LGM (MIS 2).This situation appear similar to the neighbouring Laurentide ice sheet, for which the importance of residual ice for the glacial history leading up to the LGM has been illustrated by the MIS 3 residual ice bodies in northern and eastern Canada as nucleation centres for a much more extensive MIS 2 configuration (Kleman et al., 2010).
The presence of a Skeena Mountains ice cap during most of the last glacial cycle can be explained by meteorological conditions more favourable for ice growth there than elsewhere.In fact, reanalysed atmospheric fields used to force the surface mass balance model show that high winter precipitations are mainly confined to the western slope of the Coast Mountains, except in the centre of the modelling domain where they also occur further inland than along other east-west transects (Fig. 4).In fact, along most of the northwestern coast of North America, coastal mountain ranges form a pronounced topographic barrier for westerly winds, capturing atmospheric moisture in the form of orographic precipitation and resulting in arid interior lowlands.However, near the centre of our modelling domain, this barrier is less pronounced than elsewhere, allowing westerly winds to carry moisture further inland until it is captured by the extensive Skeena Mountains in north-central British Columbia, thus resulting in a more widespread distribution of winter precipitation (Fig. 4).
The modelled sea-level potential corresponding to these persistent ice-dispersal centres attains a minimum of 1.54 m s.l.e.(42.9 ka) in the GRIP simulation and of 2.55 m s.l.e.(52.4 ka) in the EPICA simulation, corresponding to respectively 18 and 30 % of the MIS 2 ice volumes.However, these numbers should be considered with caution as our sensitivity study (Sect.4) shows that the minimum ice volume during MIS 3 is highly sensitive to ice flow parameters, with relative errors of 45 % to the range of ice rheological parameters tested and 109 % to the range of basal sliding parameters tested (Fig. 7; Table 5).

Erosional imprint on the landscape
A correlation is observed between the modelled duration of warm-based ice cover (Fig. 11) and the degree of glacial modification of the landscape (mainly in terms of the development of deep glacial valleys and troughs).We find evidence for this on the slopes of the Coast Mountains, the Columbia and Rocky mountains, the Wrangell and Saint Elias mountains, and radiating off the Skeena Mountains (Figs. 10 and 11;Kleman et al., 2010, Fig. 2).The Skeena Mountains, for example, indeed bear a strong glacial imprint that indicates ice drainage in a system of distinct glacial troughs emanating in a radial pattern from the centre of the mountain range (Kleman et al., 2010, Fig. 2) that appear to predate the LGM phase of the Cordilleran ice sheet (Stumpf et al., 2000).We suggest that persistent ice cover (Fig. 10) associated with basal ice temperatures at the pressure melting point (Figs. 11 and 12) explains the large-scale glacial erosional imprint on the landscape.A well-developed network of glacial valleys running to the north-west on the west slope of the Selwyn and Mackenzie Mountains (Kleman et al., 2010, Fig. 2;Stroeven et al., 2010, Fig. 8) is modelled to have hosted predominantly warm-based ice (Fig. 12).However, because it is only glaciated for a short fraction of the last glacial cycle in our simulations (Fig. 10), this perhaps indicates that the observed landscape pattern originates from multiple glacial cycles and witnesses an increased relative importance of the Selwyn and Mackenzie mountains ice-dispersal centre (Fig. 10, SMKM), prior to the Late Pleistocene (cf. Ward et al., 2008;Demuro et al., 2012).
The modelled distribution of warm-based ice cover (Figs.11 and 12) is inevitably affected by our assumption of a constant, 70 mW m −2 geothermal heat flux at 3 km depth (Sect.2.2).However, the Skeena Mountains and the area west of the Mackenzie Mountains experience higher-thanaverage geothermal heat flux with measured values of ca.80 and ca. 100 mW m −2 (Blackwell and Richards, 2004).We can therefore expect even longer durations of warm-based ice cover for these areas if we were to include spatially variable geothermal forcing in our Cordilleran ice sheet simulations.

Pace and patterns of deglaciation
Similarly to other glaciated regions, most glacial traces in the North American Cordillera relate to the last few millennia of glaciation, because most of the older evidence has been overprinted by warm-based ice retreat during the last deglaciation (Kleman, 1994;Kleman et al., 2010).From a numerical modelling perspective, phases of glacier retreat are more challenging than phases of growth, because they involve more rapid fluctuations of the ice margin, increased flow velocities and longitudinal stress gradients, and poorly understood hydrological processes.The latter are typically included in the models through simple parametrisations (e.g.Clason et al., 2012Clason et al., , 2014;;Bueler and van Pelt, 2015), if included at all.However, next after the mapping of maximum ice sheet extents during MIS 2 and MIS 4 (Sects.5.1.2and 5.1.3),geomorphologically based reconstructions of patterns of ice sheet retreat during the last deglaciation provide the secondbest source of evidence for the validation of our simulations.
In the North American Cordillera, the presence of lateral meltwater channels at high elevation (Margold et al., 2011(Margold et al., , 2013b(Margold et al., , 2014) ) and abundant esker systems at low elevation (Burke et al., 2012a, b;Perkins et al., 2013;Margold et al., 2013a) indicate that meltwater was produced over large portions of the ice sheet surface during deglaciation.The southern and northern margins of the Cordilleran ice sheet reached their last glacial maximum extent around 17 ka (Sect.5.1.1;Porter and Swanson, 1998;Cosma et al., 2008;Stroeven et al., 2010Stroeven et al., , 2014)), which we take as a limiting age for the onset of ice retreat.The timing of final deglaciation is less well constrained, but recent cosmogenic dates from north-central British Columbia indicate that a seizable ice cap emanating from the central Coast Mountains or the Skeena Mountains persisted into the Younger Dryas chronozone, at least until 12.4 ka (Margold et al., 2014).
In our simulations, the timing of peak ice volume during the LGM and the pacing of deglaciation depend critically on the choice of climate forcing (Table 4;Figs. 5 and 13).Adopting the EPICA climate forcing yields peak ice volume at 17.4 ka and an uninterrupted deglaciation until about 9 ka (Fig. 13, lower panel, red curves).On the contrary, the simulation driven by the GRIP palaeotemperature record yields peak ice volume at 19.1 ka and a deglaciation interrupted by two phases of regrowth until about 8 ka.The first interruption occurs between 16.6 and 14.5 ka and the second between 13.1 and 11.6 ka (Fig. 13, lower panel, blue curve).
Hence, the two model runs, while similar in overall timing compared to the runs with other climate drivers, differ in detail.On the one hand, the EPICA run depicts peak glaciation about 2 ka later than the GRIP run, in closer agreement with dated maximum extents, and shows a faster, uninterrupted deglaciation which yields sporadic ice cover more than 1 ka earlier.On the other hand, the GRIP run yields a deglaciation in three steps, compatible with marine sediment sequences offshore Vancouver Island, where the distribution of ice-rafted debris indicates an ice margin retreat from the Georgia Strait in two phases that are contemporary with warming oceanic temperatures from 17.2 to 16.5 and from 15.5 to 14.0 14 C cal ka (Taylor et al., 2014).
Modelled patterns of ice sheet retreat are relatively consistent between the two simulations (Figs. 14 and 15).The southern sector of the modelling domain, including the Puget Lowland, the Coast Mountains, the Columbia and Rocky mountains, and the Interior Plateau of British Columbia, becomes completely deglaciated by 10 ka, whereas a significant ice cover remains over the Skeena, the Selwyn and Mackenzie, and the Wrangell and Saint Elias mountains in the northern sector of the modelling domain.After 10 ka, deglaciation continues to proceed across the Liard Lowland with a radial ice margin retreat towards the surrounding mountain ranges, consistent with the regional meltwater record of the last deglaciation (Margold et al., 2013a).Remaining ice continues to decay by retreating towards the Selwyn, Mackenzie, and Skeena mountains.The last remnants of the Cordilleran ice sheet finally disappear from the Skeena Mountains at 6.7 ka in both simulations.

Late-glacial readvance
The possibility of late-glacial readvances in the North American Cordillera has been debated for some time (Luckman and Osborn, 1979;Reasoner et al., 1994;Osborn and Gerloff, 1997;Osborn et al., 1995), and locally these have been reconstructed and dated.Radiocarbon-dated end moraines in the Fraser and Squamish valleys, off the southern tip of the Coast Mountains, indicate consecutive glacier maxima or standstills while in overall retreat, one of which corresponds to the Younger Dryas chronozone (Clague et al., 1997;Friele and Clague, 2002a, b;Kovanen, 2002;Kovanen and Easterbrook, 2002).Although most of these moraines characterise independent valley glaciers that may have been disconnected from the waning Cordilleran ice sheet, the Finlay River area in the Omineca Mountains (Fig. 15 importantly, their interaction with larger, lingering remnants of the main body of the Cordilleran ice sheet in the valleys (Lakeman et al., 2008).Additional evidence for late-glacial alpine glacier readvances includes moraines in the eastern Coast Mountains and the Columbia and Rocky mountains (Reasoner et al., 1994;Osborn and Gerloff, 1997;Menounos et al., 2008).
Although further work is needed to constrain the timing of the late-glacial readvance, to assess its extents and geographical distribution and to identify the potential climatic triggers (Menounos et al., 2008), it is interesting to note that the simulation driven by the GRIP record produces a late-glacial readvance in the Coast Mountains and in the Columbia and Rocky Mountains (Fig. 15, left panel).In addition to matching the location of some local readvances, the GRIP-driven simulation shows that a large remnant of the decaying ice sheet may still have existed at the time of this late-glacial readvance.In contrast, the EPICA-driven simulation produces a nearly continuous deglaciation with only a tightly restricted late-glacial readvance on the western slopes of the Saint Elias and the Coast mountains (Fig. 15, right panel).

Deglacial flow directions
Because a general conjecture in glacial geomorphology is that the majority of landforms (lineations and eskers) are part of the deglacial envelope (terminology from Kleman et al., 2006), having been formed close inside the retreating margin of ice sheets (Boulton and Clark, 1990;Kleman et al., 1997Kleman et al., , 2010)), we present maps of basal flow directions immediately preceding deglaciation or at the time of cessation of sliding inside a cold-based retreating margin (Fig. 16).The modelled deglacial flow patterns are mostly consistent between the GRIP and EPICA simulations.They depict an active ice sheet retreat in the peripheral areas, followed by stagnant ice decay in some of the interior regions.Several parts of the modelling domain do not experience any basal sliding throughout the deglaciation phase (Fig. 16, hatched areas).This notably includes parts of the Interior Plateau in British Columbia, major portions of the Alaskan sector of the ice sheet, and a tortuous ribbon running from the Northern Rocky Mountains over the Skeena and Selwyn Mountains and into the Mackenzie Mountains.Patterns of glacial lineations formed in the northern and southern sectors of the Cordilleran ice sheet (Prest et al., 1968;Clague, 1989, Fig. 1.12;Kleman et al., 2010, Fig. 2) show similarities with the patterns of deglacial ice flow from numerical modelling (Fig. 16).In the northern half of the modelling domain, modelled deglacial flow directions depict an active downhill flow as the last remnants of the ice sheet retreat towards mountain ranges.Converging deglacial flow patterns in the Liard Lowland, for instance (Fig. 16), closely resemble the pattern indicated by glacial lineations (Margold et al., 2013a, Fig. 2).Inside the retreating western margin, modelled south-westward flow patterns emanating from the Skeena Mountains and running over the central Coast Mountains are compatible with geological evidence from this region (Stumpf et al., 2000, Fig. 12).
On the Interior Plateau of south-central British Columbia, both simulations produce a retreat of the ice margin towards the north-east (Fig. 15), a pattern which is validated by the geomorphological and stratigraphical record for ancient proglacial lakes dammed by the retreating ice sheet (Perkins and Brennand, 2015).However, the two simulations differ in the mode of retreat.The GRIP simulation yields an active retreat with basal sliding towards the ice margin to the south, whereas the EPICA simulation produces negligible basal sliding on the plateau during deglaciation (Fig. 16).Yet, the Interior Plateau also hosts an impressive set of glacial lineations which indicate a substantial eastwards flow component of the Cordilleran ice sheet (Prest et al., 1968;Kleman et al., 2010;Ferbey et al., 2013), which is not present in any of the two simulations (Fig. 16).One explanation for the incongruent results could be that the missing feedback mechanisms between ice sheet topography and regional climate resulted in a modelled ice divide of the LGM ice sheet being too far to the east (Sect.5.1.2;Fig. 8; Seguinot et al., 2014).A more westerly located LGM ice divide would certainly result in a different deglacial flow pattern over the Interior Plateau.However, a more westerly positioned LGM ice divide would certainly be associated with an LGM ice sheet less extensive to the west and therefore thinner ice on the Interior Plateau during deglaciation than modelled here.Decreased ice thickness would not promote warm-based conditions but, on the contrary, enlarge the region of negligible basal sliding (Fig. 16).Thus, a second explanation for the incongruent results could be that the Interior Plateau lineation system predates deglaciation ice flow, as perhaps indicated by some eskers that appear incompatible with these glacial lineations (Margold et al., 2013b, Fig. 9).Finally, a third ex-  The modelled deglaciation of the Interior Plateau of British Columbia consists of a rapid northwards retreat (Fig. 15) of southwards-flowing non-sliding ice lobes (Fig. 16) positioned in-between deglaciated mountain ranges (Figs. 17 and 18).This result appears compatible with the prevailing conceptual model of deglaciation of central British Columbia, in which mountain ranges emerge from the ice before the plateau (Fulton, 1991, Fig. 7).However, due to different topographic and climatic conditions, our simulations produce different deglaciation patterns in the northern half of the model domain, indicating that this conceptual model may not be applied to the entire area formerly covered by the Cordilleran ice sheet.

Conclusions
Numerical simulations of the Cordilleran ice sheet through the last glacial cycle presented in this study consistently produce two glacial maxima during MIS 4 (62.2-56.9ka, 3.6-8.7 m s.l.e.) and MIS 2 (23.2-16.9ka, 8.6-8.8 m s.l.e.), two periods corresponding to documented extensive glaciations.This result is independent of the palaeotemperature record used among the six selected for this study and thus can be regarded as a robust model output, which broadly matches geological evidence.However, the timing of the two glaciation peaks depends sensitively on which climate record is used to drive the model.The timing of the LGM is best reproduced by the EPICA and Vostok Antarctic ice core records.It occurs about 2 ka too early in the simulation forced by the GRIP ice core record and occurs even earlier in all other simulations.The mismatch is largest for the two northwestern Pacific ODP palaeotemperature records, which are affected by the weakening of the California current during the LGM.Nevertheless, the use of palaeotemperature reconstructions from remote sites in Greenland and Antarctica produces the best agreement between modelled Cordilleran ice sheet dynamics and available geological evidence, and the significant differences remaining between the two preferred simulations highlight the need for more regional palaeoclimate reconstructions of the last glacial cycle in and around the North American Cordillera.In all simulations presented here, ice cover is limited to disjoint mountain ice caps during most of the glacial cycle.The most persistent nucleation centres are located in the Coast Mountains, the Columbia and Rocky mountains, the Selwyn and Mackenzie mountains, and most importantly, in the Skeena Mountains.Throughout the modelled last glacial cycle, the Skeena Mountains host an ice cap which appears to be fed by the moisture intruding inland from the west through a topographic breach in the Coast Mountains.The Skeena ice cap acts as the main nucleation centre for the glacial readvance towards the LGM configuration.As indicated by persistent, warm-based ice in the model, this ice cap perhaps explains the distinct glacial erosional imprint observed on the landscape of the Skeena Mountains.
During deglaciation, none of the climate records used can be selected as producing an optimal agreement between the model results and the geological evidence.Although the EPICA-driven simulation yields the most realistic timing of the LGM and, therefore, start of deglaciation, only the GRIP-driven simulation produces late-glacial readvances in areas where these have been documented.Nonetheless, the patterns of ice sheet retreat are consistent between the two simulations and show a rapid deglaciation of the southern sector of the ice sheet, including a rapid northwards retreat across the Interior Plateau of central British Columbia.The GRIP-driven simulation then produces a late-glacial readvance of local ice caps and of the main body of the decaying Cordilleran ice sheet primarily in the Coast and the Columbia and Rocky Mountains.In both simulations, this is followed by an opening of the Liard Lowland and a final retreat of the remaining ice caps towards the Selwyn and, finally, the Skeena mountains, which host the last remnant of the ice sheet during the middle Holocene (6.7 ka).Our results identify the Skeena Mountains as a key area to understanding glacial dynamics of the Cordilleran ice sheet, highlighting the need for further geological investigation of this region.
The Supplement related to this article is available online at doi:10.5194/tc-10-639-2016-supplement.
Author contributions.J. Seguinot ran the simulations; I. Rogozhina guided experiment design; A. P. Stroeven, M. Margold and J. Kleman took part in the interpretation and comparison of model results against geological evidence.All authors contributed to the text.

Figure 1 .
Figure 1.Relief map of the northern American Cordillera showing cumulative last glacial maximum ice cover between 21.4 and 16.8 14 C cal ka(Dyke, 2004, red line)  and the modelling domain used in this study (black rectangle).The background map consists of ETOPO1(Amante and Eakins, 2009) and Natural Earth Data(Patterson and Kelso, 2015).

Figure 3 .
Figure3.Effective pressure, N , as a function of water content in the till, W , for the default (δ = 0.02, W max = 2 m), hard bed (δ = 0.05, W max = 5 m), and soft bed (δ = 0.01, W max = 1 m) sliding parametrisations, using a linear scale (top panel) and a logarithmic scale (bottom panel).Calculations are made for an ice thickness, h, of 1000 m. Figure made using Eq. 5 with parameters from Table3after Bueler and van Pelt (2015, Fig.1).

Figure 4 .
Figure 4. Monthly mean near-surface air temperature, precipitation, and standard deviation of daily mean temperature (PDD SD)for January and July from the North American Regional Reanalysis (NARR;Mesinger et al., 2006), used to force the surface mass balance (PDD) component of the ice sheet model.Note the strong contrasts in seasonality, timing of the precipitation peak, and temperature variability over the model domain, notably between coastal and inland regions.

Figure 5 .
Figure 5. Temperature offset time series from ice core and ocean records (Table2) used as palaeoclimate forcing for the ice sheet model (top panel), and modelled sea-level relevant ice volume (bottom panel) through the last 120 ka, expressed in metres of sea-level equivalent (m s.l.e.).Grey fields indicate marine oxygen isotope stage (MIS) boundaries for MIS 2 and MIS 4 according to a global compilation of benthic δ 18 O records(Lisiecki and Raymo, 2005).Hatched rectangles highlight the time-volume span for ice volume extremes corresponding to MIS 4 (61.9-56.5 ka), MIS 3 (53.0-41.3ka), and MIS 2 (LGM, 23.2-16.8ka).Dotted lines correspond to GRIP-and EPICA-driven 5 km resolution runs.

Figure 6 .
Figure 6.Snapshots of modelled surface topography (500 m contours) corresponding to the sea-level relevant ice volume extremes indicated on Fig. 5.An ice cap persists over the Skeena Mountains (SM) during MIS 3. Note the occurrence of spatial similarities despite large differences in timing.

Figure 7 .
Figure7.Modelled sea-level relevant ice volume through the last 120 ka in the simulation forced by the GRIP palaeoclimate record, using default parameters (black curves), different ice rheology parameters (top panel), and different basal sliding parameters (bottom panel).Grey fields indicate marine oxygen isotope stage (MIS) boundaries for MIS 2 and MIS 4 according to a global compilation of benthic δ 18 O records(Lisiecki and Raymo, 2005).

Figure 8 .
Figure 8. Modelled surface topography (200 m contours) and surface velocity (colour mapping) corresponding to the maximum ice sea-level relevant volume during MIS 2 in the GRIP and EPICA high-resolution simulations.

Figure 9 .
Figure 9. Modelled surface topography (200 m contours) and surface velocity (colour mapping) corresponding to the maximum ice sea-level relevant volume during MIS 4 in the GRIP and EPICA high-resolution simulations.

Figure 10 .
Figure 10.Modelled duration of ice cover during the last 120 ka using GRIP and EPICA climate forcing.Note the irregular colour scale.A continuous ice cover spanning from the Alaska Range (AR) to the Coast Mountains (CM) and the Columbia and Rocky mountains (CRM) exists for about 32 ka in the GRIP simulation and 26 ka in the EPICA simulation.The maximum extent of the ice sheet generally corresponds to relatively short durations of ice cover, but ice cover persists over the Skeena Mountains (SM) during most of the simulation.See Fig. 1 for a list of abbreviations.

Figure 11 .
Figure 11.Modelled duration of warm-based ice cover during the last 120 ka.Long ice cover durations combined with basal temperatures at the pressure melting point may explain the strong glacial erosional imprint of the Skeena Mountains (SM) landscape.Hatches indicate areas that were covered by cold ice only.

Figure 12 .
Figure 12.Modelled fraction of warm-based ice cover during the ice-covered period.Note the dominance of warm-based conditions on the continental shelf and major glacial troughs of the coastal ranges.Hatches indicate areas that were covered by cold ice only.

Figure 13 .
Figure 13.Temperature offset time series from the GRIP and EPICA ice core records (Table 2) (top panel) and modelled sea-level relevant ice volume during the deglaciation (bottom panel).

Figure 14 .
Figure 14.Snapshots of modelled surface topography (200 m contours) and surface velocity (colour mapping) during the last deglaciation from the GRIP (top panels) and EPICA (bottom panels) 5 km simulations.Dashed segments (A-D) indicate the location of profiles used in Figs. 17 and 18.

Figure 15 .
Figure 15.Modelled age of the last deglaciation.Areas that have been covered only before the last glacial maximum are marked in green.Hatches denote readvance of mountain-centred ice caps and the decaying ice sheet between 14 and 10 ka.Dashed segments (A-D) indicate the location of profiles used in Figs. 17 and 18.

Figure 16 .
Figure 16.Modelled deglacial basal ice velocities.Hatches indicate areas that remain non-sliding throughout deglaciation (22.0-8.0 ka), notably including parts of the Interior Plateau (IP).Note the concentric patterns of deglacial flow in the Liard Lowland (LL).Sliding grid cells were distinguished from non-sliding grid cells using a basal velocity threshold of 1 m yr −1 .

Figure 17 .
Figure 17.Modelled bedrock (black) and ice surface (blue) topography profiles during deglaciation (22.0-8.0 ka) in the GRIP 5 km simulation, corresponding to the four transects indicated in Figs. 14 and 15.

Figure 18 .
Figure 18.Modelled bedrock (black) and ice surface (red) topography profiles during deglaciation (22.0-8.0 ka) in the EPICA 5 km simulation, corresponding to the four transects indicated in Figs. 14 and 15.

Table 1 .
Default parameter values used in the ice sheet model.

Table 3 .
Parameter values used in the sensitivity test.

Table 4 .
Extremes in Cordilleran ice sheet grounded ice extent and sea-level relevant ice volume corresponding to MIS 4, 3, and 2 for each of the six low-resolution simulations (Fig.5).
2.6 Climate forcingClimate forcing driving ice sheet simulations consists of a present-day monthly climatology, {T m0 , P m0 }, where temperatures are modified by offset time series, T TS , and lapserate corrections, T LR :

Table 4
), all model runs show relatively consistent patterns of glaciation.During MIS 4, all simulations produce an extensive ice sheet, covering an area of at least half of that The Cryosphere, 10, 639-664, 2016