Coupled ice sheet–climate modeling under glacial and pre-industrial boundary conditions

In the standard Paleoclimate Modelling Intercomparison Project (PMIP) experiments, the Last Glacial Maximum (LGM) is modeled in quasi-equilibrium with atmosphere–ocean–vegetation general circulation models (AOVGCMs) with prescribed ice sheets. This can lead to inconsistencies between the modeled climate and ice sheets. One way to avoid this problem would be to model the ice sheets explicitly. Here, we present the first results from coupled ice sheet–climate simulations for the pre-industrial times and the LGM. Our setup consists of the AOVGCM ECHAM5/MPIOM/LPJ bidirectionally coupled with the Parallel Ice Sheet Model (PISM) covering the Northern Hemisphere. The results of the pre-industrial and LGM simulations agree reasonably well with reconstructions and observations. This shows that the model system adequately represents large, non-linear climate perturbations. A large part of the drainage of the ice sheets occurs in ice streams. Most modeled ice stream systems show recurring surges as internal oscillations. The Hudson Strait Ice Stream surges with an ice volume equivalent to about 5 m sea level and a recurrence interval of about 7000 yr. This is in agreement with basic expectations for Heinrich events. Under LGM boundary conditions, different ice sheet configurations imply different locations of deep water formation.


Introduction
The understanding of past climates is a major challenge.Well-documented, non-linear climate changes can be studied by modeling past climates that differ greatly from the present-day climate.This improves the understanding of the climate system (e.g., Valdes, 2011;Braconnot et al., 2012).
The last glacial climate differed greatly from the presentday climate, although the last glacial is very recent from a geological perspective.Therefore, the proxy coverage and dating accuracy are comparatively good.The climate differences between the last glacial and the present mainly stem from the reduced greenhouse gas content of the atmosphere and from the massive ice sheets covering the continents of the Northern Hemisphere (Rind, 1987;Kim, 2004;Abe-Ouchi et al., 2007;Pausata et al., 2011).The ice sheets were at their maximum volume about 21 000 yr ago, during the so-called Last Glacial Maximum (LGM) (e.g., Clark et al., 2009).In atmosphere-ocean general circulation models (AOGCMs), the LGM climate is commonly modeled in steady-state experiments (e.g., in the Paleoclimate Modelling Intercomparison Project, PMIP).
Previous AOGCM studies of the LGM climate, especially the PMIP experiments, have relied on prescribing the ice sheets from reconstructions (Braconnot et al., 2007(Braconnot et al., , 2011(Braconnot et al., , 2012)).While the original PMIP experiments prescribed sea surface temperatures (SSTs) or oceanic heat transports, and PMIP-2 added the full ocean models to the complexity, interactive ice sheets are still missing, even in PMIP-3.This can result in climate and ice sheets not being consistent with each other.To obtain ice sheets that are consistent with the modeled climate, it is essential to overcome this method and to model the ice sheets interactively.We took exactly this step.We extended the PMIP-2 setup by an interactively coupled ice sheet model (ISM), and studied the climate as well as the ice sheets in this self-consistent system with its increased complexity.By comparing the modeled ice sheets with reconstructions, we can critically assess our model system.
An interactively coupled ice sheet-climate model opens the possibility of studying the interactions between ice sheets and the climate system.In the light of the recently observed changes in the ice sheets and their potential impact on future sea level rise (e.g., Vaughan et al., 2013), such studies become increasingly relevant for predicting future climate change.The large ice masses and the possibility of comparing the model results to proxy data make the LGM an ideally suited time period for studying ice-climate interactions.
The long intrinsic timescales of ice sheets make the assumption of the ice sheets being in equilibrium with the climate problematic.These timescales also require long simulations that are challenging to perform with AOGCMs.Complicating the modeling further, most of the snowfall and melt occurs in a narrow band along the ice sheet margin.In principle, the margin zone calls for high-resolution modeling, while the long timescales prohibit this.The discrepancy between the need for high-resolution modeling and long time spans leads to a variety of strategies (see the reviews by Pollard, 2010, andVizcaíno, 2014).
At the fast and simple end of the model spectrum, fixed climate maps or analytical expressions are scaled by the use of a time series, either from an energy balance model or from ice core data.In earlier times, these models were used to study one glacial cycle (e.g., Pollard, 1982;Greve, 1997;Tarasov and Peltier, 1997).Advances in computing power now make it possible to model the evolution of the ice sheets over millions of years with such models (Pollard and DeConto, 2009).Tarasov et al. (2012) performed multiple experiments covering the last glacial.They chose the experiments that agreed best with the proxy data and created a deglacial ice sheet chronology that is consistent with proxy data as well as with ice physics.Zweck and Huybrechts (2003) forced an ISM by interpolating between modeled climate states for LGM and modern boundary conditions using Greenland ice core data as weighting coefficients.Charbit et al. (2007) compared the effect of different AGCMs on the same ISM using a similar approach.They found the Eurasian ice sheets to be more sensitive to the choice of the climate model than the Laurentide.Abe-Ouchi et al. (2007, 2013) used multiple experiments performed with the same climate model to create a statistical model that accounts for changes in greenhouse gases, orbital parameters, and ice sheets separately.They coupled this model to an ISM and performed simulations of the last glacial cycle.They found the snow-albedo feedback and the temperature lapse rate to be the most important feedback fac-tors.For the deglaciation, the delayed isostatic rebound of the lithosphere played a key role.
Another approach uses earth system models of intermediate complexity (EMICs) and optionally downscales the fields for the ISM in an intermediate energy balance model or energy-moisture balance model.Calov et al. (2002) used such a setup, and obtained binge-purge oscillations of the Laurentide Ice Sheet resembling Heinrich events as internal oscillations of the ice sheet.Wang and Mysak (2002) and Calov et al. (2005) studied the feedbacks and bifurcations during the glacial inception, or more recently, full glacial cycles (Bonelli et al., 2009;Ganopolski et al., 2010;Ganopolski and Calov, 2011).Recently, Heinemann et al. (2014) performed transient simulations starting at 80 kyr BP and obtained a much smaller LGM Fennoscandian Ice Sheet in the transient simulation than in a steady state sensitivity study.This shows that the Fennoscandian Ice Sheet was not in an equilibrium state at the LGM.
The most expensive approaches use coupled AOGCM-ISM systems.While the EMIC or energy balance model based methods can cover long time spans with reasonable effort, they strongly simplify the atmosphere and/or ocean dynamics.In contrast, AOGCMs suffer from high computational costs but provide a detailed representation of the atmosphere and ocean dynamics.To bridge the long time spans, most research groups couple the models asynchronously.In this method, the climate model is run for a period of one to fifty years, and then the climate is used to drive the ice sheet for a substantially longer period.The resulting ice sheet is fed back into the climate model and the cycle is started over (Pollard et al., 1990;Ridley et al., 2005;Mikolajewicz et al., 2007a;Gregory et al., 2012).
Several previous AOGCM-ISM studies investigated the future development of the Greenland Ice Sheet and its impact on the Atlantic meridional overturning circulation (AMOC).Early studies by Huybrechts et al. (2002) and Fichefet et al. (2003) only fed back the fresh water fluxes, but kept the Greenland topography constant in the AOGCM.The first study with full bidirectional coupling between an AOGCM and an ISM was published by Ridley et al. (2005).It was followed by studies that also include the Antarctic Ice Sheet (Mikolajewicz et al., 2007a, b;Vizcaíno et al., 2008Vizcaíno et al., , 2010)).Most of these studies agree in that they find little effect of the shrinking Greenland Ice Sheet on the large-scale climate on a century scale.For large losses of the Greenland Ice Sheet, they find changes in the large-scale atmospheric circulation.The effect of the melt water on the ocean circulation is small compared to the effect of the increased temperatures and precipitation in most studies.The experiments presented in Mikolajewicz et al. (2007b) and Vizcaíno et al. (2010) were the first atmosphere-ocean-vegetation general circulation model (AOVGCM)-ISM studies to use an energy balance model for the surface mass balance of the ice sheet, and the first to operate without anomaly maps.This improved the physical representation of the coupling.A recent study by Gregory et al. (2012) used the AOGCM FAMOUS, a fast version of HADCM3, in combination with the ISM GLIM-MER to study the last glacial inception.They found a strong dependence of the ice volume on the way the ice-albedo feedback is treated in the coupling.None of the published studies addressed LGM conditions.
We have coupled a coarse-resolution AOVGCM (ECHAM5/MPIOM/LPJ; Mikolajewicz et al., 2007b) interactively with a state-of-the-art ISM, the modified Parallel Ice Sheet Model (mPISM).The coupling is performed bidirectionally, and without flux corrections nor anomaly methods.Our intention is to apply the model system for transient long-term simulations of, e.g., the last termination.As a first step, we will validate here the ability of the model system to reproduce the climate of two time slices, which are sufficiently different that they allow us also to validate the model's sensitivity to perturbations.We chose the pre-industrial and the last glacial maximum.The only factors we varied between the LGM and the pre-industrial setup were the greenhouse gas concentrations, the orbitals and the shapes of the continents.We let the models evolve freely sufficiently long to be largely independent of the initial state.We study the behavior of the climate system and compare the model results to observations and proxy data.As the model system is new, this paper also has a strong focus on the description of the model system and the coupling between the different components.
We describe the models, the necessary modifications, the coupling and the setups in Sect.2, analyze the mean states of the experiments in Sect.3, and summarize the main findings and draw the conclusions in Sect. 4.

Model description and setups
The AOVGCM has been applied before coupled to the ISM SICOPOLIS for studies of the future evolution of the Greenland and Antarctic ice sheets (Mikolajewicz et al., 2007b;Vizcaíno et al., 2010).For this study, we have switched to PISM, and changed various aspects of the coupling.The models and the coupling will be described in the following.

ECHAM5/MPIOM/LPJ
The atmospheric component of the coupled model is ECHAM5 (Roeckner et al., 2003), a spectral atmosphere general circulation model.For long-term simulations, the triangular spectral truncation at wavenumber 31 (T31, ∼ 3.75 • ) in combination with 19 vertical hybrid-σ -levels reaching up to 10 hPa is a compromise between computational demand and accuracy.This setup is therefore employed in all experiments described in the following.The land hydrology scheme (HD model, Hagemann and Dümenil, 1998;Hagemann and Gates, 2003) operates on a 0.5 • grid.For the HD model, we use the present-day routing direc-tions in combination with the LGM land-sea-mask.The land surface properties are modeled using the Lund-Potsdam-Jena (LPJ) vegetation model (Sitch et al., 2003).Its use in combination with ECHAM is described in Schurgers et al. (2007) and Mikolajewicz et al. (2007b).
MPIOM (Marsland et al., 2003) is a primitive equation ocean model operating on a curvilinear grid with variable resolution.In the setup employed in the following, the grid has two poles, located over Greenland and Antarctica, and a nominal resolution of 3 • .This results in an increased resolution in the deep water formation areas, and a relatively coarse resolution in the equatorial areas.
The atmosphere and ocean are coupled using the OASIS coupler (Valcke et al., 2004).The performance of a higherresolution version of the coupled AOGCM is described in Jungclaus et al. (2006).In the framework of the PMIP-2 the AOGCM ECHAM5/MPIOM/LPJ has been applied to the LGM.Results can be found in the PMIP-2 database at http://pmip2.lsce.ipsl.fr.

mPISM
mPISM is based on PISM version 0.3 from the University of Alaska, Fairbanks (Bueler and Brown, 2009;the PISM authors, 2014).PISM uses the Shallow Ice Approximation (SIA) and the Shallow Shelf Approximation (SSA) to compute flow velocities.It uses an enthalpy method to handle polythermal ice (Aschwanden et al., 2012).Details about the model can be found in the literature given above.Several aspects of the model needed to be changed for the coupling to the climate model and for obtaining pulsating ice streams.In the following, we describe the modifications of the PISM physics.For the technical changes see Ziemen (2013).We use the term PISM 0.3 when referring to the base version, mPISM when referring to the modified version, and PISM when referring to aspects of the base version that also apply to the modified version.PISM-PIK (Winkelmann et al., 2011) is a branch of PISM that is developed at the Potsdam Institute for Climate Impact Research (PIK).
Following Calov et al. (2002), we use a linear sliding law with a friction coefficient of 1 m yr −1 Pa −1 that allows sliding if there is basal water and deformable sediment available to lubricate the ice sheet.The availability of the sediment is based on the data set from Laske and Masters (1997) with a cutoff value of 5 cm, and is marked in Fig. 1.The availability of basal water is a prognostic quantity of PISM.We spread out half of a grid cell's heat flux from basal friction on the cell's eight neighbors.This slightly heats the grid cells adjacent to an ice stream.
We separate the basal water into a small locally bound fraction and an advectable fraction that we advect with the basal ice velocity, and apply the same upwind scheme as used for ice and enthalpy transport.The locally bound fraction controls the stiffness of the basal till, and thereby the sliding behavior.Ice streams.Brownish/red colors mark the fraction of time that a grid cell is sliding at more than 1 m yr −1 .Blue marks areas where the ice is not permitted to slide due to the lack of sediments in the reconstruction of Laske and Masters (1997).Only the time when the grid cell is ice covered is taken into consideration; therefore, ice shelves are considered to be constantly sliding.Numbers match the ice stream numbering in Stokes and Tarasov (2010), and are as follows: (1) Mackenzie, (16) Ungava Bay, (18) Amundsen Gulf, For the bedrock deformation, we use a Local Lithosphere Relaxed Asthenosphere (LLRA) model (e.g., Le Meur and Huybrechts, 1996) with a rebound timescale of 3 kyr.We dynamically compute the change in sea level that is caused by the ice sheets and use it to adjust the sea level in the ISM and in the climate model.
Our setup uses a Cartesian coordinate system that is a Polar Stereographic projection of the Northern Hemisphere, covers all areas north of 36.7 • N, and reaches south to 19.1 • N in the corners.This grid covers all Northern Hemi-sphere regions where we can possibly expect to grow largescale ice sheets under glacial conditions (including the Himalayas).The resolution is 20 km in both directions (625 × 625 grid cells).Outside of the mPISM domain, we use the ICE-5G (Peltier, 2004) topography matching the time slice of interest, as in PMIP-2.

The coupling from the climate to the ice sheet
The coupling scheme between ECHAM5, MPIOM and mP-ISM computes the ice sheet mass balance and surface temperature from the ECHAM5 and MPIOM output, and transfers surface topography, glacier mask and mass fluxes from mPISM back to ECHAM5 and MPIOM.We will now lay out the path from the climate model to mPISM.
To compute the mass balance for the ice sheet, we have to use a high spatial resolution that resolves the temperature distribution at the ice sheet margins.We use the 20 km grid of mPISM for the mass balance calculations.To save computational time, we compute monthly averages of the atmospheric quantities and annual averages of the ocean variables before regridding them to the ISM grid.These regridded high-resolution fields are used in the ISM to determine the mass balance.The scheme we use to process temperatures and precipitation from ECHAM5 for the use in mPISM is based on the standard scheme employed in SICOPOLIS, which stems from Braithwaite and Olesen (1989).
To account for surface elevation differences between the two models, we correct the temperatures using a lapse rate of −5 K km −1 as suggested by Abe-Ouchi et al. (2007).We use the height-corrected monthly mean temperatures to partition the precipitation into solid and liquid fractions using a linear transition between −10 and +7 • C as in Marsiat (1994).The solid fraction P solid is used as accumulation, and the liquid fraction is discarded immediately as runoff.To account for the reduced precipitation at high altitudes, we apply a height desertification parametrization as in Budd and Smith (1979) at heights above 2000 m.
where h ISM is the surface height in mPISM, h GCM is the surface height in the climate model, and λ = log(2)/1000 m.We use temperatures and precipitation in a positive degree day (PDD) model to determine the surface mass balance.Despite its limitations (e.g., van de Berg et al., 2011), the PDD method still is widely used in glaciological and coupled model studies (e.g., Gregory et al., 2012).As an extension of Fausto et al. (2011), we compute monthly temperature standard deviation maps from the 6-hourly climate model output, and use them in mPISM to represent the temperature variability better.This method represents differences in the diurnal cycles and synoptic variability on the different ice sheets.It is a compromise between computing the PDDs from the full 6 hourly temperature data on the highly resolved Clim.Past, 10, 1817-1836, 2014 www.clim-past.net/10/1817/2014/ice grid (which would be more computationally demanding) and the standard approach of prescribing a fixed standard deviation in the full model domain.The PDD scheme employs the Calov-Greve integral method (Calov and Greve, 2005) to compute PDDs from monthly mean temperatures and standard deviations.Then, from these PDD, and the snow accumulation, the actual mass balance is computed using a standard PDD scheme (Reeh, 1991).We chose a snow melt rate of m snow = 3.2 mm ice equivalent K −1 day −1 , an ice melt rate of m ice = 12.9 mm ice K −1 day −1 and allow up to 60 % of the melt to refreeze in the snowpack.The ice melt rates are based on Greve et al. (1999), who chose 3 mm and 12 mm water equivalent per day and degree Celsius in Northern Hemisphere simulations.
For the basal ice shelf melt (or growth), we use the threeequation scheme described in Holland and Jenkins (1999).This scheme computes the ice shelf basal mass balance and properties of a thin boundary layer from ocean salinity and temperature as well as the ice shelf basal temperature gradient.We force this model with the annual mean temperature and salinity averaged over the top 12 layers (203 m) from the ocean model.
We employ a simple calving scheme that operates on a grid cell and its eight neighbors.If less than three of the nine cells have an ice thickness above 200 m, and the relaxed bedrock topography is below sea level, then the ice is calved.This allows shelves to grow and retreat dynamically, and prevents one grid cell wide ice tongues.

The coupling from the ice sheet to the climate
In the following, we describe how the output fields from the ISM are transferred into the atmosphere and ocean models.For the areas outside of the mPISM domain, the mPISM surface field is combined with a background map.The combined map is adjusted for sea level changes from the ice sheet.To keep the subgrid-scale orography consistent, the background map is smoothed with respect to the standard topography map used for ECHAM5, so it matches the roughness of the ice-free areas of the remapped mPISM topography.From this map, both the mean topography as well as subgrid-scale topographic variables like the maximum surface height in a grid box are calculated.Since this smoother surface would significantly change climate in comparison to the standard model setup, we upscale the subgrid surface slopes by a factor of 2. Furthermore, we decrease the resolution-dependent threshold for activating the gravity wave drag parametrization in ECHAM5 from the standard T31 value of 400 m to 300 m peak -mean elevation, so the gravity wave drag parametrization is active in about the same grid cells as in the original setup.
All ISM grid cells with an ice thickness above 10 m are treated as glaciated.We use conservative area remapping to interpolate the glacier mask to the ECHAM5 grid.ECHAM5 does not allow for fractional glaciation of a grid box, and treats all grid cells with a fractional glacier mask value above 0.5 as fully glaciated.To generate a more gradual transition between fully glaciated and non-glaciated grid cells, we modified the ECHAM5 albedo scheme.For grid cells that ECHAM5 considers as glaciated, the albedo is computed based on a fractional glacier cover with a background albedo of 0.25 for the non-glaciated parts of the grid cell.For those grid cells that have a glacier fraction below 0.5 and are thus considered as non-glaciated in ECHAM5, we reduce the forest (tree cover) fraction by the glacier fraction to obtain an effective forest cover (forest eff ) and compute a new background albedo α eff that represents the effects of the albedo of the glaciated parts: (2) where forest is the forest fraction computed by LPJ, glacier the glacier fraction from the remapped the mPISM output, α lpj the background albedo calculated from the vegetation model, and α glacier a background albedo for melting glacier ice (0.5).The reduction of the forest fraction is important because it reduces the snow-masking effect in the albedo calculation of ECHAM5.
The mass (fresh water) flux coupling preserves flux rates.We associate fluxes from iceberg calving and shelf basal melt with a negative enthalpy flux into the ocean, so that these processes lead to ocean cooling and/or increased sea ice formation.The net surface mass balance on the ice sheets is subtracted from the precipitation-evaporation field fed into the surface runoff model.Similarly, the effects of iceocean interactions are directly fed into the MPIOM's surface scheme.All grid conversions are performed using conservative remapping.In experiments with synchronous coupling, this yields perfect mass conservation.
We keep the ocean bathymetry and land sea mask constant.The freshwater fluxes change the surface height of the ocean.In ECHAM5 we reduce the land surface elevation when the sea level rises, and vice versa.
For the experiments presented in this publication, we use an asynchronous coupling scheme, where mPISM is run for 10 years while the climate model is run for 1 year between the coupling procedures.In this setup, our scheme conserves mass fluxes instead of total mass.This method has been successfully employed in previous simulations (Vizcaíno et al., 2008), and avoids overly strong freshwater forcings for the ocean model.Such forcings would arise from conserving mass and increasing the flux rates.
The coupling scheme does not conserve energy at the iceatmosphere interface.One way to obtain this and keep the models consistent would be to run an energy balance model on different height levels in the GCM and force the ISM with the modeled surface mass balance (Vizcaíno et al., 2013).We do not change the energy balance of the surface scheme of the atmosphere model and therefore do not introduce any LGM None 100 LGM-mixed LGM None 90 new energy sources or sinks.By explicitly modeling the ice sheets, our model however is more realistic than GCM-only simulations, where the ice sheet extent is prescribed.

Setups and experiments
Our experimental strategy was to perform two coupled experiments, one under glacial and one under pre-industrial conditions, and to compare the results with each other and with observations and proxy data.We additionally performed three sensitivity studies under LGM conditions: LGM-ICE-5G, LGM-mPISM-W, and LGM-mixed.They are described in Sect.3.4.We list all experiments in Table 1.
The coupling and tuning of the models required iterative improvements.To avoid costly repeated transient spinups of the LGM state, we forced the models with constant LGM conditions (Table 2).Thus, the ISM was run for more than 100 000 years under LGM conditions with increasing completeness of the coupling.The climate model was started from an existing LGM state (Arpe et al., 2011) and run for 3850 years during the coupling.The last 21 500 ISM years (2150 climate model years), the models were fully coupled.From the resulting state, we started our main LGM experiment, LGM-mPISM.We integrated this experiment for 30 000 years in the ISM (3000 years in the climate model).Unless explicitly stated otherwise, we average over the full duration of this experiment in this publication.
The second main experiment is the pre-industrial control run PI-mPISM.The ISM was initialized with the presentday Greenland Ice Sheet shape using PISM's bootstrap methods.They provide an empirical guess for the temperature profile in the interior based on the surface air temperature and the geothermal heat flux at the base.The climate model was started from an existing pre-industrial state.The coupled model was spun up for 9000 years in the ISM (900 years in the climate model) and the following 1000 (100) years were analyzed.In the coupled experiments, the coupling was performed after every year of climate model simulations (10 years of ISM simulations).

Results and discussion
In the following, we first discuss the pre-industrial experiment (PI-mPISM) and then the LGM mean state in LGM-mPISM and LGM-ICE-5G.Finally, we describe the longterm drift of the ice sheets in LGM-mPISM.The main purpose of PI-mPISM is to show that our setup yields reasonable results under pre-industrial boundary conditions and to provide a reference state to compare the LGM climate against.We therefore limit the analysis of this experiment to key climate parameters and compare those to reanalysis and observational estimates.We largely draw on data for present-day conditions.These are available at high quality, and the differences to the pre-industrial climate are minor compared to the LGM effects we focus on in this study.We first discuss the state of the atmosphere, then the ocean, and finally the ice sheets.

The pre-industrial atmosphere
A comparison of the 2 m air temperature (SAT) in PI-mPISM with the ERA INTERIM reanalysis for 1989 to 2010 (Dee et al., 2011) shows that the annual global mean SAT in PI-mPISM is lower than the present-day values by 0.6 K (see Fig. 2).This difference can be attributed to the lower greenhouse gas concentrations in the pre-industrial time.In a Coupled Model Intercomparison Project Phase 3 (CMIP3) simulation performed with a higher-resolution version of ECHAM5/MPIOM, the temperature difference between the periods 1860 to 1899 (pre-industrial) and 1989 to 2010 (ERA INTERIM period) is 0.56 K.The temperature difference over land (−0.8 K) is larger than over the ocean (−0.5 K).The difference is also larger at the high latitudes (−3.2K north of 45 • N, and −1.2 K south of 45 • S) than at the low latitudes, where, on average, the temperature in our model is higher than in the ERA INTERIM reanalysis (+0.1 K between 45 • S and 45 • N).The model shows a significant cold bias over Siberia.
The smoothed-out T31 grid of ECHAM5 is lower than the T255 grid of ERA INTERIM in most mountain areas (the creation of the grid involves spectral smoothing and is not conservative, so the mean surface elevation is lower in a T31 grid than in reality).Differences that can directly be traced back to the lower topography in the T31 resolution of the atmosphere model are the warm biases over the Andes and the Himalayas.Over the Himalayas, typical surface altitude differences between the two model setups are about 500 to 1000 m on the T31 grid, over the Andes, they reach 2000 m, and over the southern tip of Greenland, they reach 1000 m.Assuming a lapse rate of 5 K km −1 , as it is used in our coupling scheme, this corresponds to a 2.5-5 K temperature difference over the Himalayas, 5-10 K temperature difference over the Andes, and up to 5 K temperature difference over the southern tip of Greenland.
There is a cold bias over the northern Atlantic.This is a consequence of the North Atlantic Current taking too southerly a route.Eddy-resolving modeling has been shown to improve the representation of the Gulf Stream separation and the path of the North Atlantic Current (e.g., Hurlburt and Hogan, 2000;Bryan et al., 2007), but is not yet feasible for multi-millennial simulations.
There is a strong cold bias over the Alaska Range that also occurs in stand-alone simulations (not shown) and becomes stronger because of ice-sheet growth in the coupled setup (Fig. 3).This temperature bias, and thus the glaciation, can be reduced by increasing the climate model resolution (e.g., in the CMIP3 experiments that were performed with a technically identical climate model version), but this is not yet feasible for multi-millennial experiments.
A comparison of the present-day annual mean precipitation data from the GPCP data set (Adler et al., 2003) with PI-mPISM (Fig. 4) shows a good agreement over Greenland.This is also reflected in the agreement of our Greenland mean surface accumulation values with those from regional modeling (Table 3).Over the other polar areas, the values also match closely.Most of the precipitation differences occur in the tropics and at southern mit latitudes (not shown).The differences in the tropics are standard modeling artifacts, and neither region is of special interest for this study.In the global mean, the precipitation of 1.03 m yr −1 is 5 % above the GPCP estimate, while north of 45 • N, the modeled precipitation of 0.66 m yr −1 is 10 % lower than the GPCP estimate.

The pre-industrial ocean
The North Atlantic Deep Water (NADW) cell of the AMOC peaks at 17.0 Sv (1 Sv = 10 6 m 3 s −1 ) at 32.5 • N and a depth of 1020 m.This agrees with the estimate of 16 ± 2 Sv from Ganachaud (2003) and recent measurements of 18.7±2 Sv at 26.5 • N (Kanzow et al., 2010).The NADW is formed south of Greenland and in the Nordic seas (Fig. 5).The Antarctic  (Box et al., 2006), and MAR (Fettweis, 2007) data are for the period 1958-2007(RACMO2, MAR), resp. 1958-2006 (PMM5) (Koltermann et al., 2011) shows that, in our model, the deep water flowing southward is warmer than it is in reality.This explains why our model simulates less heat transport than the estimates derived from observations indicate, while the overturning strength is similar.
The sea ice maximum and minimum extent agree very well with the long-term average of the HadISST observational data set (Rayner et al., 2003) (Fig. 3).

The pre-industrial ice sheets
Figure 3 shows the ice sheets in PI-mPISM averaged over the last 1000 yr, as well as the deviations from the reference topography.We obtain a Northern Hemisphere land ice volume of 5.9 Mio km 3 , corresponding to 14.9 m of sea level equivalent (SLE) (Table 4).Of this volume, 3.65 Mio km 3 (9.2m SLE) are stored in the Greenland Ice Sheet (see Fig. 6 for the mask used in the analysis).Today's Greenland Ice sheet has a volume of 2.9 Mio km 3 (7.3m SLE).The drift in Greenland Ice sheet volume is −14 km 3 yr −1 .It has a two-dome structure; the main dome reaches a height of 3200 m a.s.l.(above sea level) (3300 m in reality, Bamber et al., 2001) and the southern dome reaches 2700 m a.s.l.(2900 m in reality).
Along most of the coasts of Greenland, the model grows too much ice.At the northern and northern east coast (north of 70 • N), this is largely due to practically zero ablation from the PDD scheme.The glacier fraction in ECHAM is above 0.5 from the start of the experiment.The grid cells are therefore treated as glaciated, and the surface temperature in these grid cells cannot rise above 0 • C.This substantially limits the SAT and thus the melt in the PDD scheme.The sea ice along the northern and eastern coast maintains cold temperatures in the coastal ocean.The temperature is interpolated bilinearly Over the ocean, dark gray with a red outline indicates areas with perennial ice cover (15 % level) in more than 50 % of the model years.Light gray with a dark blue-black outline indicates temporary ice cover in more than 50 % of the model years.The Bering Strait appears to be ice free, because it is slightly displaced to the west in the ocean model.The orange outline marks areas that have permanent sea ice cover in more than 50 % of the years 1870 to 2010 according to the HadISST sea ice data set (Rayner et al., 2003).The light blue outline marks areas that have temporary sea ice cover in more than 50 % of these years according to HadISST.(b) The difference between the modeled pre-industrial topography and the present-day topography (ETOPO1; Amante and Eakins, 2009).
when remapping from ECHAM5 to mPISM.The mean temperature between a sea ice covered grid cell and a glaciated grid cell cannot allow for substantial surface melt.Therefore, the surface mass balance is positive practically all the way to the coast, while substantial ice melt would be needed to stop the glaciers before the coast.There is substantial melt near the western coast where two ECHAM5 grid cells with a glacier fraction below 0.5 remain that are considered as non-glaciated by ECHAM5.An energy balance scheme with detailed treatment of the different heat fluxes (e.g., Vizcaíno et al., 2010) could solve the problems at the eastern coast, but would have required substantial additional resources.Another way to improve the representation of the Greenland Ice Sheet margins is to use a very fine model resolution ( 5 km) in the ISM that allows for resolving of the mountain ranges and the individual outlet glaciers.The surface velocities in the northern part of the ice sheet agree reasonably well with the observations of Joughin et al. (2010a, b) (Fig. 7).The ridge of the ice sheet can be seen as a low-velocity band, and is very well captured in the northern part.In the northern part of the eastern coast of Greenland, the model shows too much ice because of the lack of ablation described above.The ridge of the ice sheet is displaced to the east.Since, at the eastern coast, the model shows ice in areas that are not glaciated in reality, the flow velocities cannot be compared to observations in this region.In the southeast, there is a lack of observations that prohibits a comparison.In the western part, we capture the general features well, although our velocities generally are too large, and the ice sheet reaches further towards the coast than in reality (Fig. 3).
Ice caps form on several Arctic Islands: Baffin Island and Ellesmere Island, Svalbard, Franz Josef Land, Novaya Zemlya, and Severnaya Zemlya.All of these regions presently also show glaciation.There are ice caps with a total volume of 0.94 Mio km 3 (+86 km 3 yr −1 ) growing in northeastern Siberia, and there is an ice sheet with a volume of 1.1 Mio km 3 (+31 km 3 yr −1 ) in the Alaska Range and the northern Rocky Mountains.In these regions, our climate model shows a cold bias (Fig. 2, Sect.3.1).This cold bias in regions that are characterized by many glaciers in reality leads to a glaciation in mPISM that quickly grows because of the positive feedbacks of increasing altitude and albedo.The growth of an ice cap in northeastern Siberia is fostered by the general cold bias over northern Siberia, which leads to an underestimation of the summer melt.

LGM climate experiments
In the following, we discuss the mean state in our LGM model experiments.The results from the main experiment (LGM-mPISM) are averaged over the full 3 kyr for climate model data, and over the corresponding 30 kyr for ISM data.
The mean ice sheet topography in LGM-mPISM is displayed in Fig. 8.To investigate the climatic effects of the modeled ice sheets further, we compare LGM-mPISM with a climateonly experiment with prescribed ice sheets from the ICE-5G reconstruction of Peltier (2004) (Fig. 8, called LGM-ICE-5G in the following), which is consistent with the PMIP-2 protocol.It was started from the same PMIP2 LGM experiment as the spin-ups for LGM-mPISM, and spun up over 549 years.We analyze the following 100 years (climate only).
In addition, we performed two further sensitivity studies: LGM-mPISM-W and LGM-mixed.In LGM-mPISM-W, we cut the heat fluxes from ice shelf melt and iceberg calving, so the ice enters the ocean as liquid water.The comparison of the model results with LGM-mPISM shows the effects of these heat fluxes for the Arctic Ocean sea ice.It is split off from LGM-mPISM after 16 500 years in the ISM (1650 years in the climate model) and has a duration of 7000 years in the ISM (700 years in the climate model).Data are averaged over the full duration of the experiment.In LGM-mixed, we started from LGM-ICE5G, and replaced the ICE-5G topography with the topography from LGM-mPISM, keeping the ICE-5G glacier mask in place.This (inconsistent) set of boundary conditions allows us to separate the effect of the topography on the large-scale circulation from the albedo and temperature effects of the glacier mask.
The following analysis starts with the atmosphere, then continues with the ocean, and the ice sheets.To facilitate the understanding of the ice sheet effects on atmosphere and ocean dynamics, we here provide a brief overview of the main differences between the modeled ice sheets in LGM-mPISM and in the ICE-5G reconstruction that impact the climate.In ICE-5G as well as in the other reconstructions, the European part of the Fennoscandian Ice Sheet is larger than in LGM-mPISM (Fig. 8), and reaches south across the Baltic and onto the British Isles, while in LGM-mPISM it is limited to Scandinavia and the Barents Sea in the west and reaches far into Northeast Siberia.According to the reconstructions, this region showed only small-scale glaciation.In the course of LGM-mPISM, this extension of the Fennoscandian Ice Sheet connects to an extension of the Laurentide Ice Sheet that has formed over Alaska, the Bering Sea and Kamchatka, closing the gap between the ice sheets.The mod-eled Laurentide Ice Sheet is thicker in the north and does not reach as far south as in ICE-5G.The ICE-5G reconstruction features a massive North-South-ridge over western Canada that is neither shown by the model, nor by other reconstructions (Sect.3.7).One major contributor to these differences is the experimental strategy of performing steady-state experiments and spinning up the ice sheets over a long time with constant LGM climate conditions.This leads to an ice sheet state that differs from the real LGM.Except for the growing Fennoscandian Ice Sheet, the modeled ice sheets are in quasi-equilibrium with the surface accumulation.

The LGM atmosphere
In LGM-mPISM, the global annual mean SAT is reduced by 3.5 K compared with PI-mPISM (Fig. 2 b), in agreement with 4.0 ± 0.8 K from a proxy data interpolation by Annan and Hargreaves (2013).The cooling over land (5.4 K) is larger than over the oceans (2.5 K).Applying the same T31 landsea mask to the SAT reconstruction by Annan and Hargreaves (2013) yields 6.3 K over land, and 3.0 K over ocean areas.Areas that are glaciated in LGM-mPISM cool by 11 K, non-glaciated land by 3.3 K.
The main deviations from the large-scale cooling during the LGM are as follows (Fig. 2b): the northern rim of the Pacific and the northern Atlantic south of 50 • N warm in the annual mean.The northern Pacific warms because of a change in stationary eddies due to topographic changes in the East Siberian ice sheet.This enhances the advection of warm air from the south.The warming over the Atlantic results from an increase in the oceanic heat transport (Sect.3.6).Kazakhstan, the Ural Mountains, and the West Siberian Plain warm in summer for two reasons: decreased cloud cover allows more shortwave radiation to reach the ground, and overcompensates for the effect of the higher surface albedo.Reduced precipitation leads to less evaporative cooling.
In LGM-ICE-5G, the cooling is stronger than in LGM-mPISM with a mean of 5.3 K compared to PI-mPISM (Fig. 2 c).This is slightly outside the temperature envelope of 4.0±0.8K provided by Annan and Hargreaves (2013), but within the envelope of the PMIP2 ensemble published in Braconnot et al. (2007).Over the oceans, LGM-ICE-5G cools by 4.1 K, over land by 7.7 K, over the areas that are glaciated in ICE-5G by 14 K, over non-glaciated land by 5.4 K.
We further studied the ice sheet effect on the surface air temperature in LGM-mixed.Replacing the ICE-5G topography with that of LGM-mPISM while keeping the glacier mask at ICE-5G values yielded an SAT increase by 0.87 K in years 60 to 90 compared with the same years of LGM-ICE-5G (not shown).This is about half of the 1.6 K temperature difference between LGM-ICE-5G and LGM-mPISM.Most of the northern Pacific warms.Over the western part, this is largely a downwind effect of higher temperatures over Asia, while over the eastern part, the picture is less clear with strong regional differences in the partitioning of the heat  (Tarasov and Peltier, 2003;Tarasov et al., 2012).(c) The surface topography from the ICE-5G reconstruction of Peltier (2004) and sea ice from LGM-ICE-5G.The corresponding plot for PI-mPISM is shown in Fig. 3a.fluxes (not shown).The higher topography over the Canadian Arctic Archipelago and Greenland cools these regions, and following from that the air over the Nordic seas and the Arctic Ocean.
The global mean precipitation in LGM-mPISM (Fig. 4) of 0.96 m yr −1 is 7 % lower than in PI-mPISM.This difference is especially pronounced in the regions north of 45 • N, where the precipitation is reduced by 18 % to 0.54 m yr −1 .The difference in precipitation between LGM-ICE-5G and PI-mPISM is larger (−11 % globally, and −31 % north of 45 • N) because of the lower temperature and the resulting lower atmospheric moisture content.

The LGM ocean
In LGM-mPISM, the NADW is formed southeast of Iceland (Fig. 5b).The NADW cell of the overturning circulation peaks at 22.1 Sv (PI-mPISM: 17.0 Sv) at 32.5 • N at a depth of 1020 m.The NADW cell reaches approximately 65 • N. In contrast to PI-mPISM, hardly any deep water is formed in the Nordic seas or in the Labrador Sea.This is indicated by the March mixed layer depth in Fig. 5b.Because of the lack of NADW formation in the Nordic seas, the NADW cell does not extend north of the Greenland-Scotland Ridge.The strengthening of the NADW cell is in agreement with the results of Lippold et al. (2012).Our model does not capture the shoaling of the AMOC generally found in proxies (Lynch-Stieglitz et al., 2007;Lippold et al., 2012).The AABW cell strengthens and peaks at 3.6 Sv at 21.5 • S and a depth of 3570 m.
In the LGM-ICE-5G experiment, the NADW cell is weaker than in LGM-mPISM.The stream function peaks at 18 Sv at 35.5 • N and a depth of 1020 m.The NADW cell is slightly stronger than in PI-mPISM, and, as in LGM-mPISM, it does not substantially extend beyond the Greenland-Scotland Ridge.In contrast to LGM-mPISM, NADW is formed in the Labrador Sea, and in the Nordic seas (Fig. 5), as in PI-mPISM.The deep water formation zone in the Nordic seas is smaller than in PI-mPISM, while the deep water formation zone in the Labrador Sea expands.Sensitivity experiments show that it is possible to switch between the LGM-mPISM and the LGM-ICE-5G deep water formation areas by swapping the ice sheets.This indicates that the shape of the ice sheets determines the pattern of deep water formation.A switch between the two circulation states could lead to a rapid change in the regional heat budgets and temperatures, thus providing a possible mechanism for glacial climate change events.Although the maximum of the overturning stream function stays at the same depth, the boundary between NADW and AABW shoals slightly.The AABW cell is slightly stronger (3.8 Sv) in LGM-ICE-5G than in LGM-mPISM.The North Atlantic Subtropical Gyre strengthens in LGM-mPISM compared to PI-mPISM while the Subpolar Gyre weakens (Fig. 9).Increased surface wind stress from the Westerlies (not shown) shifts the front between the two Gyres northward.The barotropic circulation in the Nordic seas is weaker in both LGM experiments (LGM-mPISM and LGM-ICE-5G) than in PI-mPISM.In contrast to LGM-mPISM, in LGM-ICE-5G, there is substantial deep water formation in the Labrador Sea and thus the Subpolar Gyre strengthens even compared to PI-mPISM (Fig. 9).In LGM-ICE-5G, the Antarctic Circumpolar The differences in ocean circulation can also be seen in the Atlantic ocean heat transport.In LGM-mPISM it reaches 1.1 PW at 23 • N (0.86 PW in PI-mPISM, 1.0 PW in LGM-ICE-5G).Between 15 • S and 35 • N, the overturning component dominates the heat transport and the differences between the experiments (Fig. 10).It is highest in LGM-mPISM, where the overturning is strongest.In both LGM experiments in the Atlantic south of 35 • N the northward Atlantic heat transport is stronger than in PI-mPISM.This explains the simulated low cooling or even slight warming in the North Atlantic.Outside of this latitude band, the gyre transports become more important.Between 40 and 60 • N, the strong Subpolar Gyre dominates the transports in PI-mPISM and LGM-ICE-5G, while in LGM-mPISM, the Subpolar Gyre is weak and the overturning contributes signifi- cantly.All experiments show similar total heat transports in the North Atlantic.
The sea ice in LGM-mPISM reaches further south in the North Atlantic than in PI-mPISM (see Figs. 3 and 8).It reaches Iceland during the entire year and covers large parts of the deep water formation areas of PI-mPISM.The winter sea ice margin reaches 63 • N east of Iceland and 43 • N at the American east coast.A part of the Labrador Sea becomes ice-free during summer.In LGM-ICE-5G, the winter sea ice is similar to that of LGM-mPISM.The summer sea ice margin shifts to the north, the Norwegian Sea becomes ice-free and the summer sea ice cover in the Labrador Sea decreases (Fig. 8).The reduced summer sea ice cover in LGM-ICE-5G is a combined effect of the missing latent heat flux from glacier melt (see below), and higher wind stress pushing more ice out of the Labrador Sea.
The heat fluxes from glacier calving and shelf basal melt contribute 30 % of the energy for the sea ice formation in the Arctic Ocean (north of the Fram Strait) in LGM-mPISM (Table 5).In LGM-mPISM-W, these heat fluxes are cut.Thus, the Arctic Ocean ice cover thins (−32 % ice thickness) and the heat loss to the atmosphere increases (+30 %) due to the reduction of the insulation effect of the sea ice.This leads to a winter (DJF) SAT warming of 2.4 K in the Arctic.This warming is almost entirely restricted to the near-surface layers, thus weakening the winter-time inversion.The ice export through the Fram Strait sinks by 16 %, the ice volume in the Nordic seas by 24 %, and its extent by 8 %.The summer sea ice cover in the Labrador Sea shrinks (orange outline in Fig. 8).These changes explain, for a part of the sea ice, differences between LGM-ICE-5G and LGM-mPISM, and show the importance of including the latent heat fluxes from calving in GCM simulations.A proper treatment of these fluxes would require the inclusion of an iceberg model.
Table 5.Heat budget of the Arctic Ocean sea ice in years 1650-2349 of LGM-mPISM and the same years of LGM-mPISM-W (the comparison experiment with cut ice shelf-ocean heat fluxes).The boundary between the Arctic Ocean and the Nordic seas is drawn at the Fram Strait.

The LGM ice sheets
The total modeled land ice volume in the Northern Hemisphere is 60 Mio km 3 , corresponding to 150 m of sea level change.Table 6 lists the ice sheet volumes, Fig. 6 shows the time evolution of the volumes, and the mask used for this analysis.We compare our results to four reconstructions.The widely used ICE-5G reconstruction (Fig. 8c, Peltier, 2004), its follow-up ICE-6G as provided by the PMIP-3 project (PMIP3 Project members, 2010), the latest reconstruction by Lev Tarasov (Fig. 8b, Tarasov and Peltier, 2003;Tarasov et al., 2012, labeled as Tarasov in the following) and the reconstruction by Kurt Lambeck as provided by the PMIP-3 project (PMIP3 Project members, 2010, ANU in the following).ICE-5G consists of a high-resolution bedrock topography combined with a coarse resolution ice thickness.The surface therefore is not smoothed out by the ice and much rougher than a real ice sheet surface.Local surface elevations can therefore not be considered for comparisons.We therefore refer to the Tarasov reconstruction for small-scale features.

The LGM Greenland Ice Sheet
Of the 60 Mio km 3 Northern Hemisphere ice volume, the Greenland Ice Sheet contains 5.8 Mio km 3 (14.5 m SLE), compared to 4.3 Mio km 3 (10.8m SLE) in the ICE-5G reconstruction (Peltier, 2004).It shows a drift of −2.7 km 3 yr −1 averaged over the full 30 kyr.The Greenland Ice Sheet's two-dome structure closely resembles that in PI-mPISM.It is wider in all directions than in PI-mPISM, and covers all the continental slopes.The height of the northern dome is at 3150 m (100 m lower than in PI-mPISM, Tarasov: 3112 m).
The southern dome profits from the cold climate and the wider ice sheet, and reaches a height of 2850 m (150 m higher than in PI-mPISM, Tarasov: 3083 m).The ICE-5G as well as the Tarasov reconstruction show a wider and flatter Greenland Ice Sheet during the LGM than the present Greenland Ice Sheet, but assume less ice at the margins than mPISM.In LGM-mPISM, the annual accumulation on the Greenland Ice Sheet of 558 Gt yr −1 is exclusively balanced by direct losses to the ocean.Melting is negligible because of the lower surface temperatures.Figure 1 provides an overview of the ice streams.In the following, capital letters refer to the figure.The Northeast Greenland Ice Stream (A) surges repeatedly.Its southern neighbor in a trough off Kejser Franz Joseph Fjord (B) shows similar surge-type behavior and matches deposits described by Wilken and Mienert (2006).South of it, at Scoresby Sund (C), an ice stream is continuously active matching glacial deposits found in this area (Wilken and Mienert, 2006).In the area of Kangerdlugssuaq Glacier (D) an ice stream continuously drains into the Atlantic.Further continuously active glaciers follow on the southern east coast.On the west coast, ice streams surge repeatedly in the region of Jakobshaven Isbrae (E) and the trough north of it, matching proxy observations in Roberts et al. (2009).Further north, in the region of Kong Oscar Glacier (F), two tributaries show surge-type activity.Since the sediment mask we use in mPISM (see Sect. 2.2) does not allow for sliding in the interior of Greenland (blue areas in Fig. 1), the ice streams are limited to the continental shelf.This does not prohibit fast-flowing ice as in PI-mPISM (Fig. 7), where fast ice flow occurs in parts of the Greenland Ice Sheet where the sediment mask prohibits sliding.This fast flow is caused entirely by internal deformation.The inclusion of temperate ice in PISM allows for a very low viscosity, so the ice can reach high speeds by pure internal deformation.

The LGM Laurentide Ice Sheet
Baffin Island and Ellesmere Island are fully glaciated in LGM-mPISM, and connect the Greenland Ice Sheet to the Laurentide Ice Sheet.The Laurentide Ice Sheet covers present-day Canada, and has a mean volume of 31 Mio km 3 (ICE-5G: 36 Mio km 3 ), corresponding to 78 m of sea level (ICE-5G: 90.7 m SLE).In the west, it terminates inland of the coast, while in the east and north, it fully covers the coasts.
The southern boundary is approximately at 50 • N in the west and at 45 • N in the east.
The Laurentide Ice Sheet is split into a main part and a western Cordilleran part by the Mackenzie Ice Stream.This ice stream cuts down to below 1500 m a.s.l. and is in continuous operation with an average strength of 694 Gt yr −1 (21 mSv water equivalent).During most of the experiment, Mackenzie Ice Stream shows net surface melt in its northern part because of a foehn-effect acting on the winds from the Pacific.This area is characterized by very low surface elevations.The main part of the Laurentide Ice Sheet has two domes that are separated by the Hudson Bay.The maximum height of the eastern dome is 3200 m (3600 m of ice thickness), the maximum height of the western dome is 3150 m (also 3600 m of ice thickness).The Hudson Bay area is largely drained by the Hudson Strait Ice Stream that approximately every 7000 yr flushes ice into the Labrador Sea.Unless otherwise noted, we average over this oscillation in this chapter.Details of the oscillation and its implications in the climate system will be covered in a follow-up publication.They basically follow the binge-purge oscillator mechanism proposed by MacAyeal (1993) and first shown to work in 3-D by Calov et al. (2002).The main part of the Laurentide Ice Sheet loses ice by surface melt at its southern boundary and by calving into the ocean at the eastern and northern boundaries.
The Cordilleran part of the Laurentide Ice Sheet reaches heights of up to 3550 m, but these elevations are reached only on high mountains, so the thickness there is about 1000 m.In valleys, the thickness reaches up to 2700 m, with surface heights of up to 2800 m.The surface accumulation of the Cordilleran part of the Laurentide Ice Sheet is balanced by the Mackenzie Ice Stream on the eastern side and by surface melt on the western side.The melt on the western side is possible because of the high temperatures in the northern Pacific that locally exceed those in PI-mPISM.The modeled ice-free western coast could be a resolution artifact.The coarse resolution of the ISM could lead to an underestimation of the extent of the ice streams.The coarse resolution of the AGCM is responsible for an underestimation of the strength of the precipitation on the slope.On the other hand, the AGCM overestimates the width of this belt of strong precipitation.The Cordilleran part of the Laurentide Ice Sheet covers Alaska fully.This is likely related to the same biases that cause the glaciation of the Alaska Range in PI-mPISM.
None of the reconstructions show an ice-free American west coast.They differ in the southern boundary and in the details of the structure of the interior of the ice sheet.All four reconstructions agree with our model in putting the southwestern edge further to the north than the eastern edge, with values for the southwestern edge between 42 • N (ANU and ICE-6G) and 50 • N (Tarasov) (50 • N in our model).For the southeastern edge the values are between 35 • N (ANU and ICE-6G) and 40 • N (Tarasov and ICE-5G), while our model yields 45 • N. Considering that our climate model has a reso-lution of about 3.75 • , this match is acceptable.In the ANU reconstruction and in our model the Mackenzie Ice Stream splits the Laurentide Ice Sheet into a western and an eastern part.In ICE-5G and in the Tarasov reconstruction this ice stream is less well represented.It is hardly discernible in ICE-6G.In the central part, the ANU reconstruction shows a higher surface elevation than our model, but the structure is very similar with a rather low surface elevation in the region of the Hudson Bay and higher surface elevation south and west of it.A similar structure can be seen in the Tarasov reconstruction.West of the Hudson bay, ICE-5G, shows a massive mountain range between 90 • W and 120 • W, reaching about 4500 m a.s.l.while ICE-6G has a peak in the Hudson Bay area.
A comparison of the ice streams simulated in the model (Fig. 1) with those found in proxy records (Stokes and Tarasov, 2010) shows several ice streams, where models and reconstructions agree.In the following, numbers relate to the numbering in Stokes and Tarasov (2010) and Fig. 1.The central Laurentide Ice Sheet is drained into the Arctic Ocean by Mackenzie Ice Stream (1).There are two major ice streams in the Canadian Arctic Archipelago, the Amundsen Gulf Ice Stream (18) just to the east of Mackenzie, and M'Clure Strait (19) north of Amundsen Gulf Ice Stream.Both show surge behavior.So does Lancaster Sound Ice Stream ( 22) with its tributaries, the Admirality Inlet ( 21) and the Gulf of Boothia ice stream (20).They drain the north-eastern part of the Laurentide Ice Sheet into Baffin Bay.Further to the south, hardly represented, the Cumberland Sound Ice Stream (23) surges into the Davis Strait.South of Cumberland Sound and well represented, the Hudson Strait Ice Stream (24) drains the Hudson Bay into the Labrador Sea.The Hudson Strait is not the only possible ice stream route for draining the Hudson Bay.The sediment distribution allows for a more northerly route joining the Lancaster Sound Ice Stream and draining into the northern corner of Baffin Bay.However, this route does not become active in our experiments (Fig. 1).A repeatedly surging ice stream drains the Ungava Bay (16) into the Hudson Strait.In the Gulf of St Lawrence, a large ice stream system forms in the Laurentian Channel (25) and neighboring tributaries.

The LGM Eurasian ice sheets
The Laurentide Ice Sheet is connected via Alaska to an ice sheet in eastern Siberia that does not exist in the reconstructions.There is however evidence of Pleistocene glaciation of the East Siberian continental margin (Niessen et al., 2013), and of late Pleistocene large-scale glaciation at the Siberian Pacific coast (Bigg et al., 2008;Barr and Clark, 2012)  N. The eastern part starts with a peak height of 2600 m.During the experiment, the ice sheet expands southward, and the peak shifts to the south and grows to 3000 m.The western part starts with a peak height of 2730 m, decreases to 2600 m during the first 10 000 yr, and then stabilizes.Along the southern border and parts of the Norwegian Sea coast, there is surface melt.At the Arctic Ocean coast, all losses occur directly into the ocean.
Among the reconstructions of the Fennoscandian Ice Sheet, ICE-5G and ICE-6G show the largest low-thickness zones.Such zones are hard to obtain as a steady-state solution in a dynamical ISM, where there are positive feedbacks for ice sheet growth, until either height desertification or a nearby coast limit the ice sheet height.They are easier to obtain as a transient state.The closure of the gap between the Fennoscandian Ice Sheet and the East Siberian Ice Sheet at the end of the coupled experiment shows such a large, flat zone that is growing by surface accumulation.ICE-5G portrays the Fennoscandian Ice Sheet as reaching far to the south and staying below 1000 m in its southern parts.In the ANU reconstruction, the region between 50 and 60 • N is covered with substantially thicker ice exceeding 1500 m in large parts and even exceeding 2500 m over Norway.The surface elevations in the Tarasov reconstruction are lower over Norway and the Barents Sea than in the ANU reconstruction, but the reconstructions largely agree.Over the Barents, Kara, and Laptev seas, our model places much more ice than any of the reconstructions.There are massive ice streams between Norway and Svalbard (α), and further streams between the present-day islands at the northern margin of the ice sheet (β, γ ).They match with proxy records (Denton and Hughes, 1981).The southern margin ice streams cannot be compared to the reconstructions, since the margin is too far to the north.
Iceland is covered by an ice cap with a volume of 278 000 km 3 (0.7 m SLE, ICE-5G: 172 000 km 3 , or 0.43 m SLE) and a maximum height of 2450 m (1400 m ice thickness).

Long-term ice sheet changes
Figure 6 shows the evolution of the ice sheet volumes in LGM-mPISM.The largest changes occur on the Bering Sea shelf and between the East Siberian and Fennoscandian ice sheets.The Bering Sea shelf is flooded with ice from Alaska between years 9000 and 11 000.The ice stream surges transport vast amounts of ice into the region (Fig. 1) and, in the first years, have to compensate for strong surface melt.Over time, the ice sheet stabilizes.
The eastern part of the Fennoscandian Ice Sheet slowly expands southward.This allows the ridge to shift southward and increase in altitude.For the first 20 000 ice model years, the snow in the region between the Fennoscandian and the Siberian ice sheets melts during the summer, except for a few cold years, when it survives the summer melt.The ice sheets slowly grow into this area by lateral ice advection and start closing the gap from the sides.During the last 10 kyr, the winter snow in the gap between the ice sheets survives the summer melt, and the gap between the ice sheets is closed by glacier growth from local accumulation.
The simulated extensive glaciation of northern Siberia is likely to be an artifact of prescribing constant LGM boundary conditions.In the coupled ice-climate simulations of Heinemann et al. (2014), a transient simulation starting at 80 kyr BP shows a rather small LGM Fennoscandian Ice Sheet, while a steady-state simulation branched off from this LGM state grows a massive ice sheet covering most of Siberia.The steady-state Laurentide and Greenland ice sheets were very similar to the transient versions.This shows that the Fennoscandian Ice Sheet is far from equilibrium at the LGM.Another factor that prevented the growth of an ice sheet in Siberia is the reduction in the snow surface albedo by increased dust deposition during the last glacial (Krinner et al., 2011).This effect led to increased temperatures and melt in Siberia.It is likely that the model's cold bias over northern Siberia for the present day could lead to a similar cold bias in the LGM simulations.This could additionally contribute to the growth of the Siberian Ice Sheet.

Summary and conclusions
We studied the glacial climate with an interactively coupled AOVGCM-ISM system.Such systems are only beginning to be developed.Therefore, we have modified the state-ofthe-art ISM PISM into mPISM, a model that can be used in coupled ice sheet-climate simulations, and coupled it to the AOVGCM ECHAM5/ MPIOM/ LPJ.Both models, as well as the coupling, work without anomaly maps or flux correction.This avoids inconsistencies between flux correction terms and modeled climate changes.In comparison to simulations using EMICs, AOGCMs represent processes of the ocean circulation and atmosphere dynamics in a much more detailed way and with higher spatial and temporal resolution.In contrast to previous AOGCM-ISM simulations (e.g., Gregory et al., 2012), the ISM covers all relevant parts of the Northern Hemisphere.mPISM is a SIA-SSA hybrid model, and is thus able to model ice streams more realistically than conventional SIA-only ISMs.The ISM is bidirectionally coupled to the atmosphere as well as to the ocean model, enabling the study of the full interactions between the ice sheets and the climate system.
We validated our setup by performing steady-state experiments under pre-industrial boundary conditions (PI-mPISM).The results agree reasonably well with the observational data.The global mean SAT in PI-mPISM is below that of ERA-INTERIM, representing the lower preindustrial greenhouse gas concentrations.The NADW cell strength agrees with the estimates obtained from observations.The NADW is formed in the Nordic and Labrador seas.In PI-mPISM, an ice sheet in the Alaska Range forms.This is due to a resolution-dependent cold bias in the atmosphere model in an area characterized by individual glaciers and ice fields.
With the same setup, we performed the first fully coupled multi-millennial steady-state AOVGCM-ISM simulations under LGM boundary conditions (LGM-mPISM).Again, the results agree reasonably well with proxy data.The NADW formation shifts to southeast of Iceland.The heat fluxes from ice shelf basal melt and calving contribute 30 % of the cooling of the Arctic Ocean.Cutting them leads to thinning of the sea ice cover and an increase in the ocean-atmosphere heat flux, as well as to a reduction in the summer sea ice cover in the Nordic and Labrador seas.
During the long steady-state LGM simulations, a spurious ice sheet forms in eastern Siberia and Alaska.This is at least partly due to neglecting the albedo effect of dust on snow and ice in our model, which would increase surface ablation in this region and probably prevent ice sheet growth (Warren and Wiscombe, 1980;Krinner et al., 2011).Further advances could be made by using a sophisticated energy balance scheme for the surface mass balance (e.g., Calov et al., 2005;Vizcaíno, 2006) and a higher model resolution that can resolve the small-scale features of the glaciation in these regions.Finally, LGM-mPISM is a multi-millennial integration under constant LGM boundary conditions, while the LGM in reality was a transient state, where the ice sheets and the climate were not in an equilibrium.The ice sheet in Siberia might in part simply be the result of running the model too long under LGM boundary conditions.Modeling the last glacial as a transient process in AOGCMs is a major challenge for the years to come.
When the model is forced with the ICE-5G ice sheet reconstruction, the LGM cooling is stronger than in LGM-mPISM.In contrast to LGM-mPISM, the NADW is largely formed in the Labrador Sea.It is possible to switch between the deep water formation regions of LGM-mPISM and LGM-ICE-5G by exchanging the ice sheets.This provides a mechanism for obtaining two different ocean circulation states in glacial climate simulations.
mPISM shows strong surge behavior in the Hudson Strait, as well as in several other regions.These surges follow the Heinrich event mechanism described by MacAyeal (1993) and first modeled in 3-D by Calov et al. (2002).The response of the climate system shows the basic features of Heinrich events (Heinrich, 1988;Clement and Peterson, 2008).We currently study the effects of a transient spinup and the processes relevant for the last deglaciation with the coupled model system.

Figure 6 .
Figure 6.Development of the ice sheet volumes in LGM-mPISM.The inset shows the split of the ice sheets into different regions for the diagnosis.

FFigure 7 .
Figure 7. Greenland Ice Sheet surface velocities.Plotted on a logarithmic color scale.Contour lines show the surface elevation in steps of 500 m with thick lines at multiples of 1000 m.(a) Modeled velocities from PI-mPISM; (b) observed velocities(Joughin et al., 2010b) and smoothed-out ETOPO1 topography(Amante and Eakins, 2009).White areas inside the ice sheet indicate data gaps.

Figure 8 .
Figure 8. Surface topography and sea ice.Isolines mark the topography at 500 m intervals.Ice covered regions are colored.Over the ocean, dark gray with red outline indicates areas with perennial ice cover (15 % level) in more than 50 % of the model years.Light gray with dark blue-black outline indicates temporary ice cover in more than 50 % of the model years.(a) The surface topography in mPISM averaged over the full 30 kyr of LGM-mPISM, orange lines in the ocean mark perennial ice cover in LGM-mPISM-W, light blue lines mark temporary ice cover in LGM-mPISM-W (b) The LGM topography reconstruction provided by Lev Tarasov(Tarasov and Peltier, 2003;Tarasov et al., 2012).(c) The surface topography from the ICE-5G reconstruction ofPeltier (2004) and sea ice from LGM-ICE-5G.The corresponding plot for PI-mPISM is shown in Fig.3a.

Figure 10 .
Figure 10.Atlantic Ocean heat transports.Solid lines mark the total heat transport, dashed lines the gyre contribution.The difference is the overturning contribution.

Table 1 .
Main experiments performed and analyzed.For the boundary conditions, see Table 2. 1 : 10 coupling means 1 climate model (CM) year per 10 ISM years.The coupling is performed after each year of climate model integrations.For easy reference, ISM years in the coupled experiments are ten times the corresponding climate model years, thus climate model year 10 corresponds to ISM years 100 to 109.

Table 2 .
Boundary conditions differing between the LGM and preindustrial setups.

Table 3 .
Greenland surface mass balance in Gt yr −1 .RACMO2, Polar MM5 (PMM5) Ganachaud and Wunsch (2003)nt S1 fromEttema et al. (2009).AABW) cell in the Atlantic peaks at 2.9 Sv at 30 • S and a depth of 3570 m.The northward heat transport in the Atlantic of 0.86 PW at 23 • N is lower than the present-day estimates.Ganachaud and Wunsch (2003)obtain 1.27 ± 0.15 PW at 24 • N as a result of the World Ocean Circulation Experiment (WOCE), and Johns et al. (2011) obtain 1.33 ± 0.14 PW at 26 • N from the RAPID mooring array measurements.A comparison of the modeled temperatures in the Atlantic and data from WOCE

Table 4 .
Lemke et al. (2007))d their changes in PI-mPISM averaged over the last 1000 yr of the simulation.Greenland volume:Bamber et al. (2001).Glaciers outside of Greenland and Antarctica are estimated between 0.05 and 0.13 Mio km 3 inLemke et al. (2007).

Clim. Past, 10, 1817-1836, 2014 www
Krinner et al. (2011)berian Ice Sheet has a maximum height of 3300 m and a volume of 9.3 Mio km 3 , corresponding to 24 m SLE.The drift is +108 km 3 yr −1 .In the model, the East Siberian Ice Sheet closes the gap between the Laurentide Ice Sheet in the East, and the Fennoscandian Ice Sheet in the West.Krinner et al. (2011)concluded that two important factors for not glaciating eastern Siberia are (1) the low snow albedo that is caused by dust deposition, and (2) moisture blocking by the Fennoscandian Ice Sheet.We do not use a locally varying glacier albedo, so the first effect is not represented in our setup.The modeled Fennoscandian Ice Sheet does not reach as far south as indicated by the reconstructions and the coarse resolution of the Atmosphere model does not allow for a realistic simulation of the moisture blocking.Furthermore, we run the model under LGM boundary conditions for a long time span, so our (steady-state) response must be expected to be different from a transient state (seeSect.3.7.4).The East Siberian Ice Sheet loses mass by surface melt along all margins except for the Arctic Ocean coast, where the losses occur by calving and shelf basal melt only.Further calving and shelf basal melt occur at the Pacific coast.The Fennoscandian Ice Sheet has a volume of 11.6 Mio km 3 (29.2m SLE; ICE-5G: 8.2 Mio km 3 , 20.7 m SLE) and shows a drift of +156 km 3 yr −1 .It consists of two main parts.One part covers the Barents, Kara and Laptev sea shelves and the islands in this region up to Svalbard in the northwestern corner.The other part covers Scandinavia south to 60 .clim-past.net/10/1817/2014/