Modelling snow accumulation on Greenland in Eemian, glacial inception, and modern climates in a GCM

. Changing climate conditions on Greenland inﬂu-ence the snow accumulation rate and surface mass balance (SMB) on the ice sheet and, ultimately, its shape. This can in turn affect local climate via orography and albedo variations and, potentially, remote areas via changes in ocean circulation triggered by melt water or calving from the ice sheet. Examining these interactions in the IPSL global model requires improving the representation of snow at the ice sheet surface. In this paper, we present a new snow scheme implemented in LMDZ, the atmospheric component of the IPSL coupled model. We analyse surface climate and SMB on the Greenland ice sheet under insolation and oceanic boundary conditions for modern, but also for two different past climates, the last glacial inception (115 kyr BP) and the Eemian (126 kyr


Introduction
Among the conditions determining the long-term evolution of an ice sheet, the amount of snow and ice accumulated on its surface is of primary importance.Together with ablation processes at the ice sheet margins and bottom, the surface mass balance (SMB) constrains the volume of the ice sheet and thus, on a large scale, impacts the evolution of the sea level.It is most directly affected by perturbations of the surface climate and hence the main contributor to the observed ice sheet response to recent anthropogenic warming (Chen et al., 2006;Mote, 2007;Fettweis et al., 2011).The representations of snow accumulated on the Earth's surface used in numerical models are various and historically depend on the intended use of the model.For general circulation and regional climate models, a wide range of schemes exists (Slater et al., 2001) from simple bucket models to sophisticated multi-layer models (e.g.Manabe, 1969;Lynch-Stieglitz, 1994;Douville et al., 1995), while models of ice sheet dynamics directed towards long-term ice sheet evolution frequently use simple parameterisations to estimate the SMB based on mean temperature and precipitation, such as the positive degree-day (PDD) method (Braithwaite, 1984;Reeh, 1991).

Published by Copernicus Publications on behalf of the European Geosciences Union.
As snow transformation processes can have considerable effects on albedo and thus the radiation budget for all surface types (Brun et al., 1992;Gallée and Duynkerke, 1997;Kuipers Munneke et al., 2011), more complex models with more accurate representation of surface processes such as melting and refreezing, albedo variation, layering and percolation are in development and more and more often used in general circulation models (GCMs).The accurate representation of such often non-linear processes allows for a more physically based estimation of the SMB on ice sheets.Some specific regional climate models have already integrated such schemes and studied the surface mass balance (SMB) of ice sheets and its variations in modern climate extensively (Box et al., 2006;Fettweis, 2007;Rignot et al., 2008;Ettema et al., 2009) to estimate ice mass variations that can be compared to those derived from remote sensing observations (Wouters et al., 2008;Rignot et al., 2008;Velicogna, 2009;van den Broeke et al., 2009;Fettweis et al., 2011).In this context, regional atmospheric models driven by global circulation data based on observation (re-)analysis are generally used.The use of GCM forcing recently permitted the first use of such a model under past climate conditions (van de Berg et al., 2011).
In the context of climate model development towards Earth system models (ESMs), comprehensive process-based representations of snow physics are beginning to find their way into global climate models to allow for realistic water and energy transformation on ice sheets (e.g.Ridley et al., 2005;Vizcaíno et al., 2010).In parallel, efforts have recently been undertaken to couple an ice sheet model with SMB fields from regional climate model (Helsen et al., 2012).Still, important spatial scale differences between atmospheric and ice sheet models can require complex interpolations and a possible alternative consists in the continued use of parameterisations including the PDD method to estimate SMB.
A weakness of the latter approach is that such parameterisations contain constants which are usually calibrated to yield results that match present-day observations, typically for Greenland (e.g.van den Broeke et al., 2010), and are thus a priori valid only for this period.van de Berg et al. (2011) indeed find that PDD factors tuned for modern days fail to describe climate under different insolation accurately.They may also be specific for the shapes and locations of the present-day ice sheets.However, ice sheet evolution generally occurs on very long time scales, for which observations are rarely available, and ice sheets take different shapes and appear in different places than at present.With climate changing in function of orbital parameters, phenomena like the daily and seasonal cycle change as well, and hence, possibly, the relationships between temperature and precipitation and the SMB.
Past climate states present an extensive test field to climate and ice sheet models.The Eemian interglacial around 125 kyr BP with its higher boreal summer insolation is an example for a warmer than present climate state (e.g.Huybrechts, 2002).Despite the different origin of the climate change, it is sometimes considered an analogue to anthropogenically caused warmer climate in the future with associated shrinking of the ice sheets (Overpeck et al., 2006;Kopp et al., 2009;van de Berg et al., 2011).The following colder period, the last glacial inception, can serve as an example for a cool climate leading to ice sheet growth (Vettoretti and Peltier, 2003;Calov et al., 2005;Khodri et al., 2005).
In the past, in paleo climate studies focussing on ice sheets, Earth system models of intermediate complexity (EMICs) are generally used.They are characterised by a rather coarse representation of the atmosphere, frequently coupled to finer resolution ice sheet models to study the ice sheets evolution.In these models, the SMB is frequently estimated by the PDD method (e.g.Goosse et al., 2010;Bonelli et al., 2009), while sometimes approaches based on energy balance are used (e.g.Ganopolski et al., 2010;Gallée et al., 1991).
A few recent studies used surface energy balance schemes with albedo parameterizations depending on surface conditions such as snow depth or surface temperature in combination with downscaled EMIC results for coupling of ice sheet and climate models to study the Eemian (Robinson et al., 2011;Fyke et al., 2011).Robinson et al. (2011) obtain constraints on model parameters based on Eemian ice sheet loss, and found the melt rate parameter to be decisive for the shape of the Greenland ice sheet (GIS).Fyke et al. (2011) find large impacts of the use of interactive or constant albedo and meltwater refreezing on the shape of the ice sheet in the Eemian and late Holocene.One of the limitations of the EMIC studies is that natural interannual variability is often poorly represented, although its effect on SMB can be very important (e.g.Wake et al., 2009).
In a general circulation model setup with a snow model including interactive albedo, Otto-Bliesner et al. (2006) found significantly lower than present-day SMB for the beginning of the Eemian interglacial, important mass losses in an ice sheet model forced by these results, and a sea level rise of 2.2 to 3.4 m.However, sea level rise estimates from direct observations are significantly higher at 6-9 m (Kopp et al., 2009) -an example for many open questions in the context of past ice sheet evolution.
In a novel approach, van de Berg et al. ( 2011) use a regional model nested into a GCM and show that the change of melt for the Eemian calculated by a PDD method appears to be 30 % lower than the one in their model, questioning the validity of the PDD approach for past climate, as also concluded by Ganopolski and Robinson (2011).
We complement these studies with the following approach: a comprehensive snow scheme already in use for regional studies (Gallée et al., 2001;Fettweis, 2007;Fettweis et al., 2011) (Braconnot et al., 2008).As this modified insolation is the principal cause of ice sheet melt in the Eemian (van de Berg et al., 2011;Ganopolski and Robinson, 2011) and growth at inception, reproducing the SMB variations in these climates is a suitable test case for models.Directly modelled SMB fields from GCMs may not be sufficient in terms of resolution to drive ice sheet models, but this study is seen as a first step in the analysis of the impact of this new snow representation on the climates simulated by the GCM.We do not focus on the evolution of the ice sheets, nor on their interactions with the oceans, but perform timeslice simulations and analyse the sensitivity of the SMB to surface climate conditions with the atmospheric component of the IPSL coupled GCM.The impact of two sources of nonlinearities on surface mass balance will be token into account in this study -the snow properties treated by a new snow model and the effect of interannual variability of the lower boundary conditions over the seas.
After describing the model versions used and the simulations performed in Sect.2, we evaluate the effect of using an improved versus a simplistic snow scheme in the LMDZ GCM (Sect.3.1) and compare the SMB obtained by the model to the SMB estimated with a PDD method (Sect.3.2) for modern and paleo climates, respectively.For all experiments, coupled simulations including the same atmospheric model provide the lower boundary conditions (LBC).This also allows us to investigate the impact of using climatological mean versus variable LBCs in Sect.3.3 before drawing conclusions in Sect. 4.

The atmosphere general circulation model LMDZ and improved snow representation on ice sheets
The LMDZ atmospheric general circulation model is used in its version 4 (Hourdin et al., 2006), as currently used in the coupled atmosphere-ocean model IPSL-CM5A for the IPCC's 5th assessment report.We use a moderate resolution 96 × 95 grid (3.75 • × 1.875 • ), with 19 layers in the vertical.
The model distinguishes 4 different surface types: "Land", "Land ice", "Ocean" and "Sea ice".At the current stage, the repartition of these types is not adjusted during the run.
The introduction of the snow model just concerns the "Land ice" component.Therefore, we limit our analysis essentially to the Greenland ice sheet, where recent ice core drillings (NEEM, NGRIP) promise new insights on the Eemian and past-Eemian time periods and have already produced numerous findings on properties of ice and snow on ice sheets (e.g.Steen-Larsen et al., 2011;Kuramoto et al., 2011).
In the standard setup, the model comprises a bulk snow layer for the "Land ice" surface, with melting in the case of surface temperatures above the freezing point, T surf > T 0 , at a rate proportional to T surf − T 0 .No transformation to ice is foreseen and runoff is immediately transferred to the ocean in the coupled model version.Snow albedo is constant at 0.77 over land ice, and for both of the two wavelength intervals in LMDZ corresponding to visible and infrared light.
The new snow model is based on the CROCUS scheme (Brun et al., 1989(Brun et al., , 1992)).It is used here in the framework of the snow and ice part of the soil-ice-snowvegetation-atmosphere transfer scheme (SISVAT) (Gallée and Duynkerke, 1997;Lefebre et al., 2002Lefebre et al., , 2003) ) as applied in the regional model MAR (Fettweis, 2007;Lefebre et al., 2005).The model includes the automatic creation and merging of up to 35 snow and ice layers each with its own properties and with evolving thickness.In particular, the snow type is described by the three variables dendricity, sphericity and grain size.The transformation from fresh snow, which is dendritic, towards more or less compact older snow, occurs as a function of environmental conditions such as temperature and the vertical temperature gradient.It will typically lead to small-grained or non-spheric snow under cold conditions but larger, more spheric grains when temperature approaches the freezing point.Melt and refreezing with ice formation are also represented.Optical properties are adjusted as a function of snow properties, permitting accurate modelling of snow albedo and mass balance (Gallée and Duynkerke, 1997;Lefebre et al., 2002Lefebre et al., , 2003)).Albedo is calculated for three wavelength bands independently and then interpolated to the two LMDZ bands.We use the same configuration as Fettweis (2007).
Below the snow model layers, we keep the permanent standard LMDZ soil model, which for the surface type land ice has the physical properties of ice.At the interface to this permanent ice, at the lower boundary of the snow model, we remove one layer when the total snow and ice load in all layers exceeds 5000 kg m −2 in order to allow for long simulations and avoid too thick layers, but supply an additional ice layer when the snow load becomes less than 200 kg m −2 , assuming the presence of a sufficiently thick ice layer on all land ice points.This is consistent with the assumption of a constant ice sheet topography but will need to be adjusted when changing surface types will be allowed and, ultimately, an ice sheet model is coupled to the ESM.The choice of the threshold is a compromise between the need to avoid the disappearance of all snow model layers or excessive snow or ice load on the one hand and minimum manipulation of the snow models layering scheme on the other hand.

Setup of numerical experiments
To examine the impact of different climate conditions on the Greenland SMB, we choose to investigate Eemian and inception climates.These mainly differ from the pre-industrial in their orbital parameters, inducing an amplified seasonal cycle in the Northern Hemisphere for the Eemian but a reduced one for the inception.The same orbital parameters as in the simulations of Braconnot et al. (2008) are prescribed.Because these changes in incoming radiation at the top of the atmosphere also have an impact on the atmosphere via ocean and sea ice feedbacks, we run the atmospheric model with the corresponding sea surface temperatures (SST) and sea ice concentrations (SIC) from a 20-yr slice of the coupled atmosphere-ocean-sea ice simulations of Braconnot et al. (2008), interpolated from the slightly different 96 × 72 grid of their IPSL-CM4 model, as LBC.We compare simulations for three model setups: -the original LMDZ4 model (Hourdin et al., 2006) driven by climatological monthly mean LBC computed from the IPSL-CM4 simulations, repeated each year, -LMDZ + SISVAT driven by the same climatological LBC, -LMDZ + SISVAT driven by the variable actual monthly mean sea surface temperatures from the IPSL-CM4 simulations.
In the following, we will refer to the first setup as standard LMDZ.The simulations are summarised in Table 1.
To improve the numerical stability of the LMDZ-SISVAT model, it was necessary to improve the vertical resolution near the ground in the simulations with this model.The first of the model's hybrid sigma levels was thereby shifted from 296.5544 hPa + 0.9882 p s to 49.8882 hPa + 0.9980 p s , corresponding approximately to a pressure shift from 1004.3 to 1011.7 hPa at a grid point at sea level, or an altitude shift from 71 to 12 m.In a sense this change anticipates the evolution in GCMs towards higher vertical resolutions near the ground in order to achieve a better accuracy and representation of near-surface processes, following the example of regional models.For LMDZ, a resolution similar to ours has already been used by Krinner et al. (1997), in light of the particular importance in of the boundary layers over ice sheets for turbulent fluxes and katabatic winds.The IPCC-AR5 simulations with the IPSL model are carried out at a lower level of a mean approximate altitude of 33 m.The vertical level change alone leads to significant changes in the modelled atmospheric fields.The additional simulations L cPIr (standard levels) and LVcPIr (shifted levels, see Table 1) were thus performed with a slightly more recent version of LMDZ to allow us to separate these changes from the effects of the snow model.
Present-day surface type repartition and topography were used in all simulations.While this clearly presents a restriction of the model with regard to the accuracy of the results obtained for conditions with seasonal cycles of insolation different from today, as provided by the climates of the Eemian and the inception, our goal here was to test the sensitivity of the new version of this model to past conditions, and equilibration of surface types in the model would happen on time scales much larger than considered in this work.Greenhouse gas concentrations were also the same as at preindustrial times in all simulations.

Surface mass balance evaluation
LMDZ-SISVAT and standard LMDZ provide total snow amount at the end of each month as output.The monthly SMB in each simulation i, h i SMB, modelled , is thus directly accessible as the difference to the previous month.
For the estimation of the annual SMB h i SMB, PDD by the PDD method we use the formulation of Reeh (1991) for the number of N i PDD of positive degree-days and the same parameters as Quiquet et al. (2012): where T i Ann and P i Ann are the annual climatological mean temperature and precipitation, respectively.
The solid fraction of f snow of climatological annual mean precipitation P i Ann is obtained with the method of Marsiat (1994) as described in Quiquet et al. (2012); C snow = 5.0 mm day −1 • C −1 and C ice = 8.0 mm day −1 • C −1 are supposed to represent the relation between melt and temperature for snow and ice, respectively.All constants have been chosen by Quiquet et al. (2012) to give a relatively good match to modelled SMB for present-day climate.
We follow this approach using the same formulation in order to compare it to the directly modelled SMB. and present-day simulations of the RACMO2 (Ettema et al., 2010(Ettema et al., , 1958(Ettema et al., -2007) ) and MAR (Fettweis et al., 2011(Fettweis et al., , 1980(Fettweis et al., -1999) ) regional models, which are driven by reanalysis data, is performed in Fig. 1.Data for further stations and details on the comparison are presented in Appendix A. For this model evaluation, we assume that the difference resulting from the different time coverage of observation data and among the models, in particular between preindustrial climate and climate at the end of the 20th century, is unimportant compared to model biases and inter-model differences.We find that in almost all cases, the LMDZ-SISVAT simulations are closer to the station data and the regional models than standard LMDZ.The simulations with modified vertical layers have a notable warm bias at Humboldt and the  other stations in the northern part of the ice sheet (cf.Figs. 1  and A1) during winter, while those with standard levels are cold biased at the stations on the ice sheet and the coastal stations in the south.Simulated summer temperature matches the observations better in the run with a new snow scheme but not in the run only changing the vertical levels (LVcPI) at all stations in Figs. 1 and 10 out of 12 in Fig. A1.Annual mean bias and RMSE for all models and stations are shown in Table A1.Despite the uncertainty expressed in the modeland measurement-derived standard deviations as well as the conceptual uncertainties of this relatively coarse model-data comparison, we attribute this improvement mainly to the change in snow scheme.Overall, we conclude that the simulations with LMDZ-SISVAT are in reasonable agreement with the observed data and the deviations from observations are often of similar magnitude with regard to the two regional models.The same is generally true for precipitation (cf.Fig. A2).
The main variables determining the SMB are near-surface temperature and precipitation.The annual mean temperature will affect long-term ice sheet dynamics, for SMB the mean temperature in the summer season (June-July-August, JJA), relatively well represented by LMDZ-SISVAT compared to the station data, will be decisive, ruling the amount of melt.We thus compare the climatological mean fields of annual mean temperature, JJA temperature and annual mean precipitation in the same five model simulations as above in the form of maps (Fig. 2).Ettema et al. (2009) and Fettweis et al. ( 2011) performed extensive validation of their models and hence we consider them as references from now on.The results of these two models also reflect the considerable uncertainties that persist, in particular regarding precipitation and SMB (Hanna et al., 2011;Ettema et al., 2009Ettema et al., , 2010)).
On the scale of the Greenland ice sheet, annual mean temperatures in the standard LMDZ version show a distinct cold bias, in particular in the centre and at the borders of the southern part of the ice sheet, which disappears in the version with new snow scheme but also in the one with modified vertical layers.Mean temperatures in the two LMDZ-SISVAT simulations and LVcPIr are similar to the results from the two regional model simulations.relatively high annual mean temperatures over the northern coast of Greenland.
The picture is similar for boreal summer (JJA) mean temperatures, however the figure confirms that most of the improvement in LScPI compared to L cPI comes from the snow model (LScPI − LVcPI) and not the change of model levels (LVcPI − L cPI), in particular towards the ablation zone.Despite an obvious lack of accuracy near the coasts due to their low spatial resolution, the simulations with LMDZ including the SISVAT snow scheme represent the field quite well.
Resolution-related differences among the annual mean precipitation fields are more pronounced, but the LMDZ simulations globally fall within the range of the regional models.This is confirmed for the total annual mean ice sheet precipitation listed in Table 2. Common to all LMDZ simulations are relatively high precipitation in the north-west of Greenland and low precipitation in the north-east.Such anomalies in the distribution of precipitation as well as regional temperature biases are the effect of imperfect representation of the large scale atmospheric circulation in the model.They can be assumed to be one of the principal causes for large scale surface mass balance biases.

Analysis of the impacts of changes in vertical levels and of the new snow scheme
We now compare the simulations run with the standard and improved versions of LMDZ for preindustrial climate in more detail in order to evaluate the effects of the model changes.Here we will assume these perturbations are small enough to cause only linear, independent responses in the model.We first evaluate the difference due to the changed model levels in standard LMDZ, LVcPIr − L cPIr, in boreal summer (JJA) and winter (December-January-February, DJF) mean temperatures (Fig. 3a and b, left hand side).These fields then have to be substracted from the difference between the simulations with a new snow scheme and the one with standard LMDZ, LScPI − L cPI, to isolate the effect of the snow scheme (Fig. 3a and b, right hand side).
The change of levels leads to a mean annual warming of about 6 K over the entire ice sheet in winter, and a mean summer warming of about 1 K over Greenland, accompanied by a cooling of the same amount over the surrounding oceans.Similar changes are observed over the other continents and oceans, but a detailed discussion of these global effects is clearly outside the scope of this article.The net effect of the SISVAT snow scheme over Greenland is a cooling of 1-2 K in winter, most pronounced at higher altitudes of the southern half, and a warming of 2-5 K in summer, most pronounced in the north-east.
The modelled surface albedo, fixed in the LMDZ simulations without SISVAT, generally increases with altitude in LMDZ-SISVAT (Fig. 3c and d).The seasonal surface albedo variations can explain the temperature differences noted above: the higher albedo of the fresh snow at the surface in winter (Fig. 3c) leads to a cooling in the south, whereas in the north the effect is absent due to missing sunlight.In contrast, the albedo is lower in summer, except at the highest grid points (Fig. 3d), due to snow transformation and melting, which leads to a higher ratio of absorbed shortwave radiation, and thus a warming.The effect is strongest in the north-east of the ice sheet where precipitation is rare and the age of the surface snow is high.
Cloud cover increases in the model simulations with modified vertical layers, both in winter and summer (Fig. 3e  and f), except for the highest altitudes.This has to be attributed to different functioning of the LMDZ cloud physics schemes with the new levels.There is however a distinct  impact of the new snow scheme on cloud concentration over Greenland, which roughly correlates with the albedo changes, except for the north-east in summer.Increased cloud cover will tend to counteract the increase of radiative loss favoured by higher albedo.As a consequence of the level shift, we find lower annual precipitation in north-west Greenland, and higher precipitation in the south-east and over the Atlantic (Fig. 3g).The net effect on precipitation due to the snow scheme is insignificant over Greenland.The combined effects of level changes and snow models (LScPI − L cPI) on precipitation amount to +33 Gt yr −1 (cf.Table 2).Evaporation (Fig. 3h) increases over the oceans and coastal areas of Greenland as a consequence of the level changes.The introduction of the new snow scheme leads to a decrease of around 0.1 m i.e. yr −1 in the coastal zones, so that the evaporation difference LScPI − L cPI is near zero in the ablation zone, and +16 Gt yr −1 over the entire ice sheet.The precipitationevaporation (P − E) change is hence +17 Gt yr −1 , an increase of 3 %, located mostly in the south-east.

Snow mass balance
Average modelled annual SMB in the simulations with the regional models and the LMDZ simulations L cPI and LScPI are given in Fig. 4a and can be compared to the ones derived with the PDD method, shown in Fig. 4b.No comparable data is available for L cPIr and LVcPIr, as a snow threshold from the coupled model version is applied in these simulations, above which all snow is melted artificially.While the results obtained by the two methods are very similar for each of the two regional models, clear differences appear in the case of the LMDZ simulations.The LMDZ model does not represent adequately the coastal high SMB zone in the south-east, certainly because of the coarse resolution, and generally underestimates SMB in the southern part.The result is more satisfactory for the north.However, the size of the modelled ablation zones with negative SMB, absent in standard LMDZ, is somewhat overestimated in LMDZ-SISVAT, in particular in the north-east.These factors lead to a rather low estimate of total GIS SMB in LMDZ-SISVAT (Table 2).
The map of SMB estimated by the PDD method, unlike the one of modelled SMB, shows several grid points with values below −0.7 m i.e. yr −1 at the very margin of the ice sheet.As the grid cells are much larger than the width of the area with this high melt in reality, their contribution to the low estimates of total ice sheet PDD SMB in Table 2 is important.
We conclude that the SMB obtained with LMDZ-SISVAT is on the right order of magnitude, but limited by the imperfect representation of the general circulation and, mostly by its coarse resolution.Future higher resolution simulations are required for a meaningful and complete validation, in particular near the borders of the GIS.Direct modelling of the SMB appears to be more accurate than the PDD method in In the previous section, we have shown the impacts of the new snow scheme on the pre-industrial SMB of the Greenland ice sheet.The first question we want to address here is whether this impact is climate-dependent.We have therefore examined the differences between the SMBs computed with LMDZ-SISVAT and standard LMDZ for the inception and Eemian climates and compare them to pre-industrial climate in Fig. 5. Modelled coastal SMB (Fig. 5a) is lower with SISVAT than in standard LMDZ in all climates.The most important contribution to this SMB decrease is the albedorelated warming in the model with new snow scheme.In the south of Greenland, the effect of the warming due to the level changes is of equal magnitude, but at least partially balanced by the increase in precipitation.Despite nearly identical temperature differences of 3-5 K (see Fig. 6) for all climates, the impacts on SMB are much more prominent in Eemian climate and much weaker at glacial inception compared to preindustrial climate: the total ice sheet surface mass balance decreases by 155 Gt for the PI, 578 Gt for the Eemian, and 79 Gt for glacial inception.This reflects the high nonlinearity of the temperature-SMB relationship: a roughly similar warming has a much more dramatic effect on an ice sheet in an already warmer climate.
The SMB computed with the PDD method (Fig. 5c) shows a similar picture.For the PI climate, the PDD method shows a slightly higher SMB in the south-west, explained by the higher precipitation in LScPI, but absent in the modelled SMB.Near the coast, the SMB difference in Eemian climate is more pronounced with the PDD method.The total PDD-estimated ice sheet surface mass balance decreases by 224 Gt for the PI, 774 Gt for the Eemian, and 39 Gt for glacial inception.The difference for the inception climate is hence near neutral, despite a similar temperature increase as in the Eemian simulation.We thus conclude from this subsection that the effect of our model developments on surface mass balance is most pronounced in warm climates and that the spread of these SMB changes among climates estimated with our PDD method is higher than the one obtained by the snow model.The annual mean temperature (Fig. 6a) on Greenland and its surroundings at 126 kyr BP is warmer than in preindustrial times by about 1.8 K in standard LMDZ; the difference is about 1.5 K in the model simulation with updated snow scheme and climatological LBC.Interestingly, the GIS mean warming is only about 1.0 K in the simulation with variable LBC.In the glacial inception climate, while a significant cooling is found on large parts of the GIS in the simulations with climatological LBC, this is not the case in the simulation with variable LBC.We conclude thus that the prescription of climatological LBCs seems to lead to an overestimation of mean near surface temperature variation on Greenland in reaction to past insolation changes.

Climate
As the differences in insolation during the two analysed past epochs concern mostly boreal summer in the Northern Hemisphere (Braconnot et al., 2008), their effect on Greenland temperatures is more pronounced in the JJA seasonal means shown in Fig. 6b.While the tendency in the annual mean is approximately the same over ice sheets and oceans, there is a clear amplification of the temperature differences on the ice sheets in boreal summer, whereas the ocean acts as a heat reservoir, and the increase of its annual cycle amplitude is smaller than over land.Eemian JJA climate is warmer by about 5 K in the centre of the GIS in all three model configurations, and inception climate is cooler by about 2.0 K.The difference between simulations with climatological and variable LBC is less pronounced in JJA temperatures than in the annual mean, although the inception cooling is a bit higher in the north in the run with climatological LBC.
Precipitation, plotted in Fig. 6c, is generally amplified in the Eemian.We find a stronger increase in LMDZ-SISVAT with variable LBC, in particular in the south-east of Greenland.High precipitation events may be better represented due to the variable LBCs and hence their increase due to the milder climate reflected in the precipitation field.Climate at glacial inception is drier in the south-east than in preindustrial in LMDZ-SISVAT with climatological LBC, but not in LMDZ-SISVAT with variable LBC, where drying is found in the south and the north-west.These signals are not very strong, and given the high variability of precipitation, care must be taken in the interpretation.
We select the Eemian to illustrate the effect of using monthly time series as lower boundary conditions compared to using their climatological mean.Figure 7a-b show the climatological temperature and precipitation differences between the two experiments with LMDZ-SISVAT using the different boundary conditions.There are indeed few significant differences in temperature and precipitation among the experiments, the most notable one is a precipitation increase of approximately 10 % on the northern tip of Greenland.

Surface mass balance
In this section we illustrate the effect of the simulated climate examined in the previous section on surface mass balance, and the impact of variable LBC.We first focus on melt, and evaluate the fraction of months with a negative modelled monthly mass balance.We find an increase of this fraction by up to 15 % for Eemian climate and a reduction of up to 10 % for inception climate over large parts of Greenland (Fig. 6d).The ratio of the fractional changes in these climates thus roughly corresponds to the ratio of boreal summer temperature changes relative to PI climate.We further note that the net melt frequency variations cover a much larger part of Greenland in LMDZ-SISVAT than in standard LMDZ.This is certainly an effect of the higher mean temperature leading to higher melt extent.There are several differences between the climatological and variable LBC runs.Notably, increase of the fraction of net melt months for the Eemian is slightly weaker in Central Southern Greenland but stronger in the north in the run with variable LBC, while the inception reduction of net melt months is less pronounced in the north in this type of run.
Figure 8a and b show the differences in modelled and PDD SMB for Eemian and inception.At 126 kyr BP, the modelled SMB decreases by 0.1-0.4m i.e. yr −1 all around the ice sheet in standard LMDZ, but by up to 1.0 m i.e. yr −1 in LMDZ-SISVAT.Also, the increase of SMB in inception climate is much more pronounced in LMDZ-SISVAT.
The SMB variations obtained by the PDD method show some agreement with the directly modelled variations, at least for LMDZ-SISVAT.The PDD method predicts however an increase of SMB at several grid points in the north-west As a consequence, the total PDD-estimated Greenland SMB decrease for the Eemian compared to the PI given in Table 2 is relatively high at up to −779 Gt yr −1 for LMDZ-SISVAT with climatological LBCs, compared to a modelled total Greenland SMB decrease of up to −554 Gt yr −1 .The decrease estimates for LMDZ-SISVAT are still high compared to van de Berg et al. (2011), who find a SMB change of the modelled SMB of −405 Gt yr −1 .Both PDD-estimated and modelled SMB decreases are lower and thus closer to the results of van de Berg et al. (2011) in the LMDZ-SISVAT experiments with variable LBC compared to the climatological LBC experiments.
For the glacial inception, the PDD-predicted SMB gains are again much higher than those obtained from the model.
Despite the limited impact of using variable instead of climatological LBCs discussed in the previous section, we find visible effects on surface mass balance in Fig. 7. Modelled and PDD SMB are higher by 0.1-0.3m i.e. yr −1 in most of the southern half of the ice sheet due to higher precipitation.Slightly higher temperature and lower precipitation in the north-west lead to a slightly lower SMB when using the PDD method.The same comparison for last inception and pre-industrial model simulations does not reveal significant differences, though.We conclude that the impact of variable LBC on mean climate is not very high in our experiments, but can have significant influence on SMB, which may in particular be the case in high-ablation climates.A more comprehensive study may address whether these signals persist in longer simulations.

Conclusions
We have analysed the effects of improving the representation of surface snow on ice sheets in an AGCM for preindustrial climate and two paleoclimates.This was achieved by introducing a sophisticated snow model.While the low spatial resolution of the GCM simulations limits the accuracy of the results compared to higher resolved models or observations, the ensemble of grid points still represents a range of typical climate and SMB responses to such a change.These responses are of considerable magnitude in modern climate, and increase significantly in warmer conditions as in the Eemian.
A change of the vertical discretization of the model was required to decrease layer thickness near the ground and ensure numerical stability, indicating that typical GCM level spacings are not sufficient to cover the physical near-surface processes represented in elaborate surface models.The level change leads to additional climatic changes, in particular a warming of the ice sheet in the winter season, and changes in cloud cover and precipitation.The climatic effects attributed to the new snow scheme include the disappearance of a cold bias of the model in the summer season, which can be explained by the seasonal variations in albedo represented in the new model, and again an alteration of the cloud cover.Overall we find a significant improvement of modelled surface climate compared to the previous model version.In particular the modelled boreal summer near-surface temperature, highly relevant for SMB, is in good agreement with station data.
The pre-industrial total modelled ice sheet surface mass balance for the Greenland ice sheet is low compared to a range of regional models (see the list in Fettweis, 2007), and in particular in the light of recent studies estimating early 20th century SMB to be higher by 50-200 Gt yr −1 than at the end of the century (Wake et al., 2009;Hanna et al., 2011).This bias is primarily due to low estimates in the southeast, but could possibly be reduced at higher model resolution including better representation of the complex topography of the region; the modelled SMB could then be more suitable for direct use as input to ice sheet models.In the near future resolutions down to 20 km will be achieved in atmospheric GCMs, and hence a representation of physical processes equivalent to present-day regional models can be expected, including feedbacks in the ablation zone such as between katabatic flows and sensible heat flux with potentially high impacts on SMB.We also note that a suitable representation of the atmospheric general circulation is a prerequisite for realistic SMB simulation that is not met by all of today's GCMs, in particular in the Arctic, and has to stay an essential goal in model development.The implementation of the numerous non-linear processes at the Earth's surface in the polar regions, including those related to snow, may in turn be an important contribution in this context.Differences in SMB results between PDD estimation and direct modelling mostly concern the margins of the ice sheet, where the PDD modelled SMB changes more drastically with surface climate.The PDD method is more sensitive to the different climatic conditions during the Eemian and the following inception, most prominent in boreal summer temperatures, than direct modelling in terms of total ice sheet SMB.We also note that these climate-induced SMB changes differ most in the south of the ice sheet.
The higher sensitivity of the SMB to the climatic conditions in LMDZ-SISVAT compared to the standard version of the model as found for the cases of the Eemian and the last glacial inception is to a high extent due to the interactively modelled albedo.This illustrates the need to represent physically based snow properties in climate simulations for the past, in particular if insolation conditions deviate significantly from modern day.In the Eemian case, the higher summer insolation and warmer temperatures lead to a melt increase higher than in the earlier model study of van de Berg et al. (2011).Such discrepancies are not surprising, given the existing spread of responses to solar forcing among AOGCM models, which are used as boundary conditions.
We further showed that atmospheric GCM simulations that use a repeated annual cycle of climatological LBC can differ significantly from those using the LBC the climatology is based on in terms of SMB for a time span of 20 yr.These differences are caused by modified mean precipitation and surface temperature differences, linked to the missing representation of extreme temperature or precipitation events in the climatological runs.This effect impacts in particular the reaction of surface climate to insolation change in our simulations.Part of these differences may become insignificant for longer term simulations, as the high variability, especially of precipitation, leads to some incertitude, but we still strongly recommend the use of variable LBC forcing where available.

Appendix A Comparison to station data
The model temperatures are interpolated linearly to the station coordinates based on the four nearest grid points.The GC-Net Automatic Weather Stations (AWS) (Steffen et al., 1996) hourly measurements span the period 1995-2008 with differing coverage, and monthly means and standard deviation are given.For the DMI weather stations, we took the average of daily minimum and maximum temperatures since 1877 (Cappelen et al., 2010) to compute the monthly climatology.Note that MAR defines near-surface temperature at the 3 m level while the other models use the 2 m level.Figure A1 shows the temperature comparison as in Fig. 1    12 out of 16 stations.Figure A2 shows the same comparisons for precipitation.
1) Clim.Past, 8, 1801-1819, 2012 www.clim-past.net/8/1801/2012/where T i (t) are the climatological monthly mean temperature from January to December near the surface for the respective model and T 0 = 0 • C. The constant σ = 5.0 • C is supposed to represent the daily cycle and the variability caused by weather in temperature.The SMB is then obtained as

3
Pre-industrial state: comparison of the results of different model versions to observations3.1 Climate variablesIn order to evaluate the mean climate in the simulations, climatologies are computed from 20 yr of the model output fields.An evaluation of modelled temperature in the original LMDZ version and the two simulations with LMDZ-SISVAT for pre-industrial times compared to selected meteorological observations(Cappelen et al., 2010;Steffen and Box, 2001)

Fig. 1 .
Fig. 1.Climatological monthly mean 2-m temperature ( • C) in Greenland at selected AWS stations of the GC-net (Steffen et al., 1996), compared to interpolated model data.(a) Humboldt, (b) Summit, (c) ETH-Camp, (d) South Dome.Vertical bars show interannual standard deviation for the observation data; triangles indicate the seasonal mean for boreal summer (left panels) and the annual mean (right panels).The respective station locations are indicated by a pink dot on the inset map.

Fig. 2 .
Fig. 2. (a) Climatological annual mean 2-m temperature ( • C) in Greenland, from left to right for the regional models RACMO and MAR, the LMDZ reference simulation (L cPI), LMDZ with modified levels (LVcPIr), and LMDZ-SISVAT with the new snow scheme (LScPI).(b) as in (a), but for climatological June-July-August seasonal mean 2-m temperature ( • C).(c) as in (a), but for climatological annual mean precipitation (m i.e. yr −1 ).

Fig. 3 .
Fig. 3. (a) Climatological mean boreal winter (DJF) 2-m temperature, difference between simulations with modified vertical levels and the standard version (LVcPIr − L cPIr, left panels) and difference from introducing the snow scheme (LScPI − L cPI − LVcPIr + L cPIr, right panels).(b) as in (a), but for boreal summer (JJA).(c) as in (a), but for boreal winter mean surface albedo, (d) as in (a), but for boreal summer mean surface albedo (only right column, vertical level changes have no impact on albedo).(e) as in (a), but for boreal winter mean total cloud cover.(f) as in (a), but for boreal summer mean total cloud cover.(g) as in (a), but for annual mean precipitation (m i.e. yr −1 ).(h) as in (a), but for annual mean evaporation (m i.e. yr −1 ).

Fig. 4 .
Fig. 4. (a) Climatological annual mean surface mass balance modelled on Greenland, from left to right for the regional models RACMO and MAR, the LMDZ reference simulation (L cPI) and LMDZ-SISVAT with the new snow scheme (LScPI).(b) as in (a), but for SMB computed with the PDD method (m i.e. yr −1 ).

Figure 6
Figure 6 illustrates the surface climate anomalies for past interglacial and inception simulations compared to the respective PI simulations, in the three model configurations using standard LMDZ and LMDZ-SISVAT with or without interannual variability of the LBC.Differences are shown only where they pass Student's t-test at the 95 % level to assure statistical significance.The annual mean temperature (Fig.6a) on Greenland and its surroundings at 126 kyr BP is warmer than in preindustrial times by about 1.8 K in standard LMDZ; the difference is about 1.5 K in the model simulation with updated snow scheme and climatological LBC.Interestingly, the GIS mean warming is only about 1.0 K in the simulation with variable LBC.In the glacial inception climate, while a significant cooling is found on large parts of the GIS in the simulations with climatological LBC, this is not the case in the simulation with variable LBC.We conclude thus that the prescription of climatological LBCs seems to lead to an overestimation of mean near surface temperature variation on Greenland in reaction to past insolation changes.As the differences in insolation during the two analysed past epochs concern mostly boreal summer in the Northern Hemisphere(Braconnot et al., 2008), their effect on Greenland temperatures is more pronounced in the JJA seasonal means shown in Fig.6b.While the tendency in the annual mean is approximately the same over ice sheets and oceans, there is a clear amplification of the temperature differences on the ice sheets in boreal summer, whereas the ocean acts as a heat reservoir, and the increase of its annual cycle amplitude is smaller than over land.Eemian JJA climate is warmer by about 5 K in the centre of the GIS in all

Fig. 6 .
Fig. 6.(a) Climatological annual mean 2-m temperature anomaly vs. preindustrial simulation ( • C) in Greenland for Eemian (top rows), and glacial inception (bottom row) simulations, in standard LMDZ (left panels) and LMDZ-SISVAT with constant climatological (centre panels) and varying (right panels) lower boundary conditions.(b) as in (a), but for climatological mean boreal summer (JJA) temperature ( • C).(c) as in (a), but for climatological annual mean precipitation (m i.e. yr −1 ).(d) as in (a), but for annual mean fraction of months with negative modelled surface mass balance.
for the remaining AWS and near-coastal weather stations.Table A1 lists the temperature bias of the annual mean of the models and the monthly root mean squared error for all locations.The two simulations with LMDZ-SISVAT have the lowest bias at 10 out of 16 stations and the lowest RMSE in www.clim-past.net/8

Fig. A1 .
Fig. A1.Climatological monthly mean 2-m temperature ( • C) in Greenland as in Fig. 1, but at the GC-net automated weather stations (Steffen et al., 1996) TUNU-N, NASA-E, NASA-U, Crawford Point, SADDLE, and at the DMI weather stations at Danmarkshavn, Upernavik, Ilulissat, and Illoqqortoormiut, Tasiilaq, Nuuk and Narsarsuaq compared to interpolated model data.Vertical bars show interannual standard deviation for the observation data; triangles indicate the seasonal mean for boreal summer (left panels) and the annual mean (right panels).The respective station locations are indicated by a pink dot on the inset map.Also shown is the comparison of the mean model temperature for all of Greenland.
is integrated into the global atmospheric circulation model LMDZ.This allows to study accumulation and SMB of ice sheets under paleo climate conditions.In addition to the preindustrial climate state, we test our model for two paleo periods with contrasted climates, the Eemian, around 126 kyr BP, and the following glacial inception, around 115 kyr BP.Earth's orbital parameters during these periods leads to increased (decreased) summer insolation in the Northern Hemisphere during the Eemian (last glacial inception) Clim.Past, 8, 1801-1819, 2012 www.clim-past.net/8/1801/2012/

Table 1 .
Summary of simulations with the respective forcings.

Table 2 .
Summary of simulations, with GIS area and total GIS SMB components.Also given is the average percentage of months with negative SMB over all ice sheet grid points.

Table A1 .
Bias and RMSE of models for all station locations of Figs. 1 and A1, in K.The mean bias and RMSE for all stations are also given.Values for the best LMDZ simulation are printed in bold, respectively.