Adding a dynamical cryosphere toiLOVECLIM (version 1.0): coupling with the GRISLI ice-sheet model

We present a coupling approach to and the first results of the GRISLI ice-sheet model within the iLOVECLIM-coupled climate model. The climate component is a relatively low-resolution earth system model of intermediate complexity, well suited for long-term integrations and thus for coupled climate–cryosphere studies. We describe the coupling procedure with emphasis on the downscaling scheme and the methods to compute the snow fraction from total precipitation fields. We then present results for the Greenland ice sheet under pre-industrial climate conditions at the end of a 14 000 yr long integration. The simulated ice sheet presents too large a thickness in its central part owing to the overestimation of precipitation in the atmospheric component. We find that including downscaling procedures for temperature improves the temperature distributions over Greenland for both the summer and annual means. We also find an ice-sheet areal extent that is overestimated with respect to the observed Greenland ice sheet.


Introduction
The extensive work carried out within the framework of the last Intergovernmental Panel on Climate Change (IPCC) exercise revealed that ice sheets are likely to become the main contributor to sea-level rise by the end of the present century.In order to take up the challenge of accurately assessing the contribution of both present-day ice sheets (Greenland and Antarctica), efforts in the development of new numerical tools have been made in the past few years.For the time being, regional atmospheric climate models, including detailed snow models, are considered to provide the best estimates of the present-day Greenland ice-sheet surface mass balance (SMB), defined as the sum of snow accumulation and ablation (Fettweis et al., 2008(Fettweis et al., , 2013;;Rae et al., 2012;Ettema et al., 2009Ettema et al., , 2010a, b), b).However, these models do not account for the effect of the future evolution of the ice sheet on the climate, a feedback process that could in turn affect the surface mass balance.In parallel, low-resolution general circulation models (GCMs) coupled to 3D thermomechanical ice-sheet models (ISMs) have also been developed to study the effects of the anthropogenic perturbation on the climate-ice-sheet system in the future (Ridley et al., 2005;Mikolajewicz et al., 2007a, b;Vizcaíno et al., 2008;Gregory et al., 2012).In these studies, the ice-sheet model is forced by anomalies of temperature and precipitation and surface melting is computed with the widely used positive degreeday approach (Braithwaithe, 1984) that relates ablation to air temperature.(Vizcaíno et al., 2010) went a step further by using absolute climatic fields and an energy balance scheme to compute the surface mass balance.Due to their computational costs, even at low resolution, these models are well suited to studying the evolution of the climate-ice-sheet system over a few thousand years at most despite increasing performances of supercomputers.The study of periods expanding throughout an entire glacial-interglacial cycle such as those that punctuated the Quaternary era requires climate models that are computationally faster by several orders of Published by Copernicus Publications on behalf of the European Geosciences Union.D. M. Roche et al.: iLOVECLIM-GRISLI-coupled model magnitudes.Earth system models of intermediate complexity (EMIC) coupled to an ice-sheet model generally meet this requirement.However, there is a wide range of complexities among EMICs.The pioneering studies were based on 2-D climate models coupled to simplified ice-sheet models (Gallée et al., 1992;Berger et al., 1998).They mainly focused on the effects of both insolation and CO 2 variations on the long-term evolution of the ice sheets.Three-dimensional ISMs (Greve, 1997;Peyaud et al., 2007) were then used and coupled to the 2.5-D CLIMBER-2 climate model to investigate, for past and future periods, the mutual interactions between atmosphere, ocean, vegetation and ice sheets over timescales ranging from tens to hundreds of millennia (Charbit et al., 2005(Charbit et al., , 2008;;Beghin et al., 2014;Ganopolski et al., 2010;Calov et al., 2005;Ganopolski and Calov, 2011).With a comparable level of complexity, the UViC-coupled model, which includes an energy-moisture balance atmosphere, has also been coupled to a three-dimensional ice-sheet model (Fyke et al., 2011) and applied to past climate equilibrium simulations.The level of coupling achieved in Fyke et al. (2011) is comparable in terms of processes to a different treatment (e.g. the ice-shelf melt in UViC is parameterized to be variable in time, which is not the case for our approach), though the complexity of the atmosphere model is clearly lower, in particular with a simplified hydrological cycle.Another notable difference is the adoption of a bias-corrected approach in Fyke et al. (2011).
Within that wide possible spectrum of climate models, we aim at developing a climate model, including climate and cryosphere components, that is simple enough to be run over multi-glacial cycles and sufficiently complex to provide meaningful comparison to proxy data from the different realms.This is why we chose a model that runs approximately one millennium per 24 h of computation with standalone climate (ocean-atmosphere-vegetation), retaining an oceanic general circulation component, a simplified dynamical atmospheric model with a full hydrological cycle and a simplified biospheric component.
In the current study, we present the initial coupling procedure implementation developed for the iLOVECLIM climate model and the GRISLI ice-sheet model.We cover the downscaling procedure and present the results of sensitivity tests to show the impact of our modelling choices.

Short description of the two models prior to coupling
In the following we first summarize the general characteristics of the climatic components of the iLOVECLIM model.It is followed by a description of the version of the GRISLI ice-sheet model.

iLOVECLIM version 1.0
iLOVECLIM is a coupled climate model of intermediate complexity.It is a code fork of the LOVECLIM1.2climate model (Goosse et al., 2010), of which it retains only some of the physical climate components: the atmosphere (ECBilt), the ocean (CLIO) and the vegetation (VECODE) modules.It is thus a direct descendant of the ECBilt-CLIO-VECODEcoupled model that has successfully simulated a wide range of different climates from the last glacial maximum (Roche et al., 2007) to the future (Driesschaert et al., 2007) through the Holocene (Renssen et al., 2005(Renssen et al., , 2009) ) and the last millennia (Goosse et al., 2005).Details on the recent developments included in the present version and on its difference to the previous ones can be found in Goosse et al. (2010).In the following, we summarize the main features of the model as given in that reference.
The atmospheric component ECBilt was developed at the Dutch Royal Meteorological Institute (KNMI) (Opsteegh et al., 1998).Its dynamical core is based on the quasigeostrophic approximation with additional ageostrophic terms added to improve the representation of the Hadley cell dynamics.It is run on a spectral grid with a T21 truncation ( 5.6 • latitude/longitude in the physical space).ECBilt has three vertical layers at 800, 500 and 200 hPa.Only the first layer contains humidity as a prognostic variable (thus the integrated humidity in the first layer is the total humidity content of the atmosphere).Precipitation, which is the main concern here, is computed from the precipitable water of the first layer and falls in the form of snow if the temperature is below 0 • C. The timestep of integration of ECBilt is 4 h.
The oceanic component (CLIO) is a 3-D oceanic general circulation model (Goosse and Fichefet, 1999) based on the Navier-Stokes equations.It is discretized on an Arakawa Bgrid at approximately 3 • × 3 • resolution.The vertical discretization follows a "z coordinate" on 20 levels.It has a free surface that allows the use of real freshwater fluxes, a parameterization of downsloping currents (Campin and Goosse, 1999) and a realistic bathymetry.CLIO includes a dynamical-thermodynamical sea-ice component that is an updated version of Fichefet andMorales Maqueda (1997, 1999).
The dynamic vegetation model (VECODE) was specifically designed for long-term computation and coupling to coarse-resolution models (Brovkin et al., 1997).VECODE consists of three sub-models: (1) a model of vegetation structure (bioclimatic classification) calculates plant functional type (PFT) fractions in equilibrium with climate; (2) a biogeochemical model computes net primary productivity (NPP), allocation of NPP and carbon pool dynamics (leaves, trunks, soil carbon pools), and (3) a vegetation dynamics model.The latter computes two PFTs (trees and grass) and a dummy type (bare soil).The vegetation model is resolved on the atmospheric grid (hence at T21 resolution) and allows the fractional allocation of PFTs in the same grid cell to account for the small spatial scale needed by vegetation.
An iceberg trajectory module is also implemented (Jongma et al., 2009) but is not activated in the present study.The different modules exchange heat, stress and water.It should be noted that precipitation correction is needed to avoid the overestimation of precipitation over the Arctic and the North Atlantic in ECBilt.Removed precipitation is then applied homogeneously in the North Pacific for water conservation purposes.
For the sake of clarity, we note that the LOVECLIM1.2,as described in Goosse et al. (2010), also includes a dynamical ice-sheet model (AGISM) (Huybrechts, 2002;Goosse et al., 2010).However, this component was not publicly available, hence our motivation to develop our own coupling to a dynamical ice-sheet model (GRISLI) for iLOVECLIM.The two coupled systems are different both in the coupling method and in the ice-sheet model itself.The Greenland icesheet model used in LOVECLIM1.2(AGISM) does not account for ice shelves (Huybrechts et al., 2011), whereas the GRISLI model in iLOVECLIM does.Concerning coupling methods, the main difference is that LOVECLIM1.2(AG-ISM) uses an anomaly mode for the coupling where we use absolute fields.Finally, there are some differences in the refreezing schemes used for computing the surface mass balance (compare Fausto in Charbit et al. (2013) and the scheme in Janssens and Huybrechts (2000)).

GRISLI ice-sheet model
GRISLI is a large-scale three-dimensional thermomechanical ice-sheet model.It was first developed for the Antarctic (Ritz et al., 2001) and then adapted to the Northern Hemisphere (Peyaud et al., 2007).GRISLI is used with the same parameter set as found in Peyaud et al. (2007).The resulting Eurasian ice-sheet extent was found to be in good agreement with reconstructions of the Weichselian.The model runs at 40 km × 40 km spatial resolution on a Lambert azimuthal equal area grid.It includes three different types of ice flow: inland ice, ice streams and ice shelves.The evolution of the ice-sheet surface and geometry is a function of surface mass balance, ice flow and basal melting: where t is time, H the ice thickness, U the depth-averaged horizontal velocity, M the surface mass balance and b melt is the basal melting.The isostatic adjustment of the bedrock in response to the ice load is governed by the flow of the asthenosphere with a characteristic time constant of 3000 yr, and by the rigidity of the lithosphere.The temperature fields are computed both in the ice and in the bedrock by solving a time-dependent heat equation.
Ice flow in grounded ice-sheet areas is governed by the 0-order shallow ice approximation (Hutter, 1983;Morland, 1984).Due to the 40 km grid spacing, single ice streams are not explicitly resolved.Rather, regions of fast-flowing ice are represented using the shallow-shelf approximation (MacAyeal, 1989).This also applies to ice-shelf regions.The difference between ice streams and ice shelf is that the latter obey to the flotation criterion and have zero basal drag, except for pinning points for which a basal drag is applied (20 times lower than that of grounded ice).The location of the ice streams is determined by the basal water head, with ice stream regions corresponding to areas where the sediment layer is saturated (Peyaud et al., 2007).
Calving at the ice-shelf front occurs when two criteria are met: (a) the front grid point has a thickness below 150 m and (b) ice coming from an upstream point fails to maintain the thickness above that threshold.This method is built in present-day observations and yields ice shelves similar to observations in West Antarctica when applied to the Antarctic ice sheet (Ritz et al., 2001).

Description of the coupling procedure
As described above, GRISLI includes land ice sheet but also a floating ice-sheet (ice-shelf) component.A complete coupling of GRISLI to a climate model would therefore include the coupling of the oceanic component (CLIO) to the iceshelf model to allow an interactive computation of the basal melting rate of the ice shelves and subsequent freshwater release into the ocean.While desirable, the question of how to parameterize the melting/refreezing under the ice shelves (a very small-scale process with respect to our model grids) from an oceanic temperature is an ongoing research question of its own (Beckmann and Goosse, 2003;Alley et al., 2008) that will be the subject of future studies.For the present work, we use the crude but simple assumption that the melting rate under the ice shelves is constant at a prescribed value depending on the local water depth.We use 2 m per year where the water depth is less than 600 m and 5 m per year where water depth is more than 600 m.This has been shown (Ritz et al., 2001) to be a reasonnable assumption for present-day in Antarctica.However, since our simulations are for preindustrial conditions in the Northern Hemisphere, no significant ice-shelf areas are expected.
In the following we therefore focus on the coupling of EC-Bilt to GRISLI, that is, the exchange of precipitation and surface temperature on the one hand and of surface altitude and ice-sheet mask on the other hand.Moreover, instead of using anomaly fields with respect to the simulated presentday climate as inputs of the ice-sheet model, we use absolute fields from ECBilt for precipitation and temperature.In order to remove potential biases of the climate model, the perturbation method is sometimes used in studies based on climate models including an interactive ice-sheet component (e.g., Vizcaíno et al., 2008;Huybrechts et al., 2011).However the perturbation method relies on the strong assumption that model biases prevailing in a given climatic context are of the same order of magnitude as those in the present-day context.Also, use of perturbation or bias-correction methods makes analysis of feedbacks less robust.Using absolute fields is therefore an important requirement to be able to use the model in different climatic contexts.A detailed overview of the coupling scheme outlining the exchange of energy and mass and the timesteps of coupling is given in Fig. 4.

Coupling method: accumulation and PDD (positive degree-day method)
The upper boundary conditions for the ice-sheet model are the surface temperature and the surface mass balance (SMB).
The SMB is the sum of ice accumulation minus the surface ablation, that is, sublimation of ice and meltwater from melting ice.Both accumulation and surface ablation are computed from the state of the atmosphere overlying the ice sheet.
In our simplified model setup, accumulation is simply the sum of falling snow precipitation, converted into an ice accumulation as follows: (2) Ablation is controlled by the energy exchange between the snow layer at the surface of the ice sheet and the atmosphere.However, the spatial scales of the processes that need to be described are at least one order of magnitude smaller than the spatial resolution of the type of climate model needed for multi-millennia integration.One classical approach to overcome such limitations is to use the widespread empirical PDD as a unique surrogate for ablation.Originally introduced by Braithwaithe (1984) and further developed by Reeh (1991), the PDD represents the sum over one year of the excess of surface temperatures above the melting point.It is expressed as follows: where T m is the monthly temperature and σ the standard deviation of temperature distribution.The conversion of the given PDD to snow-and icemelt rates requires melt rate coefficients.Additionally, water melted at the surface of an ice sheet may refreeze.To take into account those mechanisms, several refinements of the original formulations have been proposed.The reader may refer to Charbit et al. (2013) for a detailed discussion of the impact of the different formulations on ice-sheet build-up as well as for the impact of the different parameters used.In the following, we chose the method of Fausto et al. (2009) that includes a temperature dependence for the ice-and snowmelt rate parameters and an altitudinal dependence for the amount of refreezing and for the σ coefficient of Eq. (3).

Interpolation of climatic variables
For the practical implementation we need to interpolate the climatic variables of ECBilt from the T21 spatial resolution to the finer GRISLI grid at 40 km × 40 km resolution.resolution, two methods were tested to obtain the temperature and precipitation used to compute the surface mass balance (SMB) at the fine resolution of the GRISLI grid.The first method is to perform a simple interpolation for each GRISLI cell, using the neighbouring ECBilt cell.Whenever this is done in our scheme, we use a bilinear interpolation considering a GRISLI grid point and the four surrounding corresponding ECBilt centre grid points.Applying this simple interpolation to both temperature and precipitation already yields reasonable results, as shown hereafter.To further improve the representation of the local surface temperature and therefore the SMB estimate, we additionally added a vertical downscaling approach that explicitly takes into account the differences in altitude between the ECBilt and the GRISLI grid.

Vertical downscaling
In some places, there is a large height difference between the ECBilt and the GRISLI surfaces (Fig. 3).This is especially true in areas where the topography is steep (i.e.varies a lot over a short distance) like on the flank of the Greenland ice sheet.By contrast, when the topography is relatively flat, like in central Greenland, the differences are smaller.Using the temperature at the altitude of ECBilt, even when it is interpolated to the GRISLI grid, does not account for the large temperature differences in both models resulting from the different altitudes( as exemplified by the number of points of GRISLI within an ECBilt cell; cf.Fig. 2).To take the altitude into account within the coupling procedure, we therefore need to include some form of vertical downscaling.This is also true for accumulation, as the shift from liquid precipitation to snow is based on temperature in the ECBilt grid.A good procedure also needs to include the downscaled temperature to convert liquid precipitation from ECBilt to snow accumulation in the GRISLI grid, as done hereafter.

Temperature downscaling
Surface temperature dependance on altitude is computed in ECBilt in a parameterized way.In fact, the model has only three vertical layers and therefore does not fully resolve the vertical profiles of temperature in the atmosphere (Opsteegh et al., 1998).In particular, ECBilt does not explicitly resolve the atmospheric boundary layer.The temperature between the near-surface and 200 hPa-level is assumed to be linear in the logarithm of pressure, the profile being forced to pass through the prognostic temperatures computed at 650 and 250 hPa.Furthermore, ECBilt assumes no heat capacity at the surface of the Earth implying a zero net heat flux between the atmosphere and the Earth's surface, allowing computation of a surface temperature.
To obtain the surface temperature at the GRISLI altitude we therefore compute the ECBilt surface temperature at its own height, but for two virtual surfaces corresponding respectively to the lowest and highest GRISLI points within the same grid cell (cf.Fig. 5).We thus obtain a total of three surface temperatures along a virtual slope, consistent with the temperature computed within the ECBilt model.surface that is then used to add a corrective term to the temperature interpolated to the GRISLI grid as follows: where T interp. is the ECBilt surface temperature on the GRISLI grid, γ is the vertical surface temperature gradient and H is the altitude difference (positive or negative) between the considered GRISLI cell and the corresponding EC-Bilt cell.It should be noted that the vertical surface temperature gradient computed with this method is different from that of the free-atmosphere vertical temperature gradient, since it is computed using only surface temperatures.Our method thus provides a different surface temperature gradient than would be computed using the free-atmosphere one.We refer to it as the "along-slope surface temperature gradient".The γ variable is computed in ECBilt every model month, on the basis of the maximum and minimum temperatures along-slope that are accumulated every four model hours.This procedure ensures that the downscaled temperature obtained (in contrast to procedures using a constant value -both in time and in space) is coherent with the internal physics of the climate model and is thus useable for any climate that ECBilt can simulate.The spatial effect of the vertical temperature downscaling is shown in Fig. 6 for the annual mean.Mountain ranges appear colder than the rest of the region (e.g. the Alps), as (2) A vertical lapse rate γ is computed from the these two extreme temperatures, using the line defined by the two temperatures and elevation extrema in addition to the EC-Bilt cell temperature.(3) Using γ , temperatures are derived for all altitudes in GRISLI if the latter are known.−10.−9.−8.−7.−6.−5.−4.−3.−2.−1.−0.5 0.5 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.

°C
Figure 6.Annual mean temperature anomaly due to the downscaling procedure, in degrees, computed as the difference between the downscaled temperature and the interpolated temperature fields.This is one example taken from one particular year.expected for higher altitudes.There are regions in Greenland where the ice-sheet surface is lower than the mean of the ECBilt cell, others where the ice sheet is higher than the corresponding ECBilt cell.Hence, the alternation of warmer and colder areas as a consequence of the downscaling over Greenland.Computed lapse rates are shown in Fig. 7 for February and July.The computed lapse rates are in good accordance with present-day observations for Greenland (Steffen and Box, 2001), indicating a higher gradient between coastal Greenland and central Greenland for winter (6 to 9 • C per km) and a reduced gradient for summer (4 to 5 • C per km).Our values of 5 to 8 • C per km in February and 4 to 5.5 • C per km in July thus fall within the range of observations.We can therefore conclude that our simple method gives reasonable results for the computation of lapse rates at the Greenland surface in the present-day climate.

Accumulation and downscaling
For the computation of the accumulation on the GRISLI grid, two methods were developed.The first (called SNOW hereafter) consists of using the snow amount calculated in ECBilt every four hours, depending on the temperature calculated in the ECBilt grid.This method has a high temporal resolution, with changes in the amount of liquid precipitation and snow every four hours, but it does not reflect the facts that more snow may occur at higher altitude on the GRISLI grid since higher altitudes are generally colder.Higher altitudes are also generally drier, but that effect is ignored here.To transfer that accumulation from the ECBilt to the GRISLI grid, a simple bilinear interpolation is performed.The second (called PRE-CIP hereafter) is the opposite: we take the total precipitation (liquid plus snow) on the ECBilt grid, accumulate it over time and transfer it to the GRISLI grid once a month (cf.Fig. 4).Since we have also downscaled the surface temperature in the GRISLI grid, it is possible to use it to compute the fraction of the total precipitation that is delivered as snow.Deriving a snow fraction directly within ECBilt would require performing the interpolation between the two grids every atmospheric timestep, that is, every four hours, whereas the downscaled temperature allows performance of this interpolation every coupling step, thus saving computing time.To convert total precipitation into accumulation for GRISLI, one needs to define a (monthly) threshold at which liquid precipitation is turned into snow.We have implemented several different solutions (not shown) with simple assumptions (a limit in temperature or a function in a temperature range).Overall, we found that the model is not sensitive to the choice made and thus we decided to use the following: we assume that the snow fraction is zero for cases above a particular threshold in www.geosci-model-dev.net/7/1377/2014/Geosci.Model Dev., 7, 1377-1394, 2014 monthly temperature and is one for cases below that threshold.In the current study, we assume that T threshold = 2 • C.
In the present version, we do not account for vertical drying out of the atmosphere in the downscaling procedure; hence the total precipitation taken in one ECBilt cell is assigned to the respective GRISLI cells without specific vertical redistribution.

Orography and ice-sheet mask
The information that needs to feed back from the GRISLI model to ECBilt is the altitude of the surface computed in GRISLI (which depends on the dynamics of the ice sheets but also on isostatic adjustment that results from ice loading on continents) and an ice mask, as ECBilt distinguishes between the different surface types.The surface albedo in ECBilt is then computed from the ice mask provided, as in the standard LOVECLIM model.The orography is averaged onto the ECBilt grid as follows: all the GRISLI cells contained in one ECBilt cell are averaged, using a simple mean.We compute an ice mask on the GRISLI grid defined as "one" when the ice thickness is greater than 50 m and "zero" when it is less than 50 m.The rationale behind such a mask is to eliminate very small areas of ice that cannot be seen correctly by EC-Bilt because of its coarse grid and cannot be adequately computed from the shallow ice approximation used in GRISLI.This ice mask is averaged onto the ECBilt grid as is done for the orography.
The coupling between ECBilt and GRISLI is performed every coupling timestep, a value that can be freely chosen and taken as one year in the following (real-time coupling); see Fig. 4. A possibility of an asynchronous coupling is present to allow the computation of more ice-sheet years than climatic years, that is, a number of ice-sheet years with a fixed climate (as, for example, in Calov et al., 2009).For example, a factor of 10 in the asynchronous coupling means that we compute 10 ECBilt model years, then use them to the GRISLI ice-sheet model, which computes 100 yr (10 × 10), and then only feed back orography and ice mask to ECBilt.The use of an asynchronous coupling allows the simulation of more ice-sheet years than climatic years and thus makes it possible to speed up the computation.

Experimental setup
All experiments are run from a present-day Greenland topography (Bamber et al., 2001) for GRISLI and start from a previous pre-industrial equilibrium climate run of the iLOVECLIM model for the climate state.
To evaluate the effect of ice-sheet coupling on the climatic fields and on the ice sheet itself, we choose to perform two types of experiments.The CTRL simulation is performed with the iLOVECLIM model without being coupled to GRISLI, using the fixed present-day ice sheet.A second set of two experiments including the integration of the GRISLI ice-sheet model is started from a present-day Greenland topography and a pre-industrial equilibrium climate run.Two runs are performed using the interactive climateice-sheet coupling: either the coupling is achieved through snow accumulation calculated by ECBilt (hereafter SNOW) or with precipitation from ECBilt converted to snow in the GRISLI grid after downscaling (hereafter PRECIP).The coupling between ECBilt and GRISLI takes place at the end of every year and the runs are performed until the ice sheet is equilibrated.As can be seen from Fig. 8, equilibrating the ice sheet under present-day conditions in terms of volume requires about 12 000 yr with our setup.We integrated a total of 14 000 yr and used the last 1000 yr for the analysis.It should be noted that integrating the coupled climate system during 14 000 yr under constant climate forcing is a classical theoretical state study (Fyke et al., 2011;Lipscomb et al., 2013, for example); indeed the ice-sheet evolution over the last 14 000 yr saw part of the last deglaciation and was therefore further from equilibrium than simulated in our setup.Our setup is therefore not quite comparable to the presentday ice sheet.To further evaluate the effect of climate forcing and of ice-sheet evolution, we performed three additional experiments.The first is an offline equilibrium of the GRISLI ice sheet with exactly the same PDD setup as in the coupled run but the forcing fields come from the observational climatology instead of iLOVECLIM.As in Charbit et al. (2007), it is built from ERA-40 reanalyses for the temperature field; the precipitation field is derived from a compilation between the Climate research unit (CRU) data set Table 1.Summary of experiments performed."Snow acc." stands for SNOW accumulation computation for the PDD; "liq.prc.acc." stands for liquid precipitation accumulation for the PDD.In the case of liq.prc.acc., the accumulation given to the SMB calculations is the snow computed from the liquid precipitation and the temperature downscaling."Coupling" refers to the ice-climate feedback."Climatology" refers to the use of observed fields instead of iLOVECLIM for the climate forcing.

Exp. label
Coupling Snow acc.Liq.prc.acc.Climatology over continents (New et al., 1999) and the Global precipitation climatology project (GPCP) data set over oceans (Adler et al., 2003).In addition, precipitation data for the arctic area comes from Serreze and Hurst (2001).This experiment is referred to as CLIMICE in the following.Furthermore, since we also want to understand the effect the ice-sheet feedback has on the atmospheric model, we performed two semicoupled climate-ice-sheet simulations where the feedback of the ice sheet on climate is removed.That is, the atmosphere model is permantly forced by the observed ice-sheet thickness and extent.There are two simulations since we need one each for the PRECIP and SNOW method.They are referred to hereafter as PRECIP-SO and SNOW-SO (SO standing for semi-offline).
In the following we analyse the general outcome and differences between the three simulations CTRL, SNOW and PRECIP and use the CLIMICE, PRECIP-SO and SNOW-SO experiments for investigating the strength of the different feedbacks.
Our rationale for comparing various simulations is as follows.Differences between the observed ice sheet and CLIMICE reflect the inaccuracies in the climatology, the use of the PDD scheme, the tuning of the GRISLI ice-sheet model and the use of a permanent climatological equilibrium instead of the climate-evolving forcing that occurs in reality.Differences between the CLIMICE and PRECIP-SO and SNOW-SO arise due to the use of the iLOVECLIM preindustrial climatology instead of the observed climatology.They are mainly due to the biases of the climate model, as seen through the PDD scheme.Finally, differences between PRECIP-SO (SNOW-SO) and PRECIP (SNOW) are due to the long-term effects of coupling back the ice sheet to the climate model.Differences between the SNOW and the PRE-CIP experiments are due to the different treatment of the accumulation, as detailed in Table 1.The CTRL experiment does not include a dynamical ice sheet.

Simulated thickness of the Greenland ice sheet
The thickness of the observed present-day ice sheet (Bamber et al., 2001) and the modelled ice sheets are shown in Fig. 9.It should be noted that once interpolated to the GRISLI grid, the initial volume of the observed ice sheet is 2.8 × 10 15 m 3 , 1 × 10 14 m 3 lower than in Bamber et al. (2001).The calculated ice-sheet thickness and extent are overestimated in both SNOW and PRECIP experiments, with an excess volume of about 1.05 × 10 15 m 3 (cf.Fig. 8), that is, one third too much with respect to the observed present-day ice sheet.The excess volume is, however, of the same order of magnitude as what is obtained in comparable studies: Fyke et al. ( 2011) obtained 3.61×10 15 m 3 , while Lipscomb et al. (2013) obtained between 3.2 and 3.9 × 10 15 m 3 .Over the transient part of the simulation, the PRECIP setup consistently yields higher ice volume than the SNOW setup.However, when the equilibrium is reached, the remaining difference is minimal.From Fig. 9d and f, differences between the SNOW and PRECIP experiments are not readily visible, indicating a relatively small impact of the different accumulation schemes on the simulated ice-sheet thickness and spatial distribution.
Comparing our results with present-day observations (Fig. 9d, f and a), our simulated ice sheet appears too extensive towards the sea.This excessive extent is particularly visible in the northeastern and southwestern parts where no ice currently exists in observations.There is also slightly too much ice over North America where a 750 m thick ice sheet is present over Devon island.
Using ice-thickness anomalies with respect to the observations for the two modelling setups (Fig. 10a and b), we observe an excess of ice of 500 m in central Greenland, reaching up to 1000 m in the northeast.The western part of the Greenland ice sheet is much more consistent with observations.
Analysing further the differences between the PRECIP and SNOW experiments (cf.Fig. 10b), we infer that the two accumulation treatments yield differences of a few hundred metres at most in ice-sheet thickness.Compared to the SNOW experiment, the PRECIP experiment produces a smaller ice-sheet thickness in northern Greenland and Baffin islands and a thicker ice sheet in southern Greenland.
To investigate further the causes of this overestimation in simulated ice-sheet volume, it is useful to turn to the offline and semi-offline ice-sheet runs.The CLIMICE experiment already shows a significant overestimation (Fig. 8) of the icesheet volume by 4×10 14 m 3 , a little less than half the discrepancy between the SNOW and PRECIP experiments and the observations.Three main factors contribute to this discrepancy: (1) the use of the simplistic PDD approach of the SMB ;(2) the fact that the CLIMICE ice sheet is equilibrated under constant climatic conditions for 50 kyrs, which is not the case for the real ice sheet, and (3) simplified dynamics in the GRISLI ice sheet with respect to reality.In any case, since the climatology of the iLOVECLIM model has biases with respect to observations, we cannot expect it to perform better than the CLIMICE results in our setup.Spatially (Fig. 9b), CLIMICE is too thick in the southern part of Greenland and too extensive in the southwestern part, while showing an underestimated extent in the northern part.
Comparing the semi-offline simulations (PRECIP-SO and SNOW-SO) to CLIMICE (Fig. 8), we again observe an additional overestimation of simulated ice volume of 4.4 × 10 14 m 3 , with no significant difference between SNOW and PRECIP.Since the only difference between PRECIP-SO and SNOW-SO on the one hand and CLIMICE on the other is the use of the iLOVECLIM climatology instead of the observations, we can thus conclude that the biases in the simulated climate in our model result in a thicker ice sheet.In other words, there is too much accumulation and insufficient ablation.The spatial shape of the ice sheet is very similar to the SNOW and PRECIP coupled simulations (e.g.Fig. 9c and  d).This indicates that the mean climate is the main factor in shaping the ice sheet in our simulations.
Finally, the dynamical ice-sheet feedbacks on climate tend to again overestimate the ice volume, by an additional 1.3 × 10 14 m 3 of ice (Fig. 8).From this analysis, we can therefore conclude that the simulated climate biases contribute 46 % to the observed differences in ice-sheet volume, the equilibration and coupling method effect 40 % and the feedback of having the simulated ice sheet affecting the climate accounts for the remaining 14 %.

Simulated surface mass balance
Since the obtained simulated ice sheets are the product of the surface mass balance from the beginning of the simulation, it is useful to deal with the simulated SMBs right at the start of the SNOW and PRECIP experiments and to compare them to high-resolution SMB derived from a regional model forced by reanalysis climate data.The latter is required since there is no direct observation of the SMB.The simulated positive SMB is obviously overestimated in central Greenland and in most of the southern part of Greenland.Though the overall pattern of more positive SMB values in the south and less in the north and negative SMB on the coast is broadly reproduced in our model setup, the positive SMB values are clearly overestimated.It is thus logical that the equilibration of the GRISLI model with such an SMB yields an overestimation in ice volume.We further analyse the origin of this overestimation by comparing the simulated accumulation and temperatures to observations.The SNOW setup shows larger areas of negative mass balance in the coastal areas, in better accordance with observations.Their effect is seen in Fig. 8, where SNOW has consistently less ice volume than the PRECIP setup until about 12 000 years.

Simulated accumulation
Differences in accumulation between the two model setups cause differences in accumulation from the very beginning of the simulations.Furthermore, these changes in shape create some further changes in the accumulation pattern.Therefore, to analyse the sole effect of the two model setups without the ice-sheet dynamical feedbacks, it is useful to compare the two accumulation fields (SNOW-and PRECIP-type) from the CTRL experiment (where the ice sheet is fixed to observed present-day conditions) to the ones obtained at the end of the PRECIP and SNOW experiments.The differences in precipitation (in %) are displayed in Fig. 12, with the experiments conducted with a fixed ice-sheet in panel a and the runs that re performedwith an interactive ice-sheet in panel b.
Comparing the accumulations for the CTRL (with an SMB from the SNOW setup) and CTRL(with an SMB from the PRECIP setup) reveals that the results over Greenland are very similar with less (overall 10 %, locally 30 %) accumulation for the PRECIP setup.Conversely, the same computation of snow from the total precipitation (PRECIP) tends to increase accumulation on the southern border of Greenland, where the topography is steep and the mean temperature close to the freezing point (see hereafter).In this region, the higher resolution of GRISLI allows for snowfall whereas ECBilt mainly computes rain.The fact that the total precipitation is taken into account by GRISLI as either snowfall or rain is seen south of 70 • latitude as the CTRL(-PRECIP) run displays up to 150 % more accumulation falling as rain, which is not included in the CTRL(-SNOW) experiment.
We further analyse the differences in accumulation between the PRECIP and SNOW experiments, presented in Fig. 12.There is more accumulation in the PRECIP experiment south of 75 • latitude and on the eastern and western sides of the Greenland ice sheet (1.25 × 10 12 versus 0.93×10 12 m 3 yr −1 ).This is expected since the downscaling of accumulation helps to take into account the height differences between ECBilt and GRISLI.Similarly, there are very small differences in central Greenland where ECBilt predicts a high ice sheet already and where the downscaling thus does not result in much information.
So far, we have concentrated on the differences between our simulations.However, since we overestimate the icesheet extent under pre-industrial conditions, it is useful to analyse the accumulation patterns with respect to a presentday climatology.Figure 13 displays the difference in accumulation between the climatology (a) and the CTRL (b and c), the SNOW (d) and the PRECIP (e) experiments.A pattern that is common to all panels of Fig. 13 is the overestimation of accumulation in central Greenland and centred on Devon Island, up to northern Baffin and southern Ellesmere islands.Moreover, all panels show an underestimation of accumulation in northwestern and in southern Greenland, except for the PRECIP experiments in the case of the latter.These common features thus originate from the ECBilt model itself and not from the coupling procedure.As noted before, the downscaling procedure for the accumulation in PRECIP helps to reduce the discrepancies in southern Greenland.All together, the overestimation of accumulation in central Greenland seen in SNOW and PRECIP is certainly one of the causes of the overestimation of the size of the simulated ice sheet.
We note from our analysis that our model is unable, due to its simplification, to reproduce the very high accumulation in southern Greenland, linked to oceanic moisture advection over the cold and high-altitude ice sheet, nor the extremely dry conditions pertaining to central Greenland.At ECBilt resolution, Greenland is never dry enough in the high-altitude regions and moist advection is too widespread in the interior of Greenland.From the accumulation pattern, it is difficult to choose between the PRECIP and SNOW experiments.

Simulated temperature fields
Temperature is an important governing factor for the surface mass balance of the ice sheet.In the CTRL configuration, iLOVECLIM does not exhibit a coherent, systematic bias over the whole Greenland area (Fig. 14a and b).The mean annual temperatures are generally within ±2 • C of the climatological value, except for specific regions.Those include regions with high topography over a small spatial extent (overestimation of temperature by 5 • C in southern Greenland) and the sides of the ice sheet where the altitude varies a lot over a short distance (underestimation of 5 • C on the western and eastern flanks).Yet, taking into account that the climatology consists of present-day data and that it was probably about 2 • C colder during the pre-industrial times that we are simulating (Kobashi et al., 2011), the model performance for Greenland is reasonable given its low spatial resolution.Biases discussed above are of the same order of magnitude as the ones obtained with low-resolution general circulation models (Quiquet et al., 2012).
In the PRECIP and SNOW experiments, there is a common pattern of cooler conditions than climatology of about 2 to 4 • C in central Greenland.In northern Greenland, the cooler bias is even more pronounced in the PRECIP experiment (Fig. 14d).These differences can be readily explained by the large overestimation of the ice-sheet thickness in the SNOW and PRECIP simulations, causing a higher elevation.The overestimation of 500-800 m, already discussed, with respect to the observed ice sheet causes an annual mean cooling of 2-4 • C by altitudinal lapse rate effect, in very good accordance with what we observe in the simulated temperature.The slightly higher elevation of the ice sheet in northern Greenland causes an additional cooling in the PRECIP experiment.
We thus may conclude that from the mean temperature perspective, the SNOW experiment is in better accordance with observations, though this result is achieved through compensation of a warm surface bias in the CTRL by an anomalously high simulated ice sheet.In south Greenland, the temperatures are overestimated in both the PRECIP and SNOW by up to 5 • C, causing more ablation and a lower icesheet thickness than observed (Fig. 10a).
The ice-sheet mass balance is very sensitive to the temperature of the melt season, as expressed by the formulation of the PDD that relies on mean annual and July temperatures.Though we use the complete computed seasonal cycle here, it is instructive to compare our simulated summer mean temperatures to the climatology.Figure 15 presents such a comparison for July.The first striking feature is the large overestimation of the temperature in the CTRL simulation over Greenland and the adjacent Greenland-Iceland-Norwegian seas with up to 15 • C differences.Over Baffin Island, the opposite pattern is observed.In the SNOW and PRECIP experiments, the overestimated altitude of the Greenland ice sheet tends to reduce this bias to approximately 2 • C, with an opposite sign in the south and in the north.The PRECIP and SNOW simulations again show a very similar pattern.Part of the latter mismatch is due to the lack of energy coupling between ECBilt and GRISLI in the present version of the coupling approach.Temperatures in the atmospheric model should be buffered by the presence of the ice sheet (through the take-up of latent head needed to melt the ice), a process absent in the current version of our coupling procedure.

Discussion
By running two fully coupled climate-ice-sheet simulations with different assumptions for the accumulation scheme, our goal was to evaluate whether a relatively coarse-resolution atmosphere model would yield a reasonable ice sheet with respect to observations (Bamber et al., 2001) and what the differences between the two schemes and the CTRL, uncoupled simulation would be.
It should be noted that running 14 000 yr under equilibrated climatic conditions is a very unlikely scenario and not fully comparable to the present-day conditions.In fact, though the climate has been relatively stable for the 7000 yr preceding the industrial era, there were still some climate changes that triggered some of the evolution of the Greenland ice sheet, with a general thinning of the ice sheet over time (Vinther et al., 2009).The ice sheet we observe today is thus not fully in equilibrium with climate and much less so since global warming has started.It is also dependent on the complex climatic history extending back to the last glacial period.Another aspect to look at is the spatial scale of the icesheet dynamics itself.The precise ice-sheet extent of presentday Greenland is shaped by very small-scale processes like fast-flowing glaciers in localized valleys and by very local effects that directly influence the SMB.We cannot expect to represent such small spatial scales, let alone in the GRISLI ice-sheet model at 40 km resolution.It is even less possible to do so in a T21-resolution atmospheric model.
Considering these shortcomings inherent to our modelling effort, what is the result of our coupling process?We have demonstrated that the use of a simple atmospheric component and a simple downscaling method without introduction of an anomaly mode in the coupling approach allows the simulation of an ice sheet in Greenland under pre-industrial conditions and even of some of the smaller ice sheets in the Baffin area.The simulated ice sheets are too thick due to a large overestimation of the accumulation, but this does not result in the start of a hemispheric-scale glaciation that could be triggered through ice-feedback processes.The basic elements of the pre-industrial ice sheet are present and further refinement of the mass balance calculations will certainly help achieve a more realistic simulation of the present-day ice sheet.In particular, the energy consistency between the atmospheric and ice-sheet models is not achieved: the bilinear interpolation of temperature, though providing the necessary heat to the SMB of the ice sheet, does not dynamically modify the energy fluxes in the atmospheric model.A coupling based on energy conservation is still to be developed and will be the subject of future studies.
From a climatic point of view, the coupling of GRISLI to ECBilt provides cooler temperatures which seem to be in better agreement with pre-industrial temperature reconstructions (Kobashi et al., 2011) and seem to fit climatology better than the results of the CTRL experiment.However, this result is a consequence of the large altitudinal bias of the simulated ice sheet that compensates for the warming observed in this region at the CTRL.
Regarding the two accumulation techniques, both versions yield too thick an ice sheet, owing to overestimated precipitation in central Greenland.The small differences induced by the downscaling of snow do not result in thickness differences of more than 200 m locally, and mostly less than 50 m.Only the southern tip and the Arctic coast of Greenland behave differently with altitudinal differences of up to 500 m between both at the end of the simulation.As mentioned previously, it is difficult to compare our simulated ice sheet to observations since the observed ice-sheet extent and thickness is the result of climate history over the last glacial cycle.Comparing our results to an offline run forced by a climatology (CLIMICE), we found that, using the coupled iLOVECLIM model, forcing fields yield a volume bias of 15 % with respect to CLIMICE, or 30 % compared to observations.While these figures cannot be claimed to be in good accordance with observations, they are comparable to those found using Climate Model Intercomparison Project, phase 3 model outputs (Lipscomb et al., 2013) and regional climate models (Quiquet et al., 2012) when integrating the ice sheet to equilibrium.

Conclusions
We coupled an intermediate-complexity climate model with an ice-sheet model, including the exchange of water, topography and ice mask for albedo, for the purpose of long-term climate studies.Results of experiments for a pre-industrial equilibrium show a strong dependence of the simulated icesheet distribution and ice thickness to the background climate produced by the atmospheric component of the climate model.Though there are substantial biases in the climatology, we prefer here to keep absolute climate fields to force the ice-sheet model.This is done in order to be able in future studies to consistently assess the response of the ice sheets to past climate conditions, without complications associated with use of an anomaly correction approach to determining temperature, precipitation and associated SMB.The model thus obtained is clearly inadequate for decadal-to centuryscale studies but can be used profitably in large-scale studies on ice-sheet-climate interactions.
Testing different coupling methods, we find that the PRE-CIP and SNOW methods yield very similar ice-sheet thickness in central Greenland where the prevailing cold conditions at both spatial resolutions (ECBilt and GRISLI) always produce snow accumulation.Conversely, the simulated distribution of ice in coastal regions can significantly differ between simulations, thereby altering the shape of the ice sheet.Most notably, the two different iLOVECLIM versions (PRE-CIP and SNOW) tested here do not agree in the ice-sheet distribution at the northern and southern tips of Greenland.Regarding accumulation, the large overestimation of accumulation in central Greenland yields too much ice thickness there.It calls for future development regarding a scheme for precipitation redistribution with altitude that seems to be missing here but also calls for improvement in the model's atmospheric physics to reduce the precipitation bias over Greenland.Finally, we note that the results are relatively insensitive to the method used to compute the fraction of snow from total www.geosci-model-dev.net/7/1377/2014/Geosci.Model Dev., 7, 1377-1394, 2014 precipitation content.The advantage of the SNOW setup is that there is a direct conservation of snow in the atmosphere model and the PDD scheme, since there is no recomputation of the snow fraction.We thus tend to favour the SNOW setup to the PRECIP setup over that reason.The SNOW setup also has a smaller bias in accumulation but this agreement may be coincidental owing to the low resolution of the atmosphere model at T21.

Figure 1 .
Figure 1.Northern Hemisphere surface topography comparison: ECBilt (a) and GRISLI (b).The lower colour scale (in m) indicates the altitude of the topography over non-glaciated areas in GRISLI and everywhere in ECBilt; the upper colour scale shows the altitude of the glaciated areas in GRISLI.The two panels show the restricted area of GRISLI in its Northern Hemisphere configuration.
Figure  1gives an overview of the two model grids that need to be coupled together.Due to the large difference in spatial

Figure 2 .
Figure 2. Number of GRISLI cells per ECBilt cells on the NorthernHemisphere GRISLI grid.The figure shows nicely that the highest number of GRISLI cells per ECBilt cell is reached close to the corners of the GRISLI grid, situated over the oceans in ECBilt under the present-day configuration.The low numbers at the very edge of the grid are due to the partial overlap of the two grids there.

Figure 3 .
Figure 3. Altitude differences on the ECBilt grid of the maximum height of the GRISLI cells contained in the same ECBilt cell for the two reference topographies shown in Fig. 1.Colour scale is given in m.

Figure 4 .
Figure 4. Scheme presenting the method used for coupling the ECBilt, CLIO and GRISLI models, including water exchange fluxes, energy fluxes and timesteps.

γFigure 5 .
Figure 5. Scheme presenting the method used for the vertical downscaling.Numbering indicates the order of processing.(1) the temperature at the highest and lowest GRISLI point (tails of the blue and red arrows respectively) is retrieved for the given ECBilt grid cell (boundary in violet).(2)A vertical lapse rate γ is computed from the these two extreme temperatures, using the line defined by the two temperatures and elevation extrema in addition to the EC-Bilt cell temperature.(3) Using γ , temperatures are derived for all altitudes in GRISLI if the latter are known.

Figure 7 .
Figure 7. Calculated lapse rates within the downscaling procedure, in degrees per km.(a) is for February, (b) for July.Since the computation of the lapse rate involves a linear regression over three temperature points, there is no lapse rate calculated when the correlation coefficient of the regression line is too low: hence the large white areas over North America in panel (a).The orange line at the border is an artifact due to the mask computation and should be ignored.

Figure 8 .
Figure 8. Transient equilibration of the equilibrium runs SNOW and PRECIP: vertical axis is ice volume, horizontal axis is simulation years.For comparison, horizontal dashed lines mark the equilibrium ice volume obtained from observations(Bamber et al., 2001), simulated in the offline simulation forced by the climatology and using the control climatic fields from iLOVECLIM (no ice-sheet feedback).

Figure 9 .
Figure 9. Ice-sheet thickness (in m) for the observations and the experiments performed.(a) Observed Greenland thickness (Bamber et al., 2001) interpolated to the GRISLI grid (b) is the CLIMICE offline simulation, (c) the SNOW-SO semi-offline, (d) the SNOW, (e) the PRECIP-SO semi-offline and (f) the PRECIP simulations.The red contour line corresponds to the observed present-day ice margin.

Figure 10 .
Figure 10.Difference between the observed and the calculated ice-sheet thickness in the SNOW experiment (a) -(SNOW-observed) in m and between the SNOW and the PRECIP experiment (b) (PRECIP-SNOW).The red contour line corresponds to the observed present-day ice margin.

DFigure 11 .Figure 12 .
Figure 11.Comparison of simulated surface mass balance for the Greenland ice sheet.(a): SNOW; (b): PRECIP (both for simulated preindustrial conditions).(c): results from the MAR regional climate model at 25 km resolution forced by reanalyses over the period 1979-1988 from Fettweis et al. (2011).

Figure 13 .
Figure 13.Differences in accumulation between the climatology (a), and the snow accumulation method in the CTRL (b), the precipitation accumulation method in the CTRL (c), the SNOW (d) and the PRECIP (e) experiments.

Figure 14 .
Figure 14.Difference in mean annual temperature between the climatology and the different experiments: CTRL (b), SNOW (c), PRE-CIP (d).The reference climatology is shown in (a).

Figure 15 .
Figure 15.Difference in July temperature between the climatology and the different experiments: CTRL (b), SNOW (c), PRECIP (d).The reference climatology is shown in (a).