Interactive comment on “ A new coupled ice sheet-climate model : description and sensitivity to model physics under Eemian , Last Glacial Maximum , late Holocene and modern climate conditions

This paper describes a significant advance in the coupled modeling of ice sheets and climate. The authors have coupled the Penn State University Ice model to the U.Vic. Earth System Climate Model and have carried out detailed control and sensitivity simulations of both the Greenland and Antarctic ice sheets for a variety of climate conditions (Eemian, LGM, late Holocene, and present-day). To my knowledge, these are the


Introduction
Ice sheets are important components of the climate system, but are among the least understood (Lemke et al., 2007).Dramatic changes to the global climate in the past are likely due in part to the evolution of large ice sheets in response to changing external forcing and internal ice sheet dynamics.Ice sheets may also exert a large influence on the climate system in the coming millennia, a time when large changes to the planet are expected in response to human forcing of climate.Recent research indicates that present and projected anthropogenic carbon emissions and corresponding temperature will be retained by the surface-Earth carbon system over millennial timescales (Montenegro et al., 2007;Eby et al., 2009).Such timescales are relevant to ice sheet dynamics and evolution.It is therefore important that the effect of such emissions on the cryosphere be explored, before thresholds on ice sheet stability, based on increases to atmospheric or ocean temperatures, are surpassed and the Earth becomes irreversibly committed to near-permanent ice sheet decay or loss and associated sea level rise (Ridley et al., 2010).A coupled Earth system modelling approach is a powerful way to explore interactions between ice sheets and climate.Coupled Earth system models incorporate realistic subcomponent models (e.g., atmosphere, ocean, or land surface) and interactive communication of fluxes between subcomponents in such a way that allows feedbacks within the climate system to be resolved.To date only a few fully coupled Earth system models with the ability to realistically simulate climate change over millennial timescales have been synchronously coupled to high resolution, dynamic ice sheets (e.g., Kageyama et al., 2004;Ridley et al., 2005;Driesschaert et al., 2007;Vizcaíno et al., 2008;Charbit et al., 2008;Vizcaíno et al., 2010;Pritchard et al., 2008a).Furthermore, to date no coupled ice-climate models have explicitly represented the large Antarctic ice shelves.
A new ice-sheet/ice-shelf/climate model consisting of the University of Victoria Earth System Climate Model (UVic ESCM) and the Pennsylvania State University ice sheet model (PSUI) is presented here.It is designed to efficiently simulate the ice sheet/climate system over 1-10 kyr timescales.First, the model and its components are described.Then, the performance of the model under equilibrium Eemian (ca.129 kyr BP), LGM (ca.21 kyr BP), late Holocene (ca.0.2 kyr BP), and transient modern climate conditions is evaluated in response to simple perturbations to model physics, including the treatment of climate model surface air temperature (SAT) bias, surface albedo, and refreezing.Finally, a transient simulation that models present-day SMB conditions is described.

Model description
Several major design requirements dictated the choice of both climate and ice sheet models, as well as the design of the coupling software, and motivated the choice of the UVic ESCM, version 2.9 and PSUI as suitable candidates for coupling.An earlier version of the UVic model had previously been coupled to an ice sheet model (Marshall and Clarke, 1997) and used to study the role of ice sheets during deglacial periods (Yoshimori et al., 2001;Schmittner et al., 2002).However, the representation of ice in these studies was simplified, in part due to very low ice sheet resolution, isothermal ice and an inability to resolve floating ice shelves.Given the recent advances in ice sheet modeling we felt it prudent to update the ice sheet component of the UVic ESCM to the full PSUI model.

UVic ESCM
UVic ESCM is an Earth System model of intermediate complexity that includes fully coupled atmospheric, oceanic and land surface components, a thermodynamic/dynamic sea ice model and a closed and coupled complex carbon cycle model.The ocean component is the full 3-D ocean general circulation model MOM2 (Pacanowski, 1995).The atmospheric component consists of a vertically integrated energymoisture balance model in which the radiative forcing associated with an increase in atmospheric CO 2 is applied as a decrease in outgoing long wave radiation.Heat and moisture are diffused and advected throughout the atmosphere, with winds being provided by present-day NCEP climatological long-term monthly means (Kalnay et al., 1996).Dynamic modification to the prescribed wind fields is obtained by a routine which generates anomalies to the monthly mean winds as a function of surface air temperature (Weaver et al., 2001).The land surface component of the model includes the MOSES-TRIFFID global vegetation model (Meissner et al., 2003).The model has been extensively used for a large number of studies of paleoclimate and future climate projections and is currently under active development and use.The UVic ESCM was recently run at ocean eddy-resolving resolution (Spence et al., 2008).While regional differences exist between this version of the UVic ESCM and the standard lower resolution version, the overall response is similar between the two.We therefore retain the standard UVic ESCM resolution of 3.6 • in longitude by 1.8 • , as this allows millennial-scale simulations to be performed on existing computers.

PSUI
PSUI is a hybrid ice-sheet/ice-shelf model that simulates grounded and floating ice.Over inland ice on stiff bedrock, the shallow-ice approximation applies.Over ice with reduced basal traction (for example, ice shelves and regions of the West Antarctic Ice Sheet), the shallow shelf equations are employed.The transition between the two is dealt with heuristically by iteratively blending the velocities of both solutions.Additionally, a grounding-line ice velocity is imposed, and shelf back-stress is included, based on the analytic treatment of Schoof (2007).Three-dimensional advection and vertical diffusion of heat, and lithospheric flexure/local bedrock isostasy are calculated.Detailed descriptions of the model are in the Supplementary Material of Pollard and De-Conto (2009), and Pollard and DeConto (2007).

Energy-moisture balance model (EMBM)
The ice sheet surface mass balance (SMB), consisting of melt, runoff, sublimation and accumulation, is calculated using a modified version of the existing UVic ESCM vertically integrated energy-moisture balance atmosphere model (EMBM) and simple snow model.We were motivated to use the existing EMBM, instead of a positive-degree-day (PDD) model, to calculate the SMB for reasons that have been recently described (e.g., Bougamont et al., 2005).PDD approaches rely solely on surface air temperature (SAT) to derive snow or ice melt rates.The PDD factor must therefore be tuned to capture the effect of changing albedo and related non-sensible surface energy fluxes.However, this implicit adjustment may break down in scenarios where, for example, insolation remains constant but SAT increases due to higher concentrations of atmospheric greenhouse gases, or ice sheet elevation changes dramatically.
Inputs to the energy balance calculation are incoming shortwave radiation, specific humidity, wind speed, surface skin temperature, albedo and roughness, and surface air temperature (elevated from sea level using a constant lapse rate, nominally 5.0 • C km −1 ).This lapse rate is comparable to the summer lapse rate over the Greenland Ice Sheet (GIS) (Fausto et al., 2009), and like Stone et al. (2010) we suggest that this low lapse rate is a reasonable choice, as it reproduces summer ablation reasonably without significantly affecting wintertime accumulation.The snow temperature is found by solving the energy balance equation: for temperature, using Newton's method.T is temperature, and Q dswr , Q dlwr , Q lt , Q s and Q ulwr are incoming shortwave, incoming longwave (both supplied by the UVic ESCM), outgoing latent, outgoing sensible and outgoing longwave heat fluxes, respectively.Snow and ice emissivity is set to 0.94.During the iteration, the surface equilibrium vapor pressure is determined and compared to the ambient specific humidity.If the equilibrium specific humidity is greater, solid-tovapor wind-enhanced sublimation occurs to equalize vapor pressures and the outgoing latent heat and moisture fluxes are determined.If the surface temperature is at the melting point after the iteration and a net flux of energy into the surface exists, it is entirely used to melt snow/ice.Sublimation of scoured, airborne snow (Box et al., 2006;Lenaerts et al., 2010) is not explicitly modelled, although wind speed is used in the calculation of the sublimation-derived moisture flux and therefore has a direct effect on surface sublimation rates.Representation of snow and runoff used here is relatively simple, in line with the intermediate-complexity character of the UVic ESCM.An arbitrary maximum snow thickness of 10 m is prescribed, above which any additional accumulation is added to the ice sheet as a positive SMB.Conversely, the entire 10 m snowpack must melt before ice ablation can occur.The snow model is greatly simplified by assuming that snow/ice heat capacity c p = 0 (i.e., the snowpack temperature is always equal to the surface temperature) for the sake of energy balance calculations.This assumption leverages the fact that heat fluxes used to change snow and ice temperatures are an order of magnitude less than the the heat fluxes required for phase changes and can be neglected.This simplification is mirrored in the approach of Robinson et al. (2009).Pfeffer et al. (1991), Bougamont et al. (2007) and others have recognized the importance of refreezing in limiting runoff from ice sheets.In order to capture the effect of refreezing within the relatively simple snow model, the parameterization described in Pfeffer et al. (1991) for representing refreezing of meltwater is adopted.An annual snow thickness (i.e., between last fall melt and first spring melt) is first accumulated.When melt is first detected, the thickness of meltwater that is stored in the snowpack is calculated based on the requirements that the annual snowpack average temperature be raised to freezing point and all the available pore space within the snowpack be filled.The scheme effectively assumes that an impermeable layer exists between annual layers in the snowpack and that the annual snow layer has a constant density (set to 330 kg m −3 ).
Over ice sheets, surface albedo varies as a function of snow/ice thickness and snow/ice temperature.Default albedo values for cold/melting snow albedo are 0.8/0.6, with a linear interpolation between cold and warm snow values occurring between −5 • C and 0 • C. Glacial ice albedo is set to 0.45.At snow thicknesses below 10 cm, snow and ice albedos are linearly blended.These values are very similar to those of Robinson et al. (2009), though they were obtained independently by setting albedo values to approximate the glacier zone-specific values of Nolin and Payne (2007).The current model does not yet account for age-related changes in snow albedo (e.g., van de Wal and Oerlemans, 1994).
Moisture is delivered to the ice sheets via advective/diffusive transport within the atmospheric model.Precipitation occurs when the relative humidity at the surface surpasses 0.85.This relative humidity is calculated as a function of the subgridded elevation-dependent SAT and the specific humidity of the overlying atmospheric grid cell.Precipitation falls as rain/snow when the sub-gridded temperature is above/below a threshold temperature, set over ice sheets to −2 • C which is approximately the midrange of the observationally-derived snow-to-rain temperature parameterization of Robinson et al. (2009).

Subgridded elevation-binned mass balance calculations
Realistic transfer of climate-model-derived SMB to high resolution ice sheet grids is a non-trivial exercise due to spatial and temporal differences in scale between ice sheet and climate models (Pollard, 2010).Ice sheet models require grids on the order of 10 0 −10 1 km to appropriately resolve ice flow and gradients in SMB (particularly ablation).However, due to the slow flow of ice, time steps can be long, on the order of months to years.Climate models, on the other hand, typically operate on grids with spatial scales on the order 10 2 km and temporal scales of minutes to days.In order to accommodate the variation in scales while maintaining computational efficiency and conservation of mass and energy, a routine was developed to calculate SMB based on the distribution of elevation within each climate grid cell.Instead of one EMBM calculation per climate grid cell, the scheme allocates N EMBM calculations: where R is the total relief (defined as the difference between the maximum and minimum ice sheet surface elevations falling within the larger UVic ESCM grid cell, Fig. 1) and D is the binning threshold of the elevations bins, set to 100 m.Within each bin, the elevation is set to the average area-weighted elevation of the ice sheet grid cells that fall within the bin elevation range, and for each bin the EMBM algorithm evolves independent albedos, Q lt (T ), Q s (T ), Q ulwr (T ), snow thicknesses, and SMB values.When remapped back onto the ice sheet grid, the resolution of ablation zones and accumulation over short-wavelength, highrelief topography is greatly improved, in agreement with the similar approach and results of Wild et al. (2003).Over topography with wavelengths greater than the overlying climate model grid resolution (e.g., the continental East Antarctic Ice Sheet, EAIS) the subgrid scheme-derived accumulation field rapidly approaches the equivalent non-subgridded distribution.
Prior to each atmospheric time step, fluxes from all subgrid ice sheet bins within a climate model cell are accumulated, after area-weighting, and passed to the atmospheric model.This implies that the atmosphere is well mixed within each climate model cell, which was the implicit assumption before subgrid binning was introduced.At initialization and after each ice sheet timestep, the elevation bins are re-calculated according to the updated ice sheet geometry and the existing snowpack is redistributed into the new elevation bins.Where the ice dynamically retreats or melts over land or ocean any residual snow, if it exists, is sent to the ocean.Where the ice sheet dynamically expands, snow from the overridden bare land is transferred to the ice sheet surface.Finally, the updated fractional ice coverage is relayed to the overlying climate model grid.Model simulations with and without the sub-gridded elevation binning scheme indicate that the subgrid-enabled scheme is much better suited to resolve narrow ablation zones (limited in width only by the underlying ice sheet grid resolution) and orographically driven precipitation over short wavelength topography than interpolationbased SMB downscaling procedures.Furthermore, it retains computational efficiency by automatically focussing computing resources over high-relief regions, not unlike the manual procedure of van de Wal and Oerlemans (1994).For the modern ice sheet distribution, the number of EMBM calculations is reduced by a factor of 14, compared to the same number of equivalent calculations over each high-resolution ice sheet grid-point.The resulting system is designed to conserve mass and energy to within machine precision, in accordance with the rest of the UVic ESCM.

SAT bias correction
Surface mass balance and evolution of ice sheet geometry are highly sensitive to variations in SAT that are of the same magnitude as SAT biases in present-generation climate models (van de Wal and Oerlemans, 1994;Pollard, 2000).For this reason, previous coupled ice sheet-climate model studies (e.g., Ridley et al., 2005Ridley et al., , 2010) ) have utilized approaches in which modelled temperature anomalies, obtained by subtracting a perturbed model state from the control value, are added to observational climatologies prior to processing by PDD models.This method reduces spurious SMB values resulting from persistent and intrinsic modelled temperature biases but is handicapped by drawbacks of the PDD model approach (Sect.2.3.1).Other studies (Bonelli et al., 2009;Pritchard et al., 2008a;Charbit et al., 2008) also utilized the PDD approach, but with "uncorrected" climate model SAT that still retained a climate model SAT bias.Recently, Calov et al. (2009) and Vizcaíno et al. (2010) employed surface energy balance models to calculate SMB over the Northern Hemisphere and Greenland, respectively, without SAT corrections.They suggested that their annual present-day modelled climatology was sufficiently similar to observations to justify use of raw model output to supply their SMB model.However, visual inspection of Fig. 1 in Vizcaíno et al. (2010) suggests that differences of ∼5-10 • C still exist in their annually averaged surface air temperature and Calov et al. (2009) do not provide an analagous difference figure of SAT.
In order to use a surface energy balance model to generate SMB, but also analyze and optionally correct the effect of climate model-derived SAT biases on ice sheet evolution, all without losing global conservation of energy, the following bias correction procedure for use in the EMBM was developed.
1.A UVic ESCM simulation forced with late Holocene orbital and CO 2 conditions was run to equilibrium.Ice sheet model elevations were held constant at presentday values.Equilibrium late Holocene monthly-mean raw model SAT fields were extracted.
2. The model was then forced with transient CO 2 concentrations to year 2000, and the monthly-mean temperature difference between year 2000 and the late Holocene steady state was obtained.
3. This difference was subtracted from the 1971-2000 ERA40 (Uppala et al., 2006) monthly-mean SAT to derive an observed SAT, extrapolated to late Holocene conditions.
4. Another late Holocene UVic ESCM simulation was carried out, in which the model-derived SAT values were replaced with the late Holocene "observed" SAT values used in the EMBM calculations over the ice sheets (the SAT used in non-ice-sheet EMBM calculations remained the value provided by the atmospheric component of the UVic ESCM).This procedure effectively forced the energy balance model over ice sheet domains to respond to SAT anomalies and generated a regionally altered pattern of atmospheric heat transport.This altered transport field perturbed the global modelled SAT field in bias-corrected simulations, but the magnitude of The effect of the bias correction procedure is largest over the AIS and GIS and diminishes with distance from the correction regions.The global effect of the bias correction over both ice sheets is to redistribute heat from the Southern Hemisphere to the Northern Hemisphere.However, the resulting temperature change is less than 0.25 • C everywhere and does not significantly affect the non-ice-sheet components of the system.
the perturbation was less than 0.25 • C anywhere (Fig. 2) and lessened with distance from the ice sheets.
5. The monthly-mean SAT difference between the two simulations was obtained.These monthly-mean bias maps were added as a correction factor to the UVic ESCM atmospheric temperature used in the calculation of the surface energy balance and the rain/snow decision.
Results of the bias-corrected models are described in Sect.3.2.

Ice-land coupling
As the simulated grounded ice sheets advance or retreat they cover or expose bare land.The dynamic vegetation model MOSES-TRIFFID is currently excluded from any potential ice sheet region (i.e., the island of Greenland and the Antarctic continent) in the present version of the model.Instead, the albedo and roughness lengths of bare land areas in Greenland and Antarctica are set to values typical of tundra (when free of snow).For GIS and AIS simulations initialized with existing ice that do not undergo drastic retreat this is likely a reasonable approach.Surface meltwater is routed instantaneously via large-scale drainage basins to the ocean model.The AIS drainage basin distribution (Fig. 3) is based roughly on present-day ice surface topography and would become unrealistic only under a drastic change to ice geometry or significant continentalscale deviation in subglacial drainage from the surface topography.The GIS was assigned one drainage basin, resulting in an even distribution of meltwater flux around the GIS coast.
The elevation field of climate model grid cells that overlay ice sheet model instances evolve according to the gridaveraged ice sheet or exposed bedrock elevations calculated by the ice sheet lithosphere model.Elevations remain static over climate grid cells that do not overlay ice sheet grids.

Ice-ocean coupling
PSUI explicitly simulates ice shelves and a dynamic grounding line.Ice shelves are represented in the UVic ESCM by fractional surface coverage of ocean grid cells.Increasing fractional ice shelf coverage progressively "shades" underlying ocean and sea ice grid cells from atmospheric heat and moisture fluxes until the climate grid cell is fully covered, at which point the ocean and sea ice models are completely insulated from all atmospheric heat and moisture fluxes.The treatment of momentum fluxes is slightly different.If atmospheric momentum fluxes are simply shaded, calculation of sea ice dynamics within the dynamic/thermodynamic sea ice component of the UVic ESCM results in artificial convergence and divergence of sea ice in partially-covered ocean cells.Therefore, sea ice is allowed to advect freely in partially-covered ocean cells (although the thermodynam-ically derived rate of change of sea ice thickness is reduced due to decreased heat and moisture fluxes).
Simulated bathymetry in the ocean model MOM2 was improved over the default UVic ESCM bathymetry to better represent modern under-shelf conditions (Lythe and Vaughan, 2001).The ocean is therefore free to circulate to the equivalent of the modern-day AIS grounding line.MOM2 is not capable of capturing the effect of a depressed upper boundary geometry.For this reason, as well as the the coarse scale of the ocean model, oceanic effects associated with ice shelf draft are currently neglected.
Modelling and observations indicate that oceanic melting of ice shelves plays an important role in ice sheet evolution (Pollard and DeConto, 2009;Payne et al., 2004).Modelling results that attempt to simulate the interaction between ice shelves and subshelf ocean circulation (Holland et al., 2008a;Olbers and Hellmer, 2009) suggest that an ocean temperature increase of only a few degrees Celcius results in a quadratic increase in ice shelf melt rates.The obvious relationship between ocean temperature and ice shelf melt rates calls for a parameterization that links ice shelf melt rates to coarse modelled ocean temperatures.However, this version of the coupled model does not currently include such a parameterization as a default option.Instead, melt rates and calving are prescribed according to the method described in the Methods section of Pollard and DeConto (2009) for either modern or LGM conditions for the AIS.This parameterization, which generates a time-and-space-evolving distribution of sub-ice-shelf melt, is based on a triplet of melt values (M p < M d < M e ) for protected, deep-ocean and exposedshelf conditions, respectively.This triplet of melt values is combined based on gridpoint-specific bathymetry, access to open ocean, and distance to the current shelf edge, to generate point-specific sub-shelf melt rates.In lieu of a universal calving law, calving is parameterized as high sub-shelf melt, which increases with decreasing straight-line distance to the shelf edge, as described above.The overall parameterization was tuned by Pollard and DeConto (2009) in order to obtain realistic late Holocene and LGM ice shelf extents when combined with realistic sea level.The associated triplet values are given in Table 1.Fluxes of moisture and heat derived from AIS sub-shelf melting are delivered to the ocean along the ice shelf edge.In contrast to the AIS, sub-shelf melt is simply set to high values for the GIS (although the model does not generate spurious GIS shelves in the present-day with the standard Pollard and DeConto (2009) parameterization values).
The ice sheet model has the ability to simulate large-scale grounding line migration (Pollard and DeConto, 2009) similar to that found in West Antarctic Ice Sheet (WAIS) paleoclimate records (Naish et al., 2009).However, the current UVic ESCM ocean component does not support lower and lateral boundary geometry (i.e., bathymetry and coastline) adjustment.For this reason, any potential climatic effects of large-scale (order 10 2 km and above) grounding line migration, such as that implied by the endmember no-WAIS case of Bougamont et al. (2007), are not captured.Here we effectively assume that large-scale oceanographic effects resulting from grounding-line-migration-induced changes to ocean basin geometry are small and recognize that the ability to make ocean boundaries track dynamic grounding lines is an obvious avenue for future ocean model development.
In general, we believe the existing model provides an excellent foundation for development of a complete ocean-ice coupling procedure.

Multiple ice sheet capability, initialization, ice acceleration and I/O
One global, high-resolution ice sheet grid within the global coupled model is not feasible for geometric and computational reasons (Rutt et al., 2009).Therefore, following Rutt et al., the model was designed to accommodate multiple ice sheet "instances", which contain variables describing the entire state of each independent ice sheet in separate derived data type (DDT) variables.The ice sheet instances are accessible via corresponding Fortran 90 pointers (Metcalf and Reid, 1999).This method permits independent regional ice sheet model instances (each with different resolutions, grid sizes, time steps, and ice sheet physics parameters) to exist within a single coupled climate/ice sheet model simulation.
Ice sheets are among the slowest of the exogenic Earth system components to respond to applied climatic forcings.For this reason we have included the ability to accelerate ice sheet evolution relative to the climate.For equilibrium simulations, annual SMB is accumulated, and then used to force the ice sheets for multiple years.Resulting moisture and heat fluxes from ice sheets to the surrounding climate are scaled proportionally (Yoshimori et al., 2001) resulting in a loss of global conservation of heat and moisture.Ice acceleration is limited to no more than 20 times (i.e., one year of climate forcing drives no more than 20 years of ice sheet evolution) for equilibrium simulations.Given the efficiency of the UVic ESCM, this still allows 100 kyr simulations to be carried out within a reasonable timeframe.Long-term (i.e., over 10 kyr) transient simulations can be carried out using a technique in which the orbital and CO 2 forcings are accelerated at the same rate as the ice sheet.Following Calov et al. (2009) we limit the acceleration in these cases to no more than 10 times due to the long dynamic and thermal response times of the ice sheet, as well as the abyssal oceans.Interannual variability is not simulated by the UVic ESCM, so we are not concerned with extraction of a longterm climatology for use in forcing the ice sheet model, despite suggestions that millennial-scale ice sheet evolution may in fact be sensitive to ENSO-timescale climate variability (Pritchard et al., 2008b).
The current model configuration contains two ice sheet instances that represent the GIS and AIS.Both operate on 20 km resolution stereographic grids with grid dimensions of 76 by 140 and 282 by 282 grid cells, respectively.At yearly intervals each ice sheet is targeted and run through the ice sheet subroutines for A i /dt ice timesteps, where dt ice is the ice sheet timestep (years) and A i is the acceleration factor.In the following experiments each ice sheet utilizes dt ice = 0.25 yr.
Both the ice sheet model and climate model have the capability to initialize from observations or the saved states of previous model simulations.Initial ice sheet observations are obtained from Bamber et al. (2001) for the GIS and Lythe and Vaughan (2001) for the AIS.As the dimension of subgridded elevation bin arrays changes interactively with ice sheet model size and geometry, restarts are checked at initialization to ensure they are consistent.If not, new subgridded arrays are constructed and the existing snowpack is redistributed.

Results and discussion
A suite of sensitivity tests was carried out in order to explore the response of the simulated AIS and GIS to model physics, and to approach the model design that produced optimal equilibrium ice sheet states for Eemian, LGM, and late Holocene climate forcings.Equilibrium ice sheet simulations imply non-changing climate forcings, which is not the case in reality.However, since the equilibrium experiments were primarily meant to isolate the effect of model physics on GIS/AIS ice sheet evolution, transient climate forcings (e.g., CO 2 , orbitally-induced changes to insolation and sea level) were intentionally set to constant values in these equilibrium simulations.To assess the simulated effect of recent climate forcings, we carried out a transient simulation to evaluate simulated SMB trends over modern AIS and GIS geometry and the approximate state of disequilibrium of the present-day ice sheets with respect to anthropogenicallyperturbed radiative forcing.
For each equilibrium experiment the model was forced with relevant orbital and CO 2 forcings.The paleoclimate equilibrium CO 2 concentrations were obtained from a composite of the Vostok (Petit et al., 1999), Taylor Dome (Indermühle et al., 2002), and Law Dome (Meure et al., 2006) ice core records.In all equilibrium simulations the ice sheets were initialized from ice-present conditions (glacial  Eemian GIS with refreezing neglected was 31% smaller than the control Eemian simulation with refreezing enabled, equilibrating at 1.51x10×10 6 km 3 .Additional ice loss, compared to the control run with refreezing enabled, occurred mainly along the northern margin, where the NMRR model generated substantially more runoff than the control simulation.The SLR generated from the NMRR model (relative to the late Holocene NMRR simulation) was 4.1 m, probably an overestimate.

Ice sheet evolution sensitivity to surface albedo
Warm snow rapidly becomes less reflective, lowering the surface albedo and absorbing additional incoming radiation in an important positive ice-albedo feedback loop (van de Wal and Oerlemans, 1994) that can be simulated within the cou-pled model.In order to test the strength of the temperaturederived surface albedo feedback within the coupled model we compared a set of simulations with cold, warm, and bare glacial ice albedos set to the cold-snow value of 0.8 (here called the 'constant albedo' or CA simulations) to the control simulations in which albedo varied as a function of surface snow/ice melt.Note that as the ice retreated, the surface land albedo reverted a value typical of tundra.This neglects the effect of boreal forest invasion of ice-free regions, which would reduce the albedo further and enhance the ice-albedo feedback.However, the runs presented here, with the possible exception of the Eemian GIS, do not retreat enough for this effect to be significant.inception was not simulated).The ice sheets were accelerated with respect to the climate by a factor of 20 to hasten the approach to equilibrium.The results described below focus on the effects of climate-model-generated SAT biases, albedo, and parameterized refreezing for each climate state.Simulations with corrected SAT biases, refreezing of meltwater, and variable surface albedo were treated as control simulations.The effect of perturbing the model physics was then analysed by comparing runs without SAT bias correction, refreezing, or surface albedo against the control runs.AIS and GIS ice volumes and areas are summarized in Tables 2 and 3, respectively.

Late Holocene
In order to simulate equilibrium ice sheet conditions representative of the late Holocene, orbital and CO 2 concentrations were set to 0.2 kyr BP.The control AIS late Holocene volume (Fig. 4a) grew to 30.9 × 10 6 km 3 , from the initial (present-day observed) volume of approximately 26 × 10 6 km 3 .This overestimation was not due to the neglect of transient forcings or spurious SMB.Rather, since excess ice grew primarily in regions that are presently observed to discharge ice through prominent ice streams (the Ross Sea sector of the WAIS, the Pine Island Glacier ice drainage basin, the Transantarctic Mountains, and the Amery/Prydz Bay region) it is very likely that the overestimation was due to the   inability of the model to fully represent ice stream features.This is a common problem among recent modelling studies (Stone et al., 2010;Vizcaíno et al., 2010) that requires a fuller representation of outlet glaciers in whole ice sheet models to alleviate.However, the elevation of the central dome of the EAIS was captured to within 100 m, and the total ice extent, including ice shelves, equilibrated at 14.1 × 10 6 km 2 , close to the initial ice area.The control GIS late Holocene volume (Fig. 5a) grew to 3.47 × 10 6 km 3 from an initial observed volume of approximately 2.9 × 10 6 km 3 .Excess GIS ice formed primarily along the coasts.In particular, ice expanded westwards along the southwest and central eastern margins, covering exposed land surface.Elsewhere around the margins, existing ice thickened.It is likely that here, too, incomplete representation of ice streams allowed the excess ice to build up instead of dynamically discharging.Much of the elevations of the central and eastern interior and the interior of the south dome were within 100 m of their initial observed elevations.The northern margin retreated slightly due to the sensitive balance between low accumulation and summer ablation.The combination of slight northern margin retreat and intermittent central and southern margin expansion resulted in a total areal ice extent of 1.91 × 10 6 km 2 , very similar to the initial observed ice area.
Further detailed discussion of SMB is in Sect.3.5, where we examine the effect of possible anthropogenic influences on recent measured/modelled surface conditions, which have not yet had time to impact large-scale ice geometry.

Last glacial maximum
In the LGM control simulation, the model was forced with orbital and CO 2 conditions corresponding to 21 kyr BP.The control LGM AIS total ice volume (Fig. 4e) grew significantly to 40.4×10 6 km 3 compared to the late Holocene simulation, largely due to an imposed 125 m sea level drop, which allowed thick grounded ice to advance into the Ross and Filchner-Ronne basins.The prescribed basal melt rates were also reduced, but as found in Pollard and DeConto (2009) sea level was the predominant control on AIS grounding lines, which advanced to the shelf edge.Interior sectors of the WAIS, particularly those that feed the present-day Ross and Filchner-Ronne ice shelves also thickened, by up to 600 m, in response to the increase in seaward ice thickness.
The EAIS response was less drastic, with ice thicknesses increasing by 0-200 m over much of the interior of the ice sheet.Slight thinning (<100 m) occurred along the coastal EAIS.The Lambert Glacier sector also thinned slightly compared to the late Holocene simulation, due to the spurious lack of ice retreat there in the Holocene simulation.The overall difference in ice elevations was in broad agreement with the reconstruction of Denton and Hughes (2002), with the exception of the Lambert Glacier region.
The control LGM GIS (Fig. 5e) also grew in comparison to the late Holocene control simulation, to 3.68 × 10 6 km 3 .The area also increased, to 2.12×10 6 km 2 , as ablation zones largely vanished and sea level allowed the ice to expand.Despite this growth in ice area, the total SMB decreased slightly as a result of decreased atmospheric moisture transport.The lower SMB resulted in a slight decrease in the elevation of the central ice sheet.

Eemian
For the Eemian control simulation, the model was forced with orbital and CO 2 conditions corresponding approximately to the midst of the Eemian interglacial, −129 kyr BP.Middle-Eemian CO 2 concentrations were very similar to late Holocene levels, leaving insolation anomalies as the only significant change in forcings between late Holocene and Eemian periods.The Northern Hemisphere experienced a strong positive spring-summer insolation anomaly of up to 80 W m −2 at 70 • N compared to the present day, while the Southern Hemisphere insolation remained similar or slightly less than present (Overpeck et al., 2006).
The control Eemian AIS (Fig. 4i) did not respond strongly to the change in forcing and equilibrated at a volume of 31.2 × 10 6 km 3 , almost identical to the late Holocene control run.This agrees with the modelling results of Overpeck et al. (2006), who found a slight AIS surface cooling under Eemian conditions that would not change the AIS volume significantly.
The control Eemian GIS, in contrast, retreated significantly along the northern margin, and to a lesser extent the western margin, as the increased Northern Hemisphere spring-summer insolation resulted in greater ablation.The final equilibrated GIS ice area (Fig. 5i) decreased by 32% compared to the late Holocene control, to 1.29 × 10 6 km 2 .Total simulated GIS ice loss was equivalent to 3.1 m of sea level rise (compared to the equivalent late Holocene simulation), very similar to the mean model-based estimate of Otto-Bliesner et al. ( 2006) of 1.9-3.0m of global mean sea level rise (SLR), but the pattern of ice loss between the two simulations differed noticeably.Otto-Bliesner et al. ( 2006) lost ice primarily from the south, resulting in an ice-free corridor that developed across the location of the Dye-3 ice core.In support of our modelling result, the validity of an icefree Eemian Dye-3 location was challenged by recent work (Willerslev et al., 2007) which tentatively limited the age of Dye-3 basal ice to no younger than 450 kyr.The sensitivity of the northern margin to retreat is also reflected in other recent GIS simulations (Stone et al., 2010;Greve et al., 2011).However, the extent of northern retreat contradicts evidence of Eemian ice at the base of the Camp Century and NEEM ice cores, suggesting either model overestimation of northern margin retreat or a real GIS that did not reach full equilibration with transient orbital and CO 2 forcings during the Eemian.Simulated ice distribution also contradicts evidence of significant vegetation over southern Greenland during the Eemian (de Vernal and Hillaire-Marcel, 2008).
These results warrant further examination that is beyond the scope of this evaluation.However, two points are worth noting: teleconnections that could transport excess heat from north to south (or vice versa) are not significant in the control Eemian simulation, resulting in AIS surface temperatures and SMB that remained relatively unchanged despite the large positive Northern Hemisphere insolation anomaly.An associated simulation in which the Eemian AIS was also forced with the GIS-derived SLR resulted in negligible AIS volume change, ruling out that as a teleconnection, and leaving local oceanic thermal forcing, which was not resolved here as a coupled connection between ice and ocean, as the remaining mechanism available for driving Eemian AIS response.The depth-averaged Southern Ocean water temperature difference (at depths roughly typical of Antarctic ice shelf drafts) between late Holocene and Eemian simulations (Fig. 6) is not typically large, but does reach +0.2 • C in the Ross Sea, suggesting that moderate increases in simulated basal melting of the ice shelf could occur there.

Ice sheet evolution sensitivity to climate model surface air temperature bias
To gauge the effect of seasonally-resolved modelled temperature biases on SMB and long term ice sheet evolution in the model, we first compared the 1970-2001 UVic ESCM SAT with the equivalent ERA40 SAT field to determine the extent of the model SAT bias, then compared the Eemian, LGM and late Holocene control (bias-corrected) simulations to similar simulations forced with raw UVic ESCM SAT (here called non-bias-corrected "NBC" simulations).Use of a constant bias corrector field assumes that the modelled SAT bias is independent of the climate or ice sheet state; this is a common assumption among modelling studies but is difficult to confirm given the lack of reliable and comprehensive gridded paleoclimate SAT records.The simulated 1970-2001 annual average SAT over the ice sheets is compared with ERA40 in Figs.7 and 8.The simulated annual average AIS SAT exhibited a warm bias in the cold interior of Wilkes Land and a cold bias over the interior of the WAIS.Around the coast, however, the absolute bias value was typically less than 5 • C. The annual average SAT over the GIS agreed well with ERA40, with absolute values below 6 • C everywhere except two local regions in the southeast.
Annual SAT maps do not resolve seasonal SAT biases, which are of greater relevance to ice sheet SMB.Figures 7 and 8 also show summer and winter SAT compared to equivalent ERA40 values.In winter, locations on the WAIS are up to 15 • C too cold, while Wilkes Land is up to 15  biases, but still experiences warm biases of up to 10 • C locally on the southeast coast.However, despite their magnitude, the wintertime biases over the AIS and GIS are unlikely to play a significant role in the SMB (Stone et al., 2010) as they are not large enough to drive the modelled SAT to the melt point.In contrast, summer-time biases have the ability to affect melt season length, melt magnitude (through both direct sensible heat flux and the impact of SAT on albedo) and the determination of whether precipitation falls as rain or snow.The general distribution of summer biases over the AIS is similar to the winter but with a much reduced range, from −6 • C over the WAIS to locally +8 • in Wilkes Land.SAT biases over the GIS exhibit a general increasing northwest trend with maximum biases of up to −6 • C over the central northern ice sheet and along the west coast.
The following sections describe the results of the NBC simulations, in comparison with the control simulations.

Late Holocene sensitivity to model bias
The total ice volume of the AIS in the late Holocene NBC simulation (Fig. 4b) decreased to 28.2 × 10 6 km 3 , closer to the observed AIS volume.Despite lack of a bias correction, no additional ablation zones appeared, as the biases over the AIS were too small to push the SAT over the freezing point in the late Holocene climate.The majority of the ice was lost over the central EAIS, indicating that the large scale atmospheric moisture transport was affected by the temperaturebased bias correction over the AIS.
In contrast to the AIS, the total ice volume of the NBC late Holocene GIS increased relative to the control late Holocene ice sheet, stabilizing at 3.59 × 10 6 km 3 (Fig. 5b).The majority of the additional ice volume formed at the northern margin, where the UVic ESCM exhibited a cold summertime bias of up to 6 • C. Neglect of this bias was sufficient to tip the region into a net accumulation zone, allowing ice to expand to the coastlines everywhere and highlighted the fine balance between the simulated accumulation and ablation along the northern GIS margin in the model.

LGM sensitivity to model bias
The NBC LGM AIS ice volume decreased compared to the control LGM AIS volume, to 38.1 × 10 6 km 3 (Fig. 4f).As with the late Holocene non-bias-corrected simulation, ice loss was concentrated over the central EAIS, again suggesting a temperature-bias-related control on large-scale atmospheric moisture transport.
The GIS ice volume in the NBC LGM simulation (Fig. 5f) equilibrated at 3.66 × 10 6 km 3 , a volume almost identical to that of the bias-corrected control simulation.As in that simulation, the margins reached the regressed coastline everywhere.In contrast to the late Holocene climate, neglect of the northern GIS cold bias did not affect ice volume as temperatures in both the control and NBC models were sufficiently low to avoid threshold melt behaviour.

Eemian sensitivity to model bias
The AIS in the NBC Eemian simulation (Fig. 4j) responded similarly to the late Holocene NBC AIS, equilibrating at a volume of 28 × 10 × 10 6 km 3 .Once more the ice volume loss was concentrated over the EAIS, implying a change to the atmospheric moisture transport field.
The NBC GIS (Fig. 5j) diverged strongly from the control Eemian GIS simulation.Whereas the northern margin of the control GIS retreated, generating 3 m of sea level rise, the northern margin of the NBC GIS remained stable and the total volume of 3.59 × 10 6 km 3 was maintained.As with the late Holocene NBC model, northern GIS accumulation in the Eemian NBC model exceeded ablation, which was artificially weakened due to the cold model bias, such that the GIS contribution to Eemian SLR in this simulation was negligible.

Ice sheet evolution sensitivity to meltwater retention and refreezing
Meltwater retention and refreezing is expected to be more important in regions of high melt over existing snowpack, and in general in simulations in which more melt is available to be refrozen.It is also likely to dampen the response of ice sheets to transient climate change forcing (Bougamont et al., 2007).In order to test the effect of refreezing on long-term ice sheet evolution within the coupled model we compared Eemian, LGM and late Holocene simulations with meltwater retention and refreezing neglected (here called the NMRR simulations) to the control simulations with the refreezing parameterization enabled.

Late Holocene
The late Holocene AIS in the NMRR simulation (Fig. 4c) displayed a small response to the removal of the refreezing process due to the lack of significant surface melting, but still exhibited perceptible changes, largely over the Antarctic Peninsula ice shelves, but also over the Ronne and to a lesser extent Ross ice shelves, where ephemeral surface runoff in the model resulted in marginally thinner shelves and correspondingly thinner sheets near the grounding lines.
In contrast, the northern and (to a lesser extent) central western GIS margins in the late Holocene NMRR simulation retreated noticeably compared to the control simulation (Fig. 5c).This change, and associated slight inland thinning, resulted in a late Holocene NMRR GIS volume of 30.9 × 10 6 km 3 and a smaller surface area.As expected, retreat relative to the control simulation occurred in regions where neglect of the refreezing process enhanced annual runoff of surface meltwater.

LGM
Neither the LGM AIS nor the GIS (Figs. 4g and 5g) were significantly affected when refreezing was not resolved, compared to the control LGM simulation.Given that the refreezing process is activated by melting, this result is expected, because melt extent and ablation were negligible over both ice sheets in both simulations.

Eemian
As with the late Holocene and LGM NMRR simulations, the AIS was not significantly affected due to a lack of sustained surface melting anywhere on the continent in the Eemian climate (Fig. 4k).
In contrast to the LGM simulations, the NMRR Eemian GIS (Fig. 5k) differed strongly from the Eemian control simulation, underlining the importance of resolving refreezing when simulating warmer climate states.The equilibrium Eemian GIS with refreezing neglected was 31% smaller than the control Eemian simulation with refreezing enabled, equilibrating at 1.51 × 10 6 km 3 .Additional ice loss, compared to the control run with refreezing enabled, occurred mainly along the northern margin, where the NMRR model generated substantially more runoff than the control simulation.The SLR generated from the NMRR model (relative to the late Holocene NMRR simulation) was 4.1 m, probably an overestimate.

Ice sheet evolution sensitivity to surface albedo
Warm snow rapidly becomes less reflective, lowering the surface albedo and absorbing additional incoming radiation in an important positive ice-albedo feedback loop (van de Wal and Oerlemans, 1994) that can be simulated within the coupled model.In order to test the strength of the temperaturederived surface albedo feedback within the coupled model we compared a set of simulations with cold, warm, and bare glacial ice albedos set to the cold-snow value of 0.8 (here called the "constant albedo" or CA simulations) to the control simulations in which albedo varied as a function of surface snow/ice melt.Note that as the ice retreated, the surface land albedo reverted to a value typical of tundra.This neglects the effect of boreal forest invasion of ice-free regions, which would reduce the albedo further and enhance the icealbedo feedback.However, the runs presented here, with the possible exception of the Eemian GIS, do not retreat enough for this effect to be significant.

Late Holocene
The AIS in the CA late Holocene simulation (Fig. 4d) did not change volume or geometry appreciably compared to the control simulation, reflecting the fact that albedo over the AIS is rarely affected by melting conditions, minimizing the effect of the temperature-dependent ice albedo feedback.However, with consistently high albedos, the CA late Holocene GIS simulation (Fig. 5d) grew dramatically compared to the control simulation with the variable albedo scheme, reaching a volume of 3.75 × 10 × 10 6 km 3 (higher than the control LGM GIS volume).This indicates that the temperature-dependent ice albedo feedback plays a strong role in determining GIS geometry in the model by enhancing melt along the coastal ablation zones and above-equilibriumline-altitude (ELA) melt regions.The difference between the two ice sheets is indicative of the threshold nature of the temperature-dependent ice-albedo effect: when temperatures are largely below freezing (e.g., the AIS), melt and melt-derived ice exposure cannot occur, but if temperatures rise even slightly above the freezing point, the ice-albedo effect rapidly becomes a dominant factor in the overall SMB.

LGM
As with the late Holocene CA simulation, the AIS (Fig. 4h) was not strongly affected by the imposition of constantly high snow/ice albedo compared to the control LGM simulation, due to the lack of significant ablation or melt zones.Similarly, the GIS volume (Fig. 5h) remained the same as in the CA late Holocene simulation (in contrast to the variation in ice volume simulated in the late Holocene/LGM control simulations).

Eemian
The CA Eemian simulation AIS (Fig. 4l) did not respond significantly to imposition of constantly high albedos in the Eemian climate.In contrast, removal of the ice-albedo feedback over the Eemian GIS (Fig. 5l) resulted in a large change in ice volume compared to the control simulation GIS.When albedos were held constant, the resulting lack of additional ice loss was sufficient to keep the GIS at a near-late Holocene volume of 3.45 × 10 6 km 3 , contained in an area only slightly smaller than the CA late Holocene ice area (1.92×10 6 km 2 ).The GIS contribution to the Eemian sea level high-stand in this simulation, compared to the late Holocene CA simulation, was 1 m.

Modern transient simulation
In order to explore the ability of the coupled model control configuration to capture present-day SMB conditions over the GIS and AIS, we carried out a 200 year transient simulation ending at the present day, initialized from the equilibrium late Holocene simulation.The simulation was forced with transient CO 2 and orbital forcings.In this simulation the bias correction procedure, variable albedo, and refreezing were all enabled.
At year 1850 SMB over the AIS in this simulation was 22.2 × 10 14 kg yr −1 .Between years 1850 and 2010, the AIS SMB increased to 23.5 × 10 14 kg yr −1 .This SMB increase was primarily a function of increased accumulation resulting from the increased moisture-carrying capacity of warmer air as a consequence of the Clausius-Clapeyron relation (Weaver et al., 1998), a result that is common to other AIS SMB modelling studies (e.g., Krinner et al., 2007).As expected, despite increases in SAT, AIS temperatures were still generally well below 0 • C at year 2010, limiting the contribution of ablation to the overall SMB.The final value for modern SMB was approximately 7% lower than the 1980-2004 average SMB of van de Berg et al. (2006) and 11% higher than the 1980-2000 average AIS SMB value of Krinner et al. (2007).Additionally, the modelled AIS SMB range and spatial distribution were similar to the ranges and distributions of van de Berg et al. (2006) and Krinner et al. (2007) (Fig. 9a).Accumulation increased from extremely low positive values over the continental EAIS to maximum values of up to 500 kg m −2 yr −1 on the coastal flanks of the EAIS and over the Antarctic Peninsula.Despite the subgridded approach to SMB calculation, relics of the overlying EMBM grid remained in the high-resolution SMB field, and due to the simplified nature of the EMBM accumulation patterns associated with local microclimates were not resolved.At year 2010 simulated ice ablation was limited to Antarctic Peninsula ice shelves and the northern extremity of the simulated Ronne Ice Shelf.Much of the simulated interior of the Ronne Ice Shelf began to experience minor surface melting, but this melt was almost completely refrozen.Overall AIS  melt extent by year 2010 was 6.99 × 10 × 10 5 km 2 , within the range of the 1987-2008 annual average melt extent derived from Special Sensor Microwave Imager radiometer data (Tedesco, 2008).Due to the present lack of interactive ocean-ice shelf coupling, recent Antarctic mass loss attributed largely to oceanic thermal forcing of ice shelves (Payne et al., 2004;Rignot et al., 2008) was not resolved.Late Holocene simulated SMB over the GIS was 4.70 × 10 14 kg yr −1 .Between years 1850 and 2010, accumulation over the GIS increased due to increased atmospheric moisture transport and ablation rose as temperatures warmed.These trends are common to other GIS SMB studies (Alley et al., 2007) and expected in an initially warming climate.The superposition of the competing processes led to an overall SMB that did not change appreciably throughout the simulation period, in contrast with observations and other modelling efforts, which show lower total GIS SMB values, and indicate that either modern GIS SMB is being primarily driven by decadal-scale climate variability, which the model does not resolve, or that the initial model response to imposed increases in CO 2 at high latitudes is too low.Our modern SMB value was very similar to the modelled 1958-2007 SMB of Ettema et al. (2009), which is itself notably higher than previous estimates (e.g., Box et al., 2004).The broad pattern and magnitude of GIS SMB over the GIS (Fig. 9b 2010), although the magnitude of the SMB gradient between the dry interior and the west coast was much less than these studies.The accumulation pattern over south-central/southeast Greenland was anomalously high/low and the south dome accumulation maximum displays a structure that is clearly a relic of the overlying EMBM atmosphere.These observations indicate that the atmospheric model has difficulty capturing the proper small-scale pattern of accumulation here, an issue also noted by Robinson et al. (2009).Simulated ablation accurately captured both the magnitude and distribution of observed GIS melt regions.In particular, high melt (up to 2300 kg m −2 yr −1 ) was simulated in the central western margin and along the north coast.A narrow ablation zone was also resolved along the northeast coast.Over the 2000-2010 period, simulated ablation zones covered 4.81% of the GIS ice area, lower than the 10% coverage modelled by Ettema et al. (2009).18.4% of the GIS surface experienced surface melting.The modelled GIS melt extent area (3.64×10×10 5 km 2 ) was roughly 27% lower than the 1978-2007 average melt extent obtained by passive microwave data (Fettweis et al., 2007) but within the range of observed values.As with the AIS, the lack of interactive ocean forcing meant that recent dynamically-driven changes to ice volume due to ocean forcing (Holland et al., 2008b) were not resolved.
Comparison of the final state of this simulation with a control equilibrium simulation forced with constant year 2000 CO 2 and orbital conditions gives an indication of the simulated disequilibrium between the present-day ice sheets and present-day climate forcing.The control modern day AIS did not differ appreciably from the final state of the transient simulation.However, the GIS ice volume in the modern equilibrium simulation diminished noticeably, resulting in an equivalent SLR of approximately 1.1 m.This result indicates that response of the modelled GIS to CO 2 increases such as those associated with current trends is significant but delayed in time, relative to observations.

Conclusions
A new coupled model ice sheet/climate model has been constructed, consisting of the University of Victoria Earth System Climate and the Pennsylvania State University Ice Sheet models.The two models are coupled across the ice-atmosphere interface using an energy/moisture balance model on subgrid elevation bins that can correct for model surface air temperature bias.Precipitation is delivered to the ice sheets by the overlying vertically integrated, advective/diffusive atmosphere model.Moisture and heat fluxes are routed to the atmosphere and ocean, respectively.Ice shelves shade the underlying ocean from heat and moisture fluxes.The ocean receives heat and moisture fluxes from basal shelf melting and calving.Basal melting is currently prescribed in lieu of a parameterized temperature-dependent coupling.After each timestep the ice sheets return an updated ice-sheet/ice-shelf geometry to the climate model.The overall system is designed to conserve heat and moisture to machine precision.Multiple independent ice sheet model grids within a global climate simulation are made possible through use of Fortran derived data types and pointers.The model can run in "synchronous", "asynchronous", or "synchronous accelerated" modes depending on the length and character of the simulation.The current coupled model is set up to simulate the Antarctic and Greenland Ice Sheets (AIS/GIS) at 20 km resolution within a global climate simulation.
The model was run to equilibrium under Eemian, LGM and late Holocene conditions to test the sensitivity of the simulated ice evolution to surface albedo feedback, refreezing, and modelled seasonal SAT biases.Sensitivity of the final equilibrium ice volumes to these aspects of model physics increased in warmer climates, as they are threshold processes that are activated when the surface temperature rises above freezing.Conversely, cold AIS and GIS climates are less sensitive to variations in treatments of surface albedo, temperature biases, and refreezing.The difference in the changes to ice sheet volumes under different forcing scenarios indicates that perturbations to the model physics can affect not only the equilibrium ice volume and ice geometry but the change in volume and geometry under transient forcing scenarios.This is primarily due to the strong non-linear behavior of ice sheets in the vicinity of the freezing point.
The control simulation is able to reasonably simulate the Eemian, LGM and late Holocene ice sheets, albeit with a overestimate of late Holocene ice volume.As the modern SMB distribution and total fluxes over both ice sheets generally agree well with recent detailed model studies, the overestimation of ice volume is perhaps due to incomplete representation of ice discharge through ice streams.The response of the coupled ice sheets to LGM and Eemian climates is similar to previous modelling and observations in terms of total ice volume and area, and suggests that the model responds reasonably to imposed forcings.The modelled modern AIS/GIS is/is not in equilibrium with boundary climatic forcings, and a substantial contribution to SLR from the GIS is obtained when the model is run to equilibrium with constant present-day forcings.
The lack of significant AIS change to SMB over all the climate scenarios strongly suggests that oceans are the prime control on AIS volume, and that past and future coupled model simulations of AIS response to climate require an interactive coupling between ice shelves and ocean.The model represents an ideal platform with which to test ocean-ice coupling parameterizations.More work (ideally, a transient simulation through the Eemian) is needed to constrain locations of GIS Eemian ice loss and identify the interglacial that resulted in an ice-free Dye-3 location.The magnitude of the GIS contribution to SLR generated by this study agrees with the equilibrium simulations of Cuffey and Marshall (2000) and Otto-Bliesner et al. (2006) but the location of ice loss disagrees sharply, with loss in the present study coming from the northern GIS ablation zone.
Total relief used to determine number of elevation bins in each climate mode rid cell (m) for the Antarctic Ice Sheet (AIS) and Greenland Ice Sheet (GIS), where tal relief is the maximum elevation minus the minimum elevation existing in the ice heet grid cells contained within the boundaries of a UVic ESCM grid cell. 39

Fig. 1 .
Fig. 1.Total relief used to determine number of elevation bins in each climate model grid cell (m) for the Antarctic Ice Sheet (AIS) and Greenland Ice Sheet (GIS), where total relief is the maximum elevation minus the minimum elevation existing in the ice sheet grid cells contained within the boundaries of a UVic ESCM grid cell.

Fig. 2 .
Fig.2.The difference ( • C) between annual mean UVic ESCM atmospheric SAT fields of late Holocene simulations with/without the bias correction applied over the ice sheets.The effect of the bias correction procedure is largest over the AIS and GIS and diminishes with distance from the correction regions.The global effect of the bias correction over both ice sheets is to redistribute heat from the Southern Hemisphere to the Northern Hemisphere.However, the resulting temperature change is less than 0.25 • C everywhere and does not significantly affect the non-ice-sheet components of the system.

Fig. 3 .
Fig. 3. Antarctic drainage basins.Moisture and heat fluxes associated with ice discharge are evenly distributed to coastlines within each drainage basin.

Fig. 6 .
Fig. 6.Eemian minus late Holocene depth averaged ocean temperatures ( • C) , from 130 m to 750 m depth.Red: warmer Eemian, blue: warmer modern.Note the asymmetry of the color bar towards colder values (white is centered on 0 • C).
) corresponded to the general pattern recognized by other modelling and observational studies Ettema et al. (2009); Burgess et al. (2010).Accumulation decreased from south to north, with the highest SMB (approximately 980 kg m −2 yr −1 ) simulated on the south dome and the lowest values (in the absence of ablation) occurring in the central northern GIS.The model also simulated an in-crease in accumulation along the northwest coast that was captured by Ettema et al. (2009); Burgess et al. (

Table 1 .
Pollard and DeConto (2009)ues used in the ice shelf melting parameterization that is employed here.M p = protected shelf melt values, M d = deepwater shelf melt values, M e = exposed shelf melt values.For more information, seePollard and DeConto (2009).Sub-shelf melt triplet values (m yr −1 ice) M p M d M e