Assessing the impact of Laurentide Ice Sheet topography on glacial climate

Abstract. Simulations of past climates require altered boundary conditions to account for known shifts in the Earth system. For the Last Glacial Maximum (LGM) and subsequent deglaciation, the existence of large Northern Hemisphere ice sheets caused profound changes in surface topography and albedo. While ice-sheet extent is fairly well known, numerous conflicting reconstructions of ice-sheet topography suggest that precision in this boundary condition is lacking. Here we use a high-resolution and oxygen-isotope-enabled fully coupled global circulation model (GCM) (GISS ModelE2-R), along with two different reconstructions of the Laurentide Ice Sheet (LIS) that provide maximum and minimum estimates of LIS elevation, to assess the range of climate variability in response to uncertainty in this boundary condition. We present this comparison at two equilibrium time slices: the LGM, when differences in ice-sheet topography are maximized, and 14 ka, when differences in maximum ice-sheet height are smaller but still exist. Overall, we find significant differences in the climate response to LIS topography, with the larger LIS resulting in enhanced Atlantic Meridional Overturning Circulation and warmer surface air temperatures, particularly over northeastern Asia and the North Pacific. These up- and downstream effects are associated with differences in the development of planetary waves in the upper atmosphere, with the larger LIS resulting in a weaker trough over northeastern Asia that leads to the warmer temperatures and decreased albedo from snow and sea-ice cover. Differences between the 14 ka simulations are similar in spatial extent but smaller in magnitude, suggesting that climate is responding primarily to the larger difference in maximum LIS elevation in the LGM simulations. These results suggest that such uncertainty in ice-sheet boundary conditions alone may significantly impact the results of paleoclimate simulations and their ability to successfully simulate past climates, with implications for estimating climate sensitivity to greenhouse gas forcing utilizing past climate states.


Abstract.
Simulations of past climates require altered boundary conditions to account for known shifts in the Earth system. For the Last Glacial Maximum (LGM) and subsequent deglaciation, the existence of large Northern Hemisphere ice sheets caused profound changes in surface topography and albedo. While ice-sheet extent is fairly well known, numerous conflicting reconstructions of ice-sheet topography suggest that precision in this boundary condition is lacking. Here we use a high-resolution and oxygen-isotopeenabled fully coupled global circulation model (GCM) (GISS ModelE2-R), along with two different reconstructions of the Laurentide Ice Sheet (LIS) that provide maximum and minimum estimates of LIS elevation, to assess the range of climate variability in response to uncertainty in this boundary condition. We present this comparison at two equilibrium time slices: the LGM, when differences in ice-sheet topography are maximized, and 14 ka, when differences in maximum ice-sheet height are smaller but still exist. Overall, we find significant differences in the climate response to LIS topography, with the larger LIS resulting in enhanced Atlantic Meridional Overturning Circulation and warmer surface air temperatures, particularly over northeastern Asia and the North Pacific. These up-and downstream effects are associated with differences in the development of planetary waves in the upper atmosphere, with the larger LIS resulting in a weaker trough over northeastern Asia that leads to the warmer temperatures and decreased albedo from snow and sea-ice cover. Differences between the 14 ka simulations are similar in spatial extent but smaller in magnitude, suggesting that climate is responding primarily to the larger difference in maximum LIS elevation in the LGM simulations. These results suggest that such uncertainty in ice-sheet boundary conditions alone may significantly impact the results of paleoclimate simulations and their ability to successfully simulate past climates, with implications for estimating climate sensitivity to greenhouse gas forcing utilizing past climate states.
However, the simulated glacial state in any one model is sensitive to the boundary conditions used as a starting point for the simulation (Manabe and Broccoli, 1985;Broccoli and Manabe, 1987;Hyde and Peltier, 1993;Justino et al., 2006;Abe-Ouchi et al., 2007;Liu et al., 2009;Pausata et al., 2011;Hofer et al., 2012). Thus, uncertainty in a particular set of glacial boundary conditions may overshadow a GCM's simulated change in climate and its ability to constrain climate sensitivity using the LGM as a starting point. Along this line, Hargreaves et al. (2012) suggested that disregarding highlatitude response at the LGM avoids bringing glacial boundary condition ambiguity into climate sensitivity work, although this omission then ignores many climate feedbacks. Some of the important boundary conditions are well documented, such as the orbital parameters that drive the magnitude and seasonality of insolation (Berger and Loutre, 1991), as well as the concentration of atmospheric greenhouse gases (Monnin et al., 2004;Jouzel et al., 2007;Lüthi et al., 2008;Lemieux-Dudon et al., 2010). Other boundary conditions of past glacial climates have proven to be more elusive. The aerosol forcing (particularly dust and black carbon) at the LGM and through the deglaciation is known only at individual locations, with limited spatial resolution (Petit et al., 1981;Thompson et al., 1995;Mahowald et al., 1999Mahowald et al., , 2006aReader et al., 1999;Power et al., 2008;Ganopolski et al., 2010). While large uncertainties still remain on the impact of aerosols on radiative forcing (Penner et al., 2004;Lohmann and Feichter, 2005;Forster et al., 2007;Chylek and Lohmann, 2008), the LGM dust loading may issue a forcing on the order of ∼ −1 W m −2 that is comparable to changes in greenhouse gases alone (∼ −2.85 W m −2 ), at least in the tropics, where ice-sheet albedo is not a factor (Harrison et al., 2001;Claquin et al., 2003;Crucifix, 2006). Likewise, reconstructions of LGM and deglacial vegetation and the related impact on surface albedo are limited to low-resolution global approximations and higher resolution estimates that are only regional in scope (Jackson et al., 2000;Prentice and Jolly, 2000;Harrison et al., 2001;Ray and Adams, 2001;Bigelow et al., 2003;Williams, 2003;Pickett et al., 2004;Williams et al., 2011). The impact of such uncertainties in boundary conditions may even be amplified through possible internal feedbacks of climate models with the inclusion of new dynamic components, such as vegetation and dust (Ridgwell and Watson, 2002;Mahowald et al., 2006a, b;O'ishi and Abe-Ouchi, 2013).
Perhaps the greatest source of known uncertainty in glacial boundary conditions relates to ice-sheet thickness. Although the geographical extents of LGM and deglacial ice-sheets are fairly well mapped Hughes, 1981, 2002;Dyke and Prest, 1987;Anderson et al., 2002;Clark and Mix, 2002;Bennike and Björck, 2002;Dyke, 2004;Gyllencreutz et al., 2007), direct observations of ice thickness and topography are limited and usually absent for the highest/thickest portions of the ice sheets (e.g., Hughes, 1981, 2002;Dyke et al., 2002;Clark, 1992;Clark et al., 2009;Goehring et al., 2008;Carlson and Clark, 2012). Therefore, ice-sheet height must be simulated through geophysical or glaciological modeling approaches. The reconstructions of ICE-4G, ICE-5G, and ICE-6G (Peltier, 1994(Peltier, , 2004Toscano et al., 2011) have proven to be useful ice-sheet boundary conditions to employ in GCMs due to their geophysically based solutions, global scope, and product availability. However, the differences between these reconstructions (Peltier, 2004;Toscano et al., 2011), as well as other ice-sheet-specific reconstructions Licciardi et al., 1998;Tarasov and Peltier, 2004;Tamisiea et al., 2007;Argus and Peltier, 2010;Lambeck et al., 2010) suggest that there is a large range in existing ice-sheet boundary conditions, particularly for the Laurentide Ice Sheet (LIS).
Here we test the impact of LIS geometry on LGM and deglacial climate using two alternate reconstructions of LIS topography that provide upper and lower bounds for this boundary condition. As an upper bound, we use ICE-5G (VM2) (Peltier, 2004;hereafter 5G), which has the highest LIS topography of any LIS reconstruction, with a dominant ice dome over the western Keewatin sector (Fig. 1). As a lower bound, we employ the dynamics-driven flowline reconstruction of Licciardi et al. (1998, hereafter L) (Fig. 1). The maximum height of this reconstruction is comparable to ICE-4G (Peltier, 1994;Clark et al., 1996), but it is ∼ 20 % lower than that of ICE-5G, particularly over the western LIS (Fig. 1). The use of these reconstructions is  Licciardi et al. (1998) reconstructions used in 21ka-L and 14ka-L. The middle column presents the ICE-5G reconstructions (Peltier, 2004) used in 21ka-5G and 14ka-5G. For comparison, the single plot in the right column shows the LIS topography used in the PMIP3 boundary conditions (http://pmip3.lsce.ipsl.fr/) a particularly relevant assessment of the range of variability in the Paleoclimate Modelling Intercomparison Project 3 (PMIP3) topographic boundary conditions, as the maximum ice-sheet heights of the three ice-sheet reconstructions used to create the averaged topography of the LIS in the PMIP3 (http://pmip3.lsce.ipsl.fr/) lie in between the upper and lower bounds used in our study. We compare the effects of these LIS boundary conditions at two time periods, 21 and 14 ka, to assess the impact of varying LIS topography under different climate forcing conditions. The 21 ka simulations provide a test of topographic variability under full LGM boundary conditions with different LIS maximum elevations. At 14 ka, both LIS topographies are reduced and differences between maximum LIS elevation are smaller (∼ 15 %). The 14 ka simulations provide a test of ice-sheet geometry differences under a different forcing scenario, with changes in precession, obliquity, and increased atmospheric greenhouse gas concentrations.

Model
The NASA Goddard Institute for Space Studies (GISS) Mod-elE2, coupled to the Russell (R; Russell et al., 2000) ocean is a fully coupled atmosphere-ocean global climate model (Schmidt et al., 2014). We use the same version of ModelE2-R as is being used for the Coupled Model Intercomparison Project Phase 5 (CMIP5) simulations (NINT physics version), with an atmosphere resolution of 2 • latitude by 2.5 • longitude with 40 vertical layers up to 0.1 mb, and an ocean resolution of 1 • latitude by 1.25 • longitude with 32 depth layers. The 21 ka simulations presented here are the GISS ModelE2-R contribution to LGM experiment of CMIP5 (ensembles r1i1p150 and r1i1p151; see Sect. 2.3). ModelE2-R includes passive water isotopologue tracers (i.e., H 18 2 O, HDO, hereafter referred to as water isotopes), removing one degree of uncertainty when comparing its diagnostics to paleoclimate reconstructions based on water isotope variability (e.g., ice cores, speleothems, ocean foraminifera) to test model skill (e.g., LeGrande et al., 2006;Carlson et al., 2008a, b;LeGrande and Schmidt, 2009).

Boundary conditions
Each simulation uses appropriate insolation for the time period owing to changes in orbital parameters (Berger and Loutre, 1991), along with appropriate glacial greenhouse gas concentrations (Joos and Spahni, 2008) (Table 1). Coastline and basin geometry were adjusted to reflect an ∼ 120 m lowering in sea level at 21 ka and a ∼ 86 m lowering at 14 ka, in accordance with the volume of water contained in icesheet boundary conditions and glacial isostatic adjustment (Peltier, 2004). The volume of water held in the two alternate LIS reconstructions is different, but in order to avoid small coastal idiosyncrasies between the two simulations for each time slice, we kept the ocean bathymetry (and the land-sea mask) the same for each pairing of LIS boundary conditions. Because we aim to study differences in ocean circulation due to changes in the LIS topography alone, bathymetry was held consistent between the simulations at each time slice so as  Table 1. Boundary conditions for experiments presented in this paper. The abbreviations for the orbital parameters refer to eccentricity (E), obliquity in degree (O), and angular precession (longitude of perihelion, or omega) in degree (P ) (Berger and Loutre, 1991 not to introduce basin differences along ocean boundaries that might impact ocean dynamics. In addition, ModelE2-R ocean uses straits -essentially pipelines for ocean tracers and mass -from present day to accommodate sub-grid-scale passages. All straits were eliminated (because of sea level lowering) except for Gibraltar and Bab Al Mandab (Red Sea), which were reduced, and no new straits were created. ModelE2 does not have a dynamic vegetation component, so terrestrial vegetation cover was assigned according to mapped conditions at the LGM (25 000-15 000 yr BP; Ray and Adams, 2001) for regions not covered by ice sheets. A separate reconstruction for 14 ka does not exist; therefore we continued to use this LGM reconstruction, modified slightly to account for changes in ice-sheet extent by expanding nearby vegetation types into these areas no longer covered by ice. Continental river routing was assigned in accordance with ice-sheet configuration and its impact on drainage basins (e.g., Licciardi et al., 1998Licciardi et al., , 1999Peltier, 2004).

Ice-sheet geometry sensitivity
To test the influence of ice-sheet topography on glacial climate, two separate simulations were conducted at each time slice using two alternate reconstructions of LIS topography: one using the geophysically based ICE-5G (5G; Peltier, 2004; CMIP5 ensemble r1i1p150 for the LGM), and the other using the same extent of ICE-5G but replacing the LIS sector topography with that of an alternate reconstruction based on a flow-line model that simulates glacier dynamics over deformable and rigid beds with varying till viscosities (L; Licciardi et al., 1998; CMIP5 ensemble r1i1p151 for the LGM). We use the maximum reconstructions from Licciardi et al. (1998) that incorporate higher specified effective till viscosities. These Licciardi et al. (1998) maximum reconstructions were concluded to provide a better representation of the LIS than the "minimum" reconstructions due to better agreement with other ice-sheet reconstructions and relative sea level records Licciardi et al., 1998). In both simulations, all other ice-sheet topographies and extents (Greenland, Scandinavian, etc.) are assigned using ICE-5G. We focus on the LIS because it was the largest of the Quaternary ice sheets with presumably the greatest effect on climate Clark et al., 1999;Clark and Mix, 2002;Carlson and Clark, 2012).
The two LIS reconstructions have significantly different topographies ( Fig. 1). At 21 ka, the 5G reconstruction provides a much taller LIS, particularly the Keewatin Dome over western Canada, with elevations across this region at or above 4000 m. The 21 ka-L reconstruction is lower in maximum height and moves more of the ice mass to the east with three distinct ice domes centered over Keewatin, southern Ontario, and central Quebec. Comparatively, the 5G-LIS has a lower topography over eastern Ontario and Quebec, but the large Keewatin Dome dominates the western topography. Additionally, the 5G reconstruction has more abrupt changes in topography, whereas the L reconstruction has smoother transitions from high to low elevations, as would be expected from a flow-line model of ice deformation over hard and soft beds and similar to inferences from the geologic record Prest, 1987, Jenson et al., 1996).
At 14 ka, both reconstructions retain the same general geometry as 21 ka, but absolute topography is diminished (Fig. 1). ICE-5G still has a dominant Keewatin dome, while the L reconstruction continues to have three prominent domes. The L reconstruction exhibits only a small diminution in ice topography over the Ontario and Quebec domes, such that a west-east dipole in topography difference exists between the two reconstructions ( Fig. 1).
Our climate model simulations were started before the PMIP3 protocol for LGM conditions had been determined, but our selection of LIS boundary conditions provide limits of topographic uncertainty that go into the PMIP3 LIS elevations. The LIS in the PMIP3 boundary conditions is similar in geometry to ICE-5G (with higher Keewatin dome over western Canada and lower topography over eastern Ontario and Quebec), but the maximum topography of the LIS in PMIP3 is more similar to the LIS reconstruction used in 21ka-L ( Fig. 1).

Equilibrium simulations -21, 14, and 0 ka
A control simulation (0 ka) with water isotopes was initiated from a 700 yr long spinup from Levitus and Boyer (1994a, b) and  and run an additional 1000 yr. The 21 and 14 ka simulations were also initiated from the same 700 yr long spinup from Levitus and Boyer (1994a, b) and , but these experiments have altered mean ocean δ 18 O, δD, and salinity (see section below), the last of which requires a much longer integration for the ocean density structure to adjust. Surface air temperature (SAT) drift is significant for the first 500 yr in each of the deglacial/glacial simulations -some of which is related to a temporary enhancement of AMOC due to the method with which enhanced salinity was applied. Fractional increases in salinity were uniformly added to each grid box, where in reality additional salinity would have likely been preferentially buried in the deep ocean. Such colder and saltier conditions are more conducive to the enhanced overturning seen in the initial spinup, which has been shown to require an extended initialization (Zhang et al., 2013). After 1500 yr of integration, SAT drift is less than 0.03 • C per century, surface ocean temperature drift less than 0.02 • C per century, and deep ocean temperature drift less than 0.03 • C per century and in each of the simulations, but overall drift is never absent. Note that these simulations run 500 yr longer than those present in the CMIP5 archive to insure that the water isotopologues have come into equilibrium (water isotopes are not a part of the CMIP5 archive). Initially, each simulation had an enhanced ocean circulation, which diminished with approach toward equilibrium. In addition, cold surface temperatures drove sea ice growth to the bottom of the ocean grid in some shallow coastal polar regions, requiring an expansion of the land mask as such "ice shelves" are defined in other regions of the boundary conditions. Here we present the final 100 yr of these simulations.

Water isotopes
ModelE2-R includes water isotopologue tracers that allow for the direct comparison of δ 18 O measurements taken from ice cores, speleothems, ocean sediments, and other paleoclimate records (Carlson et al., 2008a, b;LeGrande and Schmidt, 2009). We have selected terrestrial data that provide a proxy of the δ 18 O in precipitation (hereafter δ 18 O a ) throughout the LGM and deglacial climate (see Table S1 in the Supplement for list of δ 18 O proxy records used in this analysis). To estimate LGM and deglacial δ 18 O in the ocean (hereafter δ 18 O o ), data from ocean sediment cores with continuous coverage from the LGM to the preindustrial era were selected (Table S1 in the Supplement). We used the average of all data within 1000 yr of each time slice (i.e., 20-22 ka for 21 ka, 13-15 ka for 14 ka, 2-0 ka for control) to estimate the δ 18 O a and δ 18 O o anomaly from each record for direct comparison with the model.

Results
Our 21 and 14 ka simulations are broadly consistent with previous GCM simulations of LGM and deglacial climate (e.g., CLIMAP, 1981;Kageyama et al., 1999;Justino et al., 2004;Otto-Bliesner et al., 2006Abe-Ouchi et al., 2007;Braconnot et al., 2007a, b;Timm and Timmerman, 2007;Liu et al., 2009). However, this paper is primarily concerned with the differences that arise due to variation in the LIS topography alone. Since most of the climate differences directly over the ice sheet are due to orographic differences alone, we focus on the LIS impacts in other regions in the Northern Hemisphere where GCM anomalies are maximized. For a detailed description of modeled 21 and 14 ka anomalies relative to the GISS ModelE2-R preindustrial control simulation, please see the Supplement.

Surface air temperature
SAT anomalies are significantly different between the two simulations at 21 ka (Fig. 2a). Global area-weighted mean SAT anomalies (21 ka minus control) for the 21ka-5G and 21ka-L simulations are −5.1 ± 0.1 • C and −5.4 ± 0.1 • C, respectively (unless otherwise stated, uncertainty is expressed as standard deviation across decadal variability about the 100 yr mean). This difference in SAT is significant, developing solely in response to changes in LIS boundary conditions. Since differences in SAT between our two scenarios are exacerbated directly over the LIS due to topographic differences (i.e., vertical lapse rate of temperature and an ∼ 1000 m maximum LIS height difference), the global mean SAT anomalies with the LIS area removed are −4.7 ± 0.1 • C and −5.0 ± 0.1 • C, thus confirming that differences between the simulations are not due to differences directly over the ice sheet alone, but rather due to differences in other regions that significantly impact the global mean (Fig. 2a).
The most prominent downstream differences in SAT between the 21 ka simulations occur over northeastern Asia, Beringia, and the North Pacific, where 21ka-L is 6-9 • C colder than 21ka-5G (Fig. 2a). Proxy records from this region are limited, but one recent modern analog approach based on pollen spectra from Lake El'gygytgyn in northeastern Siberia suggests that the mean temperature of the warmest month (July) at 20 ka is ∼ 6 • C colder than present (Melles et al., 2012), whereas the model results in a cooling of July temperatures by 16.5 ± 1.1 • C (21ka-L) and 12.8 ± 0.7 • C (21ka-5G). Model results appear to be much too cold in this region, but the lake record may not extend far enough back to capture full LGM cooling, and lake ice permanence at the LGM may have strongly limited pollen transport into the lake during this time, biasing the record towards warm years (Melles et al., 2007(Melles et al., , 2012. One SAT reconstruction from fossil lipids in paleosols of eastern Asia suggests LGM summer cooling of 8-10 • C (Peterse et al., 2011), whereas our 21 ka simulations suggest summer cooling of 9.0 ± 0.7 • C (21ka-L) and 8.1 ± 0.4 • C (21ka-5G). Reconstructions from Zagoskin and Burial lakes in northwestern Alaska suggest July LGM cooling of 4 ± 1.6 • C and 5 ± 1.6 • C, respectively, which is consistent with general LGM cooling of the region (Viau et al., 2008;Kurek et al., 2009). Our 21 ka simulations suggest July LGM cooling of 8.8 ± 0.9 • C (21ka-L) and 4.6 ± 0.4 • C (21ka-5G) for Zagoskin Lake, and 12.2 ± 1.0 • C (21ka-L) and 8.5 ± 0.7 • C (21ka-5G) for Burial Lake. Other pollenbased reconstructions may actually provide some evidence for slightly warmer in Alaska during the LGM , although there may be no-analog issues in this region due to the limited number of individual records available in this synthesis. Unfortunately, no sites from northeastern Asia exist in the LGM climate reconstructions of Bartlein et al. (2011), limiting the model-data comparison of SAT, where our 21 ka simulations provide the greatest differences.

Precipitation minus evaporation
The location of the Intertropical Convergence Zone (ITCZ) is consistent between 21ka-L and 21ka-5G, but 21ka-L simulates greater precipitation north of the Equator in the tropical Pacific (Fig. 2b). The 21ka-5G simulation shows greater precipitation over the northern mid-latitude Pacific and across eastern America, south of the LIS, indicative of a southward displacement of the primary storm tracks (Justino et al., 2004).
At 14 ka, differences in P − E between simulations are small, but highlight a slight enhancement of Pacific precipitation south of the Equator in 14ka-5G (Fig. 3b). Small differences in precipitation also occur over the North Atlantic, with drier conditions in 14ka-L along the east coast of North America and the Atlantic coast of Europe and slightly wetter conditions in the interior of the North Atlantic.

δ 18 O of the atmosphere differences
Differences in δ 18 O a largely reflect the differences in SAT, with the greatest relative depletion of 21ka-L occurring in northeastern Asia and the North Pacific, where SATs show the greatest cooling relative to 21ka-5G (Fig. 2c). In the 14 ka simulations, there is a decoupling between the differences in SAT and the differences in δ 18 O a , with relatively enriched δ 18 O a in 21ka-L across northeastern Africa and into southcentral Asia (Fig. 3c). Additionally, the 21ka-L shows relatively enriched values across much of the Arctic, but depleted values immediately downwind of the LIS and over southern Greenland.

Atmospheric circulation
The Northern Hemisphere subpolar jet in 21ka-5G is stronger and more zonal across the mid-latitude Pacific and to the south of the LIS (Fig. 2e). This difference in jet location is associated with the increased precipitation in the 21ka-5G across eastern North America as the primary storm tracks are confined further to the south than in 21ka-L (see Sect. 3.1.2; Fig. 2b). The subpolar jet in 21ka-L is more wavelike and lies to the north of the 21ka-5G jet across most of Asia and the Pacific. Differences in surface wind speed reflect this jet displacement, particularly over the North Pacific, where the more northern 21ka-L jet results in greater surface winds. Across much of the North Atlantic, 21ka-5G surface winds are stronger than 21ka-L, where both simulations are already enhanced relative to 0 ka (Fig. S1d in the Supplement).
Differences in 500 mb height between the two 21 ka simulations express a shift in the location and depth of the dominant stationary wave patterns, particularly over Siberia and Beringia, where 21ka-L heights are substantially lower than those in the 21ka-5G simulation (Fig. 2f). These lower heights in 21ka-L reflect a deepening of the trough over eastern Asia and a weakening of the ridge over Beringia relative to 21ka-5G. This general reduction of 500 mb heights across this region provide a greater influence of Arctic air masses that help to drive the colder SAT (Sect. 3.1.1).
The higher LIS in 14ka-5G still induces an increase in jet speed and southward displacement of the polar and subtropical jets relative to 14ka-L (Fig. 3e). However the stronger subpolar jet in 14ka-5G now extends further to the north over northern Europe. The general differences in 500 mb heights  Fig. 2 but with the atmospheric differences between 14 ka simulations (14ka-L minus 14ka-5G). and the implied stationary wave pattern in the 21 ka simulations still exists between the 14 ka simulations, with a continued lowering of 14ka-L 500 mb heights in northeastern Asia and Beringia that is associated with the colder SAT in this region relative to 14ka-5G (Fig. 3f).

Albedo
The 21ka-L simulation presents a region of higher albedo in northeastern Asia, where there is additional snow cover over Siberia and Beringia, along with a relative expansion of sea ice over the Sea of Okhotsk (Fig. 2g). The 21ka-5G has a higher albedo in a region along the southern margin of the LIS, likely due to the southward tracking of the subpolar jet and enhanced snowfall. The disparity in the area extent of these albedo changes results in globally averaged ground albedo of 18.5 % for 21ka-L compared with 18.1 % for 21ka-5G. Planetary albedo differences mimic the ground albedo differences (Fig. 2h), although cloud cover over the region in northeastern Asia is enhanced by 8-10 % in the 21ka-L simulation (not shown).
At 14 ka, surface albedo is largely the same between the two simulations, except for the expansion of sea ice over the Nordic Seas in 14ka-5G, which drives an increase in albedo relative to 14ka-L (Fig. 3g). Globally averaged ground albedo is 16.3 % for 14ka-L and 16.4 % for 14ka-5G. Again, planetary albedo differences mimic ground albedo differences, suggesting the prominence of surface albedo changes in overall albedo (Fig. 3h), with cloud cover differences of less than 2 %.
Regional differences in SST between the 21ka-L and 21ka-5G simulations correlate with SAT differences (Figs. 2a and 4a), particularly in the North Pacific, Sea of Okhotsk, and Bering Sea, where 21ka-L simulates colder SST by up to 3.5 • C. The large differences between the simulations in the North Pacific indicate that this is a sensitive region for testing the model results. SST proxy records are sparse from the North Pacific, but limited data seem to indicate colder conditions. Alkenone and transfer function SST reconstructions from off the coast of Oregon document LGM cooling of 4-7 • C (Doose et al., 1997;Ortiz et al., 1997;Mangelsdorf et al., 2000;Herbert et al., 2001;Rosell-Melé et al., 2004), which are more consistent with the 21ka-L result of 4.0 ± 0.3 • C cooling than with the 21ka-5G cooling of 3.1 ± 0.2 • C. Alkenone records from the Sea of Okhotsk, however, suggest LGM SST cooling of only 0-1 • C (Harada et al., 2006;Seki et al., 2004Seki et al., , 2007. The 21ka-5G simulation has LGM cooling of 4.6 ± 0.2 • C in this region, and the 21ka-L simulation shows cooling of 7.4 ± 0.2 • C, suggesting that neither simulation adequately captures SST in this somewhat isolated basin. At 14 ka, globally averaged SST anomalies are equivalent between the two ice-sheet topographies (−1.5 ± 0.1 • C in 21ka-L and −1.4 ± 0.1 • C in 21ka-5G), but total global mean ocean anomalies are still colder in the 14ka-L simulation with an anomaly of −0.83 ± 0.1 • C (21ka-L) relative to −0.71 ± 0.1 • C (21ka-5G). In the deep ocean (bottom 2000 m), however, temperature is nearly consistent with 0 ka with anomalies of −0.09 ± 0.01 • C (14ka-L) and +0.05 ± 0.01 • C (14ka-5G). Differences between the 14 ka simulations are most pronounced across the North Pacific (Fig. 5a), where 14ka-L SSTs are colder, similar to SAT differences over these regions (Fig. 4a).

Sea surface salinity differences
Surface waters are generally fresher in 21ka-L, particularly in the Arctic Ocean and the Gulf of Mexico, reflecting a greater contribution of LIS melt water through the MacKenzie and Mississippi rivers (Fig. 4b). The 21ka-5G simulation presents a greater contribution of Scandinavian Ice Sheet melt water through the Ob and Yenesei rivers, leading to a localized freshening of sea surface salinity (SSS) along the coast, but in general, Arctic surface waters are more saline in 21ka-5G. Northern tropical Pacific surface waters are slightly fresher in 21ka-L in association with enhanced P − E over this region (Fig. 2b).
SSS differences at 14 ka are largely similar to those in the 21 ka simulations, with fresher waters in 14ka-L in the Arctic Ocean, the Gulf of Mexico, and along the Atlantic coast of North America (Fig. 5b). The 14ka-L simulation continues to present localized freshening due to greater melt water contributions of the LIS to the MacKenzie and Mississippi rivers. SSS differences are negligible over the Pacific, reflecting the minimal changes in P − E (Fig. 3b).

δ 18 O of the surface ocean differences
The differences in δ 18 O o at 21 ka are largely similar to the differences in SSS, with depleted values in 21ka-L relative to 21ka-5G across the Arctic Ocean and the Gulf of Mexico (Fig. 4c), reflecting the enhanced melt water contributions in 21ka-L to the MacKenzie and Mississippi rivers (Fig. 4b). In addition, 21ka-L waters are depleted along the north shore of the Scandinavian Ice Sheet, consistent with reduced SSS and enhanced melt water contribution in 21ka-5G along this margin. The Sea of Japan shows slightly depleted δ 18 O o values in 21ka-L but no change in SSS, suggesting a shift in the isotopic composition of runoff arriving to this basin but a consistent volume of this runoff between the simulations.
In the 14 ka simulations, however, this relationship between SSS and δ 18 O o becomes decoupled, with fresher arctic conditions in 14ka-L associated with more enriched values of δ 18 O o relative to 14ka-5G (Fig. 5c). However, the 14ka-L simulation continues to show depleted values in relation to the greater contribution of fresh water from the MacKenzie and Mississippi drainage outlets. In addition, the Sea of Japan continues to reflect a difference in the isotopic composition of river runoff to the basin, but here the 14ka-5G simulation presents depleted δ 18 O o values. Again, the lack of a SSS difference in this basin suggests that freshwater runoff is consistent between the simulations (Fig. 5b).

Sea ice differences
Annually averaged sea ice differences are negligible across much of the ocean basins, except in the Sea of Okhotsk and Bering Sea, where the 21ka-L has a significant increase in sea ice over that of 21ka-5G (Fig. 4d). The expansion of sea ice across the Nordic Seas in the 21 ka simulations relative to the control simulation is consistent between 21ka-L and 21ka-5G, with 21ka-L having a slightly greater sea-ice fraction (no more than 5 %).
In the 14 ka simulation, global sea ice differences are small (Fig. 5d). Differences in sea ice are minimized in the Sea of Okhotsk compared to the 21 ka simulation. The only notable differences occur in the Nordic Seas, with enhanced sea ice in the 14ka-5G simulation.

Ocean circulation differences
All simulations result in a mean AMOC that is stronger than the preindustrial control simulation (see Supplement; 0 ka control AMOC strength: 28.2 ± 0.7 Sv). The 21ka-5G simulation exhibits a mean AMOC transport that is ∼ 10 % stronger than 21ka-L (33.2 ± 0.7 Sv and 30.8 ± 0.6 Sv, respectively). The difference in AMOC stream function reflects not only the increase in AMOC transport for 21ka-5G but also a deepening of the main overturning and a slight increase in Antarctic bottom water in flow (Fig. 4e). In the Pacific, the 21ka-5G simulation shows an enhanced contribution of Antarctic intermediate and bottom waters throughout most of the basin (Fig. 4f), with the 21ka-5G having the greatest Antarctic bottom water influence (i.e., deepwater stream function more negative than 21ka-L). In addition, the 21ka-5G shows enhanced production of North Pacific intermediate water relative to 21ka-L (or less of a reduction relative to 0 ka).
Kinematic and water-mass proxy records continue to refine reconstructions of AMOC at the LGM, with overturning strength anywhere from 30 to 40 % weaker than present (McManus et al., 2004) to about the same as today (Lynch-Stieglitz et al., 2007;Ritz et al., 2013), and even the possibility of more rapid overturning at the LGM (Gherardi et al., 2009;Lippold et al., 2012). Water mass tracers suggest the shoaling of glacial North Atlantic deepwater and a greater contribution of Antarctic bottom water (Curry and Oppo, 2005), but this increased stratification may also imply a more vigorous AMOC. Such uncertainty in circulation strength and depth is mirrored in a diverse array of modeled AMOC results, some of which present stronger AMOC during the LGM , similar to our results. In the Pacific, the expansion of enriched δ 13 C throughout the deep ocean suggests the increased influence of Antarctic-sourced waters in this basin (Matsumoto et al., 2002;Herguera et al., 2010). Other North Atlantic simulations have shown that enhanced AMOC is associated with the strengthening of the deep ocean circulation in the Pacific (Chikamoto et al., 2012), which is consistent with the enhanced negative stream function in our 21 ka results (Fig. S3f in the Supplement).
The 14 ka simulations also present strong AMOC relative to the preindustrial control (see Supplement). The 14ka-5G simulation continues to have stronger AMOC transport (32.7 ± 0.7 Sv for 14ka-5G; 30.5 ± 0.6 Sv for 14ka-L), but the depth of overturning is similar to the 14ka-L simulation (Fig. 5e). The 14ka-5G results in a greater contribution of Antarctic bottom water into the deep Atlantic. In the Pacific, 14ka-5G shows an increased contribution of Antarctic intermediate water, but the differences do not extend to include Antarctic bottom water (Fig. 5f).
Proxy records suggest that 14 ka AMOC strength was similar to that of the LGM (McManus et al., 2004;Lynch-Stieglitz et al., 2007). However, other records suggest a more intermediate rate of overturning in shallower waters between the relatively elevated values of the LGM and the reduced values in the Holocene (Gherardi et al., 2009). In the Pacific, overturning at 14 ka was in the middle of the transition from glacial to modern conditions (Herguera et al., 2010), which is similar to our results.

The effects of ice-sheet topography
Uncertainty in the height of the LIS has a measurable impact on the simulated climate in GISS ModelE2-R, particularly at 21 ka, where global mean SATs are significantly different between 21ka-L and 21ka-5G. This response to changes in LIS topography alone is due to a series of differences in the 21 ka climate, initiated by the differences in atmospheric circulation that arise from the enhanced topographic barrier in the LIS of 21ka-5G. In both 21 ka simulations, the polar jet is forced to the south of the LIS (Fig. S1e in the Supplement). However, the elevated topographic barrier in the 21ka-5G LIS drives the polar jet to be more zonal, whereas the 21ka-L jet circulation has a greater meridional component, reflecting a shift in the downstream stationary waves in the Northern Hemisphere relative to 21ka-5G, as seen in 500 mb height differences (Fig. 2f). The primary impact of these stationary wave differences is downstream colder temperatures across Siberia, Beringia, and the North Pacific, leading to elevated sea ice in the Sea of Okhotsk and Bering Sea, as well as enhanced snow cover across Siberia (Fig. S1 in the Supplement and Fig. 4). Both of these impacts lead to an increase in ground albedo, thus reducing the total incoming shortwave budget in 21ka-L ( Fig. 2g and h). This snow-cover-albedo positive feedback thus leads to further global cooling.
In the 14 ka simulations, the difference in the magnitude of maximum LIS height is smaller between the two reconstructions, but there is still downstream cooling over Siberia and Beringia in 14ka-L related to a similar difference in atmospheric circulation and stationary waves ( Fig. S2e and f in the Supplement). However, differences in snow cover and surface albedo are no longer present across this region (Fig. 2g), and the global mean temperatures between 14ka-L and 14ka-5G are equivalent. This comparison with the 21 ka simulation might suggest a minimum cutoff in LIS maximum elevation difference that induces strong enough changes in atmospheric circulation to cause significant differences in global mean surface temperature through a snow-albedo feedback. However, the 14 ka simulations are also forced with elevated boreal summer insolation and greenhouse gases, such that generally warmer globally averaged temperatures at 14 ka (relative to 21 ka) may preclude the expansion of snow and sea ice in this region, eliminating the possibility of inducing this snow-albedo cooling feedback, even with the differences in planetary wave strength induced by LIS topographic differences.
Significant differences in ocean circulation also arise due changes in the LIS topography. All of our simulations show increased AMOC transport relative to the control, with southward displacement of the North Atlantic deepwater formation (Figs. S3e and S5e in the Supplement). Glacial surface winds are enhanced over the North Atlantic (relative to 0 ka) in relation to the southward shift in the polar jet (Figs. S1d and S2d in the Supplement). Such windier conditions may provide the mechanism driving enhanced overturning in our simulations (Wunsch, 2003;Montoya and Levermann, 2008). In addition, the 21ka-L simulation shows a slightly elevated freshwater balance (P − E, Fig. 2b) over the eastern North Atlantic, consistent with lower SSS, which may also contribute to the weaker overturning circulation relative to 21ka-5G.
In each of the 21 and 14 ka simulation pairings, the higher LIS (21ka-5G and 14ka-5G) results in a stronger AMOC (Figs. 4e and 5e). In either case, there is no significant change in river runoff between the simulations, but both 21ka-5G and 14ka-5G have stronger wind speeds over the North Atlantic relative to their 21ka-L and 14ka-5G counterparts, leading to enhanced wind-driven overturning. These differences in AMOC strength are transferred to the deep Pacific with enhanced contribution of Antarctic bottom and intermediate waters associated with the stronger AMOC (Chikamoto et al., 2012). Despite these changes in circulation, the 21ka-L simulation has colder ocean temperatures (relative to 21ka-5G) in both the globally averaged ocean and the deep ocean, suggesting ocean temperatures are reflecting the overall colder climate and reduced heat content of the system, as would be expected from the snow-albedo feedback mechanism driving colder SAT.

Implications for ice-sheet stability
The extent and volume of the LIS was determined by a number of glaciological factors that impact ice-sheet stability. LIS surface mass balance was paced by boreal summer insolation (Hays et al., 1976;Imbrie and Imbrie, 1980). Such continuous cycling of boreal summer insolation limits the size of the LIS through latitudinal shifts in the equilibrium line altitude moving in step with insolation (Oerlemans, 1991;Ruddiman, 2001), with variation in latitudinal extent of the LIS likely leading to ice volume hysteresis . Ice-sheet size is also limited by an ice-thicknessbasal-melting negative feedback, whereby increases in ice thickness can lead to a decrease in the pressure melting point, decoupling the ice sheet from its bed and allowing for enhanced basal sliding and ice height drawdown (Clarke, 1987;Payne, 1995;Marshall and Clark, 2002). These thermal instabilities may be exacerbated by subglacial till rheology conditions that might also limit ice-sheet height (Clark, 1994;Clark and Pollard, 1998;Licciardi et al., 1998).
Our results may suggest an additional limit on LIS height. As shown in our LGM simulations, the larger LIS in 21ka-5G leads to a warmer global mean SAT due to atmospheric circulation changes and the snow-albedo feedback in northeastern Asia. That an increase in LIS maximum height leads to global surface warming might suggest a climatically driven negative feedback on LIS surface mass balance that limits ice-sheet height above a certain threshold. Satellite-based gravity field measurements of glacial isostatic adjustment have suggested the need for revisions to ICE-5G (Peltier and Drummond, 2008;Peltier, 2009) to be included in the upcoming ICE-6G that lower maximum LIS elevation ; otherwise see PMIP3 boundary conditions, http://pmip3.lsce.ipsl.fr/). These revisions may imply that LIS elevation in ICE-5G is above the threshold of this elevation-warming feedback. In order to test this feedback mechanism, more analysis should be conducted on the impact of this SAT difference and ice-sheet topography on LIS surface mass balance.

Implications for testing LGM climate reconstructions
We present a range of climate responses to LGM external and internal forcings that develop solely due to two different reconstructions of the LIS. This range of uncertainty in this boundary condition alone is enough to provide significant differences in LGM climate. Previous sensitivity experiments testing the impact of "ice" vs. "no ice" in model boundary conditions have shown large differences in LGM climate due to changes in atmospheric circulation and its attendant influence on ocean circulation (Justino et al., 2004;Otto-Bliesner et al., 2006;Abe Ouchi et al., 2007;Pausata et al., 2011;Hofer et al., 2012;Tharammal et al., 2012). Our 21 ka model results show that even a 20 % change in LIS elevation between two LIS reconstructions has a similar impact on atmospheric and oceanic circulation and their influence on the LGM climate state. This suggests that the range of uncertainty within the existing LIS reconstructions can lead to significantly different results in simulated regional and global climate. The ability of GCMs to capture LGM AMOC is often used as a test of model skill (e.g., Timm and Timmerman, 2007;Liu et al., 2009). Despite a range of proxy estimates for an AMOC target value (Yu et al., 1996;Adkins et al., 2002;McManus et al., 2004;Curry and Oppo, 2005;Lynch-Steiglitz et al., 2007), coupled climate models simulate different strengths and depths of the AMOC Weber et al., 2007). The strength of the PMIP approach toward simulating the LGM is in the assessment of climate variability (i.e., AMOC) across a variety of model physics and design with a common fixed set of boundary conditions. However, our results suggest that AMOC strength in a single model may vary by nearly 10 % in response to changes in LIS topography alone, meaning that any uncertainties in the LIS boundary condition may be translated into the uncertainty in this important test of model skill. Thus, the uncertainty in the range of PMIP assessed AMOC strength might be expanded with inclusion of LIS topographic uncertainty.

Implications for climate sensitivity
The LGM provides the most recent large-scale change in global climate and greenhouse gas concentrations through which global climate sensitivity can be assessed (Crucifix, 2006;Hansen et al., 2008;Schmittner et al., 2011;Hargreaves et al., 2012;PALAEOSENS, 2012). However, our model results suggest that significant differences in global mean temperature arise due to variability in LIS topography alone, indicating that the use of models to constrain CO 2 sensitivity for the LGM should include some assessment of this boundary condition uncertainty in the overall range of possible sensitivity estimates. This analysis is particularly relevant in the assessment of how uncertainty in boundary conditions leads to the development of fast feedbacks (i.e., snow-albedo) that are important in driving sensitivity (PALAEOSENS, 2012). Our 14 ka simulations show that even small differences in LIS height can lead to differences in atmospheric circulation, but these differences are not enough to initiate the snow-albedo feedback. Depending on the model (or other boundary conditions), this preconditioning may be sufficient to initiate fast-feedback mechanisms that lead to significant differences in global mean SAT, while other models may not provide the same response. Changes in LIS height within the uncertainty of reconstruction estimates may have a similar effect in other models.
Given the large uncertainties in the high-latitude boundary conditions and their associated impact on high-latitude climate, Hargreaves et al. (2012) instead correlated simulated LGM SAT from 20 • S to 30 • N with the global mean SAT change from CO 2 doubling sensitivity (2 × CO 2 ) simulations of each model in a grouping of PMIP2-CMIP3 pairings, which circumvented earlier issues with correlating global LGM SAT with 2 × CO 2 simulations (Crucifix, 2006). Our results suggest an additional source of high-latitude variability in LGM SAT that arises due to uncertainty in LIS topography alone, giving credence to the former study's focus on tropical SAT at the LGM to constrain CO 2 sensitivity. However, this additional source of uncertainty from ice-sheet topography in a region where LGM SAT anomalies are already at their greatest also suggests that focusing on the tropics could underestimate the full range of uncertainty in CO 2 sensitivity estimates.

Uncertainty in other boundary conditions
Finally, we have only assessed the impact of changing LIS height in one GCM. A number of other LGM boundary conditions have a significant impact on the global mean climate state, particularly dust/aerosol forcing (Penner et al., 2004;Lohmann and Feichter, 2005;Forster et al., 2007;Chylek and Lohmann, 2008) and vegetation land cover (Jahn et al., 2005). The large range in uncertainties of the LGM and deglacial dust forcing (Petit et al., 1981;Thompson et al., 1995;Mahowald et al., 1999;Reader et al., 1999; , 2006a, b;Power et al., 2008) and vegetation dynamics (Jackson et al., 2000;Prentice and Jolly, 2000;Harrison et al., 2001;Ray and Adams, 2001;Bigelow et al., 2003;Williams, 2003;Pickett et al., 2004;Williams et al., 2011) suggests the need to test the sensitivity of the climate response to the uncertainties in these boundary conditions. In the case of GISS ModelE2-R, future sensitivity studies will have offline coupling of ModelE2 boundary conditions to the Lund-Potsdam-Jena (LPJ) dynamic vegetation model (Sitch et al., 2003;Gerten et al., 2004; to assess the impact of this land surface choice and any associated internal feedbacks in the model.

Conclusions
We attempt to assess the range of climate uncertainty that results from variability in the possible reconstructions of the LIS elevation at two different time periods, 21 and 14 ka. The simulated climate at each time slice results in significant differences in atmospheric and oceanic climate. The 21ka-L is significantly colder than the 21ka-5G simulation in both SAT and ocean temperatures, which is due to a snow-albedo feedback in northeastern Asia that reduces the shortwave contribution to the system. This enhanced snow cover over northeastern Asia in 21ka-L develops due to a diminished atmospheric ridge relative to 21ka-5G. The higher elevation LIS in 21ka-5G is a larger impediment to atmospheric circulation, impacting the location of the atmospheric jet through resulting differences in downstream planetary wave generation, as has been suggested previously with LGM sensitivity tests removing the LIS (Justino et al., 2004;Otto-Bliesner et al., 2006;Abe Ouchi et al., 2007;Pausata et al., 2011;Hofer et al., 2012;Tharammal et al., 2012). We conclude that the range of uncertainty in existing LIS reconstructions is enough to have a similar "planetary wave effect" on Northern Hemisphere atmospheric circulation, contributing to the large differences in temperature, snow cover, and surface albedo over northeastern Asia. These fully coupled AOGCM results build upon earlier lower-resolution simulations that indicate a similar dynamical adjustment with the comparison of ICE-4G and ICE-5G LIS reconstructions in an Earth system model of intermediate complexity (Justino et al., 2006), suggesting a robust impact of paleotopographic uncertainty on reconstructed climates in a variety of models. Future research should work on minimizing the uncertainty in the LIS reconstructions, thus reducing their impact on the uncertainty in simulated glacial climate. As we show in our 21 ka simulations, this LIS topographic uncertainty may significantly impact estimates of LGM climate sensitivity through the impact on global mean SAT.
Given the significant differences between the resulting climates from changing LIS topography alone, it is tempting to use climate proxy data to test which reconstruction is closer to the actual LGM LIS. However, the regions with largest model differences are the regions with the least proxy coverage and with the greatest variation between proxy reconstructions. In particular, the region of Siberia, Beringia, and the North Pacific has the largest modeled differences in SAT, δ 18 O, surface albedo, SST, and sea concentrations, all due to variation in LIS topography alone in both 21 and 14 ka simulations. In theory, this region could serve as an important model-data test location on the LIS reconstructions used in GCMs. Unfortunately, the distribution of LGM SST proxy records from the North Pacific is limited (Kucera et al., 2005;Waelbroeck et al., 2009) and most records that do exist are confined along coastal regions, where GCMs may not adequately resolve changes in the Kuroshio and California Currents, and sea level boundary conditions may also significantly influence ocean current behavior (Ortiz et al., 1997;Kao et al., 2006). Similarly, few terrestrial reconstructions of SAT and δ 18 O a exist from lake records in the region (Melles et al., 2012;Peterse et al., 2011;Viau et al., 2008;Kurek et al., 2009;Bartlein et al., 2011), but perennial ice cover, proxy uncertainties, and the possibility of no-analog environments in pollen reconstructions may also limit the ability of such reconstructions to distinguish the climate differences discussed in this paper. To better constrain the topographic uncertainty of the LIS through the comparison of simulated climates in GCMs, more data are needed from regions where resulting climate differences are the greatest. Future data collection should focus on northeastern Asia and the North Pacific to help test LIS boundary conditions.
In addition to large differences in atmospheric circulation and attendant impacts on surface conditions, the two LIS simulations at each time slice result in significant differences in the AMOC, with the larger ICE-5G ice sheet leading to stronger overturning at each time slice, suggesting that any given model's ability to capture LGM anomalies in ocean overturning may be influenced as much by the LIS boundary condition as by limitations in model physics. At present, proxy uncertainty may preclude determination of whether the LGM AMOC was different than for the modern period (Burke et al., 2011). This uncertainty in LGM ocean circulation and surface temperature, particularly in regions where model differences are maximized, limits the determination as to which LIS reconstruction is more correct in our simulations.
One metric that may help differentiate between the LIS reconstructions is the Kr / N 2 reconstruction of global mean ocean temperature (Headly and Severinghaus, 2007), as it should be indicative of the overall thermal state of the global system. Here, the 21ka-L simulation results in a global mean ocean temperature change more in line with the proxy reconstructions, whereas the 21ka-5G mean ocean is too warm (see Sect. 3.2.1). As this oceanic temperature difference is the result of the different atmospheric cooling from changes in planetary albedo, we would suggest that a lower elevation LIS topography than 5G may be closer to what existed at the LGM. However, the uncertainty in Kr / N 2 global mean ocean temperature is large relative to the reconstructed change, therefore limiting the distinction of the two LGM simulations based on this proxy. Further refinement of LGM and 14 ka climate reconstructions (particularly from Siberia and the North Pacific) will aid in the indirect determination of LIS topography, just as refinement of LIS topographic reconstructions will reduce inherent model uncertainty due to the current range in LIS topographic boundary conditions.