An efficient regional energy-moisture balance model for simulation of the Greenland Ice Sheet response to climate change

In order to explore the response of the Greenland ice sheet (GIS) to climate change on long (centennial to multi-millennial) time scales, a regional energy-moisture balance model has been developed. This model simulates seasonal variations of temperature and precipitation over Greenland and explicitly accounts for elevation and albedo feedbacks. From these fields, the annual mean surface temperature and surface mass balance can be determined and used to force an ice sheet model. The melt component of the surface mass balance is computed here using both a positive degree day approach and a more physically-based alternative that includes insolation and albedo explicitly. As a validation of the climate model, we first simulated temperature and precipitation over Greenland for the prescribed, present-day topography. Our simulated climatology compares well to observations and does not differ significantly from that of a simple parameterization used in many previous simulations. Furthermore, the calculated surface mass balance using both melt schemes falls within the range of recent regional climate model results. For a prescribed, ice-free state, the differences in simulated climatology between the regional energy-moisture balance model and the simple parameterization become significant, with our model showing much stronger summer warming. When coupled to a threedimensional ice sheet model and initialized with present-day conditions, the two melt schemes both allow realistic simulations of the present-day GIS. Correspondence to: A. Robinson (robinson@pik-potsdam.de)


Introduction
Modeling the future evolution of the Greenland Ice Sheet (GIS) has attracted considerable attention in recent years, due to a potentially significant contribution of the GIS to future sea level rise (Lemke et al., 2007).Over recent decades, significant GIS mass losses have been diagnosed by on-site measurements (Abdalati and Steffen, 2001), InSAR velocity measurements (Rignot and Kanagaratnam, 2006;Rignot et al., 2008), GRACE satellite measurements of gravity changes (Velicogna and Wahr, 2005;Ramillien et al., 2006;Velicogna, 2009) and regional modeling (Box et al., 2006).While it is expected that only a rather small portion of the GIS can melt over the 21st century (Lemke et al., 2007), modeling studies (Van de Wal and Oerlemans, 1994;Huybrechts and de Wolde, 1999;Ridley et al., 2005;Charbit et al., 2008) show that on the millennial time scale, the GIS can melt completely if temperatures stay above a certain threshold.
For the 21st century, a number of coupled general circulation model (GCM) runs for several emission scenarios are available and can be used to force ice sheet models.The use of high-resolution, regional models driven by GCMs could additionally improve the representation of climate change over Greenland (Box et al., 2006;Fettweis, 2007;Ettema et al., 2009).However, for longer time scales, coupled GCMs are not only computationally expensive, but gradual changes in the topography and ice sheet extent should also be taken into account.This requires bi-directional coupling between climate and ice sheet models, which makes these models even more computationally expensive.So far, only a few experiments of this sort have been performed, using rather coarse resolution GCMs (e.g., Ridley et al., 2005;Mikolajewicz et al., 2007;Vizcaíno et al., 2008).
Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Robinson et al.: The Greenland Ice Sheet response to climate change Most studies of the short-and long-term response of the GIS to global warming use a rather simple approach, in which a simulated temperature anomaly field or a constant temperature offset is added to the modern climatological temperatures and a simple correction for elevation change is employed (e.g., Greve, 2000;Huybrechts et al., 2004;Parizek and Alley, 2004;Gregory and Huybrechts, 2006).While such an approach is justified for short-term predictions, it becomes less applicable on longer time scales, when the GIS can change dramatically.Changes in ice sheet extent and elevation will lead to pronounced changes in temperature and precipitation patterns.In particular, the reduction in surface albedo due to a retreat of the ice sheet would cause a large temperature change that would not be captured by a simple elevation correction.These changes would also affect the spatial and seasonal distribution of precipitation, as well as the relative amount of precipitation falling as snow.
The representation of accumulation over the ice sheet using the annual observational field suffers further limitations, in that (1) estimates of the present-day annual accumulation rate are derived from a rather sparse observational network (Bales et al., 2009), and (2) real precipitation over Greenland exhibits significant seasonality.Since most precipitation over the GIS currently occurs in the form of snow, the use of an annual accumulation field has been assumed to allow a reasonable approximation for modeling the surface mass balance (SMB).Depending on the time of year, however, new snow can have a strong effect on surface albedo.In addition, precipitation over the GIS is, to a large extent, topographically controlled.Changes in the elevation (let alone a complete disappearance of the GIS) would have a pronounced effect on the distribution of precipitation and snowfall over Greenland.
To overcome some of the limitations of the conventional approach to representing the climate over the GIS, particularly when used for long-term simulations, we developed a new approach based on a regional energy and moisture balance model.Though relatively simple compared to regional climate models (RCMs), our approach accounts for most essential physical processes.It should therefore be considered as a physically-based downscaling technique, rather than a regional climate model on its own.This approach can be used to determine realistic temperature and precipitation fields over Greenland, given topographic and climatic conditions that are dramatically different from today.Importantly, it is also computationally efficient enough to permit long-term simulations of the response of the GIS to climate change.The model is evaluated here for both present-day and ice-free topographic conditions.
Furthermore, in most previous modeling studies, surface ablation has been simulated using the positive degree day (PDD) method.Besides several applications on a smaller scale (e.g., Braithwaite, 1980), the PDD method has also been utilized to calculate surface ablation in large-scale models of the GIS (e.g., Reeh, 1991;Huybrechts et al., 1991;Calov and Hutter, 1996;Ritz et al., 1997;Janssens and Huybrechts, 2000;Huybrechts et al., 2004, Ridley et al., 2005).In this method, surface melt is explicitly determined from surface air temperature alone.The effect of albedo on surface melt is accounted for implicitly, via different empirical coefficients for snow and ice.Although the method has been successfully tested against present-day empirical data, its applicability to future climate change may be compromised, since the relationship between temperature and albedo will be different under global warming induced by greenhouse gases.Van de Wal (1996) performed a comparison between the PDD method and an energy balance model for Greenland, finding that the sensitivity of the two approaches to climate change varied considerably.
In contrast to the PDD method, another simplified method for computing surface ablation explicitly includes the effects of both temperature and insolation.It has recently been employed by van den Berg et al. (2008) to simulate ice sheet changes through glacial cycles.Such a parameterization, introduced early on by Pollard (1980), found its application to the simulation of the evolution of ice sheets during the ice ages (Esch and Herterich, 1990;Deblonde and Peltier, 1992;Peltier and Marshall, 1995) and is nowadays becoming more prevalent in ice sheet modeling (e.g., Hebeler et al., 2008).However, this method is still not as widely used or accepted as the PDD method.For convenience, we will call it the insolation-temperature melt (ITM) method.ITM requires essentially the same input as the PDD method, although an additional parameterization of surface albedo is needed.It is also computationally efficient, allowing its use for longterm simulations.As it is not known a priori which melt calculation method provides more realistic ice sheet forcing, a comparison of both methods to each other and to RCM results could help quantify uncertainties in future predictions related to the choice of the surface mass balance scheme.

Model description
The model used here to compute the surface boundary conditions over the GIS consists of two parts: (1) the regional energy-moisture balance orographic model (REMBO) that computes surface air temperature and precipitation; and (2) the surface interface, which provides surface ice temperature and surface mass balance to the ice sheet model.Both REMBO and the surface interface calculate daily fields, which allow seasonal variations in surface albedo to affect the climate and melt rate.This provides an important positive feedback, since changes in planetary albedo (via changes in surface albedo) affect the computed temperature and surface mass balance.In turn, changes in topography, simulated by the ice sheet model, affect the simulated climatology via elevation and slope effects (see Sect. 2.1).
The REMBO model is coupled via the surface interface to the ice sheet model SICOPOLIS (Version 2.9, Greve, The Cryosphere, 4, 129-144, 2010 www.the-cryosphere.net/4/129/2010/1997a).SICOPOLIS is a three-dimensional polythermal ice sheet model, which is based on the shallow-ice approximation (SIA).This type of model is the standard for modeling large ice sheets, as the SIA method neglects longitudinal stress gradients, providing significant computational advantages.Its main difference to other ice sheet models is the treatment of the temperate basal ice layers, in which the total heat flux and the diffusive water flux are calculated assuming a mixture of water and ice (Greve, 1997b).The cold ice regions are treated in a similar way to other thermomechanical ice sheet models via a temperature/energy balance equation including vertical diffusion, three-dimensional advection and dissipation terms.
While REMBO calculates climate fields on a low resolution grid (100 km), surface boundary conditions are computed via the surface interface on the grid of the ice sheet model, which has a spatial resolution of 20 km.This allows the surface interface to better resolve the rather narrow ablation zone on the margin of the ice sheet.We also performed equilibrium experiments with a spatial resolution of 10 km for the ice sheet model and surface interface and found minimal difference between the results obtained for the two grids.Therefore, all simulations presented hereafter were performed using the 20 km grid.The topography and albedo computed on the 20 km grid are aggregated to the 100 km grid of REMBO.The surface temperature and precipitation are computed at this lower resolution, and are then bilinearly interpolated onto the 20 km grid to provide input for computing the surface boundary conditions for the ice sheet model.To compute the daily surface mass balance, we used a simple snowpack model combined with the one of the melt models mentioned in the introduction (PDD and ITM).In Sect.3, we will compare results obtained using these two approaches.

Regional energy-moisture balance model
The energy balance model follows a familiar form, first employed by Budyko (1969) and Sellers (1969), and reviewed by North (1981), and still found in most simplified climate models.The equation for the atmospheric moisture budget, similar to that for temperature, was later added to energy balance models to simulate precipitation (e.g., Fanning and Weaver, 1996).Unlike most energy and moisture-balance models of the climate, which are global, the model employed in this study is regional and only applied over Greenland.Compared to conventional climate parameterizations used for forcing the long-term simulation of GIS evolution, the REMBO model provides a number of important advantages, because it explicitly accounts for the ice-albedo feedback, the effect of continentality (namely, enhanced seasonal temperature variations over the central part of Greenland as compared to the coastal areas) and the orographic effect on precipitation.
REMBO is based on two-dimensional, vertically integrated equations for energy (temperature) and water content in the atmosphere.The two prognostic variables are sea level temperature and specific humidity.The temperature and moisture balance equations are only solved over Greenland.Over the boundary ocean, surface air temperature and relative humidity are prescribed, either from climatology or GCM results.The governing equations are based on a number of assumptions.First, it is assumed that the lateral exchange of energy and moisture can be described in terms of macroturbulent diffusion, which implies the dominance of synoptic-scale processes over mean horizontal advection.Second, we assume that changes over Greenland do not affect the climate outside it, i.e., we consider only uni-directional interaction.Third, vertical temperature and humidity profiles are assumed to have a universal structure (e.g., Petoukhov et al., 2000).Finally, the heat capacity of the active soil or snow/ice layer is neglected, as well as changes in cloud cover.
The vertically-integrated energy balance for the total atmospheric column is written in terms of the sea-level temperature T SL as where the first term on the right side of the equation represents the horizontal diffusion of the temperature, second -absorbed solar radiation, third -outgoing long-wave radiation, fourth to sixth -latent heat related to condensation of liquid water, snow formation and surface melting of snow/ice, respectively, and the last term -radiative forcing of CO 2 relative to the preindustrial state (set to zero in this study), c p is the air heat capacity, ρ a is the air density, H a is the atmospheric height scale, D T is the coefficient of horizontal energy diffusion, S is the insolation at the top of the atmosphere, α p is the planetary albedo, A and B are empirical coefficients in Budyko's parameterization of outgoing long-wave radiation, P w and P s are precipitation in liquid and solid form, M s,net is the net surface melt rate (including refreezing), and L w and L s are the latent heats of condensation and snow formation, respectively.The surface temperature, T , is then related to the sea-level temperature by surface elevation z s , multiplied by the free atmospheric lapse-rate, γ a , Next, the moisture balance equation is written as where Q is the surface air specific humidity, H e is the water vapor scale height, D Q is the coefficient of horizontal macroturbulent moisture diffusion and P is the total precipitation.To make the model applicable to the simulation of different climate states, outside of the modeling domain we prescribed relative humidity rather than specific humidity, since the former is less sensitive to temperature changes than the latter.Therefore, the boundary value for specific humidity is computed from the formula where r is the relative humidity and Q sat (T ) is the saturation specific humidity, which is described by the Clausius-Clayperon function of air temperature.The total amount of precipitation is computed, following an approach similar to that of Petoukhov et al. (2000) and Calov et al. (2005), as Here |∇z s | is the module of the gradient of surface elevation, τ is the water turnover time in the atmosphere (set to 5 days), and k is an empirical parameter.Note that Eq. ( 5) is also similar to that used by van den Berg et al. (2008), in that precipitation is strongly dependent on the surface gradient.
The amount of snowfall at each point is calculated as a fraction of the total precipitation, where the fraction, f , depends on the surface temperature.Below a minimum temperature, this fraction is 1 (all precipitation falls as snow), and above a maximum temperature, this fraction is 0 (no snow).The fraction follows a sine function from 1 to 0 between the minimum and maximum temperatures, which were set to −7 • C and 7 • C, respectively, as these were found to provide a reasonable ratio between total snowfall and precipitation, and follow estimates from empirical data over Greenland (Bales et al., 2009;Calanca et al., 2000).The diffusion coefficients, D Q and D T , in Eqs. ( 1) and ( 3), both decrease linearly with latitude φ (in degrees), and D T also increases linearly with surface elevation z s (in meters): where κ Q and κ T are the diffusion constants for moisture and temperature, respectively.The decrease of diffusion with latitude accounts for reduced synoptic activity from the middle to high latitudes, while the dependence on elevation assumes that wind increases with elevation.The latter dependence was necessary for the model to produce the seasonal cycle of temperature correctly over the central part of Greenland.
Outgoing long-wave radiation is parameterized as a linear function of surface air temperature.The values of parameters A and B were found using values for upward long-wave radiation and surface temperature over Greenland from the European Center for Medium Range Weather Forecasts (ECMWF) Reanalysis 40 (ERA-40) data set (Uppala et al., 2005), shown in Fig. 1a.Monthly climatological data over the entire year were used, and the best fit to these data gave parameter values close to those used by Budyko (1969).These and other important numerical parameters of the model are summarized in Table 1.
Planetary albedo is parameterized as a linear function of surface albedo (Fig. 1b).The fit was found from values of surface and planetary albedo, derived from monthly International Satellite Cloud Climatology Project (ISCCP) radiation data (Zhang et al., 2004), and is given by Only summer values (April-September) were used to obtain the fit, since at high latitudes, insolation in winter is insignificant, thus the winter albedo is not relevant.Surface albedo, α s , is calculated as a function of ground albedo (ice or bare soil) and the snow thickness, similar to that proposed by Oerlemans (1991), and more recently used by Bintanja (2002) and van den Berg et al. (2008), The maximum snow albedo, α s,max , has a value of either 0.8 or 0.6, representing either a dry-snow or wet-snow covered surface, respectively.The value is chosen based on whether any melting has occurred at that location on that day.The ground albedo, α g , has a value of 0.4 for ice and 0.2 for icefree land.If no snow is present, the surface albedo equals the ground albedo.Comparison with the ISCCP satellite data for radiation at the surface shows that this parameterization provides a quite realistic range of values of surface albedo for Greenland.
To prescribe boundary conditions for temperature and humidity, we used ERA-40 data (Uppala et al., 2005), since the ECMWF reanalysis data sets have been shown to be quite realistic in representing important climate variables for the Greenland region (Hanna and Valdes, 2001;Hanna et al., 2005).Monthly climatological fields (averaged from 1958 to 2001) of temperature and relative humidity from the 2.5 • ERA-40 grid were bi-linearly interpolated to the Cartesian 100 km grid used in REMBO.In addition, temperature fields were corrected for elevation differences (Hanna et al., 2005) between ERA-40 and REMBO using, for simplicity, the free atmospheric lapse rate γ a = 0.0065 K/m.
The equations for temperature and moisture (Eqs. 1 and 3) are solved using an alternating-direction implicit discretization scheme, which allows a larger time step than a standard explicit scheme.Still, for numerical stability reasons, the time step used to solve the energy balance equations is quite small, on the order of 1/10 of a day and, therefore, the REMBO model is more computationally demanding than the ice sheet model.This does not present a problem for shortterm (decadal to centennial time scale) simulations, but for   millennial and longer simulations, asynchronous coupling between REMBO, the surface interface and the ice sheet model was used.In the equilibrium runs described below, the surface mass balance interface was only called for every ten ice sheet model years and REMBO was called for every one hundred ice sheet model years.This calling frequency would affect the transient behavior of the GIS somewhat, but not the simulated equilibrium state reached after several thousand years.

Surface mass balance
The annual surface mass balance is computed using a simple snowpack model through equations of snow (h s ) and ice (h i ) thickness in meters water equivalent (m.w.e.), calculated daily over the year: where M s is the potential surface melt rate and r f is the refreezing fraction.Snow cover thickness is not allowed to exceed an arbitrary maximum height of h s,max = 5 m.At each time step, any excess of snow thickness above this limit is added to the ice thickness computed by Eq. ( 12), and snow thickness is reset to 5 m.The refreezing fraction is equal to zero in the absence of snow, while for 0 < h s < 1 m, it is defined following Janssens and Huybrechts (2000) as where f (T ) is the fraction of snow of the total precipitation, and r max is equivalent to the "PMAX" factor originally described by Reeh (1991), which indicates the maximum fraction of snow that is able to refreeze.For snow thickness between 1 m and 2 m, the refreezing fraction increases linearly, reaching a maximum value of 1 for snow thickness h s > 2 m, i.e., for the thick firn layer, surface melt does not contribute to runoff, but converts into ice at the bottom of the firn layer.This difference in parameterization of the refreezing from the standard PDD model does not affect the surface mass balance of the ice sheet in equilibrium but plays an important role for the transient response of the GIS surface mass balance on decadal to centennial time scales.Finally, the annual mean surface mass balance of the ice sheet is computed at the end of each year as the difference between the final and initial thickness of the ice.
The surface ice temperature used as a boundary condition for the ice sheet model is determined as T i = min(T a ,0) + 29.2h i,sup , where T a is the mean annual air temperature and h i,sup is the amount of superimposed ice (in m.w.e.) resulting from refrozen snow and rain during the year (Reeh, 1991).
To initialize the surface interface, the snow height is set to the maximum everywhere.Then REMBO and the surface interface are run interactively for 200 years until the melt variables and the snow height reach approximate equilibrium values.

Positive degree day (PDD) method
The PDD method is the conventional approach used to determine the melt potential of a given year, using calculated positive degree days from a seasonal cycle of temperature.It was initially introduced for the simulation of local glaciers by Braithwaite (1980) and was further developed by Reeh (1991).It has been described by several others and is consistently used for ice sheet model surface forcing (e.g., Ritz et al., 1997;Cuffey and Marshall, 2000;Huybrechts et al., 2004;Charbit et al., 2007).
To account for inter-and intra-annual variability, an "effective" daily temperature, T eff , is calculated from the daily temperature, T m , as The value of the standard deviation, σ , was set to 5 • C, as in many previous studies, and T eff was numerically calculated according to the method described by Calov and Greve (2005).Usually, the annual positive degree days (PDDs) are computed as the sum of the effective temperature over the year.In our case, since the model resolves the seasonal cycle, the effective temperature is used to compute daily potential melt rate in much the same way, as where the empirical coefficient b s = 0.003 m.w.e./(day K) for snow and b i = 0.008 m.w.e./(day K) for ice.

Insolation-temperature melt (ITM) method
The ITM method is based on the work of Pellicciotti et al. (2005) and van den Berg et al. (2008).In this method, the potential daily surface melt rate is determined from surface air temperature and absorbed insolation: where τ a is the transmissivity of the atmosphere (i.e., the ratio between downward shortwave radiation at the land surface and at the top of the atmosphere), L m is the latent heat of ice melting, α s is the surface albedo, S is the insolation at the top of the atmosphere, t is the day length in seconds and λ and c are empirical parameters.Unlike the PDD method, this method explicitly accounts for shortwave radiation, and the difference between snow and ice is expressed in Eq. ( 16) via different surface albedo values.
Based on the summer (April-September) ISCCP radiation data, transmissivity over Greenland was prescribed as a function of elevation, with values ranging from about 0.4 to 0.7 (Fig. 1c).The linear fit was provided by where z s is the surface elevation in meters.The winter data were again not used for the fit, because such low values of incoming radiation increase the data spread, making any trend indiscernible.While more complex radiation schemes exist (e.g., Konzelmann et al., 1994), this equation provides a reasonable range of transmissivity values without the need for additional inputs.It should be noted that, at lower elevations, where the short-wave radiation term in the melt equation is more significant, the value of transmissivity is similar to the value of 0.5 used by van den Berg et al. (2008).The parameter λ was set to 10 W/(m 2 K), equal to that used by van den Berg et al. (2008), while c was used as a free parameter.The latter can range from −40 W/m 2 to −60 W/m 2 and still produce acceptable melt values for the present-day GIS, indicating large uncertainty in the choice of this value.This is discussed further below.

Modeling results
In order to evaluate the performance of REMBO and the melt models, we first performed diagnostic simulations with the present-day topography of Greenland and modern climatological lateral boundary conditions.These results were compared with observations and a conventional parameterization of climate forcing used in the European Ice Sheet Modeling Initiative intercomparison project (EISMINT, Huybrechts et al., 1997) and many recent publications (e.g., Ritz et al., 1997;Janssens and Huybrechts, 2000;Greve, 2005), in which temperature is parameterized as a function of elevation and latitude.Ablation and the surface mass balance of the GIS were also diagnosed from these experiments and compared with existing empirical and modeling estimates.REMBO was then coupled with the ice sheet model SICOPOLIS to simulate the equilibrium ice sheet under present-day climate conditions.The equilibrium simulations were performed using both the PDD and ITM surface melt approaches.

Simulations of climatology and surface mass balance with fixed topography
For diagnostic simulations of the present-day climatology and surface mass balance of the GIS, we used the 5 km resolution gridded topography from Bamber et al. (2001), aggregated to the resolution of the ice sheet model (20 km).Temperature and accumulation fields obtained from REMBO for the present-day Greenland topography have been compared to best estimates from observational data sets (correcting for elevation differences via the free atmospheric lapse rate).Several coastal observations were obtained from Technical Report 00-18 of the Danish Meteorological Institute (Cappelen, 2001), which provides long-term means  of various climatic variables taken from automatic weather stations.Other observations were obtained from the GC-Net program for locations on the ice sheet itself (Steffen et al., 1996).Although the earliest GC-Net observations only began in 1995, they are the best resource available currently.The combination of these datasets provided mean monthly observations from 52 station locations, although due to their temporal inhomogeneity, we consider agreement with these observations only as a simple validation.
Temperatures from REMBO agree well with the observations, with an annual mean residual of −0.16±2.48• C. Temperatures obtained from the EISMINT parameterization (corrected for elevation differences via the parameterization's lapse rates) also match observational data almost perfectly in the annual mean, with a residual of 0.03±3.62• C.However, this agreement in annual mean masks some systematic seasonal biases, which are not present in REMBO simulations (Fig. 2a).REMBO temperatures around the Greenland coast are determined by the boundary ERA-40 reanalysis temperatures over the ocean, so the consistency of the REMBO temperatures with empirical data around the coast is not surprising.The EISMINT temperatures are directly based on coastal temperature data, but nonetheless show a slight warm bias at low elevations, as exemplified by the DMI station at Daneborg (id 4330), shown in Fig. 2b.
At higher elevations on the ice sheet, where REMBOsimulated temperatures have more freedom to evolve away from the boundary conditions, the observed seasonal temperature variability is reproduced, except for a small cold bias.An example GC-NET high elevation station at Summit (id 06) is given in Fig. 2b.The EISMINT parameterization matches summer temperatures reasonably well, but winter temperatures are usually underestimated (a similar conclusion was reached by van der Veen, 2002).
The actual annual and summer temperatures predicted using REMBO can be seen in Fig. 3a and b, respectively.The differences between these temperatures and those obtained via the EISMINT parameterization are shown in Fig. 3c  and d.Both the annual and summer temperatures are quite comparable, although there are notable differences in the annual mean temperature at high elevation in the South and for high latitudes.The biases in the annual EISMINT temperatures are due to the choice of latitudinal and elevation gradients optimized to improve the fit in warmer months.This results in a worse annual fit in the North and at high elevations in the South.This difference has little practical effect on ice sheet modeling for present day, given that temperatures in these regions remain well below freezing in the winter.Summer temperatures compare especially well around the coast and margin of the ice sheet, where temperature is most important for diagnosing melt.Table 2. REMBO diagnosed mass balance components for the present-day ice sheet compared to the range of RCMs, composed of results from PolarMM5 (Box et al., 2006), MAR (Fettweis, 2007) and RACMO2/GR (Ettema et al., 2009) Figure 4 shows annual accumulation and total precipitation patterns simulated by REMBO, as well as the most recent estimates obtained from station data and several ice core samples, compiled by Bales et al. (2009).The simulated fields agree reasonably well with observations in large-scale patterns, despite local discrepancies.Particularly, low accumulation in the north and on the highly elevated central part of the GIS is reproduced well, and high accumulation values can be found along the Southeast coast.In REMBO, the gradient of elevation mainly determines how much precipitation will occur, implying that it is a result of orographic uplifting.Notably, however, REMBO produces too much precipitation on the southwest coast and not enough on the southern tip of Greenland -an indication that local circulation may play a role that is not accounted for here.The annual accumulation total for the GIS simulated by REMBO is ca.560 Gt/a, which matches the best estimate from Bales et al. (2009) and is within the range of several recent regional climate model studies (see Table 2).
Given the accumulation and temperature fields, it is then possible to diagnose the surface mass balance of the ice sheet (Fig. 5a and b).Using the fixed, present-day topography, REMBO was coupled with each melt model.Using both approaches, the overall predicted area of melt is fairly consistent; however, using the ITM approach, the area of melt is rather sensitive to the choice of the free parameter c.Here we chose c = −55 W/m 2 to simulate surface melt at similar levels to the PDD model.In Fig. 6, melt obtained from the ITM approach is plotted versus melt obtained from PDDs.For low levels of melt (at high elevation), the models tend to agree, with the ITM approach producing somewhat more intense values.At low elevation, the PDD approach produces considerably higher levels of melt compared to the ITM approach, indicating that the former has a stronger dependence on elevation.Both Van de Wal (1996) and Bougamont et al. (2007) showed a similar relationship existed when comparing the PDD approach to an energy balance model.Furthermore, the ITM approach tends to show higher values of melt at high latitudes (difference shown in Fig. 5c), a tendency also shown by the energy balance model used by Van de Wal (1996).This indicates that the ITM approach likely captures the first-order behavior of an energy balance model and produces a more realistic representation of melt.
In cumulative terms, any difference between the models is difficult to discern.A summary of surface mass balance components using REMBO with both melt models is compared to results from RCMs in Table 2.For individual components, our approach is generally able to perform well within the range of RCM results.Figure 7 shows the surface mass balance versus elevation for REMBO using PDD and ITM compared to output from two RCMs: the RACMO2/GR model for 1958-2008(Ettema et al., 2009) ) and the PolarMM5 model for 1988-2004(Box et al., 2006).The four panels show results for the GIS as divided into four quadrants, with the origin near Summit (−39 • E, 72 • N).Results from the RCMs have been binned to reduce the number of points in the plot, with darker boxes indicating a higher density of points (i.e., where darker blue squares overlap with darker red squares, the RCMs agree).The trends produced using the PDD and ITM methods in all regions tend to fall in the range of the RCM results, although there are some differences.Both the PDD and ITM models produce higher maximum melt values in the North, and particularly in the Northeast, REMBO produces more accumulation than either of the RCMs.Nonetheless, these differences are minor, since accumulation values there are small.In the South, all models agree better.RACMO2/GR generally produces a much wider range of accumulation values, due to the high resolution (11 km) of topographic features (Ettema et al., 2009).The wider spread in the RCM surface mass balance also likely results from a more detailed representation of precipitation, which can vary based on regional processes, whereas precipitation in REMBO is inherently smoother.This comparison shows that for present-day conditions, both the PDD and ITM approach produce melt values that fall in the range of RCM results, and that they can also align with each other, depending on parameter choices.
As mentioned before, the surface mass balance simulated by ITM is very sensitive to the parameter c.A higher value of c shifts the snow line to higher elevations.Also, there is a gap in the near-zero negative SMB values for the ITM model.This stems from the discontinuity in surface albedo that occurs when surface melt begins (when switching from the dry snow albedo of 0.8 to the wet snow albedo of 0.6).This gap, however, has little effect on overall SMB, since it occurs only for very low melt values.Figure 6.Annual melt rate calculated using the ITM approach vs. annual melt rate calculated using the PDD approach at each point on the ice sheet.Fig. 6.Annual melt rate calculated using the ITM approach vs. annual melt rate calculated using the PDD approach at each point on the ice sheet.
To evaluate the sensitivity of REMBO to climate change, using different melt parameterizations, we performed experiments under uniform (in space and time) warming of 1 and 3 • C with present-day insolation, and warming of 3 • C with insolation corresponding to Eemian orbital parameters.In these experiments (referred to in Table 3 as "equil."),REMBO coupled to the snowpack model was run at least for 200 years, ensuring full equilibrium of the SMB was reached with the perturbed boundary conditions.When comparing two different melt schemes, it is important that they have similar present-day surface mass balance values (Table 2), because simulated anomalies strongly depend on the tuning of the models.Results presented in Table 3 show that the equilibrium response of SMB to regional tempera-ture change is rather similar for both melt schemes.In units of sea level rise, the changes in SMB for present-day conditions correspond to ca. 0.3 mm/(a • C), which is again in line with the findings of Van de Wal (1996).A similar sensitivity of Greenland SMB to temperature was found by Janssens and Huybrechts (2000) who only used the annual PDD approach.However, this number is considerably higher than in simulations using output of coupled GCMs (e.g., Huybrechts et al., 2004).This makes sense because, in transient GCM experiments, simulated warming over the ablation zone is considerably lower than the average temperature change over Greenland, while the latter is used to calculate the sensitivity of SMB to regional temperature change.
Obtaining rather similar values for SMB sensitivity to warming obtained using both the PDD and ITM approaches is in apparent contradiction to Bougamont et al. (2007), who found that the PDD scheme predicts much larger changes in Greenland SMB as compared to a physically-based energybalance model in a transient warming scenario.At least partly, this discrepancy can be explained by differences in equilibrium and transient SMB sensitivities to temperature change.The standard PDD scheme (calculated as an annual sum, as opposed to our daily scheme) does not include an evolving snowpack, which means it has no memory of previous years and, therefore, always simulates equilibrium SMB.In reality, the gradual rise of the snowline with a temperature increase will lead to slow melting of the thick snowpack at higher elevations, which will mostly refreeze and will not contribute to the mass loss of GIS.Since our simulations include a snowpack model which has memory, we can illustrate this effect via the instantaneous SMB response to the temperature rise, i.e., the change in SMB that occurs after the first year of applying the temperature anomaly (see the lower rows in Table 3, labeled "inst.").For our PDD scheme, the simulated instantaneous response of SMB to an    (Box et al., 2006) were binned to reduce data density (blue and red squares, respectively), with darker squares indicating a higher density of points.Trendlines and the slopes are shown for the negative points of the PDD (green) and ITM (black) approaches.increase in temperature is appreciably smaller than that obtained in equilibrium.For the ITM scheme, the instantaneous response is less than half of the equilibrium response.The ITM scheme reacts more slowly initially, because it is not driven by temperature changes alone, but also by decreasing albedo as the snowpack melts.The difference between the equilibrium and instantaneous response is even more pronounced for 3 • C of warming.For both melt schemes, the sensitivity of SMB to an increase in temperature in transient global warming scenarios can be expected to lie between the instantaneous and equilibrium response, and thereby should be considerably smaller than the equilibrium one.Therefore, our "equilibrium" results support the notion that the standard (annual) PDD approach does tend to overestimate the rate of GIS mass losses.
When considering the SMB response of these melt models to past climate change during the Eemian interglacial (including an increase in insolation), the ITM scheme shows a stronger equilibrium response.Since orbital variations occur on a multi-millennial timescale, one can expect that, in this case, the response of SMB to climate change can be considered to be in equilibrium.To mimic climate conditions during the Eemian, we again apply the uniform temper-Table 3. Diagnosed equilibrium (equil.)and instantaneous (inst.)change in surface mass balance for the GIS from present day, under 1 • C and 3 • C of warming with present-day insolation and 3 • C of warming with Eemian insolation (EE).All values are in Gt/a and are relative to a present day estimate of ca.300 Gt/a (see Table 2).Otto-Bliesner et al., 2006).SMB changes computed with the PDD scheme differ only slightly from the response to 3 • C warming for present day (due to additional warming over the Greenland interior simulated by REMBO).However, the SMB change simulated with the ITM scheme increases dramatically due to the increased insolation.Thus, the "Eemian" change in SMB simulated by the ITM scheme is more than 50% greater than that simulated by the PDD scheme.These results indicate that changes in insolation are of comparable importance to temperature changes on orbital timescales and, therefore, the PDD scheme likely underestimates the past ice sheet response to climate warming via insolation changes considerably.
To assess the ability of REMBO to simulate the climate under boundary conditions radically different from the present, we performed simulations with a fixed topography of ice-free Greenland with corresponding uplifted bedrock (i.e., equilibrium bedrock after isostatic rebound).This test provides insight into the sensitivity of the REMBO climate to the presence of the ice sheet and can be compared with similar GCM simulations.It is also noteworthy to compare REMBO results with the EISMINT parameterization, which clearly demonstrates the advantages of a more physicallybased approach.
Figure 8a and b shows the difference in the mean summer (June-July-August) temperature between the ice-free and ice-covered present-day Greenland simulations, obtained with the EISMINT parameterization and with REMBO, respectively.Both modeling approaches produce qualitatively similar warming patterns associated with the lowering of elevation over currently ice-covered Greenland.However, there are significant quantitative differences.REMBO shows a surface air temperature increase of 18    10 • C in winter over the central part of Greenland.These numbers agree favorably with the results of simulations performed with GCMs for ice-free Greenland (Toniazzo et al., 2004;Lunt et al., 2004).Meanwhile, the EISMINT approach produces less warming in summer, but stronger warming in winter.As a result, REMBO shows that the magnitude of seasonal temperature variation over central Greenland increases by ca. 8 • C for ice free conditions (Fig. 8c), while the EISMINT parameterization shows a decrease of 6 • C (Fig. 8d).The increase of seasonality and stronger summer warming in the REMBO model can be attributed to the additional warming caused by a lower albedo for ice free conditions during summer, while winter albedo remains practically unaffected by the removal of the ice sheet.The result is stronger warming in the summer and, therefore, an increase in seasonality.The opposite, unrealistic effect produced by the EISMINT parameterization is explained by the higher lapse rate used in winter.As a result, the decrease in elevation leads to a reduction of the temperature difference between summer and winter and a decrease in seasonality.
Due to the large increase of summer temperature, positive surface mass balance (not shown) is only diagnosed over highly elevated areas in eastern and southern Greenland.This has implications for the possible existence of two (or more) equilibrium states under current climate conditions, which will be addressed in a separate paper.

Coupled simulations of equilibrium state
For the next step of model validation, we performed equilibrium simulations of the GIS with constant (present-day) climatological conditions at the lateral boundaries of the model domain.The REMBO model and the surface interface were coupled bi-directionally and asynchronously with SICOPO-LIS, which was run for 100 000 years, ensuring all relevant characteristics reached equilibrium state.As an initial condition, we used present-day data for the GIS elevation and bedrock (Bamber el al., 2001), and the ice temperature was set to −10 • C. The geothermal heat flux was constant and set to 60 mW/m 2 .Since the longest time scale of GIS response is comparable with orbital time scales, an assumption about GIS equilibrium forced only by present-day conditions is not very accurate and, instead, a simulation over several glacial cycles would be a more appropriate procedure.However, because here we are primarily interested in understanding the sensitivity of the simulated GIS to the different methods for determining the surface boundary conditions, we prefer the simpler approach of using constant climatological forcing.
Given the good agreement between the surface mass balance partition estimates determined using PDD or ITM, it is not surprising that using either melt scheme produces quite similar results (shown in Fig. 9a and b).In both cases, the simulated equilibrium ice sheet covers almost the entire area of Greenland, which also occurs in other GIS simulations using similar approaches (Letréguilly et al., 1991;Calov and Hutter, 1996;Ritz et al., 1997;Greve, 2005).The largest discrepancies with observations appears in the southwest and the Northeast, where using the REMBO climatology results in more extended ice coverage.This is mainly due to the overestimation of accumulation in those places by REMBO and cooler temperatures compared to the EISMINT parameterization.Using the ITM approach does significantly improve the agreement with observations in the North, which follows from the stronger diagnosed melt in this region.
In our simulations, the volume of the GIS is overestimated by ca.10-15%.In other words, it was not possible to simulate an ice sheet, which has both the right geometry and the right surface mass balance components.This may be related,  not so much to the deficiencies of the model used for simulation of the surface mass balance, but rather to an intrinsic problem of ice sheet models based on the shallow ice approximation.Since such models do not properly incorporate fast ice transport by ice streams, they require too much contact between the ice sheet and the ocean to produce a considerable amount of ice calving.In reality, the areas where the GIS is in direct contact with the ocean are rather small and yet, ice calving constitutes roughly half of the total ice loss from the GIS.Furthermore, simulating the ice sheet through past glacial cycles would likely improve the present-day representation, a topic that will be addressed in future work.

Conclusions
A regional energy-moisture balance model (REMBO) has been developed, which simulates temperature and precipitation over Greenland.The model is simple and very computationally efficient.Furthermore, it is physically-based and includes an explicit representation of seasonal changes in albedo -attributes that are crucial for simulation of climate conditions considerably different from present day.
Simulated temperature fields agree well with observational data and, particularly, improve the representation of the seasonal cycle as compared to the EISMINT temperature parameterization.Moreover, for ice-free conditions, REMBO and the EISMINT parameterization predict rather different changes.REMBO simulates a large summer warming and enhanced magnitude of seasonal temperature changes, while the EISMINT parameterization shows decreased seasonality.In this respect, the results from REMBO are more consistent with GCM experiments.Simulated precipitation matches observations in large-scale patterns, as well as the annual sum.However, regional deficiencies exist that cannot be eliminated unless more processes are included in the model.
Two different melt parameterizations were evaluated: the PDD approach and the ITM approach, with the latter explicitly accounting for the effects of temperature and insolation on snow and ice melt.The melt models were used to force a simple snowpack model with a daily time step.With the appropriate choice of model parameters, both methods produce similar runoff and total ablation rates for present-day conditions, but they differ in regional details.Both schemes also exhibit rather similar equilibrium SMB sensitivity to temperature changes.However the instantaneous SMB sensitivity to temperature change is different for each model and, therefore, each can be expected to produce different SMB changes in transient global warming experiments.For climate conditions that mimic the Eemian interglacial, even the equilibrium response of Greenland SMB differs considerably between the two schemes, with the ITM model simulating a more than 50% greater change in SMB.
Equilibrium simulations of the present-day GIS with REMBO coupled to the three-dimensional, polythermal ice sheet model SICOPOLIS demonstrate that both melt models allow us to simulate a reasonably realistic GIS.However for both methods, the simulated volume and spatial extent of the GIS are overestimated.Therefore, the present-day surface mass balance and GIS extent and volume do not provide sufficient criteria for a choice between PDD and ITM methods.Nonetheless, the ITM method looks preferable for transient and paleo simulations.
Figure 1.(a) Monthly ERA-40 data of outgoing long-wave radiative flux versus temperature, shown with a linear fit to all data.Monthly ISCCP data, Apr-Sep (dark points) and Oct-Mar (light points), of (b) planetary albedo versus surface albedo with a linear fit to the months Apr-Sep, and (c) transmissivity versus elevation with a linear fit to the months Apr-Sep.All points are from data over Greenland.

Fig. 1 .
Fig. 1.(a) Monthly ERA-40 data of outgoing long-wave radiative flux versus temperature, shown with a linear fit to all data.Monthly ISCCP data, April-September (dark points) and October-March (light points), of (b) planetary albedo versus surface albedo with a linear fit to the months April-September, and (c) transmissivity versus elevation with a linear fit to the months April-September.All points are from data over Greenland.

Figure 2 .
Figure 2. (a) Average and standard deviation of monthly residuals of REMBO (blue) and EISMINT (red) temperatures compared to station data at 52 locations.(b) Monthly temperatures from REMBO (solid blue line) and the EISMINT temperature parameterization (dashed red line), orographically-corrected and compared with DMI station data (thick grey line) for one high-and one low-elevation station on Greenland.

Figure 3 .
Figure 3. REMBO model output of (a) mean annual temperature and (b) mean summer (JJA) temperature.Difference between REMBO and EISMINT for (c) mean annual temperature and (d) mean summer (JJA) temperature.Elevation contours are shown at 300 m intervals. 33

Fig. 3 .
Fig. 3. REMBO model output of (a) mean annual temperature and (b) mean summer (JJA) temperature.Difference between REMBO and EISMINT for (c) mean annual temperature and (d) mean summer (JJA) temperature.Elevation contours are shown at 300 m intervals.

Figure 4 .
Figure 4. Greenland accumulation fields from (a) present-day data, (b) REMBO for presentday topography and (c) REMBO for ice-free, uplifted bedrock topography.Total precipitation fields for the same are shown in (d), (e) and (f), respectively.Elevation contours are shown at 300 m intervals.

Fig. 4 .
Fig. 4. Greenland accumulation fields from (a) present-day data, (b) REMBO for present-day topography and (c) REMBO for ice-free, uplifted bedrock topography.Total precipitation fields for the same are shown in (d), (e) and (f), respectively.Elevation contours are shown at 300 m intervals.

Figure 5 .Fig. 5 .
Figure 5. Diagnosed surface mass balance for fixed, present-day topography using REMBO with ablation determined by (a) PDD and (b) ITM.Panel (c) shows the difference between the two.Elevation contours are shown at 300 m intervals.

Figure 7 .
Figure 7. Diagnosed surface mass balance versus elevation by region for two RCMs and REMBO using the PDD and ITM approaches.The regions are defined as quadrants on Greenland with the origin near Summit (-39 °E, 72 °N).RACMO2/GR results (Ettema et al., 2009) and PolarMM5 results(Box et al., 2006) were binned to reduce data density, with darker squares indicating a higher density of points.Trendlines and the slopes are shown for the negative points of the PDD and ITM approaches. 37

Fig. 7 .
Fig. 7. Diagnosed surface mass balance versus elevation by region for two RCMs and REMBO using the PDD and ITM approaches.The regions are defined as quadrants on Greenland with the origin near Summit (−39 • E, 72 • N).RACMO2/GR results (Ettema et al., 2009) and PolarMM5 results(Box et al., 2006) were binned to reduce data density (blue and red squares, respectively), with darker squares indicating a higher density of points.Trendlines and the slopes are shown for the negative points of the PDD (green) and ITM (black) approaches.

+1
3 • C, along with changes in insolation corresponding to Earth's orbital parameters at 126 ka BP (kiloyears before present).While Eemian temperature anomalies likely exhibit strong seasonal variations, only summer temperatures are important for the SMB simulations and the 3 • C warming is consistent with empirical and modeling estimates of summer temperature changes around the GIS (e.g.,

Figure 8 .
Figure 8. Summer (JJA) temperatures for ice-free, uplifted bedrock topography minus those of present-day conditions for (a) the EISMINT temperature parameterization and (b) REMBO.Difference in seasonality (Jun-Jul-Aug temperature minus Dec-Jan-Feb temperature) for the uplifted, ice-free topography compared to present-day conditions using (a) the EISMINT temperature parameterization and (b) REMBO. 38

Fig. 8 .
Fig. 8. Summer (JJA) temperatures for ice-free, uplifted bedrock topography minus those of present-day, ice-covered conditions for (a) the EISMINT temperature parameterization and (b) REMBO.Difference in seasonality (June-July-August temperature minus December-January-February temperature) for the uplifted, ice-free topography compared to present-day, ice-covered conditions using (c) the EISMINT temperature parameterization and (d) REMBO.

Figure 9 .
Figure 9. Equilibrium simulated GIS elevations starting using REMBO with (a) PDD ablation and (b) ITM ablation, compared with (c) the actual GIS elevation from observations.

Fig. 9 .
Fig. 9. Equilibrium simulated GIS elevations using REMBO with (a) PDD ablation and (b) ITM ablation, compared with (c) the actual GIS elevation from observations.

Table 1 .
Selected parameters used in REMBO and the two melt models, PDD and ITM.
. All values are in Gt/a.