Coupling of climate models and ice sheet models by surface mass balance gradients: application to the Greenland Ice Sheet

It is notoriously difficult to couple surface mass balance (SMB) results from climate models to the changing geometry of an ice sheet model. This problem is traditionally avoided by using only accumulation from a climate model, and parameterizing the meltwater run-off as a function of temperature, which is often related to surface elevation ( Hs). In this study, we propose a new strategy to calculate SMB, to allow a direct adjustment of SMB to a change in ice sheet topography and/or a change in climate forcing. This method is based on elevational gradients in the SMB field as computed by a regional climate model. Separate linear relations are derived for ablation and accumulation, using pairs of Hs and SMB within a minimum search radius. The continuously adjusting SMB forcing is consistent with climate model forcing fields, also for initially non-glaciated areas in the peripheral areas of an ice sheet. When applied to an asynchronous coupled ice sheet – climate model setup, this method circumvents traditional temperature lapse rate assumptions. Here we apply it to the Greenland Ice Sheet (GrIS). Experiments using both steady-state forcing and glacial-interglacial forcing result in realistic ice sheet reconstructions.


Introduction
Ice dynamical models are valuable tools to study the response of ice sheets to climate changes, and hence to constrain their contribution to observed sea level changes.Various ice sheet model experiments have been carried out for Greenland, to reconstruct ice sheet volume on time scales ranging from centennial to glacial-interglacial cycles (e.g.Huybrechts et al., 1991;Letréguilly et al., 1991;Huybrechts, 1994;Marshall and Cuffey, 2000;Lhomme et al., 2005;Otto-Bliesner et al., 2006;Graversen et al., 2010;Robinson et al., 2011;Fyke et al., 2011).How a certain climate record is translated to surface mass balance (SMB) forcing is vital for the outcome of such experiments (e.g.Letréguilly et al., 1991;Robinson et al., 2011).
Ice sheet SMB is the complex net result of accumulation in the interior, and melt, refreezing and subsequent run-off at the margins.Accumulation strongly depends on atmospheric circulation, and changes in ice sheet elevation and extent.In Greenland the present-day accumulation pattern is reasonably well constrained by measurements (Bales et al., 2009) and regional climate modelling (Box et al., 2006;Fettweis et al., 2008;Ettema et al., 2009), though uncertainties remain large in areas where measurements are sparse.Melt is a function of the surface energy balance components, which vary widely in space and time over the ice sheet (Ambach, 1977a,b;Duynkerke and Van den Broeke, 1994).Part of the melt water is retained or refreezes in the firn layer, limiting run-off (Pfeffer et al., 1991).
In an ice sheet model, assumptions and simplifications are unavoidable in the translation of a time-dependent climate record to spatially and temporally varying forcing fields of surface temperature (T s ) and SMB.The classical solution is to separate SMB into accumulation and run-off and estimate both fields separately (e.g.Letréguilly et al., 1991;Huybrechts, 1994;Ritz et al., 1997).Snow accumulation is then prescribed from either a compilation of measurements (e.g.Ohmura and Reeh, 1991;Bales et al., 2009), atmospheric reanalysis data (e.g.Uppala et al., 2005;Robinson et al., 2010), time slice products (e.g.Kiehl and Gent, 2004;Huybrechts et al., 2004;Otto-Bliesner et al., 2006) or scenario runs (e.g.Graversen et al., 2010) from global climate models.To account for changes in precipitation, sometimes a thermodynamic scaling of the accumulation is applied as a function of temperature.Melt is usually calculated from positive degreedays (PDD) (Braithwaite and Olesen, 1989;Reeh, 1991), and run-off is estimated based on assumptions on superimposed ice formation.However, the use of a PDD-model to derive Published by Copernicus Publications on behalf of the European Geosciences Union.
ablation generally leads to overestimation of the climate sensitivity ( Van de Wal, 1996), due to spatial and temporal variability of the degree-day factors (Van den Broeke et al., 2010) and not explicitly accounting for changes in e.g.lapse rates and albedo feedbacks in a transient climate (Bougamont et al., 2007).In view of the complexities that determine the spatial and temporal evolution of SMB, its inclusion in ice sheet models needs improvement.
It would be ideal to have a fully coupled ice sheet-climate model system, but such a coupling set-up is still not feasible due to the length of the required model simulations (kyr): ice sheet model experiments require at least a few millennia of evolution, whereas climate models are typically used for several decades of climate reconstruction.Given these different time scales, asynchronous coupling strategies (e.g.Charbit et al., 2002) are required.Next to that, downscaling techniques (e.g.Robinson et al., 2010;Vizcaíno et al., 2010) have been developed to translate climate model fields (often available on a lower resolution than ice sheet model grids) into useful forcing fields for ice dynamical models.However, SMB is usually not a product of climate models and hence a parameterized calculation of SMB is still required.All in all, simulations of ice sheet-climate interactions will strongly benefit from an unambiguous calculation of SMB in a coupled atmosphere-ice sheet model.
Here we introduce a novel approach to force an ice sheet model in which SMB fields of a regional climate model are used directly, which circumvents assumptions regarding the calculation of run-off.The method (Sect.2) (1) enables to apply and reproduce fields of SMB to a dynamic ice sheet topography, and (2) allows for an elevation -SMB adjustment while the ice sheet topography evolves, the latter is the main novelty of this work.For this particular application, we use SMB fields from the regional climate model RACMO2/GR, a product that realistically simulates the SMB observations (Ettema et al., 2009).In principle, this method is designed for use in asynchronous coupled ice sheet -climate simulations, but here we apply it to a set of ice sheet model experiments using a single climatology (Sect.3).In Sect. 4 we evaluate our method and conclusions are drawn in Sect. 5.

Methods
When imposing a change in SMB on an ice sheet model, the simulated ice sheet will instantly advance or retreat, thereby changing its areal extent and surface elevation (H s ).The modified ice geometry will feedback on the SMB pattern via changes in temperature, atmospheric circulation, orography, albedo, etc.Therefore, using a fixed SMB field as a forcing for an ice sheet model is not realistic for simulations longer than several decades.However, similar to e.g.surface temperature (T s ), a lapse rate could be used to correct SMB as a function of changes in H s .In the ablation zone we expect the SMB to become less negative with increasing H s , but the rate of change will depend on the partitioning of the surface energy balance during melt, i.e. the sum of net short-and net longwave radiation, and the sensible, latent and subsurface heat fluxes (e.g.Duynkerke and Van den Broeke, 1994;Van den Broeke et al., 2008, 2011).Here we assume that this variability can be accurately predicted by the local SMB gradient, instead of making it a function of 2 m temperature as in PDD models.In the accumulation zone the behavior of SMB as a function of H s is even less predictable: SMB can increase due to less ablation, but moving further into the interior SMB will decrease due to decreasing precipitation.To account for these regional differences, we use a regional H s -SMB relation that can be applied locally.This way we can account for spatial variability in the relation between H s and SMB; the method also allows us to predict SMB values at locations that become ice covered, but are currently outside the ice mask.

Spatial mass balance gradients
Here we use 1958-2007 average SMB from the regional climate model RACMO2/GR (Fig. 1, Ettema et al., 2009) forced at its lateral boundaries by 6-hourly fields from the global model of ECMWF, using ERA-40 reanalysis (Uppala et al., 2005) up to September 2002 and the operational analysis thereafter.The SMB field is calculated using a physical snow model that accounts for refreezing of meltwater and the thermodynamic evolution of the upper snow/firn/ice layers, yielding a realistic reproduction of the present-day SMB distribution with a high horizontal resolution (11 × 11 km, Ettema et al., 2009).Note that SMB in RACMO2/GR is only calculated for the area within the ice mask.
For each grid point, pairs of H s and SMB are selected within a search radius of at least 150 km.A distinction is made between accumulation and ablation area, and for each regime the search radius is extended in steps of 5 km until a minimum number of H s -SMB pairs (n = 100) is found.An example is shown in Fig. 2a, for a location in the western ablation zone.In this case, a search radius of 150 km is sufficient for the accumulation area, but for the ablation area the search radius had to be extended to 225 km.
A simple linear regression through all (ablation and accumulation) points does not lead to a useful relation of SMB as a function of H s (Fig. 2b).SMB is overestimated at high values of H s , especially for H s > 2000 m.To solve this, we split the accumulation area and the ablation area, and calculate separate linear relations for both regimes (Fig. 2c).Furthermore, we impose the constraint that the original SMB value at the specific location is reproduced at the correct elevation, so: To this end, the linear regression is forced through the reference value (green dot in Fig. 2) by adjusting the intercept of the line after the regression, without changing its slope.This method ensures a better representation of the SMB gradients compared to a regression that forces the line through the ref-  accumulation area, nor does it lead to regression lines that are a good physical representation of the actual H s -SMB pattern.This is illustrated in Fig. 3, where the dashed lines are the linear regression lines that follow from minimizing the vertical offsets between the points and the regression line.In this case, the reconstructed equilibrium line altitude (ELA) is ∼100 m too high, just as the height at which two lines cross.We improved this by minimizing the perpendicular distance between points and the regression line rather then minimizing the vertical distance.This method can only be used in scatter plots of variables with equal units.We first normalize the model data on both axes, and then perform a linear regression by minimizing the perpendicular offsets (Fig. 3).Hereafter, regression equations are rewritten to non-normalized form, such that SMB is parameterized by: where H c is the elevation of intersection between the two lines (Fig. 3).
The resulting expressions work well for small SMB adjustments following small ice geometry changes.However, when temperature perturbations are imposed (see below), positive values of b acc enforce a positive feedback between H s and SMB.This becomes unrealistic with increasing H s (Figs.2c and 3a).Negative values of b acc can also lead to problems, since these eventually lead to negative SMB values with increasing H s , which is not likely to occur in Greenland (Fig. 3c).To keep SMB within reasonable limits, tests show that the following conditions are considered most realistic for the accumulation regime: SMB pos refers to the mean positive SMB value within the search area, SMB ref is the SMB value at the particular grid point.The minimum of 25 % of the present-day value is based on reconstructions of accumulation rate from ice cores (e.g.Dahl-Jensen et al., 1993).
Results for each grid point in the domain are shown in Fig. 4. SMB gradients for the accumulation regime (b acc in Eq. ( 2), Fig. 4b) are generally small (<1000 kg m −2 yr −1 ), and mostly negative over the interior part of the ice sheet, apart from an area in the central north.b acc is positive along the western margin, implying decreasing run-off and/or an increase in accumulation with increasing elevation.High accumulation in the southeast is reflected in large positive values of a acc in this area, while b abl has high values in the south(-west) and lower values in the north.SMB gradients in the western ablation zone are in the order of We tested the sensitivity of the SMB gradients for changes in the resolution of the input data, by undersampling the 11 km data from the regional climate model (not shown).These tests revealed that the patterns in Fig. 4 are robust, and that SMB gradients are only weakly dependent on resolution.However, to apply this method, it is critical that the narrow ablation zone is resolved by the regional climate model.This will generally not be the case in coarse-resolution GCMs.
These results allow us to calculate a continuous SMB field, as a function of H s , also for areas outside the present-day  ice mask, required when the ice sheet expands.Once a grid point outside the present-day ice mask becomes ice covered, an SMB value is calculated based on H s -SMB pairs from the ice-covered area in the vicinity of this location.However, as long as such a location is not (yet) ice covered, a lower SMB value should be assigned than the value that follows from the parameterization, to account for the influence of e.g. a lower albedo of tundra compared to ice.To that end, we subtract 1000 kg m 2 yr −1 from the calculated SMB when ice thickness is below 1 m.Different values for this modification have been tested.This value best prevented an unrealistic expansion of the ice sheet for present-day SMB forcing.
Figure 5 shows the elevation of intersection (H c ) and the ELA for each individual grid point.As expected, a north-south gradient is present, reflecting ELA in the RACMO2/GR fields.The area with low ELA in the southeast is caused by high accumulation, prohibiting ablation on the 11 km ice sheet mask (in reality the ablation zone there is a few km wide).The east-west gradient over the northern part of the domain is due to the higher accumulation in the northwest.

PDD model
We compare our results with a PDD model, forced by RACMO2/GR (1958-2007) mean T s .Choices of parameters in the PDD model are made such that the ice sheet integrated value of SMB PDD agrees with the present-day SMB field as reported by Ettema et al. (2009) (Appendix A).The grey lines in Figs.2d, 3b and d show SMB -elevation relationships resulting from this PDD method.SMB PDD is calculated using different values of mean annual temperature, but to facilitate comparison with the SMB gradient method, the results are plotted as a function of elevation, using γ atm to translate T s to H s .This affects the SMB gradient: the general picture of steeper SMB gradients from the PDD method in the ablation regimes is less pronounced with a smaller value for γ atm .Figure 2d also shows that the PDD method underestimates SMB values in the higher accumulation zone, due to a too strong decrease of accumulation with increasing elevation.

Temperature adjustment by refreezing
Ice temperature influences ice velocity through the temperature-dependence of ice viscosity.The temperature of the ice is determined by ice advection, diffusion, The Cryosphere, 6, 255-272, 2012 www.the-cryosphere.net/6/255/2012/geothermal heat flux at the bottom, heat production due to ice deformation, friction of the ice at the bottom when sliding occurs, and the mean annual temperature at the surface (T s ).We use the RACMO2/GR 1958-2007 mean T s for the latter, and correct for elevation changes using the atmospheric lapse rate γ atm (Table 1).Another process that affects the ice temperature is refreezing (R) of percolating meltwater in the firn.R is often calculated as a fraction of the annual ablation, by making assumptions about the seasonal cycle of surface temperature and snowpack characteristics.
In our model set up, the effect of R on SMB is already taken into account, but we still need to take into account the thermodynamic effect of R on ice temperatures.For this, we can use the refreezing as a separate forcing field (Fig. 6), by applying a relationship between R and the associated ice temperature warming as suggested by Reeh (1991): with R in kg m −2 yr −1 .However, using a fixed field of R poses a comparable problem as for SMB, since R will likely change with a changing ice geometry.Hence, we treat this problem in a similar way as we do for SMB, by calculating the gradients of R as a function of H s , using the same set of points for each location for consistency.Figure 7 shows an example for a location in the northeast.As an additional constraint for the R H s relation, we impose a positive gradient in the ablation regime, and a negative gradient in the accumulation regime.These are the expected patterns, as R will increase with H s in the ablation zone, since the water retention capacity will increase with the thickness of the snow layer.In the accumulation zone a reverse relation between R and H s will occur, due to decreasing availability of meltwater.Reijmer et al. (2011) compare R from RACMO2/GR with different parameterization schemes.

SMB perturbations in climate change experiments
The SMB gradient method is well suited to be used in asynchronously coupled ice sheet -climate model simulations: it is called in the ice sheet model each time step to account for SMB changes as a result of changes in ice sheet elevation changes and extent.After a certain integration time of the ice sheet model, the ice sheet topography and extent should be updated in the climate model, which can then generate a new SMB field, etc.The method also allows steady-state climate experiments, by modifying the background SMB pattern only for local changes in H s , see below.
However, ice sheet model experiments are often used to simulate the effect of temperature perturbations from proxies.Therefore, we also intend to test whether the SMB gradients as calculated here can be used to translate a near-surface temperature perturbation that is applied uniformly over the ice sheet into a spatially differentiated change in SMB.To this end, we extend our method by introducing an extra term in Eq. ( 2) that accounts for a SMB change as a function of a near-surface air perturbation.Instead of using the actual ice sheet elevation (H s ) in Eq. ( 2), we use a climatic elevation (H T ) that is adjusted as a function of a near-surface temperature perturbation ( T climate ): For example using γ atm = −7.4K km −1 , a climate perturbation of +1 • C will lead to a decrease of H T of 135 m.With a typical SMB gradient of 2000 kg m −2 yr −1 km −1 in the ablation regime this leads to a drop in SMB of 270 kg m −2 yr −1 .Note that since SMB gradients differ spatially, a spatially homogeneous temperature change will lead to regionally variable SMB adjustments.

Ice sheet model
We use the 3-D thermomechanical model ANICE (e.g.Van de Wal, 1999a,b;Bintanja and Van de Wal, 2008;Van den Berg et al., 2008;Graversen et al., 2010) based on the shallow ice approximation (SIA, Hutter, 1983), and including thermodynamics to explicitly account for the temperaturedependent stiffness of the ice.Hence, ice temperature is calculated based on the 3-D advection, diffusion, friction, geothermal heat flux (G) at the bottom and annual surface temperature (T s ) adjusted for the effect of refreezing (Sect.2.3).The vertical dimension is scaled with the local ice thickness, and consists of 15 layers with increasing resolution near the bed, to accurately take into account the large gradient in ice velocity near the bed.For areas where basal temperatures reach the pressure melting point we allow the ice to slide over its bed, by using a Weertman-type sliding law (Weertman, 1964), corrected for the effect of subglacial water pressure (Bindschadler, 1983).Formation of ice shelves is not allowed; as soon as the ice starts to float, it breaks off.As such, calving by means of a floatation criterion is included, but calving physics are not incorporated explicitly, since model resolution and dynamics are not suited for a more realistic treatment of marine terminating outlet glaciers.The response of a changing ice load on bedrock elevation is taken into account using an Elastic Lithosphere-Relaxing Astenosphere (ELRA) model (Le Meur and Huybrechts, 1996).In summary, the ice sheet model is a traditional SIA model including thermodynamics and bedrock adjustment.Table 1 summarizes the values for different parameters used in the ice sheet model.

Model set up
The different model components (ice flow, thermodynamics, SMB, and bedrock response) are calculated on a rectangular domain of 141 × 77 grid points with a grid spacing of 20 km.Bedrock elevation (H b ) and initial ice thickness (H i ) are from Bamber et al. (2001b), and these fields are interpolated to our ice model grid using the mapping package www.the-cryosphere.net/6/255/2012/The Cryosphere, 6, 255-272, 2012 OBLIMAP (Reerink et al., 2010), using an oblique stereographic projection centered at 72 • N, 40 • W, with projection angle α = 7.5 • .The same mapping configuration is used to interpolate fields of SMB, T s and R from the regional climate model RACMO2/GR, and the spatially differentiated functions for SMB H s ,λ,φ and R H s ,λ,φ are interpolated likewise.
A difference exists between the areal extent of the GrIS ice thickness data as presented in Bamber et al. (2001b) and the arial extent of the ice mask in Bamber et al. (2001a), the latter also containing the spatial distribution of numerous small ice caps and glaciers along the periphery of the GrIS.The differences are especially prominent along the rugged topography of the east coast, e.g. in the area south of Scoresby Sund, where numerous glaciers and ice caps exists, without information on ice thickness.Since SMB from Ettema et al. ( 2009) is available for the ice mask from Bamber et al. (2001a), we apply a correction to the initial ice thickness field from Bamber et al. (2001b) by assigning a 10 m thick ice layer to all grid points within the ice mask but where ice thickness data is missing, and let the model freely evolve from there.
Initialization of the 3-D temperature field is done by using the Robin solution based on T s , G and SMB in the accumulation zone.Ice temperatures in the ablation zone are initialized as a linear profile between T s and the pressure melting point.

Reference experiment
We start with a steady-state run of 100 kyr using constant present-day forcing (Fig. 8a).Within 10 kyr, the ice volume initially increases from the present-day observed value of 2.90 × 10 15 m 3 to 3.20 × 10 15 m 3 .After ∼30 kyr the ice volume has reachted its steady-state value of 3.18 × 10 15 m 3 , 10 % above the present-day observed value.Only ice on the Greenland mainland has been included.Neglecting the effect of refreezing (R) on ice temperature results in a 2 % larger ice sheet (black dashed line in Fig. 8a), due to lower ice temperatures, that limit deformation rate and sliding.The dependency of T s on the amount of refreezing is controlled by the constant in Eq. ( 5).In-or decreasing this value by 20 % does not strongly influence the results, as illustrated by the red and orange curves in Fig. 8a.
To illustrate the effect of using a different climatology, we split the 1958-2007 RACMO2/GR SMB into two periods: before and after 1990.Since 1990, warming over Greenland has resulted in a significant negative trend in SMB (Ettema et al., 2009).Separate SMB gradients are calculated for each period and new steady-state runs were performed with the ice sheet model (red and blue lines in Fig. 8a).The pre-1990 SMB field results in only a small increased ice volume.However, the lower SMB of the post-1990 period results, as expected, in a smaller steady-state ice volume.
Figure 9a shows the steady-state ice sheet after 100 kyr simulation, and Fig. 9b the difference in H s compared to the present-day state.The ice sheet has advanced too far, especially in the southwest, east and north.
Figure 8b shows SMB acc and SMB abl , representing the integrated values of SMB over the accumulation area and ablation area, respectively.Note that these terms cannot be compared with the accumulation and run-off.The steady-state SMB equals 362 Gt yr −1 , and is in balance with the calving flux.This is 23 % lower than the integrated SMB value as calculated for the present-day ice sheet by RACMO2/GR (Ettema et al., 2009;Table 2).This difference can be explained by the expansion of the ablation area in the simulated ice sheet, reducing the SMB.
The too large ice sheet in the southeast is caused by high accumulation in combination with the (initial) absence of a significant ablation zone (Fig. 1).In reality, most of the ice in this area is lost by calving of fast-flowing outlet glaciers, flowing through deep, narrow fjords.This process is not well-simulated for two reasons: (1) narrow fjords are not resolved in the 20 km grid, leading to a seaward displacement of the model coastline; and (2) our SIA-type model does neither accurately describe fast flowing glaciers, nor the calving process.This results in an ice margin advancing towards the coast, increasing the calving flux.In addition it allows for the formation of an ablation zone in areas that were previously ice-free (Fig. 10).Increasing the resolution to 10 km does not improve the results; outlet glaciers in these fjords have typical widths of less than 5 km.
These problems are typical for many ice sheet modelling studies (e.g.Greve, 2005;Graversen et al., 2010;Robinson et al., 2010).The thinner ice sheet interior can be explained by the fact that we have performed a steady-state experiment, whereas the present-day GrIS is not in steady-state with the current climate; rather, it consists of colder ice originating from the glacial that deforms at a lower rate.
The resulting mass balance pattern (Fig. 10) is in reasonable agreement with the original fields (Fig. 1).Differences occur mainly along the eastern margin due to higher elevations in combination with negative SMB gradients in the accumulation area.Furthermore, the simulated ice sheet has formed a wider ablation area.Reconstructed SMB is higher than the original values in the western ablation area, where ice sheet elevations are also greater, but here SMB gradients are positive.

Temperature perturbations
Table 2 shows the sensitivity of integrated SMB to temperature perturbations as calculated using the method in Sect.2.4.SMB becomes negative only at climate perturbations of 4 K and higher, but for smaller climate perturbations the ice sheet will also shrink, as described below.
Temperature-perturbation experiments were carried out (Fig. 11).In these experiments, the steady-state ice sheet is perturbed for 100 kyr, to reach a new equilibrium state.The perturbation has a direct SMB effect and an indirect thermodynamic effect on ice volume.The effect on SMB is controlled by Eqs. ( 2) and ( 6), the dominant mechanism being that a cooler climate results in a more extensive accumulation area and smaller ablation area.However, accumulation regions with negative SMB gradients (b acc ) will receive less accumulation, which can locally result in a SMB.Insets show scatter plots of area and volume as function of temperature perturbations using the SMB gradient method (dots) and using the PDD method (crosses).
Both ice sheet extent (Fig. 11a) and ice volume (Fig. 11b) show a clear nonlinear relation with the applied temperature perturbation, with stronger effects for positive values of T climate .The ice sheet extent hardly increases at lower temperatures, since it almost entirely fills the island of Greenland.Increased ice volume is thus mainly due to thickening of the ice sheet.A slight decrease in ice volume can be identified in the experiments with a temperature perturbations in the range of T climate = −5 to 0 K.This is due to decreasing SMB values in the accumulation area, which outweighs the effect of enlargement of the accumulation area.
A large difference in ice sheet volume (89 %) occurs between the +1 and +2 K experiments.These results are highly dependent on the value of γ atm in Eq. ( 6).A sensitivity test using γ atm = −6 K km −1 in the +2 K experiment has the same effect as a +3 K experiment using γ atm = −9 K km −1 (not shown), which points out that care must be taken with the quantitative robustness of these results.However, it is expected that a threshold value exists for the SMB perturbation, above which the GrIS will eventually retreat to only a fraction of its current size.This is in agreement with e.g.Van de Wal (1999a), who performed a set of similar experiments using a different approach to estimate the SMB forcing.
This set of experiments has also been carried out for a model set-up neglecting the effect of refreezing (not shown).Although the temperature adjustment due to refreezing can be substantial (Fig. 6), the impact on the final results in terms of ice volume are mostly minor (see e.g.Fig. 8).
When the steady-state climate perturbation experiments are repeated using the PDD method (Table 2, crosses in Fig. 11), several differences can be identified.Maximum ice sheet volumes are found for smaller negative values of T climate ; a decrease of ice volume due to decreasing precipitation also occurs in the PDD method, but only becomes dominant at large negative temperature perturbations.
Results for the positive temperature perturbations are particularly different for the +2 K scenario, where the PDD forcing results in a steady-state GrIS volume of intermediate size (∼1.7 × 10 15 m 3 ), whereas the SMB gradient method results in a near-total retreat.

Simulating a full glacial cycle
In analogy with e.g.Letréguilly et al. (1991); Van de Wal (1999a); Greve (2005), we performed an experiment that aims to describe the GrIS evolution through a full glacial cycle.The climate record used as a proxy for T climate is based on the GRIP δ 18 O record (Johnsen et al., 2001), and converted into a surface temperature deviation following Johnsen et al. (1995).Prior to 105 kyr, the GRIP record is not a valid climate proxy due to ice-flow irregularities (North Greenland Ice Core Project members, 2004), so for this period we use the Vostok δD record and blend the two records in a similar way as described in Greve (2005).This temperature forcing is applied uniformly over the domain, and additionally a lapse rate correction on T s is applied (Table 1).Sea level is prescribed using the reconstructed sea level from Bintanja and Van de Wal (2008).We start our glacial cycle experiment at 128 kyr BP, i.e. in the maximum of the Eemian climate optimum, using the present-day ice thickness as initial conditions, just as in the reference experiment (Sect.3.1).
Figure 12 shows ice sheet volume through the glacial cycle.Since these results critically depend on the way Eq.( 6) links a climatological temperature record to changes in SMB, we also show results obtained with different values of γ atm , ranging from −6 to −9 K km −1 .Final ice volume does not vary strongly; the largest differences occur in the first part of the simulation (Eemian), because the positive temperature deviation is largest in this epoch.We do not show minimum Eemian ice volume here, given the strong dependence of the results on the value of γ atm , the exact starting time of the simulation, and initial ice temperatures.Eemian ice volume will be estimated in a future study, using a coupled ice sheetregional climate model simulation.Ice sheet growth during glacial conditions is largely controlled by the accumulation rate, which is constrained by the value of SMB min in Eq. ( 4).If we change the value of this parameter to either 10 % or 50 % of present-day SMB ref , ice volume in the glacial is directly affected, as illustrated by the blue lines in Fig. 12.The final present-day ice volume is not much affected, but the difference between the maximum ice volume at LGM and present-day is more than doubled when using a higher value for SMB min .Not surprisingly, applying the PDD forcing results in yet another ice sheet volume history (grey line in Fig. 13).
The simulated increase in ice volume through the glacial period culminates in a peak LGM ice volume of 3.56 × 10 15 m 3 , which is in the lower range of earlier reconstructions (e.g.Van de Wal, 1999a;Huybrechts, 2002;Robinson et al., 2010) and paleoclimatic evidence (Fleming and Lambeck, 2004), but slightly higher than the reconstruction by Greve (2005).This low LGM ice volume is a least partly caused by the lack of ice shelf dynamics in our model, prohibiting merging of the GrIS and the Ellesmere Island section of the Laurentide Ice Sheet during the last glacial, which probably did occur in reality (England, 1999;Alley et al., 2010).
The simulated deglaciation results in a present-day ice sheet volume close to the steady-state volume.Ice sheet elevation and extent are in reasonable agreement with observations (Fig. 13).Comparison with the steady-state experiment (Fig. 9) shows that a realistic climatic forcing results in an improved ice sheet elevation in the interior.Two marginal areas in the southwest and northwest stand out (arrows in Fig. 13) because they have thinner ice than presently observed, and also exhibit large ablation areas with lower SMB values than in the RACMO2/GR reconstruction (Fig. 14).This highlights the sensitivity of these areas to surface melting.

Discussion
The SMB gradient method is designed to improve asynchronous coupling between climate models and ice sheet models.As shown here, it can also be used as a stand-alone SMB forcing module.A rigorous test of the performance of the SMB gradient method would be to compare results against a regional climate model run obtained on a different ice sheet topography.We do this for a simulation of Eemian climate conditions, since the GrIS was by then clearly out of balance.In this experiment, RACMO2/GR is forced at its lateral boundaries by the ECHO-G GCM, with greenhouse gas concentrations and orbital parameters of 125 kyr BP (Van de Berg et al., 2011).The SMB gradients computed from these Eemian climatologies are used to simulate GrIS retreat.A different resolution (18 × 18 km) was used for these RACMO2/GR simulations, which requires a reduction of the minimum number of H s -SMB pairs to n = 37, to ensure an equal area-of-influence in comparison with the present-day fields.Here we show SMB calculations for two different time slices at 129 kyr and 128 kyr BP (Fig. 15a and b), to assess the performance of the SMB gradient method.The ice sheet extent through the Eemian is highly unknown; for 129 kyr BP we chose an ice sheet configuration intermediate between  estimates for glacial conditions and the present-day extent.
The ice sheet for 128 kyr BP is derived from a 1 kyr simulation using the SMB approach described in this manuscript.
Between 129 and 128 kyr BP (Fig. 15a, b and d), the decrease in ice sheet elevation and SMB is most pronounced along the southwestern margin, where the elevation difference reaches ∼1000 m.The effect of this elevation change on SMB is well captured by the SMB gradient method, as calculated for 128 kyr BP (Fig. 15c), based on the SMB racmo fields at 129 kyr BP (Fig. 15a).The residual difference (SMB gradients − SMB racmo , Fig. 15f) is mainly positive, but not systematic.Most deviations can be explained by the changes in precipitation, which is a dynamic response on large-scale topographic changes.For example, the ablation area in south Greenland where SMB is largest, experiences a significant reduction of precipitation because of the concave topography.Still, the correlation between the two SMB fields in Fig. 15b and c is very high (R 2 = 0.988) and the RMSD is only 91 kg m −2 yr −1 , in strong support of the SMB gradient method.When applying the SMB gradient method to asynchronous coupled ice sheet -climate simulations, the time step between the couplings must be chosen.This is not straightforward, as it depends on the rate of the climatic shifts, but also on the magnitude of the change in ice sheet topography.A 1000 yr interval for Eemian conditions is still acceptable, as shown in Fig. 15, but it should be noted that the climate forcing is kept constant in this experiment.
The large differences with the PDD method illustrate the importance of SMB forcing on the outcome of ice sheet model simulations.Advantages of using a climatological SMB from a regional climate model over a PDD based method are that day-to-day variability is implicitly taken into account.Moreover, a regional climate model captures more of the physics of processes important for SMB variability.Van de Berg et al. (2011) showed that a PDD approach can lead to erroneous SMB reconstructions in paleoclimate simulations.As long as direct coupling is practically unfeasible, the best approach to correctly simulate changes in SMB as a result of changes in ice sheet topography is to frequently feed the new topography from the ice sheet model into the regional climate model.At these couplings, SMB from the climate model is accurately reproduced by the SMB gradient method, which is not the case using a PDD approach (Fig. A1).The SMB gradient method is designed to be used in between couplings, and assumes that the regional SMBelevation dependence is persistent and a better predictor for SMB then an indirect correlation between temperature and SMB through elevation, as assumed in the PDD approach.
Although simulated ice sheet volume and extent in this study are within the range of known reconstructions with similar ice sheet models, a note on the excess of ice along the (predominantly eastern) margin is warranted here.It seems a persisting feature in ice sheet model reconstructions (e.g.Greve, 2005;Graversen et al., 2010;Robinson et al., 2010), which has become more prominent since improved bedrock topography (Bamber et al., 2001b) and improved climate fields (e.g.Ettema et al., 2009) have become available.The east coast of Greenland consists of rugged terrain, and receives relatively large amounts of precipitation.The ablation area is too narrow to be properly resolved and to keep the ice margin in place, inducing glacial advance in the model in the direction to the coast.Once glaciated, this area remains covered with ice.Improvements can be expected from regional model grids that have <1 km resolution, resolving the narrow fjords, in combination with higher order ice sheet models, with a better description of fast outlet glacier dynamics.

Conclusions
We presented a novel approach to prescribe SMB fields from regional climate models as a function of ice sheet elevation change, to account for the height-mass balance feedback in ice sheet model experiments.Using the spatial relation between elevation and SMB, a distributed field of SMB gradients is calculated, separately for the accumulation and the ablation regimes, such that SMB values can be retrieved as a function of elevation for each regime, and over the entire domain.It enables a dynamic SMB forcing of ice sheet models, also for initially non-glaciated areas in the peripheral areas of an ice sheet.We applied the method to the GrIS, in two different experiments: (1) using steady-state forcing and (2) using more realistic glacial-interglacial forcing.These experiments result in ice sheet reconstructions and behavior that compare favorably with present-day observations.An evaluation experiment using two different GrIS topographies results in close agreement between the SMB gradient method and the regional climate model, which supports the SMB gradient method.

Fig. 2 .
Fig. 2. Example of the SMB gradient method for a location in the ablation area.Blue (red) dots indicate positive (negative) SMB values and locations, green dot indicates reference H s -SMB value for this specific grid point and black lines denote relations between SMB and H s using different methods: (a) scatter plot of SMB as a function of H s ; (b) simple linear regression; (c) separate regressions for ablation and accumulation regimes; (d) final SMB gradient result maximized to a value (SMB max ) in the accumulation regime.H c is the elevation of intersection between the accumulation and ablation regime.Grey line represents SMB calculated by a PDD model.
Fig. 3. (a) The effect of the choice between least square linear regression by minimizing vertical offsets (dashed lines) and minimizing perpendicular offsets (solid lines); (b) SMB gradients after maximizing the relation for the accumulation regime and forcing the relation for the ablation regime through the reference H s -SMB values; (c) a negative SMB gradient in the accumulation regime illustrates the necessity of introducing a minimum and maximum SMB value (d) to avoid ablation at high altitudes.Grey line represents SMB calculated by a PDD model.

Fig. 4 .
Fig. 4. Coefficients a and b in the spatially varying equation SMB = a + b H s for (a) and (b) the accumulation regime, and (c) and (d) for the ablation regime.

Fig. 5 .
Fig. 5. (a) Elevation of intersection of H s -SMB relations for ablation and accumulation regimes (H c ) and (b) resulting ELA using the H s -SMB relations.

Fig. 8 .
Fig.8.Time series of (a) ice volume for experiments including (black solid line) and excluding (dashed line) refreezing, and using the PDD method (gray line) resulting from a constant climate forcing; using two different climate periods as a forcing results in the blue (pre-1990) and green (post-1990) curve; in-or decreasing the effect of refreezing on T s by 20 % results in the red and orange curve, respectively; (b) mass balance components in the reference experiment (solid black line in panel a).

Fig. 9 .
Fig. 9. (a) Steady-state ice sheet and (b) difference with present-day observed elevation after a 100 kyr run using present-day climate forcing and the SMB gradient method.

Fig. 10 .
Fig. 10.(a) Steady-state SMB and (b) difference with present-day reconstructed SMB from Ettema et al. (2009) after a 100 kyr experiment using present-day climate forcing and the SMB gradient method.

Fig. 11 .
Fig. 11.Time series of (a) ice sheet area and (b) volume resulting from temperature perturbation experiments in combination with the SMB gradient method.Insets show scatter plots of area and volume as function of temperature perturbations using the SMB gradient method (dots) and using the PDD method (crosses).

Fig. 12 .
Fig. 12.Ice volume over the last glacial-interglacial cycle resulting from experiments using a T climate forcing from ice core records, using the SMB gradient method (black line) and the PDD method (gray line).Colored lines show results or sensitivity experiments (see text).
Fig. 13.(a) Ice sheet elevation and (b) difference with present-day observed elevation after the glacial cycle experiment.

Fig. 15 .
Fig. 15.SMB gradients evaluation: (a, b) two different ice sheet geometries (dashed lines show 500 m contourlines) and SMB patterns (colors) as calculated by RACMO2/GR using Eemian maximum insolation as boundary conditions for 129 and 128 kyr BP; (c) predicted SMB pattern for the ice sheet elevation at 128 kyr using SMB gradients based on the RACMO2/GR SMB pattern at 129 kyr BP from (a); (d) differences in H s ; (e) differences in SMB RACMO ; (f) differences between SMB gradients and SMB RACMO at 128 kyr BP.

Table 1 .
Ice sheet model parameter values.

Table 2 .
Ice sheet integrated SMB as a function of T climate using the SMB gradient method and a PDD approach.