Ice sheet dynamics within an earth system model: downscaling, coupling and first results

We present first results from a coupled model setup, consisting of the state-of-the-art ice sheet model RIMBAY (Revised Ice Model Based on frAnk pattYn), and the community earth system model COSMOS. We show that special care has to be provided in order to ensure physical distributions of the forcings as well as numeric stability of the involved models. We demonstrate that a suitable statistical downscaling is crucial for ice sheet stability, especially for southern Greenland where surface temperatures are close to the melting point. The downscaling of net snow accumulation is based on an empirical relationship between surface slope and rainfall. The simulated ice sheet does not show dramatic loss of ice volume for pre-industrial conditions and is comparable with present-day ice orography. A sensitivity study with high CO2 level is used to demonstrate the effects of dynamic ice sheets onto climate compared to the standard setup with prescribed ice sheets.

Abstract. We present first results from a coupled model setup, consisting of the state-of-the-art ice sheet model RIMBAY (Revised Ice Model Based on frAnk pattYn), and the community earth system model COSMOS. We show that special care has to be provided in order to ensure physical distributions of the forcings as well as numeric stability of the involved models. We demonstrate that a suitable statistical downscaling is crucial for ice sheet stability, especially for southern Greenland where surface temperatures are close to the melting point. The downscaling of net snow accumulation is based on an empirical relationship between surface slope and rainfall. The simulated ice sheet does not show dramatic loss of ice volume for pre-industrial conditions and is comparable with present-day ice orography. A sensitivity study with high CO 2 level is used to demonstrate the effects of dynamic ice sheets onto climate compared to the standard setup with prescribed ice sheets.

Motivation
For several decades, earth system models (ESM) of growing complexity have been vastly utilized to study climate dynamics and especially anthropogenic climate change (for a summary and comparison of recent results, see e.g. IPCC, 2013). These setups generally include models for atmosphere and ocean circulation, sea ice and often also chemistry and vegetation modules (e.g. IPCC, 2013). On the other hand, the large ice sheets are often only represented by ruleof-thumb parameterizations, e.g. in the atmosphere model ECHAM (European Centre Hamburg Model) (Roeckner et al., 2003). More complex approaches include, for example, parameterizations for ice sheet melting (Swingedouw et al., 2009) and coupled ice sheet models (Ganopolski and Calov, 2011). The reason for this simplification is that ice sheets are, under normal conditions, evolving on timescales much longer than atmosphere or ocean, with maximum ice velocities of the order of magnitude of 1000 m a −1 and melting rates of less than 0.01 Sv (1 Sv = 10 6 m 3 s −1 ). Thus, a separation of timescales is adequate.
When the climatic evolution leads to large-scale loss of ice volume, this approximation becomes invalid. This is certainly the case for large-scale past climate changes, e.g. glacial-interglacial transitions. Deglaciation periods following ice ages can obviously not be simulated based on static ice sheets. Until now, most studies have prescribed the reconstructed ice sheets (Braconnot et al., 2007;Liu et al., 2009) and the mechanisms for deglaciation are therefore not known.
Furthermore, the contribution of dynamic ice sheets on recent and future sea level rise needs to be addressed.
The breakup of buttressing ice shelves in Antarctica, acceleration of ice streams and outlet glaciers and drainage of the warming ice sheets are intensely argued upon. Indeed, satellite data such as those from the GRACE mission (Sasgen et al., 2012) give indications of accelerated mass loss in southern and western Greenland, in almost equal parts due to increased melting and draining glaciers. For a careful comparison to model results, a coupled setup of an ESM to a 3-D thermodynamical ice sheet model (ISM) is therefore mandatory.
Special care is needed concerning the spatial scales involved. For example, while our atmosphere model ECHAM is operated at the low resolution of T31 (which corresponds Published by Copernicus Publications on behalf of the European Geosciences Union.

2004
D. Barbi et al.: Ice sheet dynamics within the earth system to 3.75 • , i.e. a resolution of 100 (respectively 400) km over the zonal (respectively meridional) directions) in long-term climate studies (e.g. Zhang et al., 2013;Stepanek and Lohmann, 2012;Knorr et al., 2011), the resolution of the ISM needs to be lower than or equal to 20 km for the basic processes of ice sheet dynamics to be resolved. Consequently, the forcing data for the ice sheet, obtained from the atmosphere model, is too coarse to represent realistic distributions of temperature, accumulation or ablation. Our work shows that a careful downscaling procedure can lead to stable ice sheet distributions, which can be non-trivial especially in southern Greenland (throughout the remainder of the article, we use the word "stable" in the sense of "with only small changes of the ice volume").

Model descriptions
We briefly introduce the basic characteristics of the models we use, as far as they matter in the performed experiments. For further information, we refer the reader to the literature listed within the following paragraphs. The coupling scheme is described in Sect. 3.

The ice sheet model RIMBAY
The 3-D thermomechanical ISM RIMBAY was first introduced by Pattyn (2003) and in recent years further developed by Thoma et al. (2010Thoma et al. ( , 2013. The model approach is based on continuum thermodynamic modelling and encompasses balance laws of mass, momentum and energy, extended with a constitutive equation. Treating ice as an incompressible fluid with constant density, the equations for conservation of mass, momentum and energy are iteratively solved. In its most complex version, RIMBAY can solve the momentum balance by means of the full Stokes equations (Pattyn, 2008). By solving the temperature field within the ice, fully coupled thermomechanical simulations are possible. As we aim for an understanding of the large-scale feedbacks of dynamic ice sheets on the climate system, we employ commonly used approximations, namely the shallow ice approximation (SIA) and the shallow shelf approximation (SSA) instead of using the model in its highest complexity. A comprehensive description of SIA and SSA is given by Greve and Blatter (2009).
An important advantage of RIMBAY, compared to other ice sheet models, is the ability to combine different levels of complexity within one model. For instance, a fast flowing stream can demand full Stokes, while in the surroundings, SIA can be adequate on its own. This feature makes RIM-BAY effective and accurate. We use RIMBAY in a horizontal resolution of 20 km × 20 km and apply in the vertical dimension 21 levels of sigma layers, scaling the vertical dimension to dimensionless coordinates. The highest resolution is given at the ice-bedrock interface and near the ice surface.
We use RIMBAY for simulations of the Greenland Ice Sheet. For the boundary as well as initial values, we use information about the geothermal heat fluxes, which we took from Shapiro and Ritzwoller (2004) as well as the bedrock topography and initial ice geometry taken from the bed topography and ice elevation data provided by Bamber et al. (2001a, b). Both compilations provide a 5 km resolution, which we interpolated to a Cartesian 20 km model grid as we found it to be sufficiently accurate for our purposes. The initial temperature distribution within the ice is calculated according to Robin (1955), based on the spatially varying surface temperature and the basal-pressure-dependent freezing point of the ice sheet base.

The earth system model COSMOS
The simulations described in this manuscript have been carried out using the COmmunity earth System MOdelS (COS-MOS) which have been mainly developed by the Max Planck Institute for Meteorology (MPI) in Hamburg. Our version of COSMOS includes the ECHAM5 atmosphere model in T31 resolution with 19 levels and the MPIOM (Max Planck Institute ocean model) in GR30 resolution with 40 levels.
The land-vegetation model Jena Scheme of Atmosphere Biosphere Coupling in Hamburg (JSBACH), an optional part of COSMOS, has not been used in the setup described in this manuscript. Our setup is identical to the COSMOS release, which has been developed in the Millennium project  and runs without any flux corrections. The atmosphere model ECHAM5 is used at T31 resolution (∼ 3.75 • ) with 19 vertical levels and is described in detail by Roeckner et al. (2003). The ocean model MPIOM is a hydrostatic, Boussinesq, free surface, primitive equation ocean and sea ice model (Marsland et al., 2003;Jungclaus et al., 2006). The model dynamics are solved on an Arakawa C-grid (Arakawa and Lamb, 1977). In our model setup, which is identical to the one used by , MPIOM is formulated on a bipolar, orthogonal, curvilinear GR30/L40 grid with poles over Greenland and Antarctica. The advantage of this setup is an increased resolution at many deep-water formation sites (around Greenland of up to 24 km), which facilitates a more realistic simulation of the physical processes operating in these regions. The formal horizontal resolution is 3.0 • ×1.8 • , with the vertical dimension being split into 40 unequally spaced z coordinate model levels. COSMOS has been intensely used for the study of very different aspects of climate sciences, notably for the production of IPCC scenario runs for AR4 and simulations covering the last millennium (e.g. Jungclaus et al., 2010). Our version of COSMOS (COSMOS 1.2.0) has been used for Holocene , glacial (Zhang et al., 2013), Pliocene (Stepanek and  and Miocene (Knorr et al., 2011) conditions. The spectral atmosphere model ECHAM5 was operated in a T31 (3.75 • ) resolution, which is certainly very coarse compared to the resolution of the ISM, but still a good resolution for palaeoclimate simulations as well as for technical testing. As an example, we will show our downscaling procedure later in this work. With ECHAM5, we applied the Hydrological Discharge (HD) hydrology sub-model (Hagemann and Dümenil, 1998) that routes runoff along fixed runoff masks which have to be prepared during preprocessing.
Major changes were already applied to the ice sheet mass balance parameterization of ECHAM/HD. ECHAM ice sheets are stable, meaning that the amount of water collected as snow/rain is exactly compensated by mass loss due to evaporation and runoff. We disabled this parameterization. Instead, we coupled accumulation to RIMBAY, reading back ablation, basal melting and calving, which are then sent from ECHAM5/HD to MPIOM. We also read in ice orography from the ice model output. We furthermore included checks ensuring that the corresponding fields in both ECHAM5/HD and RIMBAY agree on the mass balance of the ice sheets. No changes were applied to the MPIOM ocean model, as meltwater was returned from RIMBAY to ECHAM and from there distributed to MPIOM using the standard coupling routines.

Timing
Regarding the coupling of ice dynamics to the atmosphere and ocean, it is obvious that the processes under consideration are acting on very different timescales. The atmosphere, for example, is integrated using a time step of ten minutes, while an ice sheet will react to changing forcing only in the range of decades to millennia. This naturally led to the paradigm that ice sheets are at least quasi-stable, and was not needed to be simulated at all.
The separation of timescales enables us to couple the models iteratively, meaning that we run the coupled atmosphereocean setup with prescribed ice boundaries, alternating with ice sheet runs with fixed atmospheric forcings. The forcings and boundaries are updated depending on the coupling mode: -Synchronous (transient) coupling: when in synchronous mode, both models are run iteratively but in equal time steps, meaning after each period of atmosphere-ocean integration, we update the forcings to the ice sheets and simulate these for the same period, then transfer the changed ice orography and fresh water fluxes back to the global climate model. We chose the coupling period to be 1 year as, for the time being, we neglect seasonal effects of ice dynamics. The interest in smaller timesteps will rise once we include a direct coupling between ocean and ice shelf, and thus feedbacks from sea ice variations. Annual mean surface temperature Summer mean surface temperature -Asynchronous coupling: while synchronous coupling is certainly the most intuitive way of coupling two models, it is nevertheless not applicable over very long integration times as these would require computation time beyond any reasonable limit. Luckily, as the (computationally expensive) atmosphere exists on timescales much smaller than the internal timescale of an ice sheet, we can choose different integration periods for both models. A sensible scheme would include a rather long integration period for the slowly changing ice sheets (e.g. 1000 years), and a much smaller integration period for the atmosphere part, which we set to 25 years which is sufficient for the atmosphere to adapt to the changed ice orography, given that these variations are small as we start from a present-day ice sheet close to the equilibrium.
In both cases, we use the results of the last year of simulation as forcing for the next invocation of the other model. In asynchronous coupling, especially when a state of quasiequilibrium is reached, it is advisable to exchange multiannual means instead, to compensate for interannual variability.
As the only aim of the asynchronously coupled runs within the scope of this paper is to demonstrate the existence of a stable steady state for the ice sheets, we safely neglect the influence of these interannual fluctuations.

Exchanged fields
The fields exchanged within each coupling timestep are presented in Table 1. At the centre of attention lies the freshwater balance. Net accumulation and ablation are taken from the atmosphere; ice moving across the sheet boundary is defined to be calved and returned, together with melted snow, to the ocean. The coupling scheme is conserving the total mass of the freshwater, which was ensured and controlled by testing routines within each model. We also communicate annual mean and summer mean surface temperature to the ice model. Annual mean temperature is needed as boundary condition to the temperature distribution within the ice; summer mean temperature is needed during the downscaling process (see below).
Changed orography which is represented in ECHAM as a change in surface geopotential height is returned from RIM-BAY together with the albedo.

Interpolation method
The basic interpolation is performed using a Shepard algorithm; for a detailed discussion as well as an implementation example see, e.g. Press et al. (2007). The Shepard algorithm is a distance-weighting interpolation with weight functions of power law type. The value (x) at a point x is therefore calculated as the sum of all values λ(x i ) at the points x i , weighted by a function of the distance: where , λ are the exchanged field (temperature, snow melt or ablation) and the target (respectively, source) grid, and e > 0 is a permissible parameter, which we decided to set to 2.7. In fact, the choice of e has an influence on the outcome of the interpolation as it determines the range within which a given source grid point value is interpolated. When e is chosen larger than 2.7, a grid point is interpolated within a smaller range, leading to an interpolation result with a structure resembling the distance-weighting function, rather than the patchy structure, resembling the T31 grid. With an e much smaller than 2.7, even this structure is lost, leading to an almost homogenous interpolation result, which would be equally undesirable as regional information is lost. The sum is over all points of the source grid; in reality, we limit the sum to the 20 nearest neighbours which were determined in the initial phase of the coupling.

Downscaling
As the resolutions of the grids differ by orders of magnitude (20 km for the ice sheets, 3.75 • for ECHAM), we assume that the spatial integral of the exchanged fields is not the cause of the ice sheet instability, but rather their spatial distribution . To circumvent the unrealistic distribution of precipitation and surface melting, we established downscaling procedures for the coupled variables. Three fields need to be downscaled: net snow accumulation (i.e. precipitation minus evaporation), snow melt and surface temperature. As we redistribute the data from coarse to fine grid, we "generate" information, based on basic theoretical and empirical assumptions. All downscaling methods we employ have in common that mass/latent heat are conserved in the following sense: we calculate the net mass fluxes for accumulation/ablation over the ice as returned by the atmosphere model, and ensure that by multiplying the entire field by a constant factor, the two fields on the different grids agree on this value. In the same way, assuming that heat capacity is constant throughout the model domain (a necessary simplification), we made sure that the net sum of latent heat over the ice is also the same within RIMBAY and in COSMOS.
Moreover, we perform a regionalized downscaling -for the downscaled value at a given grid point, only data from cells within a defined distance (25 grid cells, i.e. 500 km from the given grid point) are considered. This is important, as we do not want to mix information concerning distinct regions of the ice sheet; we especially need to distinguish between northern and southern Greenland. The value of 25 grid cells was established by trial and error -if the value taken is too small, no redistribution of information is achieved; if too big, no regional information is kept.
The distribution of melting snow is difficult to estimate. It strongly depends on the amount of time during which the surface temperature is above the melting point. As we were able to obtain a downscaled, high-resolution distribution of surface temperature, it is straightforward to use a positive degree day (PDD) model (Clyde, 1931;Collins, 1934). In our model, we used the PDD model as implemented by Huybrechts (1993) and Reeh (1991). One assumes that the seasonal cycle varies as a sine function. From the summer and annual means, one can derive the number of days with temperatures above melting point. We furthermore assume that the number of days during snow is melting is directly proportional to the amount of snow melting at this location. From the constraint that overall ablation is conserved, we get the downscaled ablation pattern.
A suitable formulation for the downscaling of net snow accumulation is based on an empirical relationship between surface slope and rainfall (Fortuin and Oerlemans, 1990). Therefore, we redistribute snowfall proportionally to the surface gradient of elevation. Within a box of 50×50 grid points surrounding the grid point x, we search for the highest and the lowest absolute value of gradient δ min and δ max (neglecting directions), and for the highest and lowest value of accumulation α min , α max . If δ(x) is the absolute value of the gradient at x, we estimate the accumulation α(x) linearly as The result is corrected for the whole Greenland Ice Sheet by multiplying the resulting local accumulation by a constant factor over the ice sheet so that the total accumulation is conserved.
Mean annual and mean summer temperature are downscaled similarly. Instead for the surface gradient, we go for surface elevation, assuming that it is colder on higher elevations than further down. Calling the temperature T and surface elevation h, we thus use: We tried the classical lapse rate correction, but found that a linear redistribution completely analogous to the method described for accumulation gives rather better results (not shown). A combination of the lapse rate correction and our linear interpolation method might be the best option. We will implement and try this in a new version of the coupled system which is currently under development. Nonetheless, what we try to show is not that we have the optimal solution, but that a simple approach as the one demonstrated here can already lead to meaningful results. Again, we correct the result so that the global surface integral of temperature is conserved. If one assumes that the heat capacity is constant (a simplification, of course), then heat is also conserved.

Performed experiments
We conduct a series of experiments, first of all for testing the coupling scheme. As the initial ice sheet distribution is taken from present-day data, we expect the coupled system to be out of equilibrium, and to reach a steady state only after an initialisation phase.
The initial state for COSMOS of our integrations was taken from the results from a 800-year equilibration run, which was forced with pre-industrial greenhouse gas concentrations . The initial state of the COSMOS part of our setup is thus comparable to a pre-industrial climate.
In two sets of experiments, we applied two different values for the CO 2 -concentration, while keeping all other forcings identical to the pre-industrial climate. The reference experiments are performed with a CO 2 -concentration of 278 ppm. The results of these runs are analysed in order to validate the long-term stability of the climatic state. In a second phase, we use a high CO 2 -level (999 ppm) to demonstrate the effect of ice dynamics in a scenario with large-scale ice melt. We will compare the results to a standard COSMOS setup with prescribed ice sheets. An overview of the different experimental setups is presented in Table 2. Before the coupled system can be reliably operated, checks against drifts in the climate-ice sheet model are required. These tests are supposed to indicate mistakes and invalid approximations in the numerical implementation. Among other indices, we decided to analyse global mean surface temperature as an indicator for the climate state, and ice sheet mass balance (both locally and summed) which, as will be shown in the next paragraph, led to a reassessment of our coupling strategy.
All tests are performed starting from present-day ice thickness distributions and a well-established atmospheric circulation not subject to global drifts using pre-industrial conditions. To assess substantial changes in the ice sheet, we operated the setup in asynchronous coupling mode (Experiments 1a and b).
The more detailed values for ice sheet mass balance as well as surface temperature were then taken from a comparable, but transient run (Experiment 2).  Table 2) and after (b, Exp. 1b) the downscaling procedure, after 10 000 ice model years. This is compared to measurement data (c, Bamber et al., 2001a).

Impact of downscaling -Experiments 1a and b
Refining the resolution makes a large difference in the results. Figure 1a shows the resulting ice distribution of Experiment 1a without downscaling and after 10 000 model years. We compare the results for ice thickness of observational data from the Bamber et al. (2001a, b) data set. These ice elevation data were derived from ice-penetrating radar data collected in the 1970s and 1990s, and re-gridding to a resolution of 5 km. Ice mass is lost dramatically in southern Greenland compared to observations, much more than realistic under the chosen forcing.
The relevant forcing variables (net snow melt, snowfall, surface mass balance and surface temperature) are given in Figs. 2a, 3a, 4a, and 5a. The maps show the forcings as they are seen by the ISM, resulting from simple interpolation from the ECHAM-grid. The snow ablation is most problematic (Fig. 2a). The main mass loss is concentrated at the southern tip of Greenland.
In contrast, Figs. 2b, 3b, and 5b show the downscaled forcing fields. Accumulation (Fig. 3b) and surface temperature (Fig. 5b) can directly be compared to the respective observation data sets (Bamber et al., 2001a, Figs. 3c and 5c). In the case of ablation rates, detailed observational data is not available; the reader is be referred to Fig. 2c (Ettema et al., 2009) for an estimate taken from model simulations. The distribution of the ablation rates in Fig. 2c clearly differs from the result of our downscaling procedure (Fig. 2b). The application of the PDD model yields additional ablation, neglecting the predominant wind effect along the west coast of Greenland. However, we decided to keep the PDD model as including the wind directions and speeds into the downscaling would ultimately lead to the nesting of a regional model, which is beyond the scope of this study. Nevertheless, the downscaled  . Ablation pattern of the initial forcing over Greenland before (a, Exp. 1a) and with (b, Exp. 1b) the downscaling procedure. For reference, we also plot data from a regional climate model (c, Ettema et al., 2009). The scale is logarithmic; this is also the cause for the white area in plot (c) as the values are too close to zero to be represented on this scale.
ablation rates lead to better results in the ice sheet area than unscaled ones. The downscaled accumulation pattern (Fig. 3b) indicates redistribution from the north-western to northern Greenland at high altitudes. The low observed accumulation area in north-eastern Greenland is not resolved, but the patchiness of the southern Greenland high accumulation zone of ECHAM output is reduced in favour of more coastal marginal expression of the downscaled field. This holds also for the temperature field (Fig. 5b), where the separated low temperature patterns in northern Greenland are much more smoothed out into a low temperature field at high altitudes. Even if the temperature pattern has improved compared to observations, the overall temperature is still too high. This is clear, of course, as the space integral of temperature over Greenland given by ECHAM is too high and held constant throughout out downscaling procedure. If the source model has an overall bias, it cannot simply be fixed by a downscaling procedure. The downscaled accumulation pattern can be seen in Fig. 3b.
We conclude this discussion by showing the surface mass balance (SMB) (Fig. 4). The plot clearly demonstrates that the downscaling scheme changes the regional distribution of the surface mass balance, shifting some of the mass loss formerly concentrated on the southern tip of Greenland along the western margin and to the north. Nevertheless, the full ECHAM input does not vanish so that small localized regions of negative mass loss remain, instead of the observed narrow net mass loss zone along the west coast of Greenland.
Overall, Fig. 1b demonstrates that downscaling leads to a more stable and more realistic representation of the Greenland Ice Sheet; again, we present the result of 10 000 model years of ice sheet simulation (Exp. 1b). Our result is very similar to observational data (Fig. 1c), although showing more ice in most places than observed. This is also expected as we simulate under pre-industrial climate forcing conditions but compare to the present-day ice distribution; to check for a possible additional model bias, further experiments, e.g. transient simulations of the late Holocene, would be required.

Ice sheet mass balance -Experiment 1b
We present the details of the mass balance for the Greenland Ice Sheet in Table 3. All values are within a reasonable range and comparable to literature values. Net values from observations and Regional Atmospheric Climate Model (RACMO) data indicate a higher mass loss, which is not surprising as our setup is forced with pre-industrial atmospheric parameters.
Up to about 40 % of mass loss is caused by surface melting, and about 60 % by calving -values that are in the same range as the GRACE results (Sasgen et al., 2012). The net mass balance is close to zero, meaning the ice sheet is stable. As this is the difference between two almost equally sized numbers, the relative error is huge.

Surface temperature -Experiment 2
Another check for inconsistencies and drifts is to observe the global mean surface temperature. We used a synchronously coupled setup, exchanging fields every COSMOS/RIMBAY model year over 400 years (Exp. 2). The change in global mean surface temperature of the asynchronously coupled run (Exp. 2) can be seen in Fig. 6. Apart from the first 70 years at the beginning of the run (which can be interpreted as the adaptation phase to the changed boundary conditions of the global circulation model, GCM) the temperature remains on the same level, without any visible drift. We see that the change in temperature due to ice sheet coupling is negligible during the integration period and, moreover, the system seems to run into a statistically stable limit.

Sensitivity study: ice feedback on ocean (Experiments 3a and b)
In the following, we present a sensitivity study on the effect of the ISM on the ocean. All results are taken as a mean over the last 10 years out of 350 simulation years (Table 2). Table 3. Mass balance of the Greenland Ice Sheet, all values in Gt a −1 . Model values were taken from Experiment 1b as means over the last 1000 of 10 000 model years. For calving and net balance we also present the standard deviation as derived from the simulation. For accumulation and ablation these would be meaningless as the last 1000 model years are forced with the same atmospheric fields (see Table 2). Nevertheless, we can estimate the standard deviation over the whole 10 000 model years and get ±71 Gt a −1 for accumulation and ±45 Gt a −1 for ablation. For comparison, we also present observational data from GRACE (Sasgen et al., 2012), as well as surface mass balance (SMB) model data from RACMO (Ettema et al., 2009), in which discharge D from IceSAT (Rignot et al., 2006) is taken into account (values taken from Sasgen et al., 2012).

Accumulation
Calving Ablation Net GRACE (net) SMB-D (net) 1149 −820 ± 57 −481 −152 ± 57 −238 ± 29 −260 ± 53   Figure 7 also shows that in the reference run, a coupled setup using preindustrial greenhouse gas concentrations, the total ice volume is slowly growing. At first sight this may contradict the ice loss visible in Table 3, but the experiments are quite different -not only the total length of the run but also the coupling intervals are very different. Figure 9 shows the anomaly of sea surface temperatures at the end of the simulation, compared to results taken from . In the northern North Atlantic, surface water is colder when a dynamic ice sheet is present, compared to the setup with prescribed ECHAM ice sheets (except in the region of the Barents/Kara seas). This large-scale cooling is parallel to a relative freshening of surface water (Fig. 10), with values of about −1 psu in the North Atlantic (which corresponds to a freshwater inflow of about 0.6 Sv if mixing of the top 700 m is assumed), and positive values in the South Atlantic. Part of the freshening can be attributed to Greenland melting and additional freshwater release into the North Atlantic. Another part is due to less salinity transport related to a weakening of the meridional overturning circulation (Fig. 11) transporting less warm and saline water to the North Atlantic. The different contributions of the freshwater fluxes (atmosphere, ocean, ice) and their feedbacks will be the topic of a subsequent paper. In an earlier paper, we found almost equal contributions of the atmosphere and ocean to North Atlantic freshwater transport feedbacks (Lohmann, 2003). In future work, we will study such feedbacks in the coupled COSMOS-ISM model to be compared to the model results with static ice sheets and a moderate weakening of the Atlantic meridional overturning circulation (AMOC).

Discussion
We coupled the dynamic ISM RIMBAY to the global ESM COSMOS. The results clearly demonstrate that a coupled setup is possible, and potentially provides important corrections to simulations with fixed ice sheets, as shown in Sect. 4.4.
The details of the implementation are, however, far from trivial. Concerning the coupling scheme, we tried synchronous as well as asynchronous coupling. Both methods have their benefits depending on the scenario under investigation. If interested in non-equilibrium simulations, as e.g. past, present and future scenarios and projections, where the ice sheet can be far from an equilibrium state, a synchronous coupling with a short exchange period (∼ 1 year) is necessary. When it comes to long-term simulations, like stability tests, an asynchronous scheme is the only option available regarding the huge demand on resources imposed by the ESM. Regardless of the synchronization scheme, we use an offline coupling strategy, exchanging data between separate calls of the respective models. This surely has its disadvantages -an integrated version running both components simultaneously would be faster, and, in the long run, also easier to maintain. In principle, it would be desirable to implement a parallel version where the ice sheet model runs as a separate process. We decided against this option as at the time our project started, our ISM development was still rapidly evolving and we tried to keep porting issues as small as possible. We are planning to reorganize our implementation for future projects accordingly.
As the spatial scales of the ESM and the ISM are very different, we use a downscaling procedure to redistribute atmospheric forcing to the ISM. Our semi-empirical approach proves to be sufficient for the production of stable ice sheets comparable to present-day orographies. A refining of these methods would be interesting, especially when regions sensitive to climate disturbances are to be studied (outlet glaciers, shelf breakup). Furthermore, the influence of wind speed and direction have been completely neglected so far. For instance, Robinson et al. (2010) studied the use of a simple regional energy-moisture balance model (REMBO) to downscale ERA40 reanalysis data for ice sheet forcing. While their model is still very basic (the authors themselves refer to it as a downscaling procedure rather than a regional climate model), it nevertheless is considerably more complex than our method, including a set of fitting parameters we are able to avoid. The inclusion of all relevant processes would ultimately lead to the nesting of a regional climate model, or a refinement of the global ESM as is possible with some sophisticated circulation models (Li, 1999;Krinner and Genthon, 1998). In the long run, this is also one route of development we will pursue as it represents the most consistent approach, avoiding the introduction of additional assumptions and simplification. On the other hand, as we are mainly interested in long-term simulations of the global climate, as in palaeoclimate simulations, for the time being we can restrict our downscaling method to the most influential effects.
Up to now, we have not changed the land/sea mask dynamically, which would be required for long-term palaeoclimate applications. A changing sea-mask would clearly be desirable, as are a number of improvements to the described coupling scheme. Our aim here is to present the effects caused by a suitable downscaling procedure. Equally important is the implementation of moving margins, i.e. sheet-shelf boundaries (grounding line), ice-ocean boundaries (shelf front) and ice-rock boundaries (Nunataks). These play a crucial role in the study of ocean feedback on shelves, which have been shown to possibly induce rapid shelf breakups (Hellmer et al., 2012;Determann et al., 2012).
In other approaches, rather elaborate downscaling methods were used to couple ice sheet models to GCMs. Vizcaino et al. (2010) have shown that an energy balance scheme can sometimes give better results than the PDD; Ganopolski and Calov (2011) introduced an interface for the energy and mass balance in an additional snow layer. We decided to add as few assumptions to the original models as possible, keeping the downscaling rather basic and technical. The reason is that we did not want to write improvements to well-established and reliable models, which were mostly left untouched, but to enable the exchange of the most influential fields as a first step of a modular ice sheet coupling.

Conclusions
Modelling the earth system evolution with global earth system models like COSMOS strongly depends on the inclusion of ice sheet models. It is possible to couple dynamic ice sheets into the ESM COSMOS, even though time and spacial scales of the involved dynamics differ widely. We demonstrated that we can reproduce realistic and stable ice sheet orographies when forcing the coupled system with pre-industrial greenhouse gas concentrations. The simulated Greenland ice elevation is close to the observed one except at the centre of the ice sheet, where it is too thick.
The difference in spatial scales forces us to a careful downscaling procedure, as we have clearly shown. This is especially crucial in ice regions close to the melting point, where small changes in the temperature and/or ablation rates can lead to drastic mass losses, as could be seen in the example of southern Greenland. Concerning the downscaling method, the surprising result is that a very basic procedure is sufficient to provide stable ice sheets. A more elaborate downscaling, e.g. a nested regional model, might have many advantages when small regions sensitive to climate change are investigated. Within the context of long-term simulations of global climate, a semi-empirical scheme, as presented here, already leads to reliable results.
Unstable ice sheets can have a strong impact on climate. We have shown in a first sensitivity experiment that a global warming results in melting in Greenland which in turn leads to a freshening and cooling of surface waters in the North Atlantic and may influence global overturning circulation.