Greenland ice sheet surface mass balance: evaluating simulations and making projections with regional climate models

. Four high-resolution regional climate models (RCMs) have been set up for the area of Greenland, with the aim of providing future projections of Greenland ice sheet surface mass balance (SMB), and its contribution to sea level rise, with greater accuracy than is possible from coarser-resolution general circulation models (GCMs). This is the ﬁrst time an intercomparison has been carried out of RCM results for Greenland climate and SMB. Output from RCM simulations for the recent past with the four RCMs is evaluated against available observations. The evaluation highlights the importance of using a detailed snow physics scheme, especially regarding the representations of albedo and meltwater refreezing. Simulations with three of the RCMs for the 21st century using SRES scenario A1B from two GCMs produce trends of between − 5.5 and − 1.1 Gt 2 in SMB (equivalent to +0.015 and +0.003 mm sea level equivalent yr − ), with trends of smaller for scenario E1, in which are mitigated. Results from one of the RCMs whose present-day simulation is most realistic indicate that thereby been


Introduction
During the 21st century, loss of mass from the Greenland ice sheet in response to anthropogenic climate change is expected to make a substantial addition to global mean sea level (Meehl et al., 2007). The ice sheet contributes to sea-level rise through dynamical processes (ice flows from the interior to the coast, followed by iceberg calving) and surface mass balance (SMB; the net balance between accumulation via snowfall and ablation via melt and subsequent runoff). The accurate calculation of SMB requires a good representation of snowfall and melt. Ice sheets are steep at the margins and flat in the high-elevation interior. Most precipitation is orographically forced and falls at the ice sheet margins; most of the ice sheet ablation also occurs there. General circulation models (GCMs), which are the class of model used to make predictions of climate change, generally have insufficient resolution to represent the orography accurately at the margins of the ice sheets (see Gregory and Huybrechts, 2006, and references cited therein), unless they use subgrid-scale tiling to represent the range of elevations present in each grid box.
To reach the high resolution necessary to resolve the steep coastal topography of the Greenland ice sheet (GrIS), two methods have previously been used: statistical downscaling techniques that produce higher-resolution output from the low-resolution GCM fields (e.g. Huybrechts et al., 2004;Hanna et al., 2005;Gregory and Huybrechts, 2006;Vizcaíno et al., 2008Vizcaíno et al., , 2010; and dynamical downscaling with regional climate models (RCMs) at high spatial resolution, forced at the boundaries by GCMs or reanalysis products (e.g. Box et al., 2004Box et al., , 2006Fettweis et al., 2005;Lefebre et al., 2005;Ettema et al., 2009Ettema et al., , 2010Mernild et al., 2010). Hanna et al. (2011) used reanalysis data to drive a degree-day model and obtain estimates of SMB. Robinson et al. (2012) ran a high-resolution energy-moisture balance model, coupled to a high-resolution ice sheet model; they used a range of values for model parameters, and obtained a distribution of possible future trajectories of SMB.
A number of studies have highlighted the importance of the resolution in reproducing the effects of the topography on the surface mass balance of the GrIS (e.g. Box et al., 2004Box et al., , 2006Fettweis et al., 2005Fettweis et al., , 2011aLefebre et al., 2005;Ettema et al., 2009;Lucas-Picher et al., 2012;Bengtsson et al., 2011). High-resolution RCMs are therefore an important tool for making reliable projections of sea-level rise.
Another aspect of systematic uncertainty in projections arises from the formulation of the ice sheet SMB model. Most previous work has been done using temperature-index and positive-degree-day schemes (e.g. Ridley et al., 2005;Mikolajweicz et al., 2007;Graversen et al., 2011), but a surface energy balance model, while requiring more input data than a degree-day model, is more physically satisfactory (Bougamont et al., 2007;Vizcaíno et al., 2010).
Regional climate models offer the possibility of incorporating a detailed surface model, including mass and energy balance, coupled to the overlying atmosphere model in a physically consistent way. Finally, projections of ice sheet SMB also depend on the scenario of future emissions (Huybrechts et al., 2004;Meehl et al., 2007).
Here, a number of simulations conducted with RCMs for the GrIS are examined, with a view to understanding the uncertainties in model estimates of SMB. In Sect. 2, a brief description is given of the RCMs, and of the simulations performed. In Sect. 3, results from simulations of the recent past are presented, and evaluated against available observations. In Sect. 4, future projections from the RCMs with two emissions scenarios, A1B and E1, are presented, and the modelled SMB changes are related to changes in climate drivers. The paper ends with conclusions from this study in Sect. 5.

Regional climate models
Four RCMs (HadRM3P, HIRHAM5, MAR and RACMO2) are run for the recent past. Three of these (HadRM3P, HIRHAM5, MAR) are then used for the future scenarios. Identical boundary conditions were applied to all the RCMs, but the extent of the spatial domain used in each RCM was different. The domain used in HadRM3P, MAR and RACMO2 covered approximately the area shown in Fig. 1 of Ettema et al. (2010), while that used in HIRHAM5 was larger. With a larger domain, the RCM physics is likely to have a greater influence over the GrIS.

HadRM3P
HadRM3P (Jones et al., 2004), run at the Met Office Hadley Centre, is a limited-area, atmosphere-only model based on a version of the HadCM3 GCM (Gordon et al., 2000) with improved atmospheric physics and an improved (zero-layer) surface scheme (MOSES 2.2; see Essery et al., 2001) with the snow albedo treatment of Marshall (1989) (a linear function of near-surface air temperature). Meltwater percolation and refreezing is not included in the surface snow scheme, and therefore refreezing is calculated offline with the scheme used in the GLIMMER ice sheet model (Rutt et al., 2009;Hagdorn et al., 2010) as a constant multiplied by the instantaneous daily snow cover. HadRM3P uses a polar rotated grid, at a resolution of 0.22 • (equivalent to ∼ 25 km), with 19 vertical levels, and a timestep of 300 s.

HIRHAM5
HIRHAM5, run at the Danish Meteorological Institute, is a combination of two models. The atmospheric dynamics is from the High Resolution Limited Area Model (HIRLAM) (Eeorla, 2006) and the physics from the ECHAM5 global model (Roeckner et al., 2003). Simulations with HIRHAM5 over Greenland have been well validated with ice core and automatic weather station data (Dethloff et al., 2002;Box and Rinke, 2003;Kiilsholm et al., 2003;Stendel et al., 2008;Lucas-Picher et al., 2012;Mottram et al., 2012). The land surface scheme has been modified to account for melt and meltwater retention processes in snow, but analysis of the model results suggests that only a small amount of the meltwater is refrozen in this scheme (Mottram et al., 2012). The albedo of snow and ice is assumed to be a linear function of surface temperature, ranging between a minimum value (0.6) at the melting point to a maximum value (0.8) for temperatures below −5 • C (Roeckner et al., 2003). HIRHAM5 uses a polar rotated grid at a resolution of 0.25 • (equivalent to ∼ 27 km) with 31 vertical levels and a timestep of 300 s in the dynamical scheme.

MAR
The MAR (Modèle Atmosphérique Régional) RCM, run at the University of Liège, is coupled to the one-dimensional surface vegetation-atmosphere transfer scheme SISVAT (Soil Ice Snow Vegetation Atmosphere Transfer) model (Gallée and Schayes, 1994). The snow-ice component, based on the CEN (Centre d'Etudes de la Neige) snow model, CROCUS (Brun et al., 1992), is a one-dimensional multilayered model that determines the energy fluxes between the sea ice, the ice sheet surface, the snow-covered tundra, and the atmosphere. It includes snow thermodynamics, meltwater refreezing, snow metamorphism, snow/ice discretisation, and an integrated surface albedo (Gallée et al., 2001). The version of MAR used here (Fettweis et al., 2011a) is calibrated to compare best with the satellite derived melt extent over the period 1979-2009. MAR is run with a horizontal resolution of 25 km, with 31 vertical levels and a timestep of 150 s.

RACMO2
The Regional Atmospheric Climate Model version 2.1 (RACMO2), run at the University of Utrecht, is a combination of two numerical weather prediction models. The atmospheric dynamics originates from the HIRLAM model (version 5.0.6;Undén et al., 2002), and the physical processes are adopted from the global model of the European Centre for Medium-Range Weather Forecasts (ECMWF, updated cycle 23r4; White, 2004). The tuning of the model is described by Van Meijgaard et al. (2008). In Greenland simulations, RACMO2/GR is extended with a multi-layer snow model to represent the surface and sub-surface processes over ice sheets (Ettema et al., 2009(Ettema et al., , 2010). This snow model includes snow/ice melt, percolation, refreezing, snow compaction, meltwater runoff and heat diffusion; the surface snow/ice density determines the surface albedo.
Here, RACMO2 was run with a horizontal resolution of 0.1 • (equivalent to ∼ 11 km), on 31 vertical levels, with a timestep of between 240 and 360 s, depending on the wind speed.

Boundary conditions
The RCMs are driven at their domain boundaries by sixhourly winds, temperature, humidity and surface pressure provided from a lower-resolution global model. The RCM then downscales the boundary conditions over a transition zone over which the RCM adapts the boundary conditions to its own interpretation of the physics. The ocean surface is updated daily by fields of sea surface temperature and sea ice cover. Here, boundary conditions for the recent past from ERA-40 and ERA-Interim reanalyses, and from two different GCMs, have been used. The GCM boundary conditions were from HadCM3 (Gordon et al., 2000) at a resolution of 3.75 • × 2.5 • , with 19 vertical levels, and ECHAM5 (Roeckner et al., 2003) at a resolution of ∼ 3.8 • , also with 19 vertical levels. The boundary conditions used in the 21st century projections were from a HadCM3 simulation with the SRES A1B scenario (Nakicenovic et al., 2000), and two ECHAM5 simulations, one with SRES A1B, and one with the E1 mitigation scenario used in the ENSEMBLES project (Lowe et al., 2009).
Over the GrIS, the near-surface air temperatures in both HadCM3 and ECHAM5 are warmer in summer than those in ERA-40 ( Fig. 1). However, while ERA-40 may be considered to be more realistic than the GCMs, it has been found to have a tropospheric cold bias in the Arctic (Bromwich et al., 2002;Uppala et al., 2005). This cold bias was removed from 1997 onwards as an additional effect of improvements made to satellite data processing with the aim of solving a tropical precipitation bias (Bromwich and Wang, 2005;Bromwich et al., 2007). As a result of this change in 1997, an artificial positive trend in Arctic temperature has been found in ERA-40 data (Screen and Simmonds, 2011).
All simulations used present-day ice sheet surface topography at a resolution of ∼ 5 km (Bamber et al., 2001), interpolated to the appropriate RCM grid. The boundary conditions used in each RCM are summarised, with dates, in Table 1.

Comparison and evaluation of simulations for the recent past
In this section, simulations of the recent past from four RCMs are intercompared to identify common characteristics, and to assess the uncertainty in the SMB and the reliability of the simulations. For evaluation of the RCMs, we concentrate on simulations forced by the reanalyses since these are likely in general to provide more realistic boundary conditions, in particular because they used observed sea surface conditions, whereas the GCMs include their own ocean models, which have biases, and even if they were perfect would not reproduce actual history because of unforced variability in the climate system. However, reanalysis models are also imperfect in some respects, and differences observed here between GCM outputs and ERA-40 should not necessarily be interpreted as model biases in the GCMs. Moreover, it is relevant to analyse the GCM-driven simulations of the past, comparing them with the reanalysis-driven simulations, because the former provide the baseline for the projections.

Near-surface air temperature (T as )
The near-surface air temperatures (T as , i.e. air temperature at a height of 1.5 m above the ground for HadRM3P, 3 m for MAR and 2 m for HIRHAM5 and RACMO2) obtained from the RCM simulations for the recent past are assessed against observations from the Danish Meteorological Institute's (DMI) synoptic weather stations situated along the Greenland coast (mostly from Cappelen, 2011, with additional data from J. Cappelen, 2012, personal communication), and the Automatic Weather Stations (AWS) of the GC-Net network (Box and Steffen, 2000), shown on the maps in Fig. 3. Precise locations of the stations can be found in the table on page 39 of Cappelen et al. (2000), and in Table 1 of Box and Steffen (2000). For DMI coastal stations, longterm observations are available, and twenty-year means of modelled T as are evaluated against the means of the corresponding years in the observations (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) for the ERA-Interim-forced simulations; 1980-1999 for the others). For GC-Net stations, observations are only available from the mid-1990s onwards, and all simulations are evaluated against the mean observations for the period for which they are avail-able; this period depends on the period for which each station was operational (see Box and Steffen, 2000).
To perform the evaluation, the model T as field and orography were each interpolated horizontally by 2-D bilinear interpolation to the location of the observing site. The surface lapse rate was then calculated by first determining the nearest model grid box to the site, then using the model T as and surface elevations for that grid box and the eight surrounding grid boxes (neglecting ocean grid boxes). This lapse rate was used, along with the horizontally interpolated orography and the published elevation of the observing site, to apply a vertical correction (as much as 3 • C) to the horizontally interpolated T as .
HadRM3P has a cold bias at the coastal stations for all forcings for the recent past (Fig. 2a). When forced by ECHAM5 boundary conditions, HadRM3P is considerably colder on the ice-free areas around the coast (by ∼ 8 • C) than MAR ( Fig. 3; also seen in Fig. 2). For stations in the percolation zone, HadRM3P underestimates T as when the observed T as is above the freezing point ( Fig. 2b; the three GC-Net stations concerned can be seen, marked by triangles, close to the west coast in Fig. 3). This is caused by a lack of snow melt in this region, meaning that the snow cover persists, and the temperature is held at the freezing point, even in summer. It will be seen later that the snow melt in HadRM3P is underestimated because the surface albedo is too high near the coast. HIRHAM5 appears to give T as which are slightly too high compared to observations ( Fig. 2c-d); RACMO2, on the other hand, is slightly too cold around the coast (Fig. 2g), but performs better in the interior (Fig. 2h). Both HadRM3P and MAR give a lower GrIS mean T as when forced by GCM boundary conditions than when forced by ERA-40 (not shown). This is the opposite of what was seen in the GCMs and ERA-40 themselves, where T as was greater in the GCMs than in ERA-40 ( Fig. 1), although the 500 hPa temperature was greater in ERA-40 than in the GCMs. This illustrates the importance of the RCM physics and dynamics, in addition to the boundary conditions.
For both HadRM3P and MAR driven by ERA-Interim boundary conditions, the correlation between 20-yr mean modelled and observed T as is better in the interior than on the coast (correlations for ERA-Interim-forced simulations are given on the plots in Fig. 2). Both HadRM3P and HIRHAM5 are warmer than MAR in the interior, in some regions by 3-4 • C (Fig. 3b, c), but the root mean square errors (RMSEs) are similar (Fig. 2).

Melt area extent and meltwater production
Ablation is a key factor in the GrIS SMB, and a good representation of melt and meltwater production is therefore essential. Unfortunately, while some short time series of observations have recently become available, e.g. from the PROMICE network (van As et  production. Although not a measurement of melt water volume, a related quantity is the area over which surface snow melting occurs. This quantity can be approximately estimated from spaceborne passive microwave brightness temperature in the 19 GHz horizontally-polarised band, T19H. By comparison with in situ measurements from GC-Net (Box and Steffen, 2000), it is found that daily mean T as above 0 • C occurs when T19H exceeds 227.5 K (Fettweis et al., 2011a). Thus, the satellite observations may be used to map the progression of melting across the ice sheet. In reality, partial melt in snow can occur at midday even when the daily mean T as is below 0 • C, so this method detects only free water, which may percolate through the snow or run off. For evaluation of the RCMs, we use simulated melt rate, in mm w.e. day −1 ("w.e." denotes liquid water equivalent), exceeding some threshold as an indicator of the melt extent in daily mean fields. We chose the threshold that gives the best comparison (in Table 2) of the time series of simulated daily melt area with that from T19H. By this method, we cannot evaluate the absolute value of the simulated melt extent, because that has been used for calibration, but we can evaluate the daily and interannual variability of model melt extent. The relationship of RCM meltwater production to T as , and therefore the melt threshold, is sensitive to the surface albedo and eddy heat fluxes. Consequently, the melt threshold is model dependent ( Table 2).
The method has previously been applied to MAR and RACMO2 (Fettweis et al., 2011a); here, it is applied additionally to HadRM3P and HIRHAM5. For given atmospheric and snowpack conditions, a higher albedo leads to a lower meltwater threshold, and vice versa. The albedo in HIRHAM5 is lower than that in MAR and RACMO2 over melting snow (see Fig. 4); hence, HIRHAM5 has a higher melt threshold than MAR or RACMO2. HadRM3P has the highest surface albedo and therefore the lowest melt threshold of all the RCMs. Recalling that HadRM3P is generally biased cold and HIRHAM5 warm, consistent with their differences in albedo, we infer that it is likely that HIRHAM5 overestimates the meltwater production, while HadRM3P underestimates.
To facilitate the intercomparison, the satellite data and the output from all four RCMs were interpolated onto the 25 km MAR grid, and the RACMO2 ice sheet mask was used. The RMSE (relative to the equivalent satellite data) in daily melt area time series in the output from ERA-40-driven RCM simulations provides an indication of the daily variability. The variability compares well in all models (correlation > 0.9), except HIRHAM5.
The Cryosphere, 6, 1275-1294, 2012 www.the-cryosphere.net/6/1275/2012/  Both HadRM3P and HIRHAM5 show a shift in phase of the annual melt cycle to earlier in the year (Fig. 5), probably due to the snow albedo in their simple schemes decreasing too quickly at the end of spring, which enhances the melt. A multi-layer snow pack and a more physical albedo parameterisation would delay the onset of melt due to the thermal inertia of the snowpack and slower snow metamorphism. Snow albedo depends, among other factors, on snow grain size (Kuipers Munneke et al., 2011). MAR is the only RCM to use this physical dependency for the albedo parameterisation.
In RACMO2 and MAR, the maximum melt occurs in the low-elevation coastal ablation zone (Fig. 4a). This is because in the multi-layer snow schemes of these RCMs, the upper snow cover melts, exposing low-albedo bare ice in summer ( Fig. 4b-c). With a high sensitivity of albedo to temperature, HIRHAM5 has a large variation in albedo over the GrIS (Fig. 4d); on the other hand, the HadRM3P albedo scheme has a low sensitivity to surface temperature, and consequently shows little spatial variability in albedo (Fig. 4e) and little melt along the GrIS margins (not shown).
In Sect. 3.4 it will be seen that an important component of a good surface snow scheme is a representation of refreezing of some of the surface melt in the snow pack. Such refreezing releases latent heat and warms the snow pack. It is speculated that the loss of embedded snowpack heat in autumn may contribute to extending the surface melt season, while the simplified parameterisations of snowpack and albedo in HadRM3P and HIRHAM5 may explain the early dropoff in melt extent compared with MAR and RACMO2.
Despite the errors in the phase of the annual melt cycle, all the models reproduce well the observed interannual variability and upward trend in the total melt extent in ERA-Interimforced simulations (Fig. 6). The trends are all significant at the 2σ level.

Accumulation
The annual net accumulation, i.e. the total solid precipitation minus evaporation (PE), of the RCMs may be assessed against observations. Because there are very few observations in regions of the GrIS with net ablation, this assessment is limited to regions with a net accumulation. The observations are from shallow ice cores and stakes (Reeh, 1991; www.the-cryosphere.net/6/1275/2012/ The Cryosphere, 6, 1275-1294, 2012  Bales et al., 2009;Cogley, 2004;van de Wal et al., 2005). Because the observations are sparsely distributed (Fig. 7), it is not possible to obtain a purely observational estimate of the accumulation integrated over the whole ice sheet area. Therefore, our intercomparison and verification of the RCM annual mean PE follows the method of van de Berg et al. (2012). The method assumes that the modelled PE (PE M ) has both a systematic bias and a random error relative to observations, and it transforms the model field in a way which optimises the match to observations according to the assumptions. The bias-corrected modelled PE (PE BC ) is defined as a quadratic function of PE M : and the random error (σ BC ) is subsequently defined as a quadratic function of PE BC : The functions are defined for multi-year time series for each model grid box. The six constants in the above functions must be obtained by optimisation, as follows. The observations are weighted (depending on the RCM resolution) to reduce the over-representation of observationdense regions, and PE M is interpolated to the location of the weighted observations (van de Berg et al., 2006). For each PE observation (PE O ) the normalised difference between the observations (PE O ) and PE BC is determined as Assuming a Gaussian distribution of observational and model errors, the combined error is also Gaussian with a mean and standard deviation that give the difference in the bias and squared sum of their random errors, respectively. So, if these measurements are corrected for their biases and the differences are divided by the squared sum of the random errors, the distribution of differences δ n is ideally a Gaussian distribution around zero with unit standard deviation. However, if this procedure is applied to the whole dataset, the minimal solution is underdetermined, with six tuneable parameters and two control parameters (mean and standard deviation). Therefore, the data is subdivided into four groups based on PE M + PE O (low, medium-low, medium-high and high), thereby ensuring a good fit for the whole range of PE. Fit errors are determined for each subgroup and the sum of fit errors for the whole dataset and the four subsets is used to find the best estimate of the bias correction and random error.
Lack of adequate observations prevents a solution for all of Greenland (Fig. 7), so our solution is restricted mainly to regions of modest snowfall (PE≤ 1000 mm w.e. yr −1 ). In south and southeast Greenland, where model estimates of PE reach up to 5000 mm w.e. yr −1 , observations are lacking. The bias correction is prone to unrealistic behaviour outside the PE domain for which it is tuned. Therefore, the high accumulation area located south of the blue line in Fig. 7 is excluded from the analysis. This high-accumulation area covers only a small fraction of the area (see Table 1).
Nevertheless, this area contributes about one third of the total PE in most models. In Table 3, modelled PE estimates are compared with observed values for the different RCMs and boundary conditions. Since the weighting varies for each RCM, the mean PE O is not equal for all RCMs. Before bias correction, RACMO2 agrees well with reanalysis, HadRM3P underestimates, and MAR and HIRHAM5 overestimate. After bias correction, for each combination, the mean PE BC deviates by only a few percent from PE O , indicating that the model biases have been removed. The root mean square error (RMSE) decreases (in most cases) after the bias correction.
The Cryosphere, 6, 1275-1294, 2012 www.the-cryosphere.net/6/1275/2012/  The area-integrated PE M in the area included in the analysis range from 306 to 532 Gt yr −1 (column 6 of Table 4). Following the bias correction this range is reduced to 403-507 Gt yr −1 (column 9 of Table 4), a range of PE BC which is mainly determined by the RCMs. The spatial precipitation gradients are smallest in the MAR simulations and sharpest in the HIRHAM5 simulations (not shown). These gradients determine the PE maxima in the high-accumulation regions, and consequently total PE BC . The GrIS-integrated value of σ BC is determined assuming that the random errors are spatially autocorrelated over a distance of 200 km, and is typically about 29 Gt yr −1 (column 10 of Table 4) and thus, for most estimates, of lesser significance than the model bias. Except in the case of MAR, reanalysis-driven simulations compare better (smallest bias corrections) with the observations than the GCM-driven simulations. Systematic model biases are generally the largest source of error. Simulations of the recent past driven by ECHAM5 and HadCM3 are generally drier than those driven by the reanalyses, generating estimates that compare less well with observations. If accumulation is underestimated in the recent past, it is likely that its increase will also be underestimated, since models generally predict precipitation changes that are proportional to the precipitation in the baseline simulation (Gregory and Huybrechts, 2006). Since a low accumulation enhances melt, this underestimation implies that the RCM projections may overestimate the increase in melt in a warmer future climate. Both of these points mean that the change in SMB may be negatively biased in the GCM-driven projections.

Surface mass balance (SMB)
Surface mass balance is defined as the difference between accumulation and ablation, i.e. solid precipitation minus the www.the-cryosphere.net/6/1275/2012/ The Cryosphere, 6, 1275-1294, 2012 sum of runoff, sublimation and evaporation. Evaluation of modelled SMB is difficult because few long-term in situ observational records exist. Here, the SMB from the RCMs is compared with previous estimates inferred from model simulations forced with reanalysis datasets Fettweis et al., 2008;Wake et al., 2009), and with available in situ observations. The spread in the 20-yr mean SMB in the RCM simulations (Table 5) is due partly to the differences between RCMs and partly to the choice of boundary conditions. The GCMdriven simulations have lower SMB than those driven by reanalysis data, because the former have less precipitation and more runoff than the latter, due to differences in the driving boundary conditions. Similarly, the SMB estimates of Hanna et al. (2008), Fettweis et al. (2008) and Wake et al. (2009), based on reanalysis-forced simulations (Table 5), also give 20-yr mean SMBs greater than those from the GCM-forced RCM simulations presented here.
The SMB difference between reanalysis-forced and GCMforced simulations is especially pronounced at the ice sheet margins (Fig. 8). To examine the transition from ablation zone to accumulation zone, we evaluate the SMB from HadRM3P, HIRHAM5 and MAR against in situ observations from five sites on the K-transect in western Greenland. Only sites located on the ice sheet itself have been used. The pattern of underestimation at some sites and overestimation at others is strongly dependent on which RCM is used, and less dependent on the forcing (Fig. 9). HadRM3P and HIRHAM5 overestimate the SMB close to the ice sheet margin, and underestimate further away, while MAR tends to underestimate at all sites. The correlation between observed and modelled SMB was found to be similar in all models, but normalised root mean square errors of modelled relative to observed SMB indicate that, at these five sites along the Ktransect, MAR reproduces observed SMB more accurately than HadRM3P and HIRHAM5.
The HadCM3-forced HadRM3P simulation gives higher SMB than that forced by ECHAM5, because the former supplies more precipitation while both generate similar runoff ( Table 5). The HadCM3-forced MAR simulation has lower SMB than that driven by ECHAM5, because both have similar precipitation while the former generates more runoff. The different precipitation produced by HadRM3P and MAR in response to the water vapour from the driving GCMs is indicative of the different boundary layer physics in the two RCMs. The refreezing calculated offline for HadRM3P (see Sect. 2.1.1) tends to be lower, and to have less variability, than in MAR.
HIRHAM5 produces more precipitation than HadRM3P and MAR, possibly because its larger domain allows water vapour gained from the local ocean to supplement that provided by the boundary conditions. Despite this, HIRHAM5 gives lower SMB than either HadRM3P or MAR (Table 5), and indeed it was found to be negative in some years when forced by ECHAM5 (not shown here), because of its higher runoff, likely due to its low albedo and because refreezing is essentially omitted. HIRHAM5 also has the highest interannual variability in SMB (given by the standard deviation), which is likely related to its low surface albedo.

Near-surface air temperature (T as )
A warming trend is observed in all three RCMs (HadRM3P, HIRHAM5 and MAR) when forced by ECHAM5-A1B over the period 2000-2099 (Fig. 10a). However, T as , averaged over the ice sheet, is consistently 2 • C lower in MAR than in the other two RCMs, on account of the cold bias noted for MAR in the ice sheet interior (Sect. 3.1).
In all simulations, the temperature trends between 2000 and 2099 were found to be different from zero at the 2σ confidence level. For HadCM3-A1B forcing, the 2080-2099 mean anomaly relative to the 1980-1999 mean is 4.5 • C for HadRM3P, and 4.2 • C for MAR. The equivalent anomalies for the ECHAM5-A1B-forced simulations are 4.3 • C for HadRM3P, 3.4 • C for HIRHAM5 and 4.0 • C for MAR. For ECHAM5-E1 forcing, they are 2.7 • C for HadRM3P, 1.9 • C for HIRHAM5 and 2.2 • C for MAR. The anomaly is smaller in the simulations forced by scenario E1 because of the lower radiative forcing due to emissions mitigation. The ECHAM5-A1B and ECHAM5-E1 scenarios give similar results for Greenland T as up to about 2050, and diverge thereafter (MAR is shown as an example in Fig. 10b).
In all RCMs, T as increases almost everywhere on the GrIS with ECHAM5-A1B forcing (Fig. 11). This was also seen in the simulations forced by other future scenarios; again, the increase is smaller in the ECHAM5-E1-forced simulations than in the A1B-forced simulations. The more-detailed snow scheme in MAR, compared with HadRM3P and HIRHAM5, results in increased interannual variability (seen in the standard deviation; not shown), but not greater sensitivity to climate change (seen in the 2080-2099 anomaly; Fig. 11).

Melt season length and meltwater production
The algorithm used in Sect. 3.2 to detect meltwater production for the recent-past RCM simulations was applied to the output from the future simulations, and the melt season The Cryosphere, 6, 1275-1294, 2012 www.the-cryosphere.net/6/1275/2012/  2080-2099 2000-2099 2080-2099 2000-2099 2080-2099 2000-2099 2080-2099 2000-2099 2080-2099 2000-2099 2000- length was calculated. For all three RCMs, the relationship between length of melt season and total meltwater production is the same for the ECHAM5-A1B-forced simulation as in the recent-past simulations (shown for MAR in Fig. 12). This suggests that the ratio of melt season length to meltwater production is likely to be conserved in warmer climates; the main change is likely to be an increase in the number of grid boxes with a longer melt season, with a corresponding increase in meltwater production.

Components of surface mass balance
The trend in total snowfall over the ice sheet is significant at the 2σ level in all A1B-forced simulations except ECHAM5-A1B-forced HadRM3P, and in none of the E1-forced simulations ( Fig. 13; Table 6). In all cases, it is small compared to the trends in melt, refreezing and runoff. The trends in total precipitation (not shown) are significant at the 2σ level for all A1B-forced simulations, and are larger than those in snowfall.  14-24 %, compared with 7-17 % for snowfall, indicating that liquid precipitation increases by more than solid precipitation. The small trend in total snowfall can be explained by snowfall increasing in some regions and decreasing in others (Fig. 4.4). This was also found for the other scenarios (not shown in Fig. 4.4, but explaining the results in Table 6). At the beginning of the century, melt is similar in all three RCMs. However, HIRHAM5 and MAR have a higher sensitivity to temperature rise than HadRM3P ( Fig. 13b; Table 6), probably because of the latter's higher surface albedo. HadRM3P and MAR give similar refreezing (Fig. 13c), which, combined with the lower melt in HadRM3P, leads to a lower increase in runoff in HadRM3P compared to MAR (Fig. 13d). The similar melt in HIRHAM5 and MAR leads to a lower increase in runoff in MAR compared to HIRHAM5, because refreezing is underestimated in HIRHAM5. A coding bug identified in the analysis of the longer HIRHAM5 simulations also leads to enhanced runoff in some areas, due to the failure of the snowpack to reaccumulate after the surface snow in a grid box has melted away completely. The difference in thermal conductivity of ice compared with snow means that in these areas there is an initial decrease in summer melt when the snow pack has gone. The lower albedo of the ice surface compared with snow subsequently enhances the early-season melt, leading to a steady increase in melt. Because the ice sheet was initialised with a 10 m snow pack, this error only becomes apparent after several decades, and then it mainly affects limited areas in the ablation zone. , 6, 1275-1294, 2012 www.the-cryosphere.net/6/1275/2012/  However, the surface snow pack above the equilibrium line accumulates as expected so that the effect on the total GrIS mass loss is small. The effects of the albedo parameterisation, refreezing and the snow pack are examined in greater detail in Mottram et al. (2012). In all projections, refreezing increases but by less than melting, so that runoff increases. Snowfall in MAR is largely insensitive to the choice of boundary conditions (Fig. 15a), and there is little trend in any of the simulations. The increase in melt during the century is lower in the E1-forced simulation than in the A1Bforced simulations ( Fig. 15b; see also Table 6) because of the corresponding smaller increases in temperature and radiative forcing. Refreezing is similar in the ECHAM5-A1Band ECHAM5-E1-forced simulations (Fig. 15c), so that the former, which has more melt, also has more runoff (Fig. 15d).   Fig. 12. Scatter plot of the mean number of melt days versus the mean annual meltwater production (in mm w.e. yr −1 ), for each grid box and for each year of simulation, for MAR forced by ERA-Interim over 1989, ECHAM5 (recent past) over 1980-1999and ECHAM5-A1B over 2080-2099 Refreezing in the HadCM3-A1B-forced simulation is similar to the other two early in the century, but has a larger trend (see also Table 6); this increase is concurrent with increases in melt (Fig. 15b) and runoff (Fig. 15d), suggesting that the increase in refreezing is driven by an increase in meltwater production.

Surface mass balance
The SMB change in the ablation zone is dominated by changes in melt; it is greatest in MAR because of that model's depiction of the low bare ice albedo in the ablation zone (Fig. 4b). MAR gives a negative trend in SMB for all three forcings ( Fig. 16b; Table 6), but the trend is larger in the two A1B-forced simulations than in E1-forced simulation. Because of the widespread low albedo in HIRHAM5 (Fig. 4d), the SMB in that model decreases almost everywhere, including many areas well inland away from the ablation zone, especially for ECHAM5-A1B. This results in a large negative trend in SMB ( Fig. 16a; Table 6). A similar trend is seen in MAR; in both of these RCMs, the SMB becomes negative around the middle of the century. HadRM3P has neither the large decrease in ablation zone SMB seen in MAR nor the decreases elsewhere seen in HIRHAM5, so that while the negative trend in SMB in HadRM3P is significant (Table 6), it is smaller than in the other RCMs, and the SMB never becomes negative during the simulation (Fig. 16a).

Relation of SMB change to climate change
Finally, we examine the sensitivity of SMB in the RCMs to climate change in order to determine relationships which may then be used to estimate future SMB changes for other climate scenarios, in the way that Gregory and Huybrechts (2006) used pattern scaling to determine the functional dependence of SMB on T as and precipitation. To reduce interannual variability, decadal means are used. It must be borne in mind that the results of this section are only applicable to the GrIS as a whole, and not regionally or locally. For each RCM, for A1B boundary conditions, the GrIS annual total precipitation depends linearly on GrIS annual total precipitation in the driving GCM (Table 7), indicating that GrIS total precipitation in the RCMs is determined by moisture availability. There is a linear dependence of about 5 % • C −1 of GrIS annual total precipitation on GrIS annual mean T as in the RCMs (Table 7), also found by Gregory and Huybrechts (2006). SMB becomes more negative as temperature rises, giving an increasing contribution to sea level, because the increase in runoff outweighs the increase in precipitation. Runoff in turn is dominated by melting (Figs. 13, 16; Table 6), which occurs mainly in summer, and we find that for each RCM, for A1B boundary conditions, the relationship between summer (JJA) T as anomaly and annual SMB anomaly can be approximated by a linear function (Table 7). The trend is less pronounced in HadRM3P, due to the lower runoff anomaly in that model compared with HIRHAM5 and MAR (Table 6).
For all simulations, GrIS mean T as anomaly in the RCM depends strongly and linearly on that in the GCM (Table 7). The slope of this dependence is lower for HIRHAM5 than for HadRM3 or MAR, probably because the relatively low GrIS albedo in HIRHAM5 leads to melt occurring on the whole ice sheet in that model, limiting near-surface summer warming. Consequently, for each RCM, the relationship of RCM SMB to GCM summer (JJA) T as (Table 7) is also linear, being the product of the relationships of RCM summer T as to GCM summer T as (Table 7) and RCM SMB to RCM summer T as (Table 7). The functional dependence of RCM SMB on GCM annual mean T as is similar to that for GCM summer T as (Table 7).
For each RCM, this linear relationship, together with that RCM's estimate of SMB for the recent past, allows www.the-cryosphere.net/6/1275/2012/ The Cryosphere, 6, 1275-1294, 2012 estimation of the GCM GrIS JJA or annual T as change at which SMB is likely to reach zero (last two rows of Table 7). While our results are relative to the recent past, other studies have used the pre-industrial period as a reference point (e.g. Gregory and Huybrechts, 2006;Robinson et al., 2012). We have found that, in HadCM3, annual mean T as over Greenland is ∼ 1 • C warmer in the recent-past (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) simulation than in a 150-yr control simulation (equivalent to pre-industrial). Our results (last row of Table 7) therefore suggest that, relative to the pre-industrial period, the required increases in T as for SMB to reach zero are ∼ 9 • C, ∼ 2 • C, ∼ 3 • C for HadRM3P, HIRHAM5 and MAR, respectively. The values for HadRM3P and HIRHAM5 are outside the range of 3.2-6.5 • C found by Gregory and Huybrechts (2006) with pattern scaling, while the value from MAR is at the extreme lower end of this range. Meanwhile, the value for MAR lies within the range of 2.0-3.5 • C found by Robinson et al. (2012), while the value for HIRHAM5 is at the extreme lower end of this range and that for HadRM3P is beyond the upper limit. MAR, which has the most detailed surface scheme and was the best performing of these three models in the evaluation against observations (Sect. 3), is therefore consistent with these previously published estimates. The required temperature change is high for HadRM3P, probably because the surface scheme in that model responds weakly to temperature changes (Table 6); the offline refreezing calculation may also affect this result, but this effect is likely to be much smaller than that of the surface scheme because the 21st century trend in melt (which is unaffected by the offline refreezing calculation) is also small ( Table 6). The temperature change is low for HIRHAM5 because that model produces a low SMB in simulations of the recent past (Table 5).

Conclusions
Four regional climate models (RCMs) -HadRM3P, HIRHAM5, MAR and RACMO2 -have been run for the area of Greenland to produce results for the recent past, forced by common sets of boundary conditions which were obtained from reanalysis data and output from two general circulation models (GCMs). The RCM simulations of the recent past have been evaluated against available observations of near-surface air temperature, area where melting occurs, accumulation and surface mass balance. Three of the RCMs -HadRM3P, HIRHAM5 and MAR -were used to simulate the 21st century under two emissions scenarios, again with a common set of boundary conditions obtained from the GCMs. This is the first time an intercomparison of RCM results has been done for Greenland climate and surface mass balance (SMB), and it is motivated by the need for reliable simulations of the Greenland ice sheet contribution to future sea-level change. We have shown that a realistic SMB simulation depends critically on the use of satisfactory schemes for ice sheet surface energy and mass balance. In addition, it is important to resolve accurately the steep topography at the ice sheet margins. In the present study, we have achieved this by using RCMs, which have higher resolution than GCMs; an alternative approach would be to run GCMs with subgrid-scale tiling schemes to represent the range of elevations present in each grid box. The model evaluation reveals that HadRM3P consistently simulates low near-surface air temperatures (T as ) at stations near the coast, HIRHAM5 is generally too warm, and MAR is too cold in the interior of the ice sheet. These biases, especially in ablating areas, are related in part to the different treatments of albedo, which is generally too high in HadRM3P and too low in HIRHAM5. Consequently, there is probably too little melting in HadRM3P, leading to insufficient runoff, while in HIRHAM5 there is too much; The Cryosphere, 6, 1275-1294, 2012 www.the-cryosphere.net/6/1275/2012/ Table 7. Linear fit coefficients for relationships of the form Y = a 0 +a 1 X, denoted Y (X) to indicate dependence of Y on X, between decadal mean changes in GrIS summer T as (T JJA as ; • C), GrIS annual mean T as (T ann as ; • C), annual precipitation (%) and annual SMB (mm yr −1 sea level equivalent) calculated from A1B-forced simulations. Changes are with respect to the 1980-1999 mean in the appropriate simulation for the recent past. GCM GrIS mean T JJA as at which SMB = 0 according to the fitted linear function is also shown. Errors on T JJA as are calculated from the RMSEs on fitted SMB.

HadRM3P
HIRHAM5 furthermore, since HIRHAM5 does not represent meltwater refreezing, it simulates much greater runoff than MAR and HadRM3P. These biases are somewhat offset by biases in precipitation, which is generally too low in HadRM3P and too high in HIRHAM5 and MAR, but overall our assessment is that RACMO2 and MAR give the most realistic simulation of SMB, with HadRM3P biased high and HIRHAM5 biased low.
Despite the absolute biases in temperature and melting, the trend of increasing melt area in recent years, inferred from satellite measurement of microwave brightness temperature, is well reproduced by all the models. The form of the seasonal cycle is also similar in the models, but HadRM3P and HIRHAM5 both have melting beginning and ending too early. We attribute this also to the snow albedo representation in these models, which depends only on surface snow temperature without explicit consideration of snow metamorphism.
The RCMs simulate trends in GrIS mean T as during the 21st century of between ∼ 0.04 and ∼ 0.05 • C yr −1 for the SRES A1B scenario, and between ∼ 0.02 and ∼ 0.03 • C yr −1 for the E1 mitigation scenario (all significant at the 2σ level). The trend is smaller in E1 because of the mitigation of greenhouse gas emissions, but the A1B and E1 scenarios do not significantly diverge until about 2050. Trends in SMB were between ∼ −5.5 and ∼ −1.1 Gt yr −2 for the A1B scenario (all significant at the 2σ level), and between ∼ −1.9 and ∼ +0.4 Gt yr −2 for E1 (where the positive trend was not significantly different from zero). The most negative SMB trends come from HIRHAM5 and the least negative from HadRM3P. This is consistent with their respective runoff simulations for the recent past because the SMB trends are dominated by runoff increases, somewhat offset by precipitation increases of about 5 % • C −1 . In all models, there is an approximately linear relationship between summer (JJA) T as change in the driving GCM and SMB change in the RCM, with HadRM3P least sensitive, giving the smallest in-crease in sea-level contribution of 0.08 mm yr −1 • C −1 sealevel equivalent, and HIRHAM5 and MAR with similar sensitivities of 0.31 mm yr −1 • C −1 and 0.35 mm yr −1 • C −1 , respectively. Negative SMB has been suggested as a threshold beyond which the ice sheet would eventually be eliminated (Gregory and Huybrechts, 2006). By extrapolation, we find that the GCM JJA T as change for Greenland total SMB to become negative is ∼ 1 • C according to HIRHAM5, which gives a very low estimate of SMB for the recent past; ∼ 2 • C for MAR, which has a realistic snow scheme and performs well in recent-past simulations; and ∼ 8 • C for HadRM3P, whose surface scheme responds very weakly to temperature changes. The equivalent thresholds for annual mean T as change are similar to those for JJA, and when taken relative to the pre-industrial period, the threshold for MAR is consistent with previous estimates (Gregory and Huybrechts, 2006;Robinson et al., 2012); the threshold for HIRHAM5 is below the lower limit found by Gregory and Huybrechts (2006) and close to that found by Robinson et al. (2012), while the threshold for HadRM3P is beyond the upper limits found in both of these studies. We note that the sensitivities in the present study could all be overestimated because the GCMs generally give higher temperatures and less precipitation than the reanalysis data; if the reanalyses are more realistic, these biases would tend to lead to underestimated projections of accumulation increase and overestimated projections of ablation increase.
Overall, the models with more detailed snow schemes (MAR and RACMO2) give better agreement with observations than the other two models (HadRM3P and HIRHAM5). Some coarse-resolution GCMs apparently give reasonablyaccurate results for SMB (e.g. Ridley et al., 2005); however, this is because the underestimation of melt is compensated for by the lack of refreezing. Our results underline the need to use a multi-layer snow scheme that includes meltwater percolation, retention and refreezing, snow metamorphism, and albedo evolution in order to make reliable projections of the www.the-cryosphere.net/6/1275/2012/ The Cryosphere, 6, 1275-1294, 2012 Greenland SMB contribution to sea-level change. To compute the total mass balance of the ice sheet (including ice discharge), it is necessary to include an ice sheet model that simulates ice dynamics, calving of ice shelves, and fast-flowing ice streams. To make projections into the future with ice sheet models, accurate forcing from the atmosphere and ocean is required, and there may be feedbacks from the changing ice sheet topography, extent, liquid runoff and ice discharge upon the regional climate experienced by the GrIS. We therefore recommend the development of coupled atmosphere-ice sheet-ocean model systems that can simulate the interaction of the Greenland and Antarctic ice sheets with the changing climate as an element of the development of a new generation of Earth system models.