The impact of the North American glacial topography on the evolution of the Eurasian ice sheet over the last glacial cycle

Modeling studies have shown that the continental scale ice sheets in North America and Eurasia in the last glacial cycle had a large influence on the atmospheric circulation and thus yielded a climate distinctly different from the present. However, to what extent the two ice sheets influenced each others growth trajectories remains largely unexplored. In this study we investigate how an ice sheet in North America influences the downstream evolution of the Eurasian ice sheet, using 5 a thermomechanical ice-sheet model forced by climate data from atmospheric snapshot experiments of three distinctly different phases of the last glacial cycle: the Marine Isotope Stages 5b, 4 and 2 (LGM). Owing to the large uncertainty associated with glacial changes of the Atlantic meridional overturning circulation, each atmospheric snapshot experiment was conducted using two distinctly different ocean heat transport representations. Our results suggest that changes in the North Amer10 ican paleo-topography may have largely controlled the zonal distribution of the Eurasian ice sheet. In the MIS4 and LGM experiments, the Eurasian ice sheet migrates westward towards the Atlantic sector – largely consistent with geological data and contemporary ice-sheet reconstructions – due to a low wavenumber stationary wave response, which yields a cooling in Europe and a warming in northeastern Siberia. The expansion of the North American ice sheet between MIS4 and LGM am15 plifies the Siberian warm anomaly, which limits the glaciation there and may therefore help explain the progressive westward migration of the Eurasian ice sheet in this time period. The ocean heat transport has only a small influence on the stationary wave response to the North American glacial topography; however, because temperature anomalies have a smaller influence on an ice sheet’s ablation in a colder climate than in a warmer one, the impact of the North American glacial topography 20 on the Eurasian ice-sheet evolution is reduced for colder surface conditions in the North Atlantic. While the Eurasian ice sheet in the MIS4 and LGM experiments appears to be in equilibrium with the simulated climate conditions, the MIS5b climate forcing is too warm to grow an ice sheet in Eurasia. First-order sensitivity experiments suggest that the MIS5b ice sheet was established during preceding colder stages. 25

Abstract.Modeling studies have shown that the continentalscale ice sheets in North America and Eurasia in the last glacial cycle had a large influence on the atmospheric circulation and thus yielded a climate distinctly different from the present.However, to what extent the two ice sheets influenced each others' growth trajectories remains largely unexplored.In this study we investigate how an ice sheet in North America influences the downstream evolution of the Eurasian ice sheet, using a thermomechanical ice-sheet model forced by climate data from atmospheric snapshot experiments of three distinctly different phases of the last glacial cycle: the Marine Isotope Stages 5b, 4, and 2 (Last Glacial Maximum -LGM).Owing to the large uncertainty associated with glacial changes in the Atlantic meridional overturning circulation, each atmospheric snapshot experiment was conducted using two distinctly different ocean heat transport representations.Our results suggest that changes in the North American paleo-topography may have largely controlled the zonal distribution of the Eurasian ice sheet.In the MIS4 and LGM experiments, the Eurasian ice sheet migrates westward towards the Atlantic sector -largely consistent with geological data and contemporary ice-sheet reconstructions -due to a low wave number stationary wave response, which yields a cooling in Europe and a warming in northeastern Siberia.The expansion of the North American ice sheet between MIS4 and the LGM amplifies the Siberian warm anomaly, which limits the glaciation there and may therefore help explain the progressive westward migration of the Eurasian ice sheet in this time period.The ocean heat transport only has a small influence on the stationary wave response to the North American glacial topography; however, because temperature anomalies have a smaller influence on an ice sheet's ablation in a colder climate than in a warmer one, the impact of the North American glacial topography on the Eurasian ice-sheet evolution is reduced for colder surface conditions in the North Atlantic.While the Eurasian ice sheet in the MIS4 and the LGM experiments appears to be in equilibrium with the simulated climate conditions, the MIS5b climate forcing is too warm to grow an ice sheet in Eurasia.First-order sensitivity experiments suggest that the MIS5b ice sheet was established during preceding colder stages.

Introduction
The Quaternary period is characterized by the alternation between cold and warm phases -glacials and interglacialswhen massive ice sheets expand and retreat over the subpolar continents.The last glacial cycle began about 115 000 years ago (115 kyr BP) following a minimum in the boreal summer insolation (Berger and Loutre, 1991).Over the subsequent ∼ 90 kyr, paleo-records suggest that ice sheets progressively expanded in North America and Eurasia, with relatively rapid ice growth during colder phases followed by warmer periods when the global ice volume remained relatively constant (Peltier and Fairbanks, 2006;Stokes et al., 2012;Kleman et al., 2013).The colder phases are typically referred to as the Marine Isotope Stages (MISs) 5d (106-115 kyr BP), 5b (85-93 kyr BP), 4 (60-74 kyr BP), and 2 (12-24 kyr BP), where Published by Copernicus Publications on behalf of the European Geosciences Union.LGM, based on the ice-sheet reconstructions in Kleman et al. (2013).The shading represents ice sheets and the contour interval is 500 m. the latter includes the culmination of the last glacial cycle at the Last Glacial Maximum (LGM; 19-23 kyr BP).
The progressive increase in the Northern Hemisphere ice volume was dominated by the Laurentide and Cordilleran ice sheets in North America (Kleman et al., 2013).Subsequent to the ice-sheet inception in the Canadian Arctic and Quebec, the Laurentide ice sheet expanded over the eastern parts of the continent and eventually coalesced with the Cordilleran ice sheet to form a coherent and continent-wide ice sheet at the LGM (Fig. 1; Clark et al., 1993;Kleman et al., 2010Kleman et al., , 2013)).As opposed to the North American counterpart, the combined volume of the Eurasian ice sheets (Fennoscandian and Barents-Kara ice sheets) changed relatively little between the inception phase and the LGM (Fig. 1; Svendsen et al., 2004;Kleman et al., 2013).Instead, the most notable feature of the ice-sheet evolution in Eurasia is a progressive westward migration in time; in the early and intermediate stages (MIS5b and MIS4) the eastern margin of the Eurasian ice sheet was located in central Siberia (Svendsen et al., 2004;Kleman et al., 2013), whereas essentially only northern Europe and Great Britain (Bradwell et al., 2008) were ice-covered at the LGM (Fig. 1).Hence, in both North America and Eurasia, the ice sheets had strong zonal asymmetries toward the Atlantic sector over large parts of the glacial cycle.The driving mechanism of this asymmetry remains an open question as it has been difficult to capture this feature in conventional ice-sheet model experiments (Marshall et al., 2000;Zweck and Huybrechts, 2005;Charbit et al., 2007;Bonelli et al., 2009;Beghin et al., 2014).
The role of ice-sheet-atmosphere interactions has mostly been studied for the buildup of the North American ice sheet (Roe and Lindzen, 2001;Liakka et al., 2011;Löfverström et al., 2014Löfverström et al., , 2015)).These studies suggest that the east-heavy pre-LGM configuration arose from changes in the time-mean atmospheric circulation (stationary waves) forced by the ice sheet itself, possibly in combination with complex interactions with the North American Cordillera.The temporal evolution of the Eurasian ice sheet has received less attention.The orographic precipitation feedback, initially proposed by Sanberg and Oerlemans (1983), is generally considered an important feature to explain the westward migration of the ice sheet (Roe and Lindzen, 2001;Van den Berg et al., 2008;Liakka and Nilsson, 2010;Kleman et al., 2013;Löfverström et al., 2014) as surface winds from the Atlantic are forced vertically by the western and southern slopes of the ice sheet, hence leading to increased precipitation rates in those areas and ultimately to a (south)westward propagation of the ice sheet (Sanberg and Oerlemans, 1983).Although orographic precipitation is a robust feature in atmospheric general circulation models (Roe, 2005), questions regarding the timing of the westward migration of the ice sheet remain unanswered.For example, why did the Eurasian ice sheet propagate westward only in the latter stages of the glacial cycle and not immediately subsequent to the inception phase?The answer to this question is complicated by the fact that the orientation of the Atlantic storm track, which has a large impact on the European precipitation, appears to be controlled by the size of the North American ice sheet; for smaller ice sheets in North America (e.g., MIS5b and MIS4) the Atlantic storm track has a pronounced southwest-northeast tilt (similar to the modern climate; Löfverström et al., 2014;Pausata and Löfverström, 2015), whereas for large ice sheets (LGM) the storm track has a more zonal orientation (Li and Battisti, 2008;Kageyama et al., 2013;Löfverström et al., 2014;Ullman et al., 2014;Merz et al., 2015;Pausata and Löfverström, 2015).The zonalization of the Atlantic storm track typically yields drier (wetter) conditions in northern (southern) Europe (Löfverström et al., 2014).
The connection between the size of the Laurentide ice sheet and the orientation of the Atlantic storm track suggests that the North American glacial topography may have influenced the ice-sheet evolution in Eurasia.Studies investigating remote climate impacts of the North American and Eurasian paleo-topography typically used static ice sheets as forcing in comprehensive circulation models (e.g., Li and Battisti, 2008;Löfverström et al., 2014;Ullman et al., 2014) or dynamic ice-sheet models coupled to highly simplified atmospheric models (sometimes with parameterized climate anomalies; Beghin et al., 2014).In this study, we investigate the effect of the geologically constrained ice sheets in North America at MIS5b, MIS4, and the LGM (see Fig. 1) on the evolution of the Eurasian ice sheet.The atmospheric response to the North American ice sheets is evaluated using a comprehensive atmospheric circulation model with nonlinear dynamics.The atmospheric fields are used as forcing in a thermomechanical ice-sheet model in order to evaluate their impact on the Eurasian ice sheet.More information about the models and the experiments is given in Sect. 2. In Sect. 3 we show the main results from the atmospheric and ice-sheet model experiments, followed by a comprehensive discussion in Sect. 4. Finally, the conclusions are summarized in Sect. 5.

Models and experiments
The successive experimental approach used here is summarized in Fig. 2. The details are provided below.

Atmospheric simulations
We use the climate snapshot (steady-state) simulations from Löfverström et al. (2014), representative of the preindustrial (PI), MIS5b (88 kyr BP), MIS4 (66 kyr BP), and the LGM (20 kyr BP) climates.These experiments were conducted with the National Center for Atmospheric Research Community Atmospheric Model version 3 (NCAR CAM3; Collins et al., 2006) using T85 spectral resolution (approximately 1.4 • horizontal resolution) and 26 hybrid levels in the vertical.Land surface processes are handled by the Community Land Model 3 (CLM3; Oleson et al., 2004).The ocean is represented by a mixed-layer (slab) model with a prescribed depth and ocean heat transport.
For each glacial time slice, we conduct two sets of simulations with different surface topography: (i) the reconstructed glacial topography from Kleman et al. (2013) (hereafter referred to as the "fullGlacial" simulations; see Fig. 1) and (ii) the same topography as (i) except for using present-day topography in North America (hereafter referred to as the "EAonly" simulations).The impact of the North American ice sheet on the climate is evaluated as the difference between the fullGlacial and EAonly simulations.The orbital clock (Berger and Loutre, 1991) and greenhouse gas concentrations (Petit et al., 1999;Spahni et al., 2005) were adjusted to the nominal time of the ice-sheet reconstruction (Table 1).Other boundary conditions, e.g., aerosols, vegetation, and land fraction, were set to preindustrial values; the latter two were properly adjusted for glaciated regions (Löfverström et al., 2014).As reference climate, we use an equilibrated present-day simulation from the same model (Hurrell et al., 2006).Results presented below are based on climatologies over 25 years after the simulated climates have reached statistical equilibrium.

Slab ocean model and ocean heat transport representations
We use a simplified slab (mixed-layer) ocean model in order to facilitate a high number of experiments, bracketing the uncertainty range in the planetary boundary conditions (see also Löfverström et al., 2014).The slab ocean model has a prognostic sea-surface temperature (SST) and a dynamic sea-ice edge calculated from the surface energy balance and the prescribed ocean heat transport (OHT) in the mixed layer (Collins et al., 2004;Bitz et al., 2012).Thus, the slab ocean model does not account for changes in ocean dynamics but retains the thermodynamic feedback between the ocean and the atmosphere.
The westerly mean flow implies that the North Atlantic sea-ice cover has a large influence on the temperature and moisture availability in Eurasia (e.g., Smith et al., 2003).The mean position of the North Atlantic sea-ice margin is in turn largely maintained by the Atlantic meridional overturning circulation (AMOC; Bitz et al., 2005).The strength of the LGM AMOC is a topic of ongoing research as it cannot be explicitly inferred from proxy-data evidence; modeling studies with coupled atmosphere-ocean models disagree on the LGM AMOC strength, with some models suggesting that it was stronger than at present, whereas other models yield a weakening (Otto-Bliesner et al., 2007;Weber et al., 2007).Following Löfverström et al. (2014), we therefore use two end-member representations of the OHT to bracket the uncertainty range of the AMOC strength.Both OHT representations are derived from equilibrated simulations with the (NCAR) Community Climate System Model version 3 (CCSM3), which is a fully coupled model using CAM3 as atmospheric component.The first OHT representation, which represents a state of a relatively strong AMOC, stems from a preindustrial simulation (hereafter referred to as "PI OHT") and the second representation -representative of a weak AMOC state -from the LGM simulation in Brandefelt and Otto-Bliesner (2009) ("BO2009 OHT").Note Table 1.Top of the atmosphere insolation during the northern summer solstice (60 • N; Berger and Loutre, 1991) and greenhouse gas concentrations (Petit et al., 1999;Spahni et al., 2005) (Brandefelt and Otto-Bliesner, 2009).Unfortunately, the simulated LGM mixed-layer depth is not available on the CCSM3 data server; following Löfverström et al. (2014) we therefore use a modern annual mean mixed-layer depth in all simulations.Aside from areas covered by perennial sea ice, simulated changes in the LGM mixed-layer depth are, however, small compared to PI (of order 1-10 m; Brandefelt and Otto-Bliesner, 2009).Changes in the simulated LGM AMOC have a large impact on the sea-surface conditions in the North Atlantic; in CCSM3 (CAM3 with BO2009 OHT), the simulated AMOC is reduced with respect to the preindustrial and the annual mean SSTs are substantially lower than contemporary proxybased LGM SST reconstructions (MARGO Project Members et al., 2009), resulting in a zonal sea-ice margin at ∼ 40 • N (Fig. 3b; Brandefelt and Otto-Bliesner, 2009). 1 In CCSM4 (Brady et al., 2013), on the other hand, the simulated LGM SSTs are significantly warmer than in CCSM3 and the seaice margin is located farther to the north in the eastern North Atlantic (Fig. 3c; Brady et al., 2013) and is thus in better agreement with LGM sea-ice reconstructions (Paul and Schäfer-Neth, 2003;Toracinta et al., 2004;De Vernal et al., 2005, 2006;MARGO Project Members et al., 2009).The annual mean LGM SST response using the relatively strong PI OHT is in better agreement with the response in CCSM4 (and thus with the proxy) than in CCSM3 (Fig. 3a); for example, the sea-ice margin is located further north in the eastern North Atlantic (Fig. 3a).Therefore, we base the analysis on the simulations with PI OHT, whereas the simulations with BO2009 OHT are used in sensitivity experiments (Sect.4.2).

Model description
To simulate the evolution of the Eurasian ice sheet, we use the three-dimensional ice-sheet model SICOPOLIS (SImulation COde for POLythermal Ice Sheets, version 3.1), which treats ice as an incompressible, viscous, and heat-conducting fluid (Greve, 1997).The model equations are subjected to the shallow-ice approximation, which means that only the lowest-order terms are retained (Hutter, 1983).The model obeys Glen's flow law to calculate strain rates (deformation) from the applied stresses (e.g., Van der Veen, 2013), and a Weertman-type sliding scheme to calculate the basal velocities (Weertman, 1964).Ice streams are not explicitly accounted for.We run the model in the "cold-ice mode", i.e., temperatures above the pressure melting point are artificially reset to the pressure melting temperature.Expansion of marine ice is allowed if the bathymetry is less than 500 m (default value), otherwise instant calving is assumed.The bedrock and overriding ice sheet are assumed to relax to iso-static equilibrium with a timescale of 3 kyr, and the geothermal heat flux is 55 mW m −2 over the entire domain.
The surface mass balance is given by the difference between accumulation and ablation.In SICOPOLIS, accumulation is equal to precipitation and the ablation is parameterized using the positive-degree-day (PDD) approach (Braithwaite and Olesen, 1989;Reeh, 1991).The amount of PDDs in a year is given by the integrated sum of positive temperatures over that year and is evaluated using the semi-analytical solution in Calov and Greve (2005).It is assumed that the daily temperatures in a month are normally distributed about the monthly mean temperature.The standard deviation (dayto-day variability) of the temperature is 5 • C everywhere (default).We use the default values of the degree-day constants, which relate the PDDs to actual melt rates (3 mm day −1 K −1 for snow and 12 mm day −1 K −1 for ice).The melting procedure follows Reeh (1991).First, the PDDs are used to melt the annual snow fall.It is assumed that 60 % of that meltwater percolates into the ice and contributes to the formation of superimposed ice.Second, the superimposed ice is melted, after which the remaining PDDs, if any, are used to melt the glacier ice.
Following Charbit et al. (2002Charbit et al. ( , 2007)), the surface temperature (T ) and precipitation (P ) over the evolving ice sheet are modified according to a fixed atmospheric lapse rate γ : (1) where z(t) is the height of the evolving ice-sheet surface (t = time) and T 0 and P 0 are the reference temperature and precipitation on the initial ice-free topography z 0 , respectively (see Eqs. ( 3) and (4) in Sect.2.2.2).Hence, it is assumed that the temperature decreases linearly with height z at the lapse rate γ (set to the value of the standard atmosphere: γ = −6.5 × 10 −3 K m −1 ) and that precipitation decreases exponentially with the temperature change (due to elevation) times the parameter γ s , which relates the temperature anomaly to precipitation change (set to γ s = 0.05 K −1 following Charbit et al., 2002Charbit et al., , 2007)).Because the surface temperature on the ice sheet evolves over time, the relative amount of solid and liquid precipitation is parameterized; following Marsiat (1994), the fraction snowfall to the total precipitation is 1 if the monthly mean air temperature is below −10 • C and 0 if it is greater than 7 • C. For intermediate temperatures the fraction snowfall is linearly interpolated.

Experimental design and initial climate forcing
The SICOPOLIS simulations are carried out to steady state (at least 150 kyr) from an ice-free initial state using the CAM3 simulations as climate forcing.The horizontal resolution is set to 80 km, and the model domain covers most of the Northern Hemisphere.The relatively coarse horizontal resolution is motivated by the fact that we are interested in large-scale first-order changes in the Eurasian ice sheet (as a reference, Kleman et al., 2013 used a horizontal resolution of 95 km in their ice-sheet reconstructions).The vertical resolution amounts to 81 levels in the ice and 11 levels in the bedrock.We use the procedure described in Charbit et al. (2007) to deduce the initial fields of surface temperature (T 0 ) and precipitation (P 0 ) from the atmospheric model: To account for systematic biases in the atmospheric climatology, we first calculate anomalies of the glacial temperature (T paleo, CAM ) with respect to the temperature of the presentday simulation (T PD, CAM ).In doing so, we correct for the different orographies in the glacial and present-day simulations (z paleo, CAM and z PD, CAM , respectively) using the standard lapse rate.Subsequently, the anomalies are bilinearly interpolated to the SICOPOLIS grid and added to the observational data set (T PD, obs ), which is based on ERA-Interim reanalysis data (Dee et al., 2011).To calculate P 0 we use the same technique as for T 0 , but we use ratios instead of anomalies in order to omit negative precipitation (Charbit et al., 2007).

Summer temperature
The annual ablation is dominated by the summer conditions; we therefore focus on the surface temperature in boreal summer (June-August: JJA). Figure 2 shows the JJA surface temperature in the reanalysis data (a), present-day simulation (b), and the EAonly paleo-simulations (c, e, g).To highlight areas susceptible to inception, the temperatures in Fig. 4 are projected onto the present-day orography using the standard LGM EAonly 15.9 lapse rate.A summary of the average summer temperatures in the Northern Hemisphere and Eurasia is presented in Table 2.
The average Northern Hemisphere summer temperature decreases across the EAonly simulations; it drops by 3 • C between present day and MIS5b and by an additional 2 • C at the LGM (Table 2).The cooling across the glacial simulations has even larger regional variations: in Eurasia, the LGM summer temperature is about 5 • C lower than at MIS5b (Table 2).Regions with subfreezing summer temperatures are particularly interesting for glacial inception; the average position of the zero-degree summer (surface) isotherm is indicated by the green contour in Fig. 4a, b, c, e, and g.Similar to present day (Fig. 4a, b), the zero-degree isotherm at MIS5b is mainly located in the Arctic Ocean poleward of the Eurasian continent (Fig. 4c).Owing to the cooler conditions at MIS4 and the LGM, the (zonal) average location of the zero-degree isotherm is shifted equatorward by 6 to 7 • (Table 2); the largest regional changes are found in Scandinavia and eastern Siberia, where it reaches as far south as 60 • N at MIS4 and the LGM (Fig. 4e, g).
Figure 4d, f, and h shows the summer (surface) temperature anomalies induced by the North American ice sheet.These anomalies are calculated as the difference between the fullGlacial and EAonly simulations; a lapse rate correction has been applied to account for elevation differences.Due to an increased surface albedo and cold air advection by orographically forced stationary waves (Cook and Held, 1988;Roe and Lindzen, 2001;Abe-Ouchi et al., 2007;Liakka and Nilsson, 2010;Liakka, 2012;Löfverström et al., 2015), the largest cooling occurs in the vicinity of the North American ice sheet.In Eurasia, the temperature response to the North American ice sheet exhibits large regional variations.For all time slices, the North American ice sheet induces colder conditions in Europe.The response in Siberia is more complicated; at MIS5b the Siberian temperature response is almost negligible (Fig. 4d), whereas there is a warming in eastern Siberia at MIS4 and the LGM.The largest difference between the MIS4 and LGM responses is found in central Siberia, which becomes colder at MIS4 (Fig. 4f) and warmer at the LGM (Fig. 4h).

Annual precipitation
The large-scale features of the annual precipitation in the EAonly simulations are reminiscent of the modern climate, although the global precipitation rates are somewhat reduced in the glacial simulations (Fig. 5a, b, c, e, g).The largest precipitation rates in Eurasia are found in northwestern Europe where the cyclones from the Atlantic storm track make landfall.
As for the temperature, the largest precipitation response to the North American ice sheet is found over the ice sheet itself, with generally increased precipitation on the windward (westerly) slopes of the ice sheet and reduced precipitation over the leeward (easterly) slopes (Fig. 5d, f, h).In Eurasia, the North American ice sheet has a relatively small impact on the precipitation at MIS5b and MIS4 (Fig. 5d, f) but yields a significantly reduced precipitation in northwestern Europe at the LGM (Fig. 5h).As discussed in Löfverström et al. (2014), the reduced precipitation rates at the LGM are associated with a zonalization of the midlatitude Atlantic jet stream resulting from flow-topography interactions with the continent-wide North American ice sheet (Li and Battisti, 2008;Ullman et al., 2014;Merz et al., 2015;Pausata and Löfverström, 2015).Löfverström et al. (2014) found that this effect is not present for the smaller pre-LGM ice sheets (MIS5b and MIS4), as their location and spatial extent allow the mean flow to largely circumvent the topography, thus rendering the tilt of the Atlantic jet -and storm track -largely similar to the present day.

Summer stationary waves
In order to understand the temperature response in Fig. 4, we examine the stationary Rossby waves in the different climate states.Stationary waves, defined as zonal asymmetries in the climatological fields, are the result of large-scale orography and diabatic heating (e.g., Hoskins and Karoly, 1981;Held et al., 2002;Held, 1983;Kaspi and Schneider, 2011).Ice sheets constitute both orographic and diabatic forcing of stationary waves.Therefore, ice sheets expanding into the westerly mean flow can potentially influence the global stationary wave field (e.g., Cook and Held, 1988;Roe and Lindzen, 2001;Löfverström et al., 2014).
The lower-troposphere (700 hPa) geopotential height anomalies from the EAonly simulations are shown in Fig. 6c, e, and g.The stationary wave response is qualitatively similar in all glacial time slices; similar to the modern climate (Fig. 6a, b), the summer stationary wave field is characterized by anticyclonic circulation (ridges) over the subtropical Pacific and Atlantic ocean basins, and cyclonic circulation (troughs) over Asia and northeastern Canada (Fig. 6c, e, g).In addition, the ridge over the Atlantic Ocean extends over Europe and covers most of the ice sheet area, suggesting that the local ridge is excited by the Eurasian ice sheet.As noted by Löfverström et al. (2014), this indicates that the ice sheet's diabatic cooling dominates the stationary wave response.
The 700 hPa geopotential height responses to the North American ice sheets are shown as shading in Fig. 6d, f,  and h.As expected from theory, the stationary wave amplitudes increase with the size (spatial extent and height) of the North American ice sheet (Cook and Held, 1992;Ringler and Cook, 1997;Liakka and Nilsson, 2010;Liakka et al., 2011;Löfverström et al., 2014).Besides the amplitude, the stationary wave response to the North American ice sheet is qualitatively similar in all time slices.The local response is a ridge over the northwestern parts of the North American ice sheet and a trough in the southeast.This particular response is a robust feature across models using nonlinear stationary wave dynamics (Ringler andCook, 1997, 1999;Liakka et al., 2011).The remote downstream response consists of two wave trains: (i) a subtropical wave train with a northwest-southeast orientation and (ii) a low wave number polar wave train with a more zonal orientation.The polar wave train is characterized by a trough over Europe and western Asia and a ridge over Siberia.
The contours in Fig. 6d, f, and h depict the geopotential height anomalies at 300 hPa.The anomalies at this level have essentially the same spatial location as at 700 hPa, indicating that the climatological response to the North American ice sheet is largely equivalent barotropic.
In the summer season, high-latitude geopotential height anomalies are typically well correlated with surface temperature anomalies.Ridges are associated with reduced cloudiness and increased downwelling shortwave radiation, which leads to a surface warming, whereas troughs typically yield increased cloudiness and thus lower surface temperatures.This is also seen here (cf.Fig. 4d, f, h and 6d, f, h): the ridge over eastern Siberia and Alaska is associated with a surface warming and the trough in Europe with a cooling.Note that the magnitude of these temperature anomalies, in particular the Siberian warm anomaly, is not only controlled by the geopotential height anomalies but also by albedo feed- backs due to changes in the snow cover (see Fig. S1 in the Supplement).

Ice-sheet evolution
In this section we examine how the altered climate conditions -induced by the North American ice sheet -influence the spatial equilibrium extent of the Eurasian ice sheet.To evaluate our results, we compare the simulated extents of the Eurasian ice sheet with the geologically constrained reconstructions from Kleman et al. (2013).Note that we only compare the geographical distribution of ice (i.e., ice area) but not the ice thickness or ice volume.The reason is that the ice thickness in the Kleman et al. (2013) reconstructions is a model-dependent feature, whereas the spatial extents are constrained by geological evidence.
Figure 7 shows the simulated equilibrium ice thickness when using the atmospheric simulations summarized in Figs. 4 and 5 as climate forcing.Apart from some ice caps in the Scandinavian mountains, Eurasia remains ice free at MIS5b (Fig. 7a, b).This is consistent with a negative surface mass balance over essentially the entire domain (Fig. S2 in the Supplement).A comprehensive discussion on the po-tential shortcomings in the MIS5b simulations follows in Sect.4.3.
At MIS4 and the LGM, atmospheric circulation changes induced by the North American ice sheet serve to increase the total ice area in Eurasia by about 80 and 30 %, respectively (Fig. 7).This increase is mediated by an expansion of ice in western Eurasia and a reduced ice extent in eastern Siberia; apart from slightly too much ice in the Kara Sea region in the LGM simulation, the outlines of the simulated MIS4 and LGM ice sheets in Eurasia are in good agreement with the reconstructions from Svendsen et al. (2004) and Kleman et al. (2013); see Fig. 7d, f.In the absence of ice in North America, the MIS4 and LGM ice sheets are more zonally distributed along the Arctic coast (Fig. 7c, e).Hence, our simulations suggest that the North American ice sheet induces a westward migration of the Eurasian ice sheet, and consequently, the evolution of the Eurasian ice sheet between MIS4 and the LGM was to a large extent controlled by the growth of the North American ice sheet.
The ice sheet's westward migration is depicted in Fig. 8, which shows the longitude of the center of mass (λ c ) of the total ice distribution in Eurasia.In the reconstructions from Kleman et al. (2013), λ c decreases from 49 • E at MIS5b to 44  and 27 • E at MIS4 and the LGM, respectively (black bars in Fig. 8).Note that the westward migration between MIS4 and the LGM is captured only if the North American ice sheet is present (white bars in Fig. 8), otherwise λ c remains large (∼ 55-60 • E) for both stages (gray bars in Fig. 8).

Discussion
We have examined how the North American ice sheet (constrained by geological data) influences the extent of the Eurasian ice sheet in the MIS5b, MIS4, and LGM climate states.We found that the MIS4 and LGM ice sheets in North America yield a westward migration of the Eurasian ice sheet (Fig. 8), characterized by more ice in Europe and less ice in Siberia (Fig. 7).When accounting for the North American ice sheet, the spatial distributions of the simulated MIS4 and LGM ice sheets in Eurasia are in good agreement with contemporary proxy-based ice-sheet reconstructions (Fig. 7; Svendsen et al., 2004;Kleman et al., 2013); this suggests that the growth of the North American ice sheet between MIS4 and the LGM may have been vital for limiting and shifting the Eurasian ice sheet westward in time.

North American influence on the Eurasian climate
The westward migration of the Eurasian ice sheet in the MIS4 and LGM simulations (Figs. 7 and 8) is associated with changes in the summer stationary wave field.The North American ice sheet yields a cooling (less ablation) in Europe and a warming (more ablation) in northeastern Siberia (Fig. 4f, h).These temperature anomalies are associated with an equivalent barotropic cyclonic-anticyclonic anomaly in the target regions (Fig. 6f, h).
The geopotential height anomalies in Fig. 6d, f, and h result from (typically nonlinear) interactions between the atmospheric flow and the thermal and orographic forcing of the North American ice sheet.This leads to a complicated nonlinear response in the vicinity of the wave source (i.e., the North American ice sheet in our case) where the climate anomalies rotate clockwise for larger topographic barriers (Cook and Held, 1992;Ringler and Cook, 1997;Liakka et al., 2011;Liakka, 2012).Away from the wave source, however, the geopotential height anomalies share many similarities with linear wave theory (see Appendix A for details).For example, the low wave number polar wave train in Fig. 6d, f, and h is consistent with a latitudinal decrease in the (barotropic) stationary wave number due to the spherical geometry of the planet (see Appendix A and Fig. S3 in the Supplement).In addition, linear Rossby ray tracing arguments (Appendix A and Hoskins and Karoly, 1981) suggest that higher-latitude wave trains should have a more zonal orientation than wave trains at lower latitudes and are thus consistent with Fig. 6d, f, and h.
Although the remote stationary wave response is broadly consistent with linear theory, our findings are different from the coupled ice-sheet-climate model experiments in Beghin et al. (2014), who used CLIMBER-2 (CLIMate-BiosphERe model) with (linear) parameterized stationary waves to examine the interaction between the Northern Hemisphere ice sheets.They found that the North American ice sheet has a negligible impact on European summer temperatures but yields a slight cooling in Siberia, thus contradicting our results (Fig. 4).Although CLIMBER has proven to be a valuable model for studying the transient ice-sheet-climate evolution through the glacial cycles (e.g., Calov et al., 2002Calov et al., , 2005;;Bonelli et al., 2009;Ganopolski et al., 2010;Ganopolski and Calov, 2011), it has a very limited representation of the atmospheric circulation; in particular it does not account for Rossby wave dynamics.Unless explicitly corrected for (Ganopolski et al., 2010), the lack of Rossby wave dynamics in CLIMBER typically facilitates ice inception over the western rather than the eastern part of North America (Bonelli et al., 2009;Beghin et al., 2014); this presumably influenced the Eurasian climate anomalies in Beghin et al. (2014).

Sensitivity to OHT
Owing to the large uncertainty of the AMOC strength during glacial times (Otto-Bliesner et al., 2007;Weber et al., 2007), we perform sensitivity simulations of the equilibrium ice thickness using the atmospheric simulations with BO2009 OHT (Brandefelt and Otto-Bliesner, 2009) as climate forcing.The results are summarized in Fig. 9. Due to the colder conditions in the simulations with BO2009 OHT, the Eurasian ice sheet expands equatorward compared to when using PI OHT (cf.Figs. 9 and 7).However, despite the colder conditions in the North Atlantic, the model fails to simulate a large ice sheet at MIS5b (Fig. 9a, b).
Using the BO2009 OHT climate forcing, the North American ice sheet induces a westward migration of the Eurasian ice sheet (λ c is reduced by 6 • for MIS4 and by 11 • for the LGM; not shown); however, it is not as pronounced as with PI OHT (Fig. 8).Since the stationary wave response to the North American ice sheet is qualitatively similar in the BO2009 OHT (Figs.S4-S6 in the Supplement) and the PI OHT simulations (Figs. 4,5 and 6), the reduced westward migration in the BO2009 OHT simulations is most likely attributed to the colder climate; in the PDD model, cold background conditions (temperatures below freezing) reduces the effect of temperature anomalies on the ablation.

What prevents ice-sheet growth at MIS5b?
The vexing issue of this study is that we fail to simulate an MIS5b ice sheet of comparable size to the data-based reconstructions (Figs. 7 and 9).The lack of ice growth at MIS5b is associated with a negative surface mass balance across the entire Eurasian continent (Fig. S2) due to relatively high summer temperatures (Fig. 4 and Table 2).The relatively warm conditions at MIS5b compared to MIS4 and the LGM are attributed to both a higher insolation and higher concentrations of greenhouse gases (Table 1).It is possible that allowing for certain feedbacks, such as vegetation changes (Colleoni et al., 2009;Liakka et al., 2014), would cool the summer climate and thus support ice inception at MIS5b.However, since the MIS4 and LGM extents of the Eurasian ice sheet are in good agreement with the reconstructions when omitting these feedbacks, it seems unlikely that systematic biases in the climate forcing are the primary cause of the lack of ice growth at MIS5b.
In addition, we fail to find multiple equilibrium states (e.g., Calov and Ganopolski, 2005;Abe-Ouchi et al., 2013) of the simulated Eurasian ice sheets (Fig. S7 in the Supplement); initializing the ice-sheet simulations using the Kleman et al. (2013) reconstructions leads to very similar equilibrium extents as in Fig. 7.This suggests that preceding configurations of the Eurasian ice sheet were not crucial for maintaining the ice sheet at MIS5b.Instead, it is more likely that the MIS5b ice sheet was not in equilibrium with the prevailing climate.The successful glacial inception and good agreement between the equilibrated MIS4 and LGM ice sheets and the reconstructions suggest that the climate was locally cold enough to support glacial inception and the resulting ice sheets were in equilibrium with the prevailing climate; this is, however, not necessarily true for MIS5b.Instead, it is plausible that the MIS5b climate was too warm to support glacial inception and the ice sheet was a remnant of ice growth in preceding colder periods.In this context it is interesting to note that the Eurasian ice sheet reached a size comparable to MIS5b already at ∼ 105 kyr BP, subsequent to a relative minimum in the high-latitude boreal summer insolation (Kleman et al., 2013;Löfverström et al., 2014).
Since we do not have access to any atmospheric simulations of a colder stage prior to MIS5b, we use a crude approach by imposing a cooling of the JJA temperature artificially in SICOPOLIS.To estimate the magnitude of the cooling, we employ the parameterization of the surface temperature to changing insolation proposed by Abe-Ouchi et al. (2007, 2013); based on sensitivity experiments with a coupled atmosphere-ocean model, they obtained a linear relationship between changes in the high-latitude temperature ( T insol ) and insolation ( Q): T insol = 3.25× Q/40.The insolation at the youngest minimum preceding 88 kyr BP (at ∼ 95 kyr BP) was about 40 W m −2 lower than at 88 kyr BP (Berger and Loutre, 1991); this yields T insol ≈ −3 • C. Using the colder "minimum insolation" conditions, the extent of the Eurasian ice sheet agrees well with the MIS5b reconstruction in Scandinavia and the Barents Sea region (Fig. 10) -in particular when the North American ice sheet is included (Fig. 10b) -whereas the Kara Sea region continually remains ice free.Hence, in contrast to MIS4 and the LGM, our firstorder sensitivity analysis suggests that the MIS5b extent of the Eurasian ice sheet results from a memory of preceding colder stages rather than the prevailing climate.

Conclusions
We have examined the impact of the geologically constrained MIS5b, MIS4, and LGM ice sheets in North America on the spatial extent of the Eurasian ice sheet.The conclusions are summarized as follows: -The North American ice sheet yields cooler summer temperatures in Europe and warmer temperatures in northeastern Siberia in all time slices.The amplitude of these anomalies and the westward extent of the Siberian warming increase with the size of the North American ice sheet (Fig. 4).
-The temperature anomalies are associated with an equivalent barotropic cyclonic and anticyclonic anomaly in Europe and Siberia, respectively (Fig. 6).
The structure of the circulation anomalies away from the wave forcing is qualitatively consistent with linear barotropic stationary wave theory (see Appendix A).

1237
-Owing to its impact on the Eurasian summer temperatures, the North American ice sheet controls the westward migration of the Eurasian ice sheet; in the presence of the North American ice sheet, the spatial extents of the simulated Eurasian ice sheets at MIS4 and the LGM are consistent with contemporary ice-sheet reconstructions (Svendsen et al., 2004;Kleman et al., 2013).However, if the North American ice sheet is omitted, the Eurasian ice sheet becomes more zonally distributed with a more eastward located center of mass (Figs. 7, 8).
-The stationary wave response to the North American glacial topography is not that sensitive to changes in the ocean heat transport (compare Figs. 4 and 6 with Figs.S4 and S6 in the Supplement).Nevertheless, a weakening of AMOC reduces the influence of the North American glacial topography on the Eurasian ice-sheet evolution by imposing cooler background conditions in Eurasia (Fig. 9).
-Although the spatial extents of the MIS4 and LGM ice sheets are well captured by SICOPOLIS, Eurasia remains essentially ice free for MIS5b.Unlike MIS4 and the LGM, first-order sensitivity analysis reveals that the MIS5b ice sheet was not in equilibrium with the prevailing climate but more likely a result of preceding colder climate conditions.
-Our study suggests that the westward migration of the Eurasian ice sheet between MIS4 and the LGM was induced by the expansion of the North American ice sheet.Furthermore, our results are consistent with the notion that the east-heavy Eurasian ice sheet at the late Saalian Maximum (∼ 140 kyr BP) was accompanied by a relatively small ice sheet in North America (Svendsen et al., 2004;Colleoni et al., 2014).
www Due to the complexity of the atmospheric model, it is useful to resort to a simpler linear framework to obtain a conceptual understanding of the stationary wave field.Linear models have been shown to qualitatively capture the largescale features of the stationary waves in the present-day atmosphere (Charney and Eliassen, 1949;Held, 1983;Held et al., 2002).However, many features omitted in linear models (e.g., zonal variations in the background state and nonlinear interactions between different forcing agents) can significantly alter the stationary wave response (e.g., Cook and Held, 1992;Hoskins and Ambrizzi, 1993;Ringler and Cook, 1997).Therefore, results from linear models should only be considered as a qualitative first-order estimate of the total wave response.The equivalent barotropic structure in Fig. 6d, f, and h suggests that the wave field is dominated by orographic rather than thermal forcing; the latter has been shown to yield stationary waves with a more baroclinic structure (geopotential height anomalies tilt westward with altitude; Hoskins and Karoly, 1981;Ting, 1994;Ringler and Cook, 1999).Therefore, we use the orographically forced linear barotropic model (this is the simplest model that can be used to study meridional dispersion of stationary waves; Held, 1983).
In models linearized about a zonal mean basic state, the horizontal scale of the stationary waves is given by the "stationary wave number" K s , which is a function of the atmospheric background state.In a barotropic model, K s is given by (Held, 1983) where β and a −1 ∂[ζ ]/∂φ are the meridional gradients of planetary and (zonal mean) relative vorticity, [u] is thezonal mean background flow, φ is the latitude, and k and l denote zonal and meridional wave numbers, respectively.In the present-day atmosphere (Hoskins and Karoly, 1981;Held, 1983) as well as in our simulations (Fig. S3), K s monotonically decreases with latitude (as β ∼ cos(φ) → 0 toward the pole).This implies that stationary waves at high latitudes typically have lower zonal wave numbers than those propagating at lower latitudes.Hence, the low wave number response (small K s ) at high latitudes in Fig. 6d, f, and h (> 60 • N) is essentially a result of the spherical geometry of the planet.Following Hoskins and Karoly (1981), the propagation direction of stationary waves is given by the direction of the local group velocity: c g = (c gx i, c gy j ), where i and j are the unit vectors in the zonal and meridional direction, respectively.Because c gy is identical to c gx except for a factor l instead of k (Hoskins and Karoly, 1981;Vallis, 2006),2 the inclination (α) of the ray path (propagation direction) is given by Here, k is constant along a ray; hence, as K s decreases with latitude (Fig. S3), |l| must decrease to satisfy Eq. (A1).This implies that waves at high latitudes propagate along more zonal paths than waves at lower latitudes (Eq.A2); this is also seen in Fig. 6d, f, and h, where the polar wave train is more zonally oriented than the subtropical wave train.Hence, despite the high complexity of the atmospheric circulation model used here, the key features (wave number and orientation) of the polar wave train in Fig. 6d, f, and h -that are associated with the westward migration of the Eurasian ice sheet -are consistent with linear barotropic theory.

1239
The Supplement related to this article is available online at doi:10.5194/cp-12-1225-2016-supplement.

Figure 1 .
Figure 1.Northern Hemisphere topography representative of (a) present day and PI, (b) MIS5b, (c) MIS4, and (d) LGM, based on the ice-sheet reconstructions in Kleman et al. (2013).The shading represents ice sheets and the contour interval is 500 m.

Figure 2 .
Figure 2. Schematic flow diagram (from top to bottom) of the experimental setup.Step 1: three glacial time slice experiments are carried out with the atmospheric model, using two different representations of ocean heat transport and glacial topography.Step 2: monthly climatologies of surface temperature (T s ) and precipitation (Prec) are used to create the climate forcing fields which are interpolated to the ice-sheet model grid.Step 3: steady-state icesheet model simulations are carried out to obtain the ice thickness in Eurasia.A more detailed description of the experimental setup is provided in the main text.

Figure 3 .
Figure 3.The colored shading illustrates the simulated annual mean SST anomalies (LGM-PI; in • C) in the North Atlantic from the CAM3 simulations using (a) PI OHT and (b) BO2009 OHT (Brandefelt and Otto-Bliesner, 2009) as well as (c) the LGM simulation from Brady et al. (2013) using the Community Climate System Model version 4 (CCSM4).The equatorward location of the annual mean sea-ice margin in the respective LGM simulation is depicted by the thick black contours in each panel.The colored markers in each panel show the (annual mean) LGM SST anomaly (LGM-PI) from the MARGO (Multiproxy Approach for the Reconstruction of the Glacial Ocean surface) SST reconstruction (MARGO Project Members et al., 2009).

Figure 4 .
Figure 4. Boreal summer (JJA) surface temperature (in • C) from (a) the ERA-Interim climatology (Dee et al., 2011), (b) the present-day simulation, and the EAonly simulations (with PI OHT) of (c) MIS5b, (e) MIS4, and (g) LGM; green contour denotes the zero-degree surface isotherm.Panels (d, f, h) show the JJA surface temperature anomalies induced by the North American ice sheet (the difference between the fullGlacial and EAonly simulations; in • C) for (d) MIS5b, (f) MIS4, and (h) LGM.The temperature in the glacial simulations (c) to (h) has been projected to the present-day orography using the standard lapse rate (γ = −6.5 × 10 −3 K m −1 ).The dashed black contours depict the outlines of the Kleman et al. (2013) ice-sheet reconstructions in Eurasia and North America.Only significant changes at 95 % (based on Student's t test) are shown in (d, f, h) (nonsignificant changes are displayed in gray).

Figure 5 .
Figure 5. Annual precipitation (in m) from (a) the ERA-Interim climatology (Dee et al., 2011), (b) the present-day simulation, and the EAonly simulations (with PI OHT) of (c) MIS5b, (e) MIS4, and (g) LGM.The annual precipitation anomalies induced by the North American ice sheet (the difference between the fullGlacial and EAonly simulations; in m) are shown in (d, f, h) for (d) MIS5b, (f) MIS4, and (h) LGM.The dashed black contours depict the outlines of the Kleman et al. (2013) ice-sheet reconstructions in Eurasia and North America.Only significant changes at 95 % (based on Student's t test) are shown in (d, f, h) (nonsignificant changes are displayed in gray).

Figure 6 .
Figure 6.Same as Fig. 5 but for the JJA geopotential height anomalies (in m; zonal mean subtracted) at 700 hPa (shading) and 300 hPa (black contours in (d), (f), (h); contour interval is 30 m, and negative values are dashed).Positive anomalies refer to an anticyclonic circulation anomaly and negative anomalies to a cyclonic circulation anomaly.Only significant changes at 95 % (based on Student's t test) are shown in (d, f, h) (nonsignificant changes are displayed in gray).

Figure 7 .
Figure 7. Simulated equilibrium ice thickness in Eurasia (shading; in km) using the PI OHT climate forcing from the EAonly (a, c, e) and fullGlacial (b, d, f) simulations for MIS5b (a, b), MIS4 (c, d), and the LGM (e, f).The dashed black contours depict the outlines of the Kleman et al. (2013) ice-sheet reconstructions.The land area in the simulations is indicated by the brown color and the present-day coastline by the thin black contour.The total Eurasian ice-sheet area in each simulation is indicated in the panel titles.

Figure 8 .
Figure 8.The longitude of the center of mass of the total ice distribution in Eurasia (λ c ) in the Kleman et al. (2013) reconstructions (black bars), EAonly simulations (gray bars), and fullGlacial simulations (white bars).

Figure 9 .
Figure 9. Same as Fig. 7 but using the climate forcing from the atmospheric simulations with BO2009 OHT.

Figure 10 .
Figure 10.Simulated equilibrium ice thickness in Eurasia (shading; in km) using the MIS5b (with PI OHT) climate forcing from the EAonly (a) and fullGlacial (b) simulations in addition to an imposed 3 • C cooling of the JJA surface temperature across the entire domain.
in the time slice simulations.

Table 2 .
Average summer (JJA) temperature in the Northern Hemisphere T NH , in Eurasia T EA (average within the area 20 • W, 180 • E, 45 and 90 • N), and the average latitude of the zero-degree isotherm φ 0 in the PI OHT simulations.