Evolution of the large-scale atmospheric circulation in response to changing ice sheets over the last glacial cycle

. We present modelling results of the atmospheric circulation at the cold periods of marine isotope stage 5b (MIS 5b), MIS 4 and the Last Glacial Maximum (LGM), as well as the interglacial. The palaeosimulations are forced by ice-sheet reconstructions consistent with geological evidence and by appropriate insolation and greenhouse gas concentrations. The results suggest that the large-scale atmospheric winter circulation remained largely similar to the interglacial for a signiﬁcant part of the glacial cycle. The proposed explanation is that the ice sheets were located in areas where their interaction with the mean ﬂow is limited. However, the LGM Laurentide Ice Sheet induces a much larger planetary wave that leads to a zonalisation of the Atlantic jet. In summer, the ice-sheet topography dynamically induces warm temperatures in Alaska and central Asia that inhibits the expansion of the ice sheets into these regions. The warm temperatures may also serve as an explanation for westward propagation of the Eurasian Ice Sheet from MIS 4 to the LGM.


Introduction
Over the last 2.6 million years (Gibbard and Kolfschoten, 2004) Earth's climate has fluctuated between cold and warm periods -glacials and interglacials -characterised by the growth of major ice sheets over the Northern Hemisphere continents during cold periods and their retreat during warm periods. The most recent glacial cycle began some 115 000 years ago (115 kyr BP) following a relative minimum in the Northern Hemisphere summer insolation (Berger and Loutre, 2004). Reconstructions of ice-sheet development through the subsequent 90 kyr (Peltier and Fairbanks, 2006;Stokes et al., 2012;Kleman et al., 2013) show global ice volume increasing in a step-wise fashion, with rapid growth bursts followed by longer periods of stagnation and culminating in the Last Glacial Maximum (LGM) spanning 19-23 kyr BP.
Initially, distinct ice sheets developed in the central Canadian Arctic, Quebec, Scandinavia and the Barents-Kara seas (Kleman et al., 2013), with possible smaller ice masses also in the Canadian Cordillera. An amalgamation process subsequently took place whereby these smaller and spatially separated ice sheets successively coalesced to finally form two massive and coherent ice sheets at the LGM, the Laurentide-Cordilleran ice sheet in North America and the Fennoscandian-Barents Sea ice sheet in Eurasia. At certain times, the mostly independent British-Irish Ice Sheet also formed part of the Eurasian Ice Sheet (Bradwell et al., 2008). The Laurentide was by far the larger of these LGM ice sheets, filling the northern part of North America from east to west and reaching southwards to approximately 42 • N.
The topography of the Northern Hemisphere thus evolved in a rather complex way over the past glacial, starting from the interglacial configuration characterised by two high-elevation areas (Rocky Mountains and the Himalayas), over states with up to four substantial and separated highelevation areas, to the LGM state with two continental-scale ice sheets. Large-scale topography affects climate by exciting stationary Rossby waves, manifested as planetary-scale meanders in the time-mean atmospheric zonal flow (see the review by Held et al., 2002). In the interglacial regime, the Rossby wave train forced by flow over the Rocky Mountains Published by Copernicus Publications on behalf of the European Geosciences Union.
(in conjunction with thermal forcing over the Gulf Stream region, Kaspi and Schneider, 2011) induces northwesterly flow over central and eastern North America, yielding harsh winters there. The topographically driven wave is also largely responsible for the northeastward tilt of the Atlantic jet (Brayshaw et al., 2009), contributing to the milder winters of Europe compared with similar latitudes in North America (Seager et al., 2002).
Much effort has been devoted to understanding the circulation and climate of the LGM, including comprehensive proxy data syntheses (e.g. CLIMAP, 1981;MARGO, 2009, and QUEEN 1 ), combined data-model reconstructions (Dail and Wunsch, 2014) and a large range of modelling studies, notably within the Paleoclimate Modelling Intercomparison Projects 2 (PMIP 1, 2, and 3) (Braconnot et al., 2007). Though there are appreciable model-to-model (Braconnot et al., 2007;Li and Battisti, 2008;Otto-Bliesner et al., 2009;Kageyama et al., 2013a) and model-data discrepancies (Kageyama et al., 2006(Kageyama et al., , 2013bOtto-Bliesner et al., 2009), these studies generally depict an LGM climate substantially different from present. This is especially true in the Atlantic sector, which exhibits pronounced cooling of the northern North Atlantic Ocean, southward displacement of the sea-ice margin and southward-shifted, and, in some studies, a nearly zonally oriented atmospheric jet stream and storm track (e.g. Li and Battisti, 2008;Kageyama et al., 2013a;Ullman et al., 2014). A recent study by Ullman et al. (2014) attributed the massive mechanical forcing of the Laurentide Ice Sheet (in particular the ICE-5G reconstruction used in PMIP2; Peltier, 2004) as a key factor for the zonalisation of the jet. Similarly, in a model-based decomposition of various factors involved in creating the LGM climate, Pausata et al. (2011) ascribed the largest circulation change in the Atlantic region to the mechanical forcing of the Laurentide, rather than to increased albedo or reduced CO 2 .
The circulation during the long build-up phase to the LGM has received much less attention, despite the importance of this time-wise dominant period for understanding how the ocean, atmosphere and cryosphere reorganised as the world transitioned from an interglacial to a full-glacial state. Today, it is well established that important terrestrial glaciation traces can only be understood in the context of glacial configurations predating the LGM and of less than full-glacial size (Ljungner, 1949;Kleman, 1992;Fredin, 2002;Kleman et al., 2008). The significance of this long but less than full-glacial time period was recognised by Porter (1989), who coined the term "average glacial" conditions.
An important unanswered question about the pre-LGM climate is whether the atmospheric circulation characteristics were more similar to those in the LGM or in the inter-glacial. Though smaller than they would eventually become at the LGM, the North American ice sheets nonetheless presented a Rocky Mountains-sized topographic obstacle over eastern North America even in the early stages of the last glacial (Kleman et al., 2013). It is thus conceivable that they could have affected the circulation significantly, bringing it closer to the full LGM regime.
Another key issue is the extent to which atmospheric perturbations induced by pre-LGM ice sheets helped shape the evolution of the ice sheets themselves (Beghin et al., 2014). Studies using idealised coupled atmosphere-ice-sheet models focusing on the dynamics of a single continental-scale ice sheet on a flat continent (Roe and Lindzen, 2001;Liakka et al., 2011;Liakka, 2012) show that stationary waveice-sheet interaction can strongly influence both the spatial form and the temporal evolution of the ice sheet. Remote interactions between widely separated ice sheets mediated by stationary Rossby waves have received less attention (Lindeman and Oerlemans, 1987;Kageyama and Valdes, 2000;Beghin et al., 2014). However there have been suggestions that the North American ice sheets excited a stationary wave train acting to warm northwestern Europe, suppressing ice growth there and potentially explaining why the Eurasian Ice Sheet was considerably smaller than the Laurentide in the latter stages of the glaciation (Roe and Lindzen, 2001). Similar questions arise regarding the documented absence of major ice sheets in places where they could be expected, such as Alaska and eastern central Siberia (Svendsen et al., 2004;Kleman et al., 2013), and the reasons for the general southwestward migration of the Eurasian Ice Sheet through the last glacial cycle (Sanberg and Oerlemans, 1983). Ice core analysis has revealed that the atmospheric dust concentration was considerably higher over the glacial cycle than in the present atmosphere (Mahowald et al., 1999). Theories have therefore been put forth suggesting that dust deposition may have contributed to the absence of major ice sheets in Alaska and eastern central Siberia (Calov et al., 2005a, b;Krinner et al., 2006;Colleoni et al., 2009). Other studies show that changes in the the vegetation cover (Claussen et al., 2006;Colleoni et al., 2009) and SST (sea surface temperature) distributions (Colleoni et al., 2011) may also have contributed to hinder the development of ice sheets in these areas.
Here, we explore these questions using a general circulation model (GCM) to simulate the atmospheric circulation's response to changing topography -taken from the geologically constrained ice-sheet reconstruction of Kleman et al. (2013) -at four representative time slices spanning the last glacial. The modelling strategy is described in Sect. 2. The main results of the study are presented in Sect. 3 and discussed in Sect. 4. Section 5 summarises our conclusions.

Model
We employ the National Center for Atmospheric Research Community Atmospheric Model version 3 (CAM3) (Collins et al., 2004(Collins et al., , 2006, using a spectral dynamical core with T85 (approximately 1.4 • ) horizontal resolution and 26 hybrid sigma-pressure levels in the vertical. Continental surfacesincluding prescribed ice sheets -are represented by the Community Land Model version 3 (CLM3) (Oleson et al., 2004). The ocean is represented by a motionless slab of fixed heat capacity, with ocean heat transport (OHT) represented by a prescribed climatological seasonally varying energy convergence field. The slab ocean also contains a thermodynamic sea-ice model. Further details on the prescription of the ice sheets and OHT are given below.

Ice sheets
Continental ice sheets over North America and Eurasia are prescribed from the recent reconstruction described by Kleman et al. (2013), to which the reader is referred for full details. Briefly, the reconstruction spans the last glacial and employs a numerical ice-sheet model (specifically, the University of Maine Ice Sheet Model (UMISM); Fastook and Chapman, 1989;Fastook, 1993) constrained by geological and geomorphological data. While reasonably reliable data constraints for ice-sheet margins during pre-LGM times are available in some locations, other regions are less well constrained, and there are no constraints at all on ice-sheet elevations and cross-sectional profiles. The strategy employed by Kleman et al. (2013) was to tune the ice-sheet model's mass balance forcing using global model parameters to match well-constrained outline segments as closely as possible where such constraints exist while relying on the model physics to capture ice-sheet elevations and outlines in unconstrained regions. This procedure bears some analogies to data assimilation, where a physically based model is used as an optimal extrapolator or predictor for regions or fields not directly constrained by data.
The resulting evolution of global ice volume is shown in Fig. 1, and is a good qualitative match to previous reconstructions of global ice volume and inferences from sea-level data (Peltier and Fairbanks, 2006;Stokes et al., 2012). Snapshots of the spatial distribution of the ice sheets are shown in Fig. 2 and discussed further below.
UMISM outputs the net surface elevation, combining the topographic height with the ice thickness and isostatic depression of the bedrock due to ice loading, on a rectangular grid with 100 × 100 points, covering from approximately 45 • N to the North Pole. This dataset is interpolated to T85 resolution and combined with the present day topography to create a palaeotopography with global coverage for use in the present simulations. We assume full glaciation of grid cells where the palaeodata is at least 250 m higher than the present day topography.

Ocean heat transport (OHT)
Since a considerable number of runs are required to elucidate the evolution of the circulation over the past glacial, use of a fully coupled ocean-atmosphere model is unfeasible for the present study. Resorting to the more computationally efficient slab ocean model comes at the price of uncertainty regarding the impact of changes in ocean circulation. If the atmospheric circulation is very sensitive to OHT, results obtained with an arbitrarily prescribed OHT may be of little relevance to the real world. To address this problem we use two end-member OHT prescriptions, corresponding to interglacial conditions and to the LGM, respectively, which presumably bracket the range of OHTs over the last glacial. The two OHT specifications are derived from corresponding equilibrated simulations of the fully coupled version of the NCAR model (the Community Climate System Model version 3) under, respectively pre-industrial and LGM conditions (Brandefelt and Otto-Bliesner, 2009). All the simulations, to be described below, were carried out twice, once with interglacial and once with LGM OHT (with modern-day annual-mean mixed layer depths used to specify the slab's heat capacity in all cases). Comparison between these twin simulations permits an assessment of the sensitivity to OHT in each case. As it turns out, there are quantitative differences between the two sets of simulations, but the overall qualitative conclusions of the study are robust to changes in OHT.

Experiments
As shown in Fig. 1, global ice volume during the last glacial grew in three main surges centred around ∼ 110, ∼ 70, and ∼ 25 kyr, while remaining roughly constant for extended periods in between. To capture the main features of this evolution, we focus on four time-slices representative of the main stages of the glacial cycle: -Interglacial (Fig. 2a). This case represents conditions prior to inception of the last glacial and during modern, pre-industrial conditions. It employs modern topog-raphy and pre-industrial greenhouse gas concentrations and orbital parameters.
-Marine isotope stage 5b (MIS 5b, ∼ 88 000 kyr BP; Fig. 2b). Already at this relatively early stage of the glacial cycle, the ice sheets are of a significant size. In Eurasia they cover the Arctic coast from Scandinavia to central Siberia and stand about 2000 m high.
In North America there are two ice sheets, both in the east of the continent -the Keewatin Dome in the eastern Canadian archipelago and the Quebec Dome on the Labrador Peninsula -that are beginning to merge and form a single entity as high and wide as the Rockies. At this stage, the Eurasian and North American ice sheets contain about the same ice volume.
-MIS 4 (∼ 66 000 kyr BP; Fig. 2c). Approximately halfway into the glacial cycle, the ice volume in North America is now roughly twice as large as in Eurasia and consists of the two-domed proto-Laurentide Ice Sheet and a number of smaller freestanding ice sheets in the Cordilleran region. The Keewatin Dome covers most of the Canadian archipelago and northern central parts of the Canadian mainland. In the north, a small gap to the Rocky Mountains still exists, widening to approximately 1000 km further south. On the northeastern and southeastern fringes, towards Baffin Bay and the North Atlantic, ice everywhere reaches the sea. The Quebec Dome is the larger of the two (the figure suggests the opposite but that is due to the map projection) and extend almost as far south as 42 • N over the eastern continent. The Eurasian Ice Sheet is here at its maximum zonal extent. Both the Eurasian and proto-Laurentide ice sheets are in this case about 2500 m high.
-LGM (MIS 2, ∼ 20 000 kyr BP; Fig. 2d). Under full glacial conditions, the Laurentide Ice Sheet is a singledomed continent-wide entity that dwarfs the Rockies. The centre of mass is in the middle of the continent, as in the ICE-5G reconstruction (Peltier, 2004) used in PMIP2. Our reconstruction of the Laurentide ice dome is lower, however, about 3500 m as compared to about 4500 m in ICE-5G. In Eurasia the ice margin has retreated in the east and shifted the centre of mass to Scandinavia, but achieves its greatest meridional extent and even covers most of the British Isles. The total ice volume in the Northern Hemisphere corresponds to about 100 m of sea-level equivalent, of which the Laurentide and Eurasian ice sheets make up roughly 80 and 20 %, respectively. These numbers agree well with the PMIP3 ice-sheet reconstructions.
For each glacial simulation, we set the orbital clock to the nominal time of the ice-sheet reconstruction. We also adjust the concentrations of long-lived greenhouse gases to typical values estimated from the EPICA (European Project for Ice Coring in Antarctica) Dome C ice core (Petit et al., 1999;Spahni et al., 2005); see Table 1. These values are rounded off to the nearest 5 ppmv (parts per million by volume) for CO 2 and 5 ppbv (parts per billion by volume) for CH 4 and N 2 O. The concentrations of chlorofluorocarbons (CFC 11 and CFC 12 ) are identically zero in all simulations and we use the pre-industrial aerosol distribution. All simulations use modern pre-industrial vegetation cover adjusted by the glacier mask.
To test the sensitivity of the results to orbital parameters and greenhouse gas concentrations, we conduct an additional set of simulations in which topography is held fixed at its interglacial distribution while orbital parameters and greenhouse gases take on glacial values. For the MIS 4 case we also run two additional simulations where the Eurasian and North American ice sheets are removed in turn to further in-

Results
This section presents the main results from the simulations of the four time slices discussed above. For brevity we present only results from the interglacial and MIS 5b simulations with interglacial OHT and MIS 4 and LGM simulation with LGM OHT; as noted in Sect. 2.3, changing OHT does not affect our qualitative conclusions. The full set of simulation results, showing the sensitivity to OHT, is presented in the Supplement.

Surface temperature
As shown in Table 3, global-mean surface temperatures continuously decrease across the simulations, dropping by about 5 • C from the interglacial to the LGM irrespective of season. To put this number in perspective with the rest of the modelling community, Braconnot et al. (2007) and Otto-Bliesner et al. (2009) found that the annual global-mean surface temperature in the fully coupled PMIP2 models is between 3.1 and 5.8 • C colder at the LGM compared to the interglacial, where the number for CCSM3 (Community Climate System Model version 3), which uses the same atmospheric component as in this study, is 4.5 • C . More recent generations of models with updated boundary conditions report values between 4.5 and 5.5 • C (see, e.g. Brady et al., 2013;Kageyama et al., 2013a;Ullman et al., 2014).
The annual global-mean surface temperatures for the preindustrial range about 12.2-14.0 • C in the studies mentioned above and the surface temperature in the, alleged, true LGM equilibrium state in CCSM3 (Brandefelt and Otto-Bliesner, 2009, on which we base our LGM ocean Table 2. Configuration of boundary conditions used in the study (marked by ×). The oceanic heat fluxes for interglacial and LGM conditions are denoted by Q IG and Q LGM and the interglacial and reconstructed palaeotopography by IG and Palaeo , respectively. The sensitivity tests of the palaeotopography with only the North American or Eurasian ice sheets present are denoted by NA and EA . heat flux convergence) is 7.9 • C. The corresponding number for IPSL_CM5 with PMIP3 boundary conditions is 7.7 • C (Kageyama et al., 2013a). The values for the pre-industrial and LGM presented in Table 3 are thus well within the uncertainty range of the PMIP models.
In the sensitivity simulations in which the ice sheets are eliminated (indicated by asterisks in the table), the temperature drop is only about 3.4 • C; this implies that the ice sheets themselves account for about one-third of the global temperature decrease from interglacial to LGM. This differs from the situation in Antarctica, where sensitivity studies with the same model used here show no appreciable temperature change when the Antarctic ice sheet is eliminated (Goldner et al., 2013).
The mean Northern Hemisphere Equator-pole temperature gradient generally increases as temperatures cool, in line with expectations of polar-amplified cooling (Singarayer and Valdes, 2010), though interestingly the LGM has a slightly smaller gradient than MIS 4. The gradient grows by approximately 25 % from the interglacial to the LGM in winter and almost 40 % in summer; this leads to substantial strengthening of the midlatitude westerlies with important consequences for the response to topography, as discussed below.
The somewhat larger meridional temperature gradient in MIS 4 compared to the LGM is related to the difference in Earth's orbital parameters at the nominal time of the simulations. The Northern Hemisphere high latitudes in MIS 4 receive more insolation in spring compared to the LGM but the summer and fall insolation is less (not shown), thus rendering the Arctic generally colder over a large part of the year. When introducing the ice sheets this effect is intensified as their high albedo helps to cool regions with a seasonal snow cover in the sensitivity simulations. The annual insolation in the tropics is also slightly higher in MIS 4 compared to the LGM, which acts to further strengthen the meridional temperature gradient. Figure 3 shows the climatological 300 hPa eddy geopotential height field (an eddy is defined here as a deviation from the zonal mean) in boreal winter (December-February, DJF) for the four time slices as well as the ERA-Interim reanalysis (Dee et al., 2011) for a modern day comparison. The corresponding zonal winds are shown in Fig. 4. The interglacial eddy heights in the North American/Atlantic sector (Fig. 3b) show a pattern familiar from present-day observations ( Fig. 3a), with a wave train arcing from western North America to North Africa; this pattern is generally interpreted as a stationary Rossby wave excited by the midlatitude westerlies flowing over the Rockies and then refracting into the subtropics where it is absorbed (Held et al., 2002). The Atlantic midlatitude jet in this simulation (Fig. 4b) is somewhat stronger than in present-day observations (Fig. 4a)   meridionally tilted precipitation field (see, e.g. Hoskins and Valdes, 1990;Lee, 2000;Chang et al., 2002;Held et al., 2002;Vallis and Gerber, 2008;Brayshaw et al., 2009Brayshaw et al., , 2011 for discussions on the existence of stormtracks and how they are influenced by various forcing agents).

Evolution of the winter circulation
The results from the glacial simulations suggest that the interglacial disposition of the stationary waves and zonal jet was maintained for a significant portion of the glacial cycle, despite the massive changes in the boundary conditions. At MIS 5b one can only see minor circulation changes and even though the amplitudes of the stationary height anomalies are generally greater at MIS 4, the geopotential height gradients maintain a similar basic disposition. The Atlantic jet is also stronger in this case, but the axial tilt remains largely unchanged. The LGM, however, presents a distinctly different circulation, especially in the Atlantic region. The trough LGM boundary conditions following the PMIP2 protocol. However, simulations with the updated PMIP3 boundary conditions yield a slightly more tilted Atlantic jet (e.g. Kageyama et al., 2013a;Ullman et al., 2014).
The sea-ice margin (Fig. 2) moves equatorward in both ocean basins as the glacial progresses. However, the North Atlantic remains largely ice-free through MIS 5b, and even in MIS 4 the ice margin has a clear northeastward tilt parallel to the prevailing winds, presumably because warm advection by the winds prevents sea-ice formation off western Europe. It is only at the LGM, when the winds become perfectly zonal, that the ice line reaches south to Iberia.
Note, however, that the zonalisation of the Atlantic jet axis is not principally determined by the southward expanding sea-ice margin. Supplement Fig. S5 shows the glacial evolution of the winter sea-ice and upper tropospheric zonal wind using the two different representations of the OHT. The interglacial OHT supports open waters in the eastern Atlantic also at the LGM, while at the same time the jet has a conspicuous zonal orientation in this region. Note that the sea-ice extent is largely the same as in the MIS 4 simulations with the same OHT where the tilt of the jet axis is maintained over the entire ocean basin. The spatial extent of the zonalised region is however less than in the simulation with the LGM OHT. This suggests that a more equatorward sea-ice margin may help to expand the zonalised region but it is primarily the high LGM ice-sheet topography that causes the shift of the jet axis.
The sea-ice cover obtained in our LGM simulation resembles the CLIMAP (1981) reconstruction which is more extensive than in more recent proxy-data reconstructions, especially in the eastern Atlantic where the perennial sea-ice line is found to not venture past 55 • N (Kucera et al., 2005;De Vernal et al., 2005, 2006MARGO, 2009). Recent modelling studies also show a more moderate Atlantic sea-ice cover than obtained here; see, e.g. Braconnot et al. (2007); Li and Battisti (2008); Pausata et al. (2011). The reason for the very extensive Atlantic sea-ice cover obtained in this study is that our OHT is based on the second LGM equilibrium state discussed by Brandefelt and Otto-Bliesner (2009), which is believed to be close to the model's true steady state. In this state the Atlantic Meridional Overturning Circulation (AMOC) is strongly reduced and sea-ice is therefore able to form over a larger area of the ocean (cf. Fig. 2 in Brandefelt and Otto-Bliesner (2009)).
Consistent with the discussion above, the storm track undergoes a similar transition as the jet stream and retains much of its meridional tilt in MIS 5b and MIS 4 ( Fig. 4). At the LGM, however, it becomes almost completely zonalised over the Atlantic Ocean.
The LGM storm track is also more meridionally confined and the precipitation does not reach as far inland over Eurasia as in the previous cases, likely because the simulated more extensive sea-ice cover in the North Atlantic limits evaporation and yields drier air masses moving into Europe. Neither case shows a significant amount of winter precipitation on the ice sheets themselves. There is some precipitation falling on the southwestern and western edges but almost nothing in the interior. In Eurasia this is due to the ice sheets' northerly location, where the bulk of the precipitation falls south of the main ice sheet. The North American continent is shielded from the Pacific cyclones by the Rockies, making the air in the interior of the continent relatively dry and less prone to precipitate. We will see later that the picture is different in the summer season.
We can make use of the low-level winds in Fig. 2 to give a rationale for the evolution of the planetary waves. The topography of the Rockies yields an upstream ridge over the eastern Pacific as part of the Rossby wave response, adding anticyclonic curvature of the mean flow: as the air flows onto the topography the effective column depth decreases and the mean flow deflects southwards to conserve potential vorticity. On the eastern side of the mountain range the situation is the opposite and the mean flow follows a cyclonic path as a response to the increasing fluid depth (see Appendix A for a discussion on potential vorticity conservation). A consequence of the meridional deflection of the mean flow is that the strongest low-level winds are channelled between the Rockies and the incipient Laurentide Ice Sheet, with little flow normal to the ice-sheet topography. Consequently, the mechanical stationary wave forcing by the ice sheet is small, and the outline of the planetary waves remains largely similar to the present day.
This explanation is applicable to MIS 5b and partly also to MIS 4, as the lowland corridor in the interior of the North American continent is present in both cases. However, in the latter case the westward expansion of the Keewatin Dome in northern Canada partially pinches off the channel, forcing flow over the ice sheet in the north. Further south, the flow crosses the Rockies in a similar fashion as in the interglacial, but overall there is a more substantial interaction between the ice sheets and the mean flow which is reflected in an amplified planetary wave response and an eastward shift of the lee-side trough from the Hudson Bay to the Labrador Sea.
The topographic profile of the LGM Laurentide Ice Sheet is structurally very different from the earlier stages of the glacial cycle. The ice volume is nearly twice as large as at MIS 4 (Fig. 1), and the lowland region east of the Rockies is entirely filled up and instead constitutes the highest part of the continent. Consequently, the topographic influence is much more substantial than earlier, and the low-level wind field over the eastern Pacific (Fig. 2d) is subject to an even stronger meridional splitting. In fact, the entire flow normal to the ice sheet is obstructed and forced on a poleward track. A similar meridional splitting was also found by Manabe and Broccoli (1985) and Cook and Held (1988). Note that the anticyclone over the northeastern Pacific is strong enough to even force low-level easterlies over parts of Alaska and Siberia. This is also hinted in Fig. 3e where the upper tropospheric anticyclone is significantly stronger than earlier and protrudes far into eastern Siberia.
The Eurasian Ice Sheet topography is located too far to the north to be able to interact with the strong midlatitude zonal-mean flow and thereby influence the stationary wave field. In both MIS 5b and MIS 4 the strongest westerly winds are found over central Europe, whereas the southern ice margin is located over Scandinavia. Even though the ice sheet expands equatorward when moving into the LGM, the zonalisation of the Atlantic storm track shifts the westerly winds southwards and the flow-topography interaction therefore remains small.

Evolution of the summer circulation
The reduced Equator-to-pole temperature gradient in the summer season yields a more quiescent troposphere and zonal asymmetries in the circulation are driven to a larger extent by radiative heating contrasts than by mechanical forcing. The outline and vertical structure of the summer stationary waves is therefore quite different from those in winter. The pressure anomalies over the continents resemble the summer monsoon, with thermally induced cyclones near the surface and high pressure anomalies aloft caused by divergent flow at higher levels. The opposite polarity of the eddy geopotential is found over the oceans; see White (1982) and  Ting (1994) for comprehensive discussions on the summer stationary waves. Figure 5 shows the upper tropospheric eddy geopotential averaged over the summer months. One can readily see a progressive development of the wave field over the course of the glacial cycle. Already at MIS 5b there is some amplification of the height anomalies over the northern parts of North America and one can also see an additional wave train propagating over the Atlantic Ocean into Europe. The response does resemble the linear topographic wave, or low mountain wave, discussed by, e.g. Valdes and Hoskins (1991), Cook and Held (1992), and Ringler and Cook (1997). The amplitude of these waves is clearly linked to the size of the ice topography as they are reinforced but stay in approximately the same location when the ice sheets grow larger. The anticyclone associated with the divergent flow in the upper troposphere over the North American continent is gradually being weakened and shifted southward as the glacial progresses. At the same time a secondary high anomaly develops at higher latitudes due to forcing of the ice sheet. This development is particularly evident when transitioning from MIS 5b to MIS 4 where the upper tropospheric anticyclone is partitioned into two distinct centres, and, at the LGM, there is even a band of relatively low geopotential in between, rendering the separation even more apparent. A systematic development of the wave field can also be seen in Eurasia. The upper tropospheric anticyclone in eastern Siberia is both strengthened and shifted westwards in time and becomes a substantial feature in the region. The relative importance of topographic and diabatic forcing of stationary waves is controlled by the magnitude of the low-level winds (Held and Ting, 1990), where a weak mean flow yields a dominance of diabatic over mechanical forcing. A weak flow is also unfavourable for wave propagation as the wave energy tends to be dissipated closer to the source region. The high pressure ridge over the ice sheet is, however, strongly amplified by the diabatic cooling (Ringler and Cook, 1999;Liakka, 2012). To better understand the influence of the ice sheets on the planetary waves the reader is referred to the supplementary material where we show the development of the eddy geopotential field both for the fully forced cases as well as the sensitivity experiments with eliminated ice sheets. Figure 6 shows the spatial distribution of the summer precipitation. The tropics dominate the picture but one can clearly see that the midlatitude storm tracks are weaker and less organised than in winter and the precipitation is therefore dispersed over a larger area. The warm air transports water vapour far inland and precipitates over the ice sheets. One can also see a precipitation maximum on the windward side of the ice sheets that encourages their westward expansion in time. There is also an additional maximum along the southern edge of the Eurasian Ice Sheet that is a re-sult of the anticyclonic circulation in Fig. 5. Likewise, the drier conditions east of the Eurasian Ice Sheet at the LGM comes from advection of cold air in the southward branch of the circulation cell.
To maintain the ice sheets over the warm season the total mass loss from melting and ablation has to be balanced by the accumulation of snow over the year. The question is therefore how the changes in the atmospheric circulation influence the surface temperature field. To investigate, we use the parallel sequence of simulations with identical forcing as in the three glacial simulations but with the interglacial topography; by subtracting the surface temperature in the no-ice-sheet cases from the glacial simulations, we obtain a qualitative idea of how the climate responds to the presence of the large-scale ice sheets and therefore how the ice sheets affect their own maintenance. Note that two forcing agents -topography and surface albedo -change between the no-ice-sheet and glacial simulations. We make no attempt to decompose this further.
Results from the sensitivity experiments are shown in Fig. 7. It is apparent that the presence of the ice sheets has a large influence on both the regional temperature and cloud cover. The Atlantic sector features a progressively expanding area of cooling while there are relatively warm anomalies in Alaska and central Asia (note that similar results are obtained regardless of the OHT; see supplementary material). A similar warm anomaly southeast of the Eurasian Ice Sheet has been found in previous studies of the LGM; see e.g. Manabe and Broccoli (1985), Rind (1987), Felzer et al. (1996), Feltzer et al. (1998. Warmer temperatures in Alaska were also reported by e.g. Otto-Bliesner et al. (2006) and Abe-Ouchi et al. (2007), where the latter compared the LGM and interglacial temperatures from several different models.
The increasingly dominant anticyclonic circulation over Alaska and Siberia implies subsiding motion there and therefore generally reduced cloudiness as seen in Fig. 7. The land surface thus absorbs a larger amount of downwelling shortwave radiation, which leads to warmer temperatures. The anticyclonic circulation also implies warm air advection into the target areas that helps to increase the temperatures even further.
The strongly diminished temperature anomaly in Eurasia in the LGM simulation (Fig. 7e) is somewhat at odds with this picture, however, since the anticyclone is at its strongest. This weakening of the warm anomaly is likely the result of cold westerly temperature advection from the Atlantic (which remains ice-covered through the summer; Fig. 2d) and counteracts the radiative heating in Siberia. This does not happen to the same extent in MIS 4 as the anticyclone is located farther to the east and the temperatures in Europe are not as low as at the LGM.
The colder surface temperature in the Atlantic region is a result of the extensive sea-ice cover. Figure 2 shows that in winter, the North American ice sheets force strong northerly winds that advect cold polar air over the eastern parts of the continent and also the western Atlantic. The cold air advection thus helps to shift the sea-ice line farther equatorward than it would if the ice sheet was not present. In the warm season the extensive sea-ice limits the evaporation over the northern North Atlantic and results in a net reduction of the cloudiness as seen in the right column in Fig. 7. A reduced cloudiness implies more down-welling of shortwave radiation but the melting process of snow and sea-ice is slow and inefficient due to the high surface albedo. The colder air over the ocean is also advected by the westerly winds (Fig. 2) to Europe where it holds down the surface temperature and helps sustain and build the western part of the Eurasian Ice Sheet over the warm season (note that similar temperature anomalies are obtained regardless of the OHT, see supplementary figures). Note that the LGM simulation features an area of increased cloudiness south of the sea-ice margin that is not reflected in the precipitation pattern. A decomposition (not shown) reveals that this is related to a larger fraction of low-level clouds. We make no attempt to evaluate the processes involved in generating this signal.
To further establish the connection between the temperature anomalies and the different ice sheets we have conducted an additional set of simulations for the MIS 4 case where only one of the Laurentide and Eurasian ice sheets are present. The corresponding temperature difference is shown in Fig. 8. It is obvious that most of the temperature response in the western hemisphere and Europe is linked to the presence of the Laurentide Ice Sheet, whereas the relative temperature anomaly in Asia is a result of the circulation change induced by the Eurasian Ice Sheet.

Winter stationary waves
The results presented in Figs. 3 and 4 show some very substantial changes in the winter stationary waves, Atlantic jet stream and storm tracks between the interglacial and LGM. These changes do not come about as a smooth progression through time, however. Our MIS 5b simulation, which is representative of the first 30-40 kyr of the glacial (see Fig. 1), shows a circulation essentially unchanged from the interglacial, despite the considerable ice sheets developing in eastern North America and Eurasia. Stationary wave amplitudes increase during the subsequent 30-40 kyr, represented by our MIS 4 simulation, and the jet strengthens considerably, but the overall structure of the stationary waves remains very similar to that in the interglacial -retaining, in particular, the characteristic southwesterly tilt of the Atlantic jet. It is arguably only in the LGM simulation, representing some 5-10 kyr around the peak of the last glacial, that qualitative changes can be seen, with a full zonalisation of the jet and a perceptible shift towards more zonally elongated stationary wave structures: note that throughout the interglacial, MIS 5b and MIS 4 simulations there are three peaks and troughs in the eddy height field along latitude lines in the midlatitudes (Fig. 3b-d), while there are clearly only two in the LGM simulation.
The shift toward longer stationary waves is consistent with the predictions of simple linear barotropic theory (see Appendix A). A key result of the theory is that the stationary wave number is the zonal-mean zonal wind and β is the meridional gradient of the Corio- lis parameter. As discussed in Sect. 3.1, the Equator-to-pole temperature gradient increases as the LGM is approached resulting in stronger winds (the zonal-mean wind averaged between 800 and 150 hPa increases by about 30 % from 18 m s −1 in the interglacial to 24 m s −1 in the LGM simulation) and thus smaller K s , or longer stationary waves. In addition, the zonal width of the North American orographywhich forces the stationary wave train over the Atlantic -increases very substantially from the interglacial to the LGM, so that the forcing of the stationary waves also shifts towards longer scales.
The muted response to the ice sheets in the early stages of the glaciation is a nonlinear effect. The North American ice sheets begin their development in the trough downwind of the Rockies, with the ice margin roughly parallel to the mean flow, where they cannot efficiently excite stationary waves (as this requires flow normal to the topography). This is a nonlinear -or, more precisely, finite-amplitudeeffect since it is the superposition of the wave generated by the Rockies onto the zonal-mean basic flow that yields the muted response; in the linear approximation, waves are excited only by the basic flow. Note that it is not by chance that the ice sheets develop where they do; it is precisely the northwesterly advection by the pre-existing flow that creates favourably cold conditions over northeastern North America. We therefore expect this to be a robust effect that should be replicated also in other models as long as they faithfully reproduce the large-scale flow over the Rockies.
A similar effect occurs for the Eurasian Ice Sheet, which also fails to excite stationary waves in the early stages of the glacial. In this case, warm southwesterly advection by the Atlantic jet presumably limits the southward extension of the Fennoscandian Ice Sheet, automatically confining this ice sheet to a region of weak climatological westerlies which are unable to excite large stationary waves.
Although the simple linear model captures the main features of the LGM stationary waves with surprising quantitative accuracy, this is likely a fortuitous result from the choice of parameter values and in particular of the model's channel geometry (Held, 1983). GCM simulations by Cook and Held (1992) show that substantial deviations from linearity are possible as topographic height increases, and these deviations are complex and model-dependent (Held et al., 2002). In fact, model intercomparison exercises show a wide range of responses by different models all driven by the same LGM topography (Li and Battisti, 2008); among these, the model used here (CAM3) has a larger-than-average response, with most other models showing an LGM circulation that departs less dramatically from the interglacial (but note that Li and Battisti, 2008, found good agreement between their LGM results using CAM3 and available proxy data). Overall, the conclusion that can be drawn from this discussion is that the atmospheric circulation's response to the Northern Hemisphere ice sheets during winter may have been muted through large parts of the last glacial cycle, and possibly even at the LGM.

Self-induced temperature anomalies in summer
A key factor controlling the evolution of the ice sheets is the summer surface temperature, as it is in the warm season that significant ablation occurs. The temperature at the top of the ice sheets is generally below freezing and the summer ablation is therefore restricted to the marginal ablation zones at lower elevation. Figure 7 shows that the presence of the ice sheets induces warm surface temperature anomalies in Alaska and Siberia. A similar temperature response has been obtained in previous modelling studies of the sensitivity to the LGM ice sheets; e.g. in Manabe and Broccoli (1985), Rind (1987), Felzer et al. (1996), Feltzer et al. (1998, Abe-Ouchi et al. (2007), and Otto-Bliesner et al. (2006). That the model simulates warmer temperatures in these particular regions is intriguing, as we know from geological evidence that Alaska was left largely ice-free over the entire glacial cycle (Clague, 1989). At the same time, the Eurasian ice complex was limited to the northernmost parts of the continent and its centre of mass shifted southwestwards as the glacial proceeded (Svendsen et al., 2004;Kleman et al., 2013). In the simulations the summer surface temperature in Alaska stays approximately the same as in the interglacial climate throughout the entire glacial cycle (the pollen-based temper-ature reconstruction of the LGM by Bartlein et al., 2011, suggests that the surface conditions in Alaska were indeed comparable, if not even slightly warmer than in the present climate) and the region south of the Eurasian Ice Sheet at MIS 4 is even slightly warmer than in the interglacial simulation (not shown). It is thus conceivable that the self-induced warm anomaly in central Siberia may have contributed to limit the glaciation in this region, possibly in conjunction with albedo effects from dust deposition (Mahowald et al., 1999;Calov et al., 2005a, b;Krinner et al., 2006;Colleoni et al., 2009) and the vegetation cover (Claussen et al., 2006;Colleoni et al., 2009). This warm anomaly may also help to explain the westward shift of the ice sheet when moving from MIS 4 to the LGM together with the upslope precipitation in the west as discussed by Sanberg and Oerlemans (1983).

Conclusions
We summarise our conclusions as follows: -During winter, the pre-existing stationary wave response to the Rocky Mountains limits the flow's interaction with the ice sheets during MIS 5b and 4, and the influence of these ice sheets on the tropospheric circulation is therefore surprisingly small despite their very substantial size.
-During MIS 5b and 4, the general outline of the winter stationary waves is largely unchanged from the interglacial and the jet axis in the Atlantic sector retains the characteristic southwest-northeast tilt.
-Only when the Laurentide Ice Sheet is fully developed are the planetary waves strongly influenced and the Atlantic jet axis is zonalised as a result. A highly simplified linear stationary wave model suggests that the increased mean wind speed during the LGM is the key parameter controlling the structural changes in the circulation.
-In the summer season, on the contrary, the ice sheets strongly influence the eddy geopotential field and an anticyclonic circulation develops in Alaska and central Asia as a result. Anticyclonic circulation implies less cloudiness and more shortwave radiation can therefore reach the surface. It is thus plausible that Alaska and central Asia were left largely ice-free due to selfinduced warm anomalies at the surface.
-The North American MIS 5b and 4 ice sheets have a fairly weak remote influence on the conditions over Eurasian Ice Sheet in the present simulations. The most pronounced "downstream response" is cooling over northwestern Europe and warming in northern Siberia during summer. Thus, this remote influence does not straightforwardly explain the relative smallness of Eurasian Ice Sheet, but may have contributed slightly to the westward shift of the ice-sheet-mass centre.
-Based on the results presented in this study we propose that even though the LGM is the most spectacular phase of the last glacial cycle, it only makes up a small fraction of its total duration and is not the best representative phase of the average glacial climate. Most of the last glacial cycle was spent in what can be described as a perturbed interglacial circulation regime.
To quantify this premise we compute the standard deviation of the wave solution in Eq. (A2) (σ = [ψ 2 0 ]) as a function of the zonal mean wind speed when using the interglacial and LGM topographic profiles at midlatitudes (Fig. 9c, d).
We apply the original model settings used by Charney and Eliassen (1949): f 0 = f (45 • N); H = 8 km; and a meridional wave number with the half-wavelength of 35 • in latitude but a weaker damping, and r −1 = 20 days, to make the resonant peaks more distinct. These results are presented in Fig. 9e, clearly showing how the model resonates for lower stationary wave numbers as the wind speed increases. It is apparent that the larger topography at the LGM primarily influence the amplitude of the stationary waves. For wind speeds of approximately 20 m s −1 , the LGM topography yields amplitudes of almost a factor of two larger than the interglacial topography. Figure 6.1 in Held (1983) shows an example of how a similar curve changes for various strengths of the damping timescale. Figure 9a and b show a comparison of the 500 hPa eddy geopotential height from the atmospheric circulation model in the interglacial and LGM simulations with the midtropospheric eddy field computed by the linear model. We here use the average tropospheric zonal mean wind speeds in winter (pressure weighted between 800 and 150 hPa) and a damping timescale of r −1 = 5 days as in Charney and Eliassen (1949). The zonal mean wind speed in the interglacial and LGM simulation is about 18 and 24 m s −1 , respectively. Both the GCM and the simple linear model yield a dominant wave number 3 anomaly in the interglacial case and a wavenumber 2 pattern in the LGM. The amplitudes of the waves are also comparable, suggesting that the Charney-Eliassen model has some predictive skills. However, as discussed by Held (1983), one has to be cautious when interpreting the results from the linear model. The solution is restricted to a one-dimensional barotropic channel flow with no meridional propagation of waves and it only accounts for the topographic forcing. Stationary waves in the real atmosphere tend to follow great circles rather than latitude circles (Hoskins and Karoly, 1981) and additional forcing mechanism and dissipative effects are important for the stationary wave field (e.g. Held et al., 2002).
The Supplement related to this article is available online at doi:10.5194/cp-10-1453-2014-supplement.