Interactive comment on “ Modelling large-scale icesheet – climate interactions following glacial inception ”

The manuscript by Gregory et al. for the first time presents results of simulations of the last glacial inception using a fully coupled AOGCM-ice sheet model. This is an important step towards the study of glacial cycles with comprehensive Earth system models and I would recommend publication of the manuscript in CP after some improvements in discussion of results and limitations of the modeling approach. However, before discussing more or less technical issues, I want to share with the authors my serious concern about the methodology used in their study. As it follows from the last paragraph of the manuscript, the authors do have an intention to improve their modeling approach in the future and I would strongly support this intention because serious flaws in the methodology undermine potential importance of this pioneering work.


Introduction
The repeated formation, advance, retreat and disappearance of extensive Northern Hemisphere (NH) ice-sheets is the defining characteristic of the glacial cycles of the Pleistocene.Ice-sheets are very sensitive to climate change 40 through its effect on their surface mass balance (SMB) i.e. accumulation (mostly snowfall) minus surface ablation (mostly melting followed by runoff).They :

Simulation of recent and incipient glacial climates
We employ the FAMOUS AOGCM (version name AD-TAN, www.famous.ac.uk,Jones et al. 2005), which is a low-150 resolution version of the extensively used HadCM3 AOGCM (Gordon et al., 2000).The HadCM3 simulation of presentday climate is of comparable realism to many more recently developed AOGCMs (Reichler and Kim, 2008).FAMOUS has a grid-spacing of 3.75 • × 2.5 • and 20 levels in the ocean, 155 and 7.5 • × 5 • × 11 levels in the atmosphere.The atmosphere timestep is 1 hour and the ocean timestep 12 hours.FAMOUS is structurally almost identical to HadCM3, and produces climate and climate-change simulations which are reasonably similar to HadCM3, but runs about ten times 160 faster.It is therefore a suitable model for climate experiments involving long or many integrations and it is continuing to be improved.Smith and Gregory (2012) have recently used a later version of FAMOUS (Smith et al., 2008) to carry out transient simulations of the whole of the last 165 glacial cycle but with time-dependent ice-sheets prescribed as a boundary condition, rather than interactively simulated.We carry out simulations for two climates: the Holocene at the start of the industrial revolution (which we will refer to as "recent"), and the climate during which ice-sheets were first growing (called "incipient glacial" or "IG").For  to 115 ka BP, when summer top-of-atmosphere (TOA) insolation at the latitudes of northern hemisphere ice-sheets was at a minimum; in June it was 33 W m −2 less than in the recent climate averaged over 50-75 • N. Reduced summer insolation is favourable for icesheet growth, and studies of glacial inception (such as cited above) have likewise chosen 115 or 116 ka BP.Although our model integrations are many thousands of years long, we fix the orbital parameters throughout, because our aim is only to study ice-sheet-climate interaction, not to produce a simulation of the last glacial cycle.
Following this aim, all our simulations also ppm).We also fix the Greenland and Antarctic ice-sheets with their present-day extent and topography as normally represented in FAMOUS.Thus, the only difference of boundary conditions between the two times of interest is the orbital forcing.We  2004) indicates that SAT was −0.9 K lower than present at 120 ka BP, while by 110 ka BP it was −5.5 K lower than present, equivalent to a rate of cooling of about 0.5 K ka −1 during the glacial inception.The simulated climate drift is tolerably 220 small compared with this.We thus consider both FAMOUS climates to be steady after 300 years.
The climatology of FAMOUS version ADTAN for recent conditions is evaluated by Jones et al. (2005).Its most notable error is a cold bias in annual-mean SAT and associ-225 ated excessive sea-ice cover in high northern latitudes (Figure 3 of Jones et al., 2005); HadCM3 has a qualitatively similar but smaller bias.Obviously this bias is relevant to glacial inception at these latitudes (cf.Vettoretti and Peltier, 2003), but it is most pronounced in winter, and it is smaller 230 in summer, the season when SAT will affect ablation.We discuss the effect of this bias in Section

could
: be counteracted by targeted adjustment of model parameters but this could be deleterious to the simulation of the incipient glacial climate (Smith et al., 2008; 235 Gregoire et al., 2011).Similarly, flux adjustment of the AOGCM or anomaly coupling to the ice-sheet model (e.g.Lunt et al., 2004;Ridley et al., 2005;Vizcaíno et al., 2008) could be used to compensate for the bias, but these techniques depend on the assumption that climate change is small 240 enough to be considered a linear perturbation, which will not be true for the large changes in local climate that develop as a result of growth of ice-sheets (Section :::::: section 5).Global-mean SAT is slightly warmer in the incipient glacial climate than in the recent climate: for years 301-245 400 it is 287.26 ± 0.03 K and 286.95 ± 0.03 K respectively.This is consistent with the positive global-mean orbital forcing in the incipient glacial climate and is a simple demonstration that global-mean quantities For a further comparison, we have produced a simulation of the incipient glacial climate using the HadSM3 climate model, which comprises the HadCM3 atmosphere coupled 260 to a "slab" ocean of 50 m depth with the same grid-spacing as the atmosphere, and heat convergence calibrated for the Table 1.List of FAMOUS experiments discussed in this work.The notation 1:n indicates an asychronous integration with n Glimmer icesheet years for each FAMOUS climate year; in this notation, the synchronous integrations are 1:1."Recent" and "IG" ("incipient glacial") identify the climate conditions, and for each FAMOUS-Glimmer experiment the range of Glimmer years is shown.The two FAMOUS-only experiments were 400 years long, preceded by 300 years of spin-up.

Model configuration
Glimmer years Recent IG FAMOUS-only none none FAMOUS-Glimmer synchronous 1-500 1-500 FAMOUS-Glimmer 1:10 1-35000 1-25000 FAMOUS-Glimmer 1:100 -10000-60000 FAMOUS-Glimmer 1:10 no albedo feedback -1-10000 FAMOUS-Glimmer 1:10 no topography feedback -1-12000 recent climate.This model has been used in previous studies of the climate of the Last Glacial Maximum (Hewitt et al., 2003) and the Holocene (Brayshaw et al., 2010).HadSM3 has higher atmosphere resolution than FAMOUS and smaller biases in its simulation of recent climate, on account of the slab ocean heat convergence.Nonetheless HadSM3 and FA-MOUS are largely similar in the patterns of the simulated differences between the incipient glacial and recent climates.
Both show low-latitude warming in the annual mean, and pronounced cooling in mid-to high northern latitudes in summer (Figure 2).
There is a small increase in annual-mean snow depth in FAMOUS at almost all longitudes north of 60N.However, there is no perennial accumulation of snow at any location except Greenland and Antarctica in the climate.

The ice-sheet model
We use the Glimmer ice-sheet model (version 1.0.4), which 280 was based on the 3D thermomechanical model of Payne (1999), was developed at the Universities of Bristol and Edinburgh (Rutt et al., 2009), and is now an open-source project (developer.berlios.de/projects/glimmer-cism).We use the shallow ice approximation for ice-sheet dynamics.This 285 is not adequate for simulation of ice-streams but that is not a serious concern because we expect the ice-sheets to be frozen at the bed during the early stages of glaciation.In the model , grounded ice cannot advance beyond the present-day coastline; hence it cannot cross the channels 290 between the Canadian islands and the mainland or invade the Baltic Sea  Glimmer can run in several discontiguous domains simultaneously.Making use of this facility, we define two separate Glimmer grids, with square grid-boxes of area 20 × 20 km 2 on Lambert azimuthal equal-area projection planes, to cover 305 two commonly accepted regions of ice-sheet inception in the north-east of Canada and in Scandinavia, which we call the "Laurentia" and "Fennoscandia" domains respectively.The former is centred over Baffin Island at 278 • E and 71.5 • N, with 150 grid-boxes in x and 160 in y, the latter in Swe-310 den at 17.5 • E and 63.5 • N, with 130 grid-boxes in x and 125 in y.We ::: We : run Glimmer with 11 σ-levels, using a timestep of 0.5 year.In both domains, we initialise Glimmer as ice-free with present-day (bedrock) topography from the ETOPO2 dataset (www.ngdc.noaa.gov/mgg/fliers/01mgg04.html, edition of 2001), which is on a 2 ′ grid, and isostatic adjustment is switched off.

Simulation of surface mass balance
We calculate SMB from FAMOUS monthly-mean SAT and precipitation using an annual degree-day ::::::::::::::::: positive-degree-day :::::::: (PDD) : scheme, which is a convenient and well-established empirical approach .This scheme does not explicitly consider energy balance or  (e.g.Reeh, 1989;Huybrechts et al., 1991;Braithwaite, 1995;Herrington and Poulsen, 2012;van den Broeke et al., 2010;Huybrechts et al., 2011) : .:::::     : potential ablation (the melting which would occur if there were no limit on the amount of snow and ice remaining to be melted) is f D, where f is a degree-day factor (8.0 mm d −1 K −1 liquid water equivalent for ice and 3.0 for snow) and D is the time-integral of the excess of T over freezing point, D = max(0,T (t) − T freezing )dt.The Glimmer scheme makes a simple allowance for meltwater refreezing in the snowpack (60% of the annual accumulation refreezes before runoff is produced).
Because we require sufficient spatial resolution to represent the pronounced gradients of SMB across the  steep ice-sheet margins, which have a width of tens of kilometres, SMB is computed on the Glimmer grid, with bilinear interpolation of the temperature and precipitation data from the lower-resolution FA-MOUS atmosphere grid, and adjustment of the temperature to the Glimmer topography using a lapse rate of 8 K km −1 , which lies between the free-tropospheric dry and saturated adiabatic lapse rates.In reality the lapse rate varies in space and time, and the surface-air lapse rate may differ from the free-troposphere lapse rate (Marshall et al., 2007;Gardner and Sharp, 2009).Precipitation is not adjusted for elevation, and all precipitation is assumed to fall as snow in Glimmer, even when it is rainfall in FAMOUS.
In the annual degree-day ::::: PDD scheme, T is assumed to follow a sinusoidal annual cycle, parametrised from monthly means, with random variations from the cycle, normally distributed with standard deviation σ dd , representing the effects of variability on all subannual timescales.To choose an appropriate σ dd , we use 100 years of SAT from the FAMOUS incipient glacial climate, and evaluate D at each Glimmer gridpoint both from hourly data and from the annual cycle as 370 parametrised by the annual scheme.We adopt σ dd = 2.8 K and 3.1 K for Fennoscandia and Laurentia respectively, since these values give the smallest differences in the domain averages between the two calculations of D. These are close to the summer values, because this is when ablation predominantly occurs, although the variability is more than twice as large in winter, as expected in mid-latitudes.
The simplifications of the SMB scheme and the choice of empirical parameters introduce systematic uncertainties.SMB is strongly sensitive to σ dd , with ±100% variation over 380 the range 3-8 K, and variations of tens of percent result from varying the degree-day factors, refreezing fraction and lapse rate within plausible ranges.
In the coupled model, FAMOUS and Glimmer have independent simulations of surface mass balance and runoff.The 385 Glimmer simulation, described in this section, is an input to the evolution of the ice-sheets only.The FAMOUS land surface scheme separately computes the snow cover which is an input to the radiation and boundary layer schemes, indirectly affected by Glimmer through the modification of land 390 surface characteristics described in the next section.The  The FAMOUS ocean receives runoff from the FAMOUS land surface scheme, not from Glimmer.In particular, solid discharge from the ice-sheet (calving) is not applied to the ocean.This excludes an ice-sheet-climate interaction which could be important, for instance in Heinrich events, although 400 not in the incipient glacial period.

Evolution of topography and land surface characteristics
As the area and volume of the ice-sheets evolve, they modify the topography and other characteristics of the ::: surface.

405
This is an essential part of ice-sheet-climate interaction.
In coupled FAMOUS-Glimmer we use similar methods to those of Ridley et al. (2005), whereby the FAMOUS surface is updated once per year to make it consistent with the Glimmer ice-sheets.The coupling code calculates various statis-410 tics of the Glimmer topography within each FAMOUS atmosphere gridbox (mean, standard deviation and mean magnitude of horizontal derivatives), used by atmosphere dynamics, boundary layer and gravity wave drag schemes.These properties are updated for all FAMOUS gridboxes which lie 415 completely within the Glimmer domain.
For each such FAMOUS gridbox the coupling code also evaluates the fraction of the land area having ice cover according to Glimmer.(We count a Glimmer gridpoint as ice- covered if it has >0.2 kg m −2 of ice.)If its ice-covered fraction is not zero, the coupling code deems it to be an icesheet gridbox in FAMOUS; if the fraction is zero, it is an ice-free gridbox.Various properties of the land surface are changed if a gridbox is converted from ice-free to ice-sheet or vice-versa, the most important being that ice-sheet gridboxes have a snow-free albedo of 0.75, while gridboxes from which ice has retreated are assigned a snow-free albedo of 0.15, appropriate to a dark coarse-grained soil.Changing the snow-free albedo when ice cover appears or disappears in a gridbox gives a positive feedback on advance or retreat of the ice-sheet, because higher albedo tends to cool the gridbox and lower albedo to warm it.The coupling code maximises this feedback, because the development of ice cover at a single Glimmer gridpoint is sufficient to convert an entire FAMOUS gridbox to an ice-sheet gridbox.We made this choice by design in these experiments in order to promote ice-sheet growth, and we examine its effect in later sections  thus expands to the size it maintains for many centuries, but because of interannual fluctuation in regional climate, some GCM gridboxes around the margin of 455 the affected area alternate continually between land-ice and ice-free states, giving rise to variability in the ice-sheet area.
In both climates, ice builds up at an almost steady rate over the first 500 years (see synchronous experiments in Figure 3), implying that the net ice-sheet-climate feedback on 460 SMB changes little during this period.In Fennoscandia the rate of growth of ice is about twice as large in the incipient glacial climate as in the recent climate, and in Laurentia about three times as large (compare the magenta lines with the cyan lines), clearly demonstrating the sensitivity of SMB 465 to orbital forcing.In both the Glimmer domains, both FA-MOUS and HadSM3 show summer cooling which is statistically significant in a 5% two-tailed t-test (Figure 4), assuming years to be independent, using the FAMOUS estimate of variance; in contrast, the change in precipitation is not significant.Otieno and Bromwich (2009) also find that summer temperature is more important than precipitation for inception in Laurentia.

Asynchronous coupling
The rate of growth of ice mass (Table 2) suggests that the build-up of ice equivalent to tens of metres of sea level would require tens of millennia.Although fast for an AOGCM, FA-MOUS would take 100 days to simulate 10 ka.To make it practical to carry out several experiments of such a length, we need to accelerate the ice-sheet evolution.Since Glimmer is only a small fraction of the computational cost of FAMOUS-Glimmer, we can conveniently achieve an acceleration by a factor of n by using each year of annual-mean SMB computed from FAMOUS as input for simulating n consecutive years of ice-sheet evolution.That is, the experiment has n ice-sheet years for each climate year.We call this "asynchronous 1:n coupling".Following Ridley et al. (2005), who successfully used this technique in a coupled simulation of the Greenland ice-sheet with HadCM3, we choose n = 10 for the initial growth of the ice-sheets, and n = 100 to reach the final steady state.Provided the ice-sheet topography does not change enough within a coupling interval to make a difference to the climate which would be significant compared with unforced variability, we do not expect this asynchronous coupling substantially to affect ice-sheet-climate interaction.Support for this is given by Calov et al. (2009), who show relatively small effects for n = 10 and n = 100.

:
In the synchronous experiment, the mass M s of either icesheet after time t is M s = t t ′ =1 S s (t ′ ), where S s (t) is the SMB for the ice-sheet for year t and the sum is over years.(Ice discharge into the sea is negligible in the early centuries, Section :::::: section 5.2.)Although the time-mean S s is positive and determines the trend in M s , unforced climate variability means that S s varies slightly from year to year, so M s (t) is not a perfectly straight line.In the asynchronous experiment, the mass after ice-sheet year t is M a = t/n t ′ =1 nS a (t ′ ) where the sum is over climate years.The factor n multiplying each S a magnifies the effect of unforced climate variability.
The rate of growth of ice appears to be different in the asynchronous from the synchronous experiments (compare magenta with black and cyan with green in Figure 3, Table 2).However, there are ten times fewer degrees of free- dom in the asynchronous timeseries.To decide whether the 520 evolution is significantly different in a statistical sense, we have to compare the distributions of the 500 S s values with the 50 S a values.We assume the values within each set to be independent, because the lag-1 autocorrelation in the S timeseries is insignificant at the 5% level.We compare the S s 525 and S a distributions for the four cases (two domains and two climates) using the Kolmogorov-Smirnov test, and in each case we accept the null hypothesis that the distributions are the same at the 5% significance level, for which the critical value (not to be exceeded) is 1.36 (Table 2).Hence we do 530 not have statistically significant evidence that the evolution of M during the first 500 years is different with synchronous and asynchronous coupling.On this basis, we consider it is acceptable to use asynchronous coupling to enable longer experiments.

Steady-state ice-sheets for the recent climate
When the asynchronous experiment for the recent climate is continued, Glimmer reaches a steady state in the Laurentide domain after about 30 ka and in the Fennoscandian after 20 ka (green in Figure 5), with ice volumes of 540 1.01 × 10 6 and 0.30 × 10 6 km 3 respectively (sea-level equivalent shown in Table 2).These are orders of magnitude greater than the real present-day masses of about 80 × 10 3 and 200 km 3 in the Canadian Arctic and Scandinavia respectively (Radić and Hock, 2010).In the FAMOUS-Glimmer 545 simulation, an ice-cap with a maximum thickness of 2200 m occupies southern Norway and Sweden, and most of the Canadian Arctic islands are covered with ice, with the greatest thickness, of 2000 m, on Baffin Island.
Table 2. Ice-sheet mass M and rates of growth dM/dt at various times through the experiments for the recent and incipient glacial ("IG") climates."L" and "F" stand for the Laurentia and Fennoscandia Glimmer domains.The notation 1:n indicates an asychronous integration with n ice-sheet years for each climate year; in this notation, the synchronous integrations are 1:1.The Kolmogorov statistic, which is dimensionless, is used to compare the SMB in the corresponding synchronous and asynchronous integrations, as explained in the text.To convert ice mass in units of metres of sea-level equivalent (m SLE) to ice volume in units of 10

Approach to steady state in the incipient glacial climate
Ice builds up initially on the islands in Laurentia, which are the highest regions and where ice-caps are found at the present day, especially along the east coast of Baffin Island

585
(Figure 6).The mainland remains ice-free, except for the Keewatin region, on the north-west coast of Hudson Bay.
The areas of accumulation have positive SMB (Figure 7a) because they are relatively cold, and the pattern is similar to that simulated by Otieno and Bromwich (2009)   found by Born et al. (2010).This area has positive SMB (Figure 7b) due to high precipitation.In both domains, the distribution of ice thickness is very similar with the two coupling methods (Figure 6).
The total ice mass in both domains continues to grow at approximately its initial rate of 0.7 m ka −1 during the first few millennia of the incipient glacial asynchronous experiment, after which the rate of growth declines (Table 2, Figure 5).Caputo (2007)  First, estimates of sea-level fall derived from foraminiferal oxygen isotope data using a fixed conversion factor could be exaggerated in the early part of the glacial cycle, because the effect of temperature on fractionation may be relatively more important then .Second, ice-sheets were probably 615 growing in other regions as well, such as Alaska, the Rockies and Siberia .(e.g.Kageyama et al., 2004;Bonelli et al., 2009;Ganopolski et al., 2010) : .::::::::  occupied by ice-sheets is about 5 m ka −1 SLE (section 5.2), which gives an upper limit on the accumulation of ice mass that exceeds the highest reconstructed rate.Within the first 10 ka, the ice-caps on the eastern Canadian Arctic islands reach a steady state, and ice growth continues only in the Keewatin region and on the Ungava Peninsula (Figure 8a,c).The most rapid growth in Fennoscandia initially occurs in the south of Norway and Sweden, but after about 10 ka the southern ice cap approaches a steady state (Figure 8b), while the northern ice cap starts to grow more quickly.
At 10 ka an asychronously coupled experiment with an acceleration of 100 for ice-sheet time was begun; it shows greater variability, but a similar trend, to the 1:10 experiment, which was continued to year 25 ka.After 20 ka the northern ice-cap of Fennoscandia becomes joined to the southern icecap and begins to extend a lobe eastwards as well (Figure 8d).After 30 ka this lobe develops into a separate eastern dome, which expands over Finland to become the largest part of the ice-sheet, with a maximum thickness of 2700 m (Figure 8f).After 40 ka, the Keewatin ice-sheet in Laurentia has attained its maximum extent, but it continues to thicken slowly for some further millennia, reaching a maximum thickness of 3200 m (Figure 8e).The ice-sheets reach a steady state after about 50 ka (Figure 5, Table 2).
The changes over time in the geographical pattern of icesheet growth suggest the influence of local climate feedbacks, which are the subject of the next section.In the same way, it is likely that further ice growth in Laurentia and Fennoscandia would be promoted by the local climate change caused by the growth of other ice-sheets outside these domains.In addition to local feedbacks, the fall in the atmospheric CO 2 concentration and the increase in surface albedo due to ice-sheet expansion during the glacial period tended to cool the global climate (shown for FAMOUS by Smith and Gregory, 2012).

Climate feedback on incipient glacial ice-sheets
There are various mechanisms whereby ice-sheets affect their regional climate.Firstly, the high surface albedo of ice and snow tends to reduces the absorbed solar radiation and surface melting.Secondly, thickening of the ice-sheet increases its surface altitude.At higher altitude, surface air temperature is usually lower, and this will tend to reduce surface melting by sensible heating.(This is less important than absorption of solar radiation.)Cloudiness and precipitation also depend on altitude, but the relationships are geographically dependent and may have either sign.Thirdly, an ice-sheet may cool the air blowing over it and thus cool surrounding ice-free areas.Fourthly, the emergence of new topography as the ice-sheet grows may influence the large-scale atmospheric circulation, with consequent local and regional effects on many aspects of climate.All of these mechanisms are represented in FAMOUS-Glimmer, to the extent permitted by the approximations of spatial resolution 680 and physical schemes.As noted above, feedbacks due to changes in runoff are not included in FAMOUS-Glimmer as used in this work.

Regional climate change
As discussed above, substantial face type produces a marked increase in surface albedo, by an average of 0.16 over the affected area, in the difference in the time-means for the first century between the synchronous FAMOUS-Glimmer experiment and the corresponding FAMOUS-only experiment.The albedo increase 700 causes a reduction in absorbed solar radiation at the surface, by 14 W m −2 in the annual average in the affected area, which is comparable with the orbital forcing in summer.This in turn produces an annual average cooling of 2.7 K (Figure 9) within the affected area, and the cooling 705 also extends outside that area i.e. the incipient ice-sheets influence the climate of the surroundings.The albedo feedback thus provides a mechanism for rapid extension of the ice-sheet area (Calov et al., 2009) : , ::::::

710
The small warming over the North Atlantic, Greenland and the Barents Sea (Figure 9) could be caused by the strengthening of 1-2 Sv which occurs in the Atlantic meridional overturning circulation.This in turn could be a response to colder air blowing over the north Atlantic from Lauren-715 tia (Hewitt et al., 2001;Smith and Gregory, 2012).In our model, unlike in that of Born et al. (2010), a decline in oceanic heat transport into the Nordic Seas is not a precondition for ice-sheet growth in Scandinavia.
In the final steady state there is a more complex picture of changes.Where the ice thickness is greatest (Keewatin, southern Norway, Finland, Figure 8e,f), there is further cooling (Figure 10a) on account of their surface elevation, and cloudiness and precipitation increase (Figure 10b), all of which feedbacks tend to promote ice-sheet growth.However, there is warming relative to the first century in various other parts of the northern hemisphere, mitigating some of the initial cooling.(Note that all areas occupied by the ice-sheets are still cooler than in the incipient glacial FAMOUS-only climate.)This is because the appearance of the ice-sheets causes a substantial perturbation to the tropospheric circulation (Figure 11).The north Atlantic is affected by a large anticyclonic anomaly centred over Baffin Island with a secondary maximum over northern Scandinavia.This feature weak-735 ens the zonal flow and the Icelandic low, and surrounding it there is anomalous descent in the mid-troposphere (positive ω amomaly) along the southern and western margins of the Laurentide ice-sheet and over the sea to the northwest of Scandinavia; the prevailing westerly flow produces 740 anomalous descent to the east of the Fennoscandian ice-sheet as well.All these regions experience reduced cloudiness.The reduction in cloud gives increased solar radiation absorbed at the surface (Figure 10c), outweighing increases in surface albedo at the ice-sheet margins, accounting for 745 the warming adjacent to the Laurentide ice-sheet, and associated with a reduction in precipitation.Thus, the circulation changes induced by the Laurentide ice-sheet in particular tend to oppose its growth by warming and drying the climate in its lower portions and around its margins.

:
The region to the east of the Fennoscandian ice-sheet has 755 strongly reduced precipitation, being in its rain-shadow, and becomes cooler, despite increased insolation, presumably because of cold air advection.There is an anomalous cyclonic circulation over southern Europe, which brings increased cloudiness, precipitation and cooling.The areas of warming 14 J. M. Gregory et al.: Modelling large-scale ice-sheet-climate interactions following glacial inception in the north-west Atlantic and the Norwegian Sea coincide with reductions of sea-ice.There is no significant change in the strength of the Atlantic meridional overturning circulation after the early decades of the experiment.

Ice-sheet mass balance
Cooling over the highest parts of the ice-sheets is enhanced in the Glimmer simulation by their greater elevation than the FAMOUS grid-box mean, and precipitation generally increases in these parts (Figure 12).The consequences of these regional climate changes can be seen in the simulated changes in SMB (Figure 7c-f).In some regions that initially have positive SMB, such as south-east Baffin Island, the north-east of the Keewatin region, and northern Norway, the SMB increases considerably.In many adjacent regions that initially have zero time-mean SMB, such as the Ungava Peninsula, most of the Keewatin region, Finland and southern Sweden, the SMB becomes positive and eventually large.Areas of negative SMB, where ice is converging dynamically, appear as narrow strips around the margins of the ice-sheets, as for instance on Greenland in the modern world.SMB increases are dominated by the effect of cooling in reducing ablation, which in general outweighs reduction in precipitation.As an example of this, it is notable that the growth of the ice-sheet in Norway and Sweden reduces the precipitation in Finland (Figure 12b), but this does not prevent a higher dome of the Fennoscandian ice-sheet from developing there, causing precipitation to increase again and the rain-shadow to move east (Figure 12d).
We quantify the effect of the developing ice-sheets on their own mass balance by evaluating timeseries of components (Figure 13) integrated over the areas covered by ice in the final steady state these timeseries apply to constant areas, which are initially mostly ice-free.The area-integral precipitation increases by about 35in Laurentia and 25in Fennoscandia, meaning that the build-up of the ice-sheet has promoted precipitation locally.However, in both domainsthe decrease in area-integral ablation is larger, both proportionately and absolutely, than the increase in precipitation .Declining ablation is associated with falling surface air temperature (Figure 12), which is caused by both decreasing absorbed solar radiation and rising elevation, as discussed above (Section 5.1).Increasing precipitation and decreasing ablation both contribute to the increase in SMB in both domains throughout the growth of the ice-sheets (Figure 13).Note that, until the steady state is attained, some of the precipitation falls on ice-free areas.In these areas, by definition, ablation cancels precipitation, so they make no net contribution to the area-integral SMB.
Ice initially accumulates inland and there is no calving until the ice-sheets reach the coast, which happens in both domains after a couple of millennia(Figure 13).Thereafter calving becomes an important term in the budget.It increases more rapidly in Laurentia because the majority of the ice 815 in the early stage is in ice-caps on islands; the rate of increase is slower later while the ice-sheet is growing on the mainland.In both domains, the area-integral SMB reaches a steady value earlier than the calving flux.This suggests that it is the increase in calving, caused by continuing thickening 820 of the ice-sheets, which finally limits their size.
A complementary picture (Figure 14) is given by the components of mass balance integrated over the ice-sheets excluding the ice-free areas, and shown as a function of icesheet area.In Fennoscandia, accumulation increases lin-825 early with area.In Laurentia, there are different linear relationships in the early (island) and later (mainland) stages.This linearity indicates that local precipitation feedback (Figure 10b) has a rather small influence.Ablation also increases with area, but not linearly: in Fennoscandia it de-830 celerates, while in Laurentia it accelerates.This suggests that local temperature feedback is important in the marginal zone, where ablation mainly occurs, with opposite sign in the two domains.
In Fennoscandia, calving and ice-sheet mass both rise 835 fairly linearly with ice-sheet area.In Laurentia, there is again more than one stage.Initially the area expands rapidly, owing to the surface albedo change, but the ice is thin.Next there is a thickening with little expansion, as the ice-caps grow on the islands.In the final stage, while the mainland ice-sheet is 840 growing, mass and calving both increase quite steadily with area, as in Fennoscandia.Theoretical considerations suggest a power-law relationship M = A γ between mass and area A for steady-state ice-caps, with γ = 1.25 (Radić and Hock, 2010).Excluding the early stages, our Glimmer simulations 845 give γ = 1.8 for Laurentia and γ = 1.2 for Fennoscandia.

Sensitivity experiments for climate feedback on SMB
We separate the effects of albedo and topography on SMB by carrying out two further experiments in which they are 850 separately suppressed.These integrations are asychronously coupled (1:10) and begin from the same starting state as the other incipient glacial simulations.
In one experiment, we suppress the albedo feedback on ice-sheet growth by not updating the snow-free albedo in the 855 FAMOUS-Glimmer coupling.The initial dM/dt is much smaller in this experiment than in the fully coupled experiment (Figure 5 and Table 2), because the area of positive SMB is much less from the outset, and expands little.A steady state is reached in both domains in less than 10 ka, 860 with M about 20 times smaller in each domain than in the fully coupled steady state (Table 2), and with ice in both cases mostly restricted to its initial areas of growth (  more rapidly than by ice dynamics; indeed, it appears that the ice-sheets cannot spread far by dynamics alone. In the other experiment, we suppress elevation feedbacks by (i) keeping the FAMOUS surface topography unchanged as the Glimmer ice-sheet surface topography evolves; (ii) using the Glimmer bedrock topography i.e. its initial surface topography, instead of the Glimmer actual surface topography, for the lapse-rate adjustment in the temperature interpolation.Thus, as far as climate and SMB are concerned, surface topography does not change as the ice-sheets grow.With this design, we cannot distinguish the effect of elevation on local climate from the effect of topography on regional climate.In reality, these two are not clearly separable, but one might think that we could suppress the latter and keep the former by doing (i) but not (ii).However, an increasingly large adjustment using the assumed lapse-rate would then be applied to convert FAMOUS surface temperature to Glimmer surface temperature as the two surfaces became further apart.The re- sults would be very sensitive to the lapse rate and we do not 885 think they would be physically reliable.
In this experiment, dM/dt is initially similar to the fully coupled experiment (Figure 5 and Table 2).It is not surprising that the elevation feedback has no effect initially, because it takes time for ice to build up.Within the first mil-890 lennium, however, M falls behind without elevation feedback, and in less than 10 ka reaches a steady state 3 times smaller in Laurentia than in the fully coupled steady state, and 9 times smaller in Fennoscandia (Table 2).Thus it appears that the topography feedback is eventually important in 895 our experiments, in contrast to the results of Kageyama et al. (2004).Without the topography feedback (Figure 15c,d), the ice-sheets remain rather thin, and they do not expand so far into the Keewatin region and Finland, suggesting that ice dynamics plays a role here.

900
Neither of the experiments achieve ice-thickness as great as in the fully coupled experiment, even in areas that initially

Conclusions
We have coupled the FAMOUS global AOGCM to the Glim-915 mer thermomechanical ice-sheet model, in order to simulate the coevolution of climate and ice-sheets following glacial inception in two domains: north-east America (Laurentia) and north-west Europe (Fennoscandia).This is the first use of a coupled AOGCM-ice-sheet model for a study on long 920 palaeoclimate timescales. :::: schemes. ::::::::

Fig. 1 .
Fig. 1.Difference between 115 ka BP (incipient glacial) and recent conditions (IG minus recent) as a function of latitude and time of year in climatological surface air temperature in years 301-330 of the FAMOUS-only experiments.
composition fixed as it was early in the industrial revolution, as normally represented in FAMOUS (CO 2 290 ppm by volume, CH 4 793 ppb, N 2 O 285 ppb;

:
carried out a FAMOUS-only integrationfor each climate, without Glimmer and without Laurentide and Fennoscandian ice-sheets.(For reference, Table 1 lists all the FAMOUS experiments we analyse in this paper.)They :::: The ::: two ::::::::::: integrations : were initialised from the same FAMOUS state and :::: were : run for 700 years.Since this starting state was obtained from a different version of FA-MOUS, both integrations exhibit fairly rapid initial climate drift, which diminishes to a small value during the first 300 years.During the remaining 400 years, the TOA net radia-tion, which indicates the energy imbalance of the climate system, is within 0.1 W m −2 of zero in both integrations, which is small compared with the global-mean orbital forcing dif-210 ference, due to reduced orbital eccentricity, of 0.74 W m −2 between incipient glacial and recent.The rate of change of global-mean SAT (surface-air temperature) is fairly constant in each integration, at +0.10 and −0.15 K ka −1 respectively for recent and incipient glacial.The Antarctic Dome C deu-215 terium record of EPICA Community Members ( inception.The relevant indicators warming in April, extending over all longitudes but most pronounced over land.

Fig. 2 .
Fig. 2. Zonal-mean surface air temperature difference between the annual and summer (JJA) means of the incipient glacial and recent climates of FAMOUS and HadCM3.
scheme,::::: with influence of variables such as solar radiation and wind .The scheme,::: the

Fig. 3 .
Fig. 3. Ice-sheet mass as a function of time for the first 500 years of the synchronous and asychronous (1:10) experiments.Glimmer domains are indicated by linestyle, and experiments by colour.Values are annual means.

Fig. 4 .
Fig. 4. Area-average difference in (a) Laurentia and (b) Fennoscandia of surface air temperature between incipient glacial and recent monthly climatology as simulated by FAMOUS (blue) and HadSM3 (orange) (IG minus recent).Hatching indicates statistically insignificant differences.

Fig. 5 .
Fig. 5. Ice-sheet mass as a function of time through the whole length of the asynchronous (1:10 and 1:100) experiments.Glimmer domains are indicated by linestyle, and experiments by colour.Values are 100-year means. 580

Fig. 6 .
Fig. 6.Ice-sheet time-mean thickness (m) in years 401-500 of the (a,b) synchronous and (c,d) asynchronous simulations for the incipient glacial climate.The letters in (a,b) indicate the locations of the Keewatin region, Baffin Island, the Ungava Peninsula, Norway, Sweden and Finland.
compares 12 reconstructions of global-mean sea-level change covering the last several glacial cycles.They indicate a sea-level fall of roughly 20-70 m during 120-100 ka BP i.e. 1.0-3.5 m ka −1 .There are a number of possible reasons why the simulated rate falls substantially short of the lowest reconstructed rate.

Fig. 8 .
Fig. 8. Ice-sheet time-mean thickness (m) (a,b) after 10 ka, (c,d) after 25 ka and (e,f) in the steady state for the incipient glacial climate.

Fig. 9 .Fig. 10 .
Fig.9.Difference in the time-means for years 1-100 of surface air temperature (K) between the synchronous FAMOUS-Glimmer experiment and the FAMOUS-only experiment under incipient glacial ("with Glimmer" minus "without Glimmer").The green line shows the extent of the ice-sheet in FAMOUS during these years, as defined in the text.

Fig. 11 .
Fig. 11.Difference in 500 hPa geopotential height (m) between the final steady state and the time-mean of years 1-100 of the synchronously coupled experiment for the incipient glacial climate (final minus initial).

Fig. 12 .
Fig. 12. Difference in the Glimmer surface air temperature (K) between the time-means indicated and the time-mean of years 1-100 of the synchronously coupled experiment for the incipient glacial climate, shown by colours.Contour lines are superimposed showing the different in precipitation (m a −1 liquid water equivalent).

Fig. 13 .
Fig. 13.Components of ice-sheet mass balance in the incipient glacial experiments shown as timeseries of area-integrals over the final ice-sheet area in each domain.Values ::: The ::: first ::: 25 :: ka ::: of ::: the

Fig. 15 .
Fig. 15.Ice-sheet time-mean thickness (m) in the steady state for the incipient glacial climate (a,b) without albedo feedback, (c,d) without topography feedback.
model, ice-sheet surface mass balance (SMB) is derived from the surface climate, and the surface topography and other surface characteristics of the cli- 6 km 3 , divide by 2.53.
Born et al. (2010)h (2009)e that the excessive ice volume simulated by FAMOUS-Glimmer for the recent climate is due to biases both from the SMB scheme and from the simulated climate.The latter may provide the additional cooling whichOtieno and Bromwich (2009)found to be necessary for inception in Laurentia andBorn et al. (2010)in Fennoscandia.In view of these biases, our simulations of the incipient glacial very likely also have biases.The approximations made by the present versions of FAMOUS and Glim-mer mean that we probably cannot accurately reproduce the real evolution of the ice-sheets at the start of the last glacial ::::::::