This study examines the evolution of a continental-scale ice sheet on a triangular representation of North America, with and without the influence of the Cordilleran region. Simulations are conducted using a comprehensive atmospheric general circulation model asynchronously coupled to a three-dimensional thermomechanical ice-sheet model. The atmospheric state is updated for every 2 × 106 km3 increase in ice volume, and the coupled model is integrated to steady state. In the first experiment a flat continent with no background topography is used. The ice sheet evolves fairly zonally symmetric, and the equilibrium state is continent-wide and has the highest point in the center of the continent. This equilibrium ice sheet forces an anticyclonic circulation that results in relatively warmer (cooler) summer surface temperatures in the northwest (southeast), owing to warm (cold) air advection and radiative heating due to reduced cloudiness. The second experiment includes a simplified representation of the Cordilleran region. The ice sheet’s equilibrium state is here structurally different from the flat continent case; the center of mass is strongly shifted to the east and the interior of the continent remains ice free—an outline broadly resembling the geologically determined ice margin in Marine Isotope Stage 4. The limited glaciation in the continental interior is the result of warm summer surface temperatures primarily due to stationary waves and radiative feedbacks.
The Quaternary period is characterized by the cyclic expansion and retreat of massive ice sheets in the mid- and high-latitude continents (Gibbard and Kolfschoten 2004). The last glacial inception occurred about 115 000 years before present (115 kyr BP), when ice sheets started to form in the central Canadian Arctic; Quebec, Canada; Scandinavia; and in the coastal regions of the Barents and Kara Seas (Svendsen et al. 2004; Kleman et al. 2002, 2013; Stokes et al. 2012). The data records suggest that the ice volume increased in a stepwise fashion over the subsequent 90 kyr, with rapid growth bursts in stadials—cold periods—centered around ~110, ~70, and ~25 kyr BP that were followed by interstadials—warmer periods—when the ice volume remained constant or even decreased slightly (Peltier and Fairbanks 2006; Stokes et al. 2012; Kleman et al. 2013). The stadials are commonly referred to as the Marine Isotope Stages 5b, 4, and 2 (MIS5b, MIS4, and MIS2, respectively), where the latter is the culmination of the glacial cycle at the Last Glacial Maximum (LGM, about 23–19 kyr BP).
The incipient Eurasian ice sheet (EIS) is believed to have been longitudinally extended along the Arctic coast, from Scandinavia in the west to central Siberia in the east, and the whole footprint of the ice sheet appears to have migrated southwestward in time (Fig. 1); see discussion in Svendsen et al. (2004) and Kleman et al. (2013). Reasons for this migration remain unresolved but the Atlantic storm track was likely the primary moisture source for the ice sheet, hence naturally shifting the center of mass toward the European side where the cyclones make landfall. Also, Löfverström et al. (2014) showed modeling results, suggesting that the ice-sheet topography induced warm summer temperatures in Siberia that possibly could help explain this development. The North American ice sheet is believed to have had a rather different growth trajectory. It initially formed in the northeastern corner of the continent and the center of mass appears to have remained in this region from the inception through MIS4 (Fig. 1a)—roughly the first 75–85 kyr of the glacial cycle, though the exact timing is debated (Kleman et al. 2002, 2010; Stokes et al. 2012). The geological data further suggests that the continental interior and the highland regions in the Cordillera were largely ice free over this time period (Kleman et al. 2002; Clague and James 2002; Clague et al. 2005), and the massive LGM Laurentide ice sheet (Fig. 1b), covering the North American continent from coast to coast, is believed to have been present for a mere 5–15 kyr before the subsequent glacial collapse (Kleman et al. 2002, 2010, 2013).
The strong asymmetry of the ice sheets toward the North Atlantic sector is enigmatic (Svendsen et al. 2004; Kleman et al. 2010) and this configuration has also been shown difficult to capture in conventional numerical ice-sheet modeling experiments (Huybrechts and T’siobbel 1995; Marshall et al. 2000; Zweck and Huybrechts 2005), and in experiments with coupled atmosphere–ice-sheet models (Bonelli et al. 2009; Beghin et al. 2014). Hence, important questions regarding the ice build-up patterns are prompted by observations and experiments within both the atmospheric circulation modeling and the glacial geology research communities:
Why was the ice invasion in the western Laurentide area (prairies and interior plains) slow and late compared to the rapid and repeated expansion of the Quebec Dome in the east (Kleman et al. 2002; Stokes et al. 2012)?
What caused Alaska and northeastern Siberia to remain largely ice free over the last glacial cycle (Clague 1989; Svendsen et al. 2004), but at the same time allowed ice-sheet expansion to 40°N over the eastern North American continent?
To what extent did the evolution of the ice sheets influence the atmospheric circulation and induce changing mass-balance patterns? Simpler put, did the ice sheets create their own growing conditions (Sanberg and Oerlemans 1983; Lindeman and Oerlemans 1987; Roe and Lindzen 2001b; Liakka et al. 2012; Liakka 2012)?
The zonally asymmetric ice-sheet development in North America, in particular the southeastward location of the major mass center (Quebec Dome) and the late ice invasion of the prairies, points to an atmospheric stationary wave influence on the buildup of the Laurentide ice sheet, that is, large-scale zonal asymmetries in the time-mean circulation that are forced by heterogeneities in the planetary boundary conditions (Held 1983; Held et al. 2002; Kaspi and Schneider 2011; Ting 1994; White 1982). The present-day winter circulation is characterized by a northwesterly flow of Arctic air over North America that gives rise to temperature differences of several tens of degrees Celsius between the “cold” eastern and “warm” western sides of the continent (Löfverström 2014). This circulation pattern is primarily induced by the flow interaction with the Cordillera, in conjunction with the thermal forcing over the northwestern Atlantic Ocean (Held 1983; Held et al. 2002; Kaspi and Schneider 2011). It is tempting to suggest that these temperature asymmetries can explain the northeastern location of the glacial inception and the asymmetric ice-sheet development. However, it is really the summer temperatures that control where an ice sheet can form, as it is in the warm season that the ablation is most significant. The recipe for a glacial inception is simple: the surface mass balance has to be positive over a long period of time (order of a thousand years), which requires cool summers and abundant precipitation; though straightforward in theory, the triggering mechanism for the last glacial inception is debated and several theories have been put forth.
One of the most widely accepted theories is that proposed by Milankovitch, suggesting that glacial cycles are controlled by variability in Earth’s orbital parameters—eccentricity, obliquity, and axial precession—and the last glacial inception occurred in a period with relatively low summer insolation in the Northern Hemisphere (Loutre 2003; Berger and Loutre 2004). It has also been shown that the “Milankovitch signal” is amplified by changes in the regional surface albedo when tundra- and cold-resistant vegetation types migrate equatorward in colder climates (e.g., Desprat et al. 2005; Kageyama et al. 2004; Colleoni et al. 2014; Vettoretti and Peltier 2003). Other theories involve atmospheric circulation anomalies induced by local and remote topographic and thermal sources (Roe 1999; Huybers and Molnar 2007). Model experiments have shown that they too can induce cold summer temperatures in the inception region and therefore could be part of the explanation.
Numerous modeling experiments, notably within the Paleoclimate Modelling Intercomparison Projects (PMIP 1–3; Braconnot et al. 2007, 2012), suggest that LGM climate was quite different from the present, with larger stationary waves and more zonally oriented jet streams and storm tracks (e.g., Cook and Held 1988; Kageyama and Valdes 2000; Otto-Bliesner et al. 2006; Li and Battisti 2008; Rivière et al. 2010; Kageyama et al. 2013; Brady et al. 2013; Löfverström et al. 2014). These changes are commonly attributed to the massive Laurentide ice sheet (LIS) in North America (Manabe and Broccoli 1985; Pausata et al. 2011; Ullman et al. 2014; Merz et al. 2015), although both the local and global climate conditions are also influenced by the sea surface temperatures (SSTs), the distribution of sea ice, the surface albedo (snow, ice, and vegetation), and the insolation and greenhouse gas concentrations (e.g., Yin and Battisti 2001; Huybers and Molnar 2007; Pausata et al. 2011; Claussen et al. 2006; Colleoni et al. 2009). However, in a recent paper, Löfverström et al. (2014) studied the interactions between first-order paleotopography (only considering the general outline and elevation of the paleotopography) and the evolution of the planetary-scale atmospheric circulation over a sequence of snapshots of the buildup of the last glacial. They argued that the extensively studied LGM climate only reflects a rather short-lived glacial extreme, not necessarily typical for most of the glacial cycle (Porter 1989). Hence, because ice buildup is a slow precipitation-limited cumulative process, the LGM ice configuration with a continent-wide LIS may potentially have been the result of long periods of “non-LGM-like” circulation conditions (Löfverström et al. 2014).
For a more complete understanding of the interaction between topography, stationary waves, and the ice-sheet evolution, the coupled atmosphere–ice-sheet system must be studied over a sequence of conditions that realistically mimics the whole build-up phase. To date, relatively few studies exist that investigate glacial states prior to the LGM: Jochum et al. (2012) studied the climate conditions of the last glacial inception in a fully coupled atmosphere–ocean model and Gregory et al. (2012) simulated the inception and early developments using a circulation model coupled to an ice-sheet model. However, because of the length of a glacial cycle (typical time scales involved are on the order 40–100 kyr; Lisiecki and Raymo 2005), transient simulations over the whole build-up phase from inception to glacial maximum are rare and the modeling strategy often involves a three-dimensional ice-sheet model forced by climate data from a simplified atmospheric model (Roe and Lindzen 2001a), an energy balance model (Tarasov and Peltier 1999), or time-interpolated data between snapshot simulations of two glacial extremes (typically the interglacial and the glacial maximum; Charbit et al. 2007). Some studies have also utilized an ice-sheet model coupled to a simplified but interactive atmospheric circulation model (Roe and Lindzen 2001b; Liakka et al. 2012; Liakka 2012) or a low-resolution intermediate complexity model in which almost all atmospheric processes are parameterized (Bonelli et al. 2009; Beghin et al. 2014).
The studies by Roe and Lindzen (2001a,b), Liakka et al. (2012), and Liakka (2012) explored the mutual interactions between atmospheric stationary waves and the buildup of the Laurentide ice sheet but did not explicitly consider the topographic influence of the Cordillera. However, in Roe and Lindzen’s (2001b) simulations the (linear mechanical) stationary wave feedback alone was capable of generating a southeast-heavy ice sheet, to first order similar to the pre-LGM configuration that is documented by geological and geomorphological data (Kleman et al. 2002, 2010). Here we address the same problem using a comprehensive atmospheric general circulation model coupled to an ice-sheet model. We resort to an idealized zonally symmetric aquaplanet setup with a simple triangular continent representing North America. Although this setup takes a step back from reality, it is deemed a necessary compromise in order to isolate and understand how asymmetries in the ice-sheet development may have been induced by local interactions over the continent. We run two transient simulations with the coupled atmosphere–ice-sheet model. In the first experiment we let the ice sheet evolve on an initially flat continent where all structural asymmetries are self-induced by the ice-sheet topography—this setup is similar in spirit to Roe and Lindzen (2001b), Liakka et al. (2012), and Liakka (2012). In the second experiment we make the setup slightly more representative for North America and also include the Cordilleran range along the continent’s west coast. Both experiments are carried out with an asynchronous coupling between the atmosphere and ice-sheet models, and extend from the glacial inception to the ice sheet’s equilibrium state.
The North American Cordillera consists of several individual mountain ranges partitioned into three main belts: the Pacific Coast Ranges in the west, the central Nevadan belt, and the Laramide belt in the east, where the Rocky Mountains are part of the latter. However, for a simpler nomenclature we refer to the entire highland region as the “Rockies” in the rest of this paper. This is also the typical naming convention used in meteorological literature.
2. Models and experiment design
a. Atmospheric model
We employ the National Center for Atmospheric Research (NCAR) Community Atmosphere Model, version 3 (CAM3.0; Collins et al. 2004, 2006), using a spectral dynamical core with T42 (~2.8° × 2.8°) horizontal resolution and 26 hybrid sigma-pressure levels in the vertical. Land surface processes are represented by NCAR’s Community Land Model, version 3 (CLM3; Oleson et al. 2004).
The lower boundary is a simplified aquaplanet with a single triangular-shaped continent representing North America. This idealized model geometry is used to isolate the local interactions between the atmospheric circulation, the background topography, and the evolving ice sheets. The land surface is prescribed with bare ground conditions, where the soil temperature and moisture are allowed to evolve freely. The ocean is represented by a fixed zonally symmetric SST distribution, using Eq. (1) in Neale and Hoskins (2000), but modified to include a seasonal cycle:
where ϕ and λ denote the latitude and longitude, respectively, and
describes seasonal changes in the SST field. Here the warmest (coldest) Northern (Southern) Hemisphere SSTs occur in August (). A representation of a seasonal cycle is necessary to yield a more realistic climate for the ice-sheet development, as the annual surface mass balance is strongly dependent on the summer ablation. The seasonal shift of the latitude of maximum SST () is represented by the parameter , which yields a maximum summer-to-winter difference of approximately 10°C. The parameter represents the latitude of the sea ice edge for annual-mean conditions. Ocean grid cells poleward of the sea ice edge are here assumed to be covered by a 2-m uniformly thick ice cover. Note that the perpetual annual-mean conditions in Neale and Hoskins (2000) is given by and . The more equatorward sea ice extension used here is motivated by the glacial context of the study. Unlike the prescribed SSTs, the surface temperature over land and sea ice is dynamic and governed by the surface energy balance.
The sea surface climate is somewhat similar to the CLIMAP (1976) LGM reconstruction (used in the initial phase of PMIP1), which has a largely zonally symmetric SST field and a sea ice cover extending well into the midlatitudes in both the Pacific and Atlantic basins; the CLIMAP Atlantic sea ice edge is at approximately 40°N in the winter season [December–February (DJF)].
The orbital parameters and greenhouse gas concentrations are specified for preindustrial conditions: CO2 = 280 ppmv (parts per million by volume), CH4 = 760 ppbv (parts per billion by volume), and N2O = 270 ppbv throughout our simulations. Note that the preindustrial (recent past) summer insolation is slightly lower than the glacial mean (see Fig. 1 in Löfverström et al. 2014), thus making it a decent representative for the “average” glacial insolation. However, the insolation (orbital forcing) is known to vary strongly over time scales of a glacial cycle and the last glacial inception occurred in a period with relatively low summer insolation at high latitudes in the Northern Hemisphere (Loutre 2003; Berger and Loutre 2004).
b. Ice-sheet model
We use the three-dimensional land-based ice-sheet model Simulation Code for Polythermal Ice Sheets (SICOPOLIS), version 2.9, which treats ice as an incompressible, viscous, and heat-conducting fluid (Greve 1997). The model uses a ⅓° horizontal resolution and 81 levels in the ice and 11 levels in the bedrock. SICOPOLIS solves the time-dependent equations governing the thickness, flow velocity, and temperature for grounded ice in response to external forcing, which is given by the surface mass balance, the eustatic sea level, and the geothermal heat flux. The model equations are subject to the shallow-ice approximation, which means that only the lowest-order terms are accounted for (Hutter 1983). It uses Glen’s flow law relating the strain rates to the third power of the applied stresses with the same flow parameters as in Greve et al. (1998). A Weertman-type sliding scheme is used to calculate the basal flow velocities (Weertman 1957, 1964). It is also assumed that the bedrock and ice sheet relax toward isostatic equilibrium with a time scale of 3 kyr and instant calving (zero ice thickness) is assumed at the continental margins.
We use the default parameter settings (see Greve 1997) with the following exceptions: (i) ice is simulated in the “cold ice mode,” which means that if the basal melting temperature exceeds the pressure melting point, it is artificially reset to the pressure melting point; (ii) for simplicity—similar to Liakka et al. (2012) and Liakka (2012)—we use a spatially uniform geothermal heat flux of 55 mW m–2; and (iii) the sea level is assumed to be constant.
The ice sheet’s surface mass balance is defined as the difference between accumulation (precipitation) and ablation (melting). The ablation is parameterized using the positive degree-day (PDD) approach developed by Braithwaite (1985), Reeh (1991), and Calov and Greve (2005). In their approach, the sum of all PDDs in a year (integrated sum of all temperatures above freezing) is calculated semianalytically using the monthly mean temperatures and a standard deviation of temperature representing day-to-day variability. The standard deviation of temperature is set to 4.5°C, and the degree-day constants that relate the PDDs to actual melt are set to 3 mm day–1 K–1 for snow and 12 mm day–1 K–1 for ice. Following Reeh (1991), the available PDDs are used to melt snow and ablate ice in the following order: as a first step, the PDDs are used to melt the snow layer on top of the ice sheet. The meltwater (also liquid precipitation) percolates down into the snowpack and refreezes as superimposed ice. However, the saturation factor for the snowpack is set to 60%, and surplus liquid water is assumed to run off the ice sheet immediately. The remaining PDDs (if any) are then used to ablate the superimposed ice and, as a last step, the glacier ice.
The ice-sheet model uses a similar but smaller continent as the atmospheric model that is truncated to not include the westernmost part where the Rockies are located. Hence, the Pacific west coast in the ice-sheet model coincides with the easternmost points of the Rockies in the atmospheric model. The ice-sheet model thus uses exactly the same continental domain in both experiments regardless of whether or not we include the Rockies (see description of the experiment design in section 2c). This is important because we want to isolate the effect of the mountain range on the atmospheric flow without altering other feedbacks affecting the ice growth. Allowing the ice sheet to freely evolve on the Rockies, it would initially grow much faster on the mountain range than elsewhere owing to the high altitude (Bonelli et al. 2009). Also, since the Rockies remained largely ice free during most of the build-up phase of the last glacial (Clague 1989; Clague and James 2002; Clague et al. 2005), ice growth on the mountain range would potentially yield an unrealistic forcing of the atmospheric circulation. To localize the glacial inception, we prescribe a 500-m-high Gaussian hill in the northeastern corner of the continent.
c. Coupling procedure and experimental design
To examine the influence of the Rockies on the downstream ice-sheet evolution, we conduct two experiments with the coupled atmosphere–ice-sheet model. The first experiment uses a flat triangular representation of North America, an experimental setup similar in spirit to Roe and Lindzen (2001b), Liakka et al. (2012), and Liakka (2012), though they all used a flat rectangular continent. The second experiment setup is identical except for a prescribed 2500-m-high and 15°-wide Gaussian mountain range along the continent’s west coast in the atmospheric model, which represents the Rockies (more precisely the North American Cordillera).
Since the growth rate is faster when the ice sheet is far from equilibrium, we update the atmospheric state based on the ice volume instead of a fixed time interval. We motivate the use of ice volume as a proxy for the ice sheet’s size because the climatological atmospheric temperature anomalies forced by the ice sheet are essentially proportional to the total ice volume, at least for small- to intermediate-sized ice sheets (Liakka and Nilsson 2010). The atmospheric fields are updated for every 2 × 106 km3 change in ice volume, implying that the coupling frequency varies from approximately 500–1000 yr–1 for the smaller ice sheets to about 5000 yr–1 when the ice sheet is close to equilibrium (see Fig. 2).
At each coupling event we compute an atmospheric monthly climatology (based on 10 years of data) using the ice-sheet topography and glacial mask from the last time step in the ice-sheet model integration. The fields are bilinearly interpolated between the model grids.
a. Ice-sheet evolution
Figure 2 shows the temporal evolution of the ice volume in the two experiments. The ice sheet grows considerably faster on the flat continent and the equilibrium ice volume is greater by almost a factor of 2, approximately 84 × 106 km3 compared to 44 × 106 km3 in the simulation with the Rockies.
Figure 3 shows a selection of snapshots of the ice sheet’s spatial evolution in the two experiments (the specific time slices are indicated in Fig. 2). The initial phase is similar in both experiments: the glacial inception is localized by the Gaussian hill in the northeast and the embryonic ice sheet rapidly expands westward and becomes continent-wide after only a few centuries. The almost instantaneous glaciation of the high latitudes is the result of low ablation rates (Figs. 4a,b) due to subfreezing (<0°C) surface temperatures in the summer season [June–August (JJA); see Figs. 5a, 6a], which allow a perennial snow cover to form along the Arctic coast (not shown).
However, subsequent to the glacial inception, the ice sheets’ growth trajectories begin to diverge significantly. On the flat continent (left column in Fig. 3), the incipient/early ice sheet evolves fairly zonally symmetric (note the largely zonal ice ablation in Fig. 4a) but with a slightly east-heavy disposition. The early and intermediate stages shown in Figs. 3c,e bear some structural resemblance to the equilibrium ice sheet in Roe and Lindzen (2001b; cf. their Fig. 13). As time progresses, a discernible double-dome structure emerges with focal points on either side of the continent, and the southern ice margin gradually becomes more zonally asymmetric. In the latter stages of the simulation, the center of mass is shifted to a more central part of the continent and the ice sheet equilibrates as a largely symmetric and continent-wide monodome. The symmetric structure of the equilibrium ice sheet is somewhat similar to the nonlinear results in Liakka et al. (2012) and Liakka (2012). Note, however, that the southern ice margin obtained here reaches much farther to the south than in those simulations, where the ice margin is located at approximately 40°N.
The right column in Fig. 3 shows the structural evolution of the ice sheet when accounting for the atmospheric flow interactions with the Rockies. The southern margin of the incipient/early ice sheet has a more pronounced northwest–southeast tilt than in the flat continent case (cf. Figs. 3c,d) and the center of mass is strongly shifted to the east. This configuration remains throughout the entire simulation as the ablation rate is high in the continental interior (Figs. 4b,d,f) and the ice sheet equilibrates with a high dome in the east, and a large ice-free region in the lee of the Rockies. This outline is structurally similar to the geologically determined margin of the pre-LGM ice sheet in MIS4 (outlined in Fig. 1a). In section 4 we discuss mechanisms and processes that are omitted in these experiments that may be of importance for building the full continent-wide Laurentide ice complex from the MIS4 ice sheet when also accounting for the North American topography.
b. Atmospheric response
1) Flat continent
The left column in Fig. 5 shows the climatological state of a number of meteorological fields on the ice-free continent. We focus primarily on the climate conditions in boreal summer (JJA) as it is in the warm season that most ablation occurs. However, we show the annual-mean precipitation (Fig. 5d) since precipitation accumulates and (superimposed) ice forms continuously over the year in the interior of the ice sheet (away from the marginal ablation zone seen in Fig. 4). The symmetric sea surface conditions (SST and sea ice cover) and the absence of topography on the ice-free continent implies that the climatological atmospheric state is essentially zonally symmetric with only small zonal deviations over the continent. This is particularly true for the surface temperature (Fig. 5a), the surface pressure (Fig. 5g), and the precipitation (Fig. 5d), although the latter has a more pronounced longitudinal asymmetry with higher values in the west as the precipitation from the Pacific storm track gradually decreases away from the ocean source. The summer stationary waves are typically weak and largely localized to their source region (Fig. 5j). In the absence of topography, the only source of stationary Rossby waves is the thermal contrast between land and ocean due to radiative heating of the land surface (Fig. 5m). This forcing is presumably weak on the ice-free continent as there are relatively small zonal variations in the surface temperature field (Fig. 5a).
The middle column of Fig. 5 shows deviations of the meteorological fields with respect to the initial state when the ice sheet has grown to an intermediate size (~50% of the equilibrium ice volume). To assess the temperature response due to circulation changes, we have eliminated the elevation effect (projected the surface temperature at sea level) by assuming a lapse rate of 6.5 × 10–3 K m−1. The temperature response to the ice sheet is a relative warming in the west and a cooling in the east (Fig. 5b). Note that the actual surface temperature is below freezing in the warm region due to the elevation effect. A similar disposition of the temperature anomalies over the North American ice-sheet topography has been found in previous studies, including both linear and nonlinear atmosphere circulation models (e.g., Cook and Held 1988; Roe and Lindzen 2001b; Liakka et al. 2012; Liakka 2012).
The temperature response is associated with an anticyclonic circulation (both at the surface and in the midtroposphere; Figs. 5h,k) with a slightly baroclinic structure, that is, tilting westward with height. The disposition of the stationary wave anomalies over the ice sheet has a different pattern than the streamfunction response obtained in linear orographic stationary wave models; such models typically yield a relatively weak anticyclone over the western part of the ice sheet and a stronger downstream cyclone (e.g., Roe and Lindzen 2001b), and the streamfunction anomalies have an equivalent barotropic structure (i.e., proportional to the temperature anomalies; Held 1983; Held et al. 2002).
Mechanical (orographic) forcing of stationary waves requires a westerly mean flow normal to the height contours (see, e.g., Held 1983). Here the low-level mean winds over the ice-free continent (Fig. 5g) are westerly on the equatorward side and easterly on the poleward side of the 0°C surface isotherm, which is indicated by the thicker contour in Fig. 5a (note that this wind pattern is somewhat artificial and is likely a result of the zonally symmetric boundary conditions; see discussion in section 4). For the intermediate-sized ice sheet, this implies that only the southernmost part of the ice-sheet topography is exposed to a westerly mean flow; hence, it is plausible that the presence of easterly winds weakens the mechanical stationary wave forcing in the early part of the simulation. The wave response in Fig. 5k shares many structural similarities with the results in Ringler and Cook (1999) and Liakka (2012). They argued that diabatic cooling and mechanical forcing couples nonlinearly and enhances the anticyclonic response over the ice sheet; for example, Liakka (2012) found that the perturbation winds induced by diabatic cooling amplify the mechanically forced stationary waves by increasing orographically induced vertical velocities. The contribution from thermal forcing also helps explain why the circulation anomalies have a baroclinic structure in our simulations (see Hoskins and Karoly 1981) rather than an equivalent barotropic structure, which typically results from mechanical forcing in isolation (Held 1983; Held et al. 2002).
The anticyclonic circulation over the ice sheet also influences the cloud cover and thus the surface radiation. We find that the amount of low- and medium-height clouds decreases, whereas the amount of high-level clouds increases with respect to the ice-free conditions. This means that the downwelling shortwave radiation at the surface increases, which has a net warming effect over the ice-free part of the continent (west of the ice sheet; Figs. 5b,n). Note that the amount of downwelling shortwave radiation also increases over the ice sheet but the net shortwave radiation (quantity shown) is reduced due to the high surface albedo.
The right column of Fig. 5 shows the climate response to the equilibrium ice sheet (100% ice volume). The warm anomaly on the western side of the continent is amplified and the warm temperatures also extend out over the sea ice in the Pacific Ocean. As the ice sheet approaches equilibrium, the temperature anomalies over the ice sheet obtain a more northwest–southeast alignment—that is, a clockwise shift relative to the intermediate-sized ice sheet (cf. Figs. 5b,c). This rotation is consistent with the results in Cook and Held (1992), Ringler and Cook (1997, 1999), Liakka et al. (2012), and Liakka (2012). They all ascribed this to an increased influence of nonlinear eddy advection that becomes important when the magnitude of the eddy winds approach, or even exceed, the zonal-mean wind. The stronger meridional deflection of the low-level winds over the ice-sheet topography (cf. Figs. 5h,i) suggests that this may be the case also in our simulation. The strength of nonlinear eddy advection in the stationary wave response is controlled by many different factors—for example, surface friction, the aspect ratio of the topography (latitude vs longitude extension), the height of the topography, the meridional temperature gradient, the strength and vertical structure of the mean wind. For a detailed discussion on this topic, see Ringler and Cook (1997).
The northerly location of the ice sheet relative to the mean flow also implies that most precipitation falls on the southern and southwestern slopes (Fig. 5e). This contributes to the development of the fairly zonally oriented southern ice margin in the early stages of the simulation (Fig. 3c). However, when a strong anticyclonic circulation cell is established over the ice-sheet topography, the poleward flow branch yields a precipitation maximum on the western and southwestern slopes. The zonal-mean flow also advects moist air south of the ice margin; hence, a significant amount of precipitation falls on the (southwestern) slopes on the ice sheet’s equatorward outcrop on the eastern side of the continent (Figs. 5e,h). This, together with the imposed calving at the continental boundary, contributes to building the double-dome structure seen in Figs. 3e,g. Note that the equilibrium ice sheet expands into a region where the ablation rate on the ice-free continent is of the order several tens of meters per year (Fig. 4a).
2) Continent with Rockies
The left panels in Fig. 6 show the climatological state of the same meteorological fields that were presented in the left column in Fig. 5, but for the initial ice-free state in the Rockies simulation.
The presence of the mountain range has a noticeable influence on the continental climate: the summer (JJA) surface temperature field (Fig. 6a) is considerably more zonally asymmetric compared to the flat continent case (cf. Fig. 5a), with warm temperatures extending poleward along the ridgeline and over the continental interior. This temperature pattern is primarily attributed to an increased net shortwave radiation at the surface due to more arid conditions (less clouds) in the lee of the mountain range (cf. Figs. 6d,m and 5d,m, clouds are not shown). The Rockies also influence the circulation in the free atmosphere: the midtropospheric eddy streamfunction shows a meridionally oriented anticyclone–cyclone dipole sitting over a highland region—a circulation feature known from the real world that is generally more pronounced in winter than in summer (Ringler and Cook 1999; Held et al. 2002). It is conceivable that the symmetric ocean boundary condition and the meridional SST gradient in the Pacific Ocean yields an anomalously strong low-level mean flow impinging on the Rockies, and, consequently, a stronger mechanical stationary wave forcing than in the present-day summer. The prescribed sea surface conditions may also yield a diabatic heating distribution that is more similar to the present-day winter, which further strengthens these circulation anomalies (Held et al. 2002). A discussion on the shortcomings of the idealized background climate follows in sections 3c and 4.
Figures 6b,c show the surface temperature (extended to sea level) relative to the ice-free state (Fig. 6a) for an intermediate-sized ice sheet (~50% of the equilibrium ice volume) and the equilibrium state (100% ice volume), respectively. Both panels exhibit a significant warming in the central parts of the continent, with temperatures up to 10°C higher than on the ice-free continent; the average summer temperature is about 15°–20°C compared to 10°–15°C in the ice-free state and the annual ablation is of the order 10–30 m yr–1 (Figs. 4d,f). The warm anomalies coincide with the region that remains ice free over the entire simulation and thereby prevent the ice sheet’s westward expansion in time, despite the abundant precipitation falling on the ice dome on the eastern side of the continent (Figs. 6e,f).
Thermal forcing of stationary waves is important for the regional climate over the North American continent in summer: radiative heating of the land surface creates a convergent flow in the lower troposphere as warm air is rising (low-level cyclone). The ascending air is slightly divergent away from the surface, implying that a thermal surface cyclone is overlaid by an anticyclone in the mid- and upper troposphere (White 1982; Ting 1994). However, in our ice-free simulation the mid- and upper-tropospheric circulation is dominated by mechanically forced waves induced by the westerly flow impinging on the Rockies; see Fig. 6j. In the glaciated cases there is a low-level cyclonic response in the ice-free region of the continental interior and an anticyclonic response over the easternmost parts of the ice sheet (in the midlatitudes). At upper levels the local response is everywhere dominated by an anticyclone.
The low-level diabatic heat field is intricate—a surface warming in the ice-free region and a low-level cooling over the ice surface (Figs. 6n,o)—and the wave response is further complicated by mechanical forcing induced by the westerly flow impinging on the ice sheet’s southwestern slopes [shown by Ringler and Cook (1999) and Liakka (2012) to yield a nonlinear interaction]. The surface warming of the ice-free part of the continent is associated with a 30%–50% reduction of the low- and medium-height clouds compared to the ice-free state, whereas the fraction of high clouds is increased by almost a factor of two (not shown). A reduction of low- and medium-height clouds, and an increase of high-level clouds, have a net warming effect on the (local) summer climate (Pierrehumbert 2010). This contributes to strengthening the thermally driven circulation in the ice-free part of the continent (Figs. 6h,i,k,l).
A consequence of the complicated structure of the stationary wave forcing fields is that the wave response is somewhat different from Hoskins and Karoly (1981) and Ringler and Cook (1999), but the baroclinic vertical structure is similar. To investigate the relative roles of thermal and mechanical forcing (and the nonlinear interactions between those) requires a hierarchy of linear and nonlinear stationary wave models that allows for an examination of the relative contribution from these wave sources in isolation. However, such investigation is out of the scope of the present study.
c. Sensitivity experiments
The idealized boundary conditions used in these simulations implies that many asymmetries in the time-mean circulation are missing. For the purpose of these experiments, perhaps the most important difference from the real world is the absence of a strong surface anticyclone over the North Atlantic in summer, which gives rise to a low-level flow of warm and moist air from the Gulf of Mexico over the interior of the continent. To investigate to what extent the simplified boundary conditions influence our main conclusions, we compare our climate response with that from the more realistic MIS4 experiment presented in Löfverström et al. (2014). They used a slab-ocean version of NCAR-CAM3.0 (same atmospheric model as used here) with insolation and greenhouse gas concentrations appropriate for the MIS4 climate state and the geologically constrained (static) surface topography shown in Fig. 1a (Kleman et al. 2013).
Figure 7 shows the difference in summer climate between the full MIS4 simulation and a sensitivity simulation with identical forcing but modern-day topography (i.e., no ice sheets). This comparison shows the climate response to the ice-sheet topography and is thus similar to that presented in the right column in Fig. 6. An important difference from our simulations is that the ocean surface in Löfverström et al. (2014) is allowed to respond to changes in the atmospheric forcing; hence, the sea surface conditions are not necessarily identical in the glaciated and nonglaciated states. This difference is small, however, and is not likely a primary driver of the climate response.
The first-order patterns are largely similar to our simplified simulation (cf. Fig. 7 and the right column in Fig. 6): there is a warm surface temperature anomaly in the region between the ice sheet and the Rockies (though shifted somewhat toward the northwest) that is driven by an increased net shortwave radiation at the surface resulting from changes in the cloud cover—an increase of high-level clouds and a reduction of low- and medium-height clouds. The time-mean atmospheric circulation also has a similar structure with a surface anticyclone over the ice-sheet topography and a westward-tilting anticyclone driven in part by diabatic cooling and mechanical forcing of stationary waves.
Though the radiation signal is of the same order of magnitude as in the idealized simulation (several tens of watts per square meter), the implied temperature signal is somewhat weaker. This is attributed to a slightly stronger northerly component of the low-level flow in the corridor between the Rockies and the ice sheet (in the full MIS4 simulation, not shown), which implies cold air advection that counteracts the radiative heating. Reasons for the changes in the low-level flow remain unknown but could possibly be the result of both local and remote flow interactions with topographic and thermal sources. We make no attempt to investigate this further.
Despite the muted temperature signal, the important message is that the climate response with a realistic model configuration is structurally similar to what is obtained in the idealized setup. This strengthens the credibility of the idealized experiment and further suggests that we capture the first-order effects that may have been responsible for the asymmetric ice-sheet development. Similar experiments, however, should be carried out in different circulation models to validate our results.
4. Summary and discussion
We have used an asynchronously coupled atmosphere–ice-sheet model to examine the influence of the Rockies (Cordilleran region) on the evolution of a continental-scale ice sheet in North America. Two experiments were conducted using an aquaplanet setup with a simplified triangular continent representing North America with (i) no background topography (similar to Roe and Lindzen 2001b; Liakka et al. 2012; Liakka 2012) and (ii) a semirealistic representation of the Rockies consisting of a 2500-m-high Gaussian mountain range along the continent’s west coast. Our experiments present a number of improvements over the aformentioned studies, most notably the use of a comprehensive atmospheric general circulation model, including moisture dynamics and a prognostic cloud water parameterization (see Collins et al. 2004, 2006, and references therein). We also account for the ice sheet’s diabatic cooling of the atmosphere, and we use a slightly more realistic continent with background topography.
a. Flat continent
The early ice-sheet development on the flat continent is qualitatively similar to Roe and Lindzen’s (2001b) equilibrium state, with a more equatorward ice expansion on the continental margins but with a slightly east-heavy disposition (cf. Figs. 3c,e with Fig. 13 in Roe and Lindzen 2001b). As time progresses the center of mass gradually migrates westward and the equilibrium state is a largely symmetric monodome, covering virtually the entire continent. The thermal and mechanical forcing of stationary waves excite an anticyclonic circulation that yields warm air advection over the western part and cold air advection over the eastern part of the ice sheet. This circulation pattern is strengthened when the ice sheet expands and it also rotates clockwise, suggesting that nonlinear interactions become more important (Cook and Held 1992; Ringler and Cook 1997, 1999; Liakka et al. 2012).
The simplified and largely zonally symmetric boundary conditions give rise to a well-defined baroclinic zone over the oceans—a low-pressure region in midlatitudes—that has no direct counterpart in the modern-day summer hemisphere. Consequently, the balance between the Coriolis force and the meridional pressure gradient force (directed from the higher pressure in the polar region) yields an artificial easterly geostrophic background flow that extends into midlatitudes (Fig. 5g). This flow and pressure configuration resembles that over Antarctica in the austral summer (December–February; see, e.g., Hartmann 1994). A similar background flow was reported by Williams and Bryan (2006), who studied glacial winds in an axisymmetric aquaplanet model. They found that low Arctic temperatures and sea ice expanding into the midlatitudes yields an artificially strong eddy-driven jet (baroclinicity) and polar easterlies migrating equatorward from the Arctic region.
Even though the flat continent experiment has little relevance for the real glacial world—the ice sheet is expanding too far to the south, we have surface easterlies in high latitudes, and we omit the Rockies that shield the continent from the Pacific cyclones—it nonetheless presents a similar local climate response as discussed in previous studies using both simplified and realistic LGM boundary conditions: that is, an anticyclone over the Laurentide ice sheet and a surface warming over the northwestern parts of the continent (e.g., Manabe and Broccoli 1985; Otto-Bliesner et al. 2006; Abe-Ouchi et al. 2007; Löfverström et al. 2014; Roe 1999; Roe and Lindzen 2001b). This is intriguing as it is known from geological data that Alaska remained largely ice free over the last glacial cycle (Clague 1989), whereas the rest of the North American continent poleward of about 40°N experienced a massive glaciation. Pollen-based reconstructions suggest that the average temperatures in Alaska were comparable or perhaps even slightly higher at the LGM than in the present climate (Bartlein et al. 2011). The simplified and zonally symmetric boundary conditions used here suggest that the warm glacial Alaska may have been the result of first-order flow–topography interactions and the ice sheet’s diabatic cooling of the lower atmosphere (Roe 1999; Roe and Lindzen 2001b).
b. Continent with Rockies
The ice sheet’s spatial evolution is strongly influenced by the presence of the Rockies; ice readily forms over the northern and northeastern parts of the continent but a large region in the lee of the Rockies remains ice free over the entire simulation. Consequently, the ice volume in the equilibrium state is only about half as large as on the flat continent. Interesting is that the equilibrium ice sheet resembles the MIS4 glacial state that is also characterized by an east-heavy configuration and largely ice-free conditions in the lee of the Rockies (Fig. 1a). This general disposition is believed to have been present over the larger part of the last glacial cycle (perhaps for as much as 75–85 kyr), and the continent-wide LGM Laurentide ice sheet is thought to have existed for a comparatively short period of time (5–15 kyr) before the glacial collapse (Kleman et al. 2013).
The Rockies yield a warm and dry summer climate in the continental interior (cf. left columns in Fig. 5 and Fig. 6) and the east-heavy ice sheet is found to enhance these climate characteristics and thus prevent its own westward expansion (Fig. 6). The ice sheet’s diabatic and mechanical forcing of stationary waves induce a westward-tilting anticyclone (baroclinic structure) that extends through the troposphere. This yields a reduction of low- and medium-height clouds but an increase of high-level clouds, which in turn contributes to a net warming of the surface climate in the interior of the continent. A sensitivity simulation with a more realistic MIS4 boundary condition shows that the warming signal induced by the ice sheet is a robust feature in this model, despite the idealized zonally symmetric boundary conditions. However, the magnitude of the warming signal is muted in the realistic MIS4 simulation, suggesting that there might be feedback loops with importance for the real world that are omitted in our idealized experiment.
We can only speculate on how different factors may have contributed to building the LGM ice sheet from the MIS4 state; however, the problem is in essence similar to that of the glacial inception. As noted earlier, a positive mass balance (low summer ablation and an abundant precipitation) is the key to building an ice sheet. Our results stress the importance of radiative feedbacks in the continental interior; hence, a lower summer insolation and possibly also lower concentrations of greenhouse gases may have been essential for breaking the MIS4 configuration and building the continent-wide LGM ice sheet. Ice core data (Petit et al. 1999; Spahni et al. 2005) have revealed that the atmospheric greenhouse gas concentrations varied considerably over the last glacial cycle, with CO2 concentrations bracketing values from approximately 280 ppmv in the interglacial to 185 ppmv at the LGM. The last significant drop of CO2 prior to the LGM (from ~210 to 185 ppmv) occurred between ~29 and 27 kyr BP, which coincides with a period of decreasing Northern Hemisphere summer insolation (from about 500 to 465 W m2 between 30 and 20 kyr BP; see Fig. 1 in Löfverström et al. 2014). It is in this time period that large and persistent ice sheets developed in the Cordillera and the Laurentide expanded westward over the prairies and interior plains (Clague and James 2002; Clague et al. 2005; Kleman et al. 2013).
Furthermore, both the paleodata and the modeling community have shown that the LGM sea surface climate was quite different from the present, especially in the North Atlantic region (e.g., Waelbroeck et al. 2009; Kucera et al. 2005a,b; Braconnot et al. 2007, 2012), and there are also indications of significant differences in the tropical ocean (Arbuszewski et al. 2013). The interaction between the ice sheet and atmosphere/ocean circulation is complex, where relatively small changes in the North American ice sheet have been shown to induce large changes in the atmosphere and ocean circulation (Jackson 2000; Zhang et al. 2014) that may feed back onto the ice-sheet development (Yin and Battisti 2001; Huybers and Molnar 2007). However, the same is not necessarily true for the Eurasian ice sheet (EIS). Although it was a large topographic obstacle, it is generally considered to have been located too far to the north to significantly influence the large-scale circulation other than regionally in Eurasia (Löfverström et al. 2014).
Our main results and conclusions are summarized as follows:
We find that the ice-sheet development on a flat representation of the North American continent is fairly zonally symmetric. In this simulation the incipient/early ice sheet develops a double-dome structure due to an anticyclonic circulation that forces zonal asymmetries in the precipitation field. In the latter half of the simulation, the ice sheet’s mass center is shifted toward the central parts of the continent.
The anticyclone over the ice sheet yields relatively warmer temperatures over the northwestern continent and cooler conditions in the southeast. This temperature configuration may possibly help explain why Alaska remained ice free over the last glacial cycle, whereas ice expanded as far south as 40°N in the east, already at an early stage of the last glacial cycle.
The ice-sheet evolution in the Rockies experiment is structurally different from the flat continent case. The presence of the mountain range yields more arid conditions in the continental interior. This shifts the ice sheet’s center of mass to the northeast, and the interior of the continent remains ice free throughout the entire simulation. This configuration resembles the geologically determined outline of the (pre-LGM) MIS4 glacial state.
We find that the ice sheet forces an anticyclone that maintains the asymmetric southwestern ice margin by inducing significantly warmer summer surface temperatures in the central parts of the continent. The warmer conditions are also amplified by the arid conditions in the lee of the Rockies. This may possibly serve as an explanation for the documented late ice invasion in the western Laurentide area (prairies and interior plains). A similar response is also obtained in sensitivity simulations using more realistic MIS4 boundary conditions.
We thank Rodrigo Caballero, Johan Nilsson, Francesco Pausata, David Battisti, Masa Kageyama, and Kerim Nisancioglu for the stimulating and useful discussions. We also thank three anonymous reviewers for their comments, which greatly helped to improve this manuscript. This work was financially supported by the Swedish Research Council, the Bolin Centre for Climate Research, and LOEWE of Hesse’s Ministry of Higher Education. The simulations were conducted on resources provided by the LOEWE Frankfurt Center for Scientific Computing and the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputing Centre (NSC) in Linköping.
Current affiliation: Climate Change Research Section, Climate and Global Dynamics Division, NCAR, Boulder, Colorado.
Current affiliation: Biodiversity and Climate Research Centre (BiK-F), and Senckenberg Gesellschaft für Naturforschung, Frankfurt, Germany.
Current affiliation: Department of Physical Geography and Quaternary Geology, and Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden.