The Cryosphere Climate of the Greenland ice sheet using a high-resolution climate model – Part 1 : Evaluation

A simulation of 51 years (1957–2008) has been performed over Greenland using the regional atmospheric climate model (RACMO2/GR) at a horizontal grid spacing of 11 km and forced by ECMWF re-analysis products. To better represent processes affecting ice sheet surface mass balance, such as meltwater refreezing and penetration, an additional snow/ice surface module has been developed and implemented into the surface part of the climate model. The temporal evolution and climatology of the model is evaluated with in situ coastal and ice sheet atmospheric measurements of near-surface variables and surface energy balance components. The bias for the near-surface air temperature (−0.8C), specific humidity (0.1 g kg −1), wind speed (0.3 m s−1) as well as for radiative (2.5 W m −2 for net radiation) and turbulent heat fluxes shows that the model is in good accordance with available observations on and around the ice sheet. The modelled surface energy budget underestimates the downward longwave radiation and overestimates the sensible heat flux. Due to their compensating effect, the averaged 2 m temperature bias is small and the katabatic wind circulation well captured by the model.


Introduction
The Greenland ice sheet (GrIS) plays a pivotal role in global climate, not only because of its high reflectivity, high elevation and large area but also because of the volume of fresh Correspondence to: M. R. van den Broeke (m.r.vandenbroeke@uu.nl)water stored in the ice mass, which is equivalent with 7 m global sea level rise.Variations in the surface mass balance (SMB) of the GrIS are determined by the balance between incoming (mass gain) and outgoing (mass loss) terms at the surface.The underlying processes are strongly controlled by atmospheric factors.Therefore, understanding the presentday climate of Greenland is important for the interpretation of the current state and prediction of the future state of the ice sheet.
Via multiple feedback mechanisms, changes in ice/snow cover can potentially influence the overlying atmosphere and, therefore, modify the local climate on the ice sheet.To quantify these strong nonlinear interactions, extensive observation campaigns were carried out on and around the GrIS (Heinemann, 1999;Oerlemans and Vugts, 1993).In 1996, the climate network GC-net was established with automatic weather stations (AWSs) to measure the near-surface atmospheric and surface conditions continuously at locations across the ice sheet (Steffen and Box, 2001).
Whereas these meteorological measurements are limited in space and time, regional climate models have the potential to be used as smart interpolators, yielding useful data for a wide range of times and locations not covered by in situ observations.Further, numerical models provide an ideal environment for testing the importance of critical processes in a controlled fashion.
In this study we used the regional atmospheric climate model (RACMO2, Van Meijgaard et al., 2008) adapted specially for the Greenland ice sheet (RACMO2/GR).RACMO2 has been successful in simulating surface heat exchange processes and accumulation in Antarctica (Van Lipzig et al., 1999;Van de Berg et al., 2006).For Greenland, RACMO2/GR showed that considerably more mass accumulates (up to 63% for the period 1958-2007) than previously thought, due to the higher horizontal resolution (11 km) and the ice sheet mask that was used (Ettema et al., 2009).The modelled SMB agrees very well with the 265 in situ observations that match the modelled period (R = 0.95).Neither the SMB nor the annual precipitation bias show a spatially coherent pattern, making post-calibration unnecessary (Ettema et al., 2009).
Here, we present a detailed description of the performance of RACMO2/GR in the lower atmosphere and at the surface.As we want to assess the quality of our model, a comparison with in situ observations is made rather than a comparison with other models, coarser re-analysis datasets or existing parameterizations.The modelled 51-year climatology of the surface and near-surface parameters is presented in Part 2 Ettema et al. (2010).First we describe the model modifications, followed by a description of the model setup and initialization.In Sect.3, we present the in situ observations used for model evaluation.In Sect.4, we assess and discuss the performance of the model, primarily in relation to near-surface and surface conditions using available in situ observations.Concluding remarks are made in Sect. 5.

Model description
In this study, the Regional Atmospheric Climate Model version 2.1 (RACMO2) of the Royal Netherlands Meteorological Institute (KNMI) is used to simulate the present-day climate of Greenland.RACMO2 is a combination of two numerical weather prediction models: the atmospheric dynamics originate from the High Resolution Limited Area Model (HIRLAM, version 5.0.6,Undén et al., 2002), while the description of the physical processes is adopted from the global model of the European Centre for Medium-Range Weather Forecasts (ECMWF, updated cycle 23r4, White, 2004).
At the lateral boundaries, ECMWF Re-Analysis (ERA-40) prognostic atmospheric fields force the model every 6 h.The underlying ECMWF model for ERA-40 has the same physical parameterizations as RACMO2/GR, except for the adjustments described below.The interior of the domain is allowed to evolve freely.In the pre-satellite era, the analyses for the Northern Hemisphere benefit from the wide extent of data available from land-based meteorological stations and ocean weather ships.Therefore, the atmospheric forcing for the Arctic area should be sufficiently well-constrained to start the model simulation in September 1957 (Sterl, 2004;Uppala et al., 2005).After August 2002, operational analyses of the ECMWF have been used to complete the model simulation up to January 2009.In the absence of an integrated ocean or sea ice model, the open sea surface temperature and sea ice fraction are prescribed from ERA-40.In the sea ice data field no distinction is made between one-year sea ice or multi-year sea ice.The minimum/maximum model time step is 240/360 s depending on the maximum wind speed in the domain, to ensure numerical stability.The 51-year simulation took approximately 100 days to run on 60 processors of the ECMWF supercomputer.
RACMO2 has 40 atmospheric hybrid-levels in the vertical, of which the lowest is about 10 m above the surface.Hybrid levels follow the topography close to the surface and pressure levels at higher altitudes.The air temperature and humidity at a standard observational height (2 m above the surface) are computed using an interpolation technique based on the similarity theory applied to the lowest atmospheric model layers (e.g.Dyer, 1974).
The model domain encompasses the Greenland ice sheet, Iceland, Svalbard and their neighbouring seas (Fig. 1).The domain includes 312 × 256 model grid points at a horizontal resolution of about 11 km (0.10 latitudinal degree).This high spatial resolution allows us to resolve much of the narrow ice sheet ablation and percolation zones, as well as the steep climate gradients in the coastal zones.For accurate topographic representation of the GrIS, elevation data and ice mask from the digital elevation model of Bamber et al. (2001) are used, which are kept constant during the model simulation.The model surface area of the ice sheet is 1.711 × 10 6 km 2 , excluding peripheral ice caps (Fig. 1).This is 1% more than previous studies (Box et al., 2006;Fettweis, 2007;Hanna et al., 2008).Sources of uncertainty include the treatment of changing shelf ice and compacted multiyear sea ice area.The underlying vegetation map is based on the ECOCLIMAP dataset (Masson et al., 2003) and has been manually corrected; the original dataset showed too little tundra and too much bare soil along the east coast of Greenland.

Atmospheric model adjustments
General adjustments to the original dynamical and physical schemes in RACMO2 are described in detail by Van Meijgaard et al. (2008).Here we only describe the adjustments to the original model formulation that have been made to better represent the melting snow conditions in the Arctic region (RACMO2/GR).
RACMO2/GR calculates the surface turbulent heat fluxes from Monin-Obukhov similarity theory using transfer coefficients based on the Louis (1979) expressions.An effective surface roughness length is used to account for the effect of small scale surface elements on turbulent transport.Originally, the roughness lengths for momentum, heat and humidity (z 0m , z 0h , z 0q ) included the effect of enclosing vegetation, urbanization and orography.This approach gave too large values over the Antarctic ice sheet (Reijmer et al., 2004).Therefore, we limited z 0m to 100 mm for tundra without snow and to 1 mm for snow-covered tundra.The value for z 0m at the snow covered ice sheet is set to 1 mm, while z 0m is set to 5 mm if bare glacier ice is at the surface.The roughness lengths for heat and humidity over snow surfaces are computed according to Andreas (1987).Based on his theory, ln(z 0h /z 0m ) or ln(z 0q /z 0m ) are calculated as a function of the roughness Reynolds number, R * = u * z 0 /υ, where u * is the friction velocity, z 0 the roughness length and υ the kinematic viscosity of air.
Simulations with RACMO2 for the Antarctic region have shown that the original model configuration overestimates liquid precipitation at the expense of solid precipitation (Van de Berg et al., 2006).We imposed that clouds with temperatures below −7 • C form snow only, so that the solid precipitation flux increases, leaving the total precipitation sum unchanged.Due to the much lower air temperatures at the higher elevations, this correction only affects the lowest areas of the ice sheet.

Snow model
The original ECMWF surface scheme (TESSEL; Tiled ECMWF Surface Scheme for Exchanges over Land) does not make a distinction between the snow cover on an ice sheet and seasonal snow cover on the tundra.In TESSEL, snow cover is treated as a single layer on top of the soil or vegetation, which is in thermal contact with the underlying soil.This is acceptable for a transient snow layer over the tundra, but not for the semi-permanent ice sheet firn layer.Snow/firn processes such as meltwater percolation, retention and refreezing are not included, while these are especially important to realistically simulate the SMB of an ice sheet with extensive summertime melting and refreezing (Genthon, 2001).
For a better representation of the processes affecting the SMB in RACMO2/GR, we introduced an additional surface tile "ice sheet" in the land surface scheme TESSEL to describe the interaction at the snow/firn/ice-atmosphere interface (Fig. 2).As the ice temperature at the bottom of the ice/firn/snow pack is kept constant, no heat flux is assumed through the lower boundary.The subsurface processes are parameterized for at least the upper 30 m with a multi-layer snow/firn/ice model (1-D) composed of a maximum of 100 layers, but of 40 layers on average.The meltwater formed at the surface is allowed to penetrate to deeper layers, where it may refreeze (internal accumulation) or runs off as described by Bougamont et al. (2005).
The optimal thickness of a snow/firn/ice layer increases linearly from 6.5 cm for the uppermost layer to 4 m for the lowermost layer.The layer thickness is continuously changing due to snow accumulation, sublimation/deposition, melting, internal accumulation and firn densification.The www.the-cryosphere.net/4/511/2010/The Cryosphere, 4, 511-527, 2010 vertical grid is adjusted by layer splitting when the layer thickness becomes more than 1.3 times its optimal thickness, or layer fusion when a layer is less than half of its optimal thickness, except for layers consisting of ice lenses in the firn.Snow/firn density ρ continually changes in time due to refreezing of capillary water (rain and meltwater) and the settling and packing of dry snow according to the empirical formulation by Herron and Langway (1980): with where a is the annual mean accumulation rate, R the universal gas constant and T the firn/snow temperature in K.The annual accumulation rate used in this formula is the spatially distributed accumulation averaged over the period 1989-2005 based on a 16-year integration with RACMO2/GR.The snow/firn/ice column is thermally coupled to the atmospheric part of RACMO2/GR through a surface skin layer formulation of the surface energy balance (SEB) and the surface albedo, α, which is also applicable to the other surface tiles, such as tundra, sea-ice and open ocean.The skin temperature is introduced for modelling purposes and is defined as the temperature of the skin layer at the surface-atmosphere interface that is infinitely thin, has no heat capacity and responds instantaneously to SEB changes.The skin temperature T s is solved by SEB closure (e.g.Brutsaert, 1982): where M is the melt energy, SW net , SW ↓ , SW ↑ , LW net , LW ↓ , LW ↑ the net, downward and upward directed fluxes of shortwave and longwave radiation, α the broadband surface albedo, the surface emissivity for longwave radiation ( = 0.98 in RACMO2/GR for the ice sheet), σ the Stefan-Boltzmann's constant, LHF and SHF the turbulent fluxes for latent and sensible heat and G s the subsurface conductive heat flux evaluated at the surface.All terms are defined as positive when directed towards the surface-atmosphere interface.
The skin temperature serves as a boundary condition to the englacial module, which treats the vertical conduction of heat as follows: where ρ is the density of the snow/firn/ice layer, c p the specific heat capacity of ice (2009 J kg −1 K −1 ), ∂T /∂t the rate of temperature change within one model time step, k the effective conductivity, z the vertical coordinate and Q the heat released by refreezing of meltwater.The term ∂G/∂z accounts for the heat diffusion driven by the vertical temperature gradient.The snow/firn/ice conductivity follows the densitydependent approach of Van Dusen (1929), which ensures the correct value for k if ice density is attained.Temperature dependence of k is neglected: Knowing the conductivity of the snow/firn/ice layers, the vertical snow/ice temperature profiles can be computed.If T s is larger than 0 • C, it is reset to the melting point of ice and the excess of energy is used for melting.Meltwater and rain are allowed to percolate into the firn until they refreeze or run off.The maximum retention capacity due to capillary forces is set to a low value of 2% of available pore space, to obtain a realistic densification rate by refreezing of capillary water (Greuell and Konzelman, 1994).If an ice surface is encountered, the remaining water runs off at the surface, or deep in the firn pack at the snow/ice transition, without delay.The snow/firn/ice albedo α follows the snow density (ρ) and cloudiness (n) dependent linear formulation of Greuell and Konzelman (1994) for the uppermost 5 cm of the snow/firn/ice pack.
where the subscript i denotes ice and subscript s denotes snow.This parameterization is based on the notion that density reflects the metamorphosis state of the snow, i.e., it represents mostly the effects of grain size on albedo.Fresh snow is characterised by a surface α of 0.825 and a density of 300 kg m −3 .Glacier ice has an albedo of 0.5 and a density of 900 kg m −3 .Refrozen meltwater or rain may increase the density of the firn pack to the ice density, but the surface albedo is limited to a minimum value of 0.7 for refrozen water (Stroeve et al., 2005).This limitation will mainly affect areas south of 70 • N, where daytime melt and nighttime refreezing occur regularly throughout the melt season.

Model initialization
The atmospheric profiles of temperature, specific humidity, wind speed and surface pressure are initialized from ERA-40 at the beginning of the integration.By starting the simulation at the end of the melting season, the tundra could realistically be prescribed as snow free.Over the ice sheet, it is important to initialize the snow/ice temperature and snow/firn density with fairly realistic profiles, since typical timescales for changes in the snow/firn/ice pack are large, in the order of decades.During the 51-year simulation, no model parameters were re-initialized.
In the dry-snow zone, where melting is rare, the mean air temperature is a reasonable approximation (within 2 • C) for the climatological deep snow and ice temperature.For this reason, the snow/firn/ice temperature is initialized vertically uniform with the climatological surface temperature as described by the empiral function of Reeh (1991), who presented a snow/ice temperature parameterization as a function of elevation and latitude based on air temperature data from Danish meteorological stations at the periphery of the ice sheet for the 1951-1961 period: with where T is the climatological ice temperature in • C, T 2 m the 2 m air temperature in • C that depends on elevation z in m and latitude φ in • N, and δ T a perturbation due to the amount of superimposed ice formed, SIF.For SIF, the melt rate is averaged over the period 1989-2005 based on a 16year integration of RACMO2.For the percolation and ablation zones, a temperature correction δ T due to refreezing energy is included in line with Reeh (1991), and the ice temperature is limited to 0 • C. The resulting deep ice temperature serves as a boundary condition for the lowest firn/ice layer, so no heat flux is allowed in the underlying ice or soil.
For the 51-year model simulation, the initial temperature and density profiles of the snow/firn/ice column were obtained by rerunning the first model year (1 September 1957 to 31 August 1958) three times to reduce spin-off effects.Analysis of the three spin-up years and the first years of the simulations shows that the initial snow pack is in a state of near-balance before the present-day climate run is started.

Observational data
A proper assessment of RACMO2/GR output is essential before its data can be used as a tool for studying the climate of Greenland and the recent changes.Moreover, identification of model deficiencies may help to improve the model formulation for future climate simulations.To verify the model results for the near-surface conditions, we use: (i) nearsurface air temperature and wind speed data from automatic weather stations (AWSs) on the ice sheet (GC-net; Steffen and Box, 2001 and K-transect;Oerlemans and Vugts, 1993) and from climate stations of the Danish Meteorological Institute (DMI) on the surrounding tundra, (ii) data of surface radiation and heat exchange processes from three K-transect AWSs ( Van den Broeke et al., 2008a,b).
Statistical procedures were applied to all observational datasets to remove occasional spurious data values.For model evaluation of monthly means, we require that at least 80% of the observations are available during one month.The length of an observational record does not influence the evaluation, since every separate month is compared independently with the same month from the model output.The elevation of model grid points closest to all observational sites is within 100 m of the observed elevation, suggesting that no height correction is needed for temperature.

GC-net
The Greenland Climate Network (GC-net) was started in 1995 and consisted of 15 AWSs until 2001 (indicated as squares in Fig. 1) near or above the 2000 m elevation contour.Station coordinates and detailed information on the measurements are given in Steffen and Box (2001).We obtained a complete and quality controlled dataset for the period 1998-2001.For this period, the biases were removed and necessary corrections were applied.As the quality of the observations for the more recent years could not be guaranteed, this dataset nor the dataset from the DMI stations, are extended.
Four parameters derived from direct observations are compared with the RACMO2/GR output: 2 m air temperature, 10 m wind speed, net shortwave radiation and net radiation, as they are described by Box and Rinke (2003).The air temperature at 2 m is calculated by using the observed temperatures at 2 levels, heights of the instruments (median heights are 1.4 and 2.6 m) and linear interpolation.A logarithmic wind profile with a roughness length of 0.5 mm is assumed to estimate the 10 m wind speed.Due to riming of the sensors, net shortwave radiation data are omitted for the springtime months March and April.Most of the available net radiation observations are excluded in this study, because these unventilated measurements often suffer from large errors due to riming inside and outside the polyethylene domes.Only the net radiation records of the sites Swiss Camp and JAR1 are believed to be reliable throughout the year.

K-transect
As part of GC-net, UU/IMAU installed three AWSs along the Kangerlussuaq transect (K-transect) in southwest Greenland in August 2003 ( Van den Broeke et al., 2008b,c) (indicated as circles in Fig. 1).Measurements have been compared to model output for the period August 2003 to August 2007.The AWSs at S5 (490 m a.s.l.), S6 (1020 m a.s.l.) and S9 (1520 m a.s.l.) are located in the ablation and percolation zone (Fig. 3).The surface at S5 is very irregular with 2-3 m high ice hummocks usually covered with a thin layer of drift snow during wintertime, while at S9 the surface is much smoother, covered by a layer of wet snow for most or all of the summer.The changing surface conditions throughout the year make this dataset valuable for a thorough model evaluation on a daily basis.
For brevity, detailed daily evaluation is only shown for S6.Monthly and seasonal means of all three sites are used to www.the-cryosphere.net/4/511/2010/The Cryosphere, 4, 511-527, 2010 assess the model performance for the seasonal cycle.The comparison of daily values is focused on the year 2004, an average year within the 51-year simulation.
The accuracy of the measured temperature and wind speed at approximately 2 and 6 m is 0.3 • C and 0.3 m s −1 , respectively, as stated by Van den Broeke et al. (2008c).As the transformation to the 2 m temperature is only done when both measurements were available and by applying the bulk method, errors in the transformation are small.Further information on the sensor specifications and data quality is described in Van den Broeke et al. (2008a).
The observed surface radiation balance, surface characteristics, cloud properties and surface energy fluxes are derived from the AWS data with a melt model as described by Van den Broeke et al. (2008a,b).The observed (corrected) net shortwave radiation and the incoming longwave radiation fluxes serve as direct input for this melt model.The measurements of wind speed, temperature and humidity at two levels (approx. 2 and 6 m) serve as input for the bulk method to calculate the sensible and latent turbulent heat fluxes (Deardorff, 1968;Van den Broeke, 1996).

DMI climate stations
DMI climate stations are operated around the Greenland periphery (indicated as triangles in Fig. 1) and provide daily records of wind speed, air temperature and precipitation (Cappelen et al., 2001).For the model evaluation we used the dataset as described by Yang et al. (2005), which comprises of measurements during the period 1 January 1973 to 1 February 2005.Data from 51 stations is compared with model output for the nearest grid point that is considered as land in RACMO2/GR.As a result, some stations on small islands or narrow peninsulas are excluded from the analyses.
Model evaluation is limited to annual and climatological means because of the inability of the 11 km model grid to resolve local complex terrain surrounding the land stations.We computed monthly means of the wind speed and temperature, and averaged them over a year or over the measuring period to obtain an annual mean or climatological value for each site for comparison with RACMO2/GR output.

Model evaluation
The comparison of model values that represent averages for a model grid cell with a typical area of 121 km 2 , with local point observations must be done carefully.The model grid box closest to the observational site does not necessarily have the same surface type, elevation, surface roughness or surface albedo.In the interior of the ice sheet, these discrepancies are smaller since the surface is more homogeneous and the climate gradients less steep.
Model evaluation is performed based on daily, monthly and climatological averages at several sites on and across the ice sheet.RACMO2/GR data are saved at 6 hourly intervals.This 6 h resolution of the model output does not allow a thorough assessment of the modelled daily cycle.For this analysis, the model output has not been post-calibrated.The model elevation bias (modelled minus observed values) at almost all measurement sites is smaller than 100 m, and as a result no elevation-based correction is applied to the model output.Evaluation of the temporal evolution on a daily basis means that the weather conditions become critical, small differences in, for example, cloudiness or surface conditions may introduce large discrepancies in the lower atmosphere.As the year 2004 was not an exceptional year within the 51year simulation, the comparison of daily model output with observations is focused on this year.Monthly averages are used for evaluation of the seasonal cycle and yearly averages for verification of the model temporal evolution and climatological values.As most observations are only available for the most recent years, the model evaluation is focused on the end of the 51-year simulation.

Temperature at 2 m
The near-surface or 2 m temperature (T 2 m ) is an important climate variable, and one of the primary variables used in climate change reports as it is measured at many sites across the globe.Moreover, the near-surface saturation specific humidity, and consequently also sublimation/deposition at the surface, all strongly depend on the near-surface temperature.Typical for the interior of the ice sheet is a surface temperature inversion, driven by surface radiative cooling and in part compensated by the downward (air-to-surface) transport of sensible heat (SHF).This temperature deficit drives a persistent katabatic wind circulation over the ice sheet (Steffen and Box, 2001).
Figure 4a shows that for the entire ice sheet (green and red dots) and the surrounding tundra (black dots), the simulated climatological values of T 2 m are in close agreement with the observations (R = 0.97) with an averaged bias of −0.8 • C. The model tends to slightly underestimate/overestimate the near-surface temperature on the tundra/ice sheet.The averaged land bias is −1.5 • C (R = 0.96), whereas the ice sheet bias is +0.9 • C (R = 0.99).Only at some of the locations along the coastline of Greenland, does RACMO2/GR deviate more than 4 • C from the observations.The largest model bias (−9.8 • C) is found for DMI station 43800, located along the southeast coast near Tingmiarmiut.Disregarding this station reduces the root mean square error (RMSE) of 2.3 • C to 2.0 • C when taking all locations into account, and from 2.1 to 1.7 • C for only the land sites.
The temperature bias is uncorrelated to the elevation bias and does not show coherent regional patterns because of the irregular distribution of the stations over Greenland, but seems to be correlated to the land surface type.In RACMO2/GR, tundra and ice sheet are considered as different surface tiles with specific characteristics, such as albedo, thermal skin conductivity and vegetation type.The calculation of the surface fluxes is done separately for these different surfaces, leading to different solutions for the SEB equation and skin temperature even if the overlying atmosphere would be identical.A similar inland warm bias has been identified in ERA-40 data (Hanna et al., 2005), in part ascribed to positive bias in downward longwave radiation from the Rapid Radiative Transfer Model (RRTM) scheme, which is also used in RACMO2/GR.
Figure 5 shows the observed and modelled 2 m temperature deviations from their annual mean value  for 4 long-term DMI stations at various locations around the ice sheet.The model closely follows the observed temperature over the measurement period, also over the most recent years when warming has been reported Hanna et al. (2008); Box et al. (2009).Comparison of the long-term measurements at all climate stations with the model output indicates that the land bias (ranging from −4.4 to 0.8 • C) is stable in time, so that the interannual variability is well captured by RACMO2/GR.The standard deviation of the observations is for 3 out of 4 shown stations larger than the modelled standard deviation, which is valid for the whole climate stations dataset.This points towards a systematic underestimation of the interannual variability by RACMO2/GR for the land stations, rather than an increasing model drift due to incorrect initializations.
To assess the seasonal cycle over the ice sheet ablation zone, Fig. 4b shows the differences between the monthly modelled and observed temperatures along the K-transect over the period September 2003-August 2007.Additionally, Table 1 shows the seasonal biases and observed standard deviation based on daily values for all three K-transect locations.During summer, the standard deviation is considerably smaller, because the surface temperature is limited to the melting point, reducing the seasonal variability.For two sites along the K-transect, S6 and S9, the mean monthly bias is 1.1 and 0.5 • C and the RMSE 0.5 and 0.7 • C, respectively.These biases and RMSE are considerable smaller than one standard deviation, which indicates that RACMO2/GR is capable of simulating the temporal variability.The warm bias is stable through the year (Table 1), except for winter (DJF) at S9 (−0.2 • C), indicating that the seasonal cycle is well captured.A similar realistic seasonal cycle in T 2 m is found for the low-elevation sites of GC-net, Swiss Camp and JAR1 (not shown).On a daily basis, Fig. 6a shows that for site S6 the difference between the observed values and RACMO2/GR is generally low for the year 2004 (RMSE = 1.9 • C).The model follows the observed temporal evolution closely throughout the year.The large day-to-day fluctuations of over 10 • C during the winter are well represented in the model output, indicating that RACMO2/GR is capable of simulating the variability in weather and the related changing atmospheric conditions over the ice sheet.The largest model biases are found in the transition months April and September, which is associated with an underestimation of the surface albedo leading to more net shortwave radiation absorption (see Sect. 4.4.1).Similar results are found for the other years.
At the lowest site S5, RACMO2/GR shows a pronounced cold monthly bias of up to 4 • C, especially in wintertime (Table 1).Here, the mean monthly bias is −2.6 • C. Compared to S6 and S9, the surroundings of S5 are more complex.S5 The Cryosphere, 4, 511-527, 2010 www.the-cryosphere.net/4/511/2010/ is located at only 6 km from the ice sheet margin on an ice tongue (Russell Glacier) that protrudes from the ice sheet onto the tundra.Its closest model grid point is classified as ice sheet, while some of its neighbouring grid points are classified as tundra.The 1 • C summer cold bias at S5 may be caused by too much nocturnal cooling of the surface in the model, whereas the ice surface is observed to be at melting point day and night.In winter, it is well known that temperatures over flat tundra are considerably lower than over the adjacent ice sheet, where katabatic winds prevent the formation of a strong temperature inversion (e.g.Van den Broeke et al., 1994).Therefore, winter temperature biases at S5 are thought to result from insufficient downward longwave radiation and/or overestimation of cold air pooling over the tundra.

Wind speed and direction at 10 m
To assess the model performance for wind over the whole ice sheet, we compare RACMO2/GR with in situ observations averaged over matching time periods (Fig. 7a).Both low and high wind speeds are well represented with a mean difference of only 0.3 m s −1 (RMSE = 1.9 m s −1 ).This suggests that the surface friction is adequately accounted for in the model and that the vertical resolution of the model with its lowest layer at about 10 m above the surface is sufficient for simulating the near-surface katabatic wind profile, as found by Reijmer et al. (2005) for Antarctica.The monthly mean observed standard deviation (2.9 m s −1 ) is considerably larger than the mean bias and RSME, which implies that RACMO2/GR is capable of simulating the near-surface wind speed variability.
The correlation between the model output and the observations is high (R = 0.74), considering that the measured wind speed may be affected by local topography.Furthermore, a considerable uncertainty exists in both the in situ and model wind speed at 10 m owing to poorly defined stability corrections in very stable surface layers, which regularly occur over the interior of the ice sheet.In situ sensors also occasionally accumulate rime, which could be expected to introduce a negative wind bias.Because the AWSs are un-attended, it is impossible to quantify how large this error is.
The seasonal cycle of wind speed is largely controlled by the strength of the katabatic forcing, which is largest in winter (Van de Wal et al., 2005).Along the K-transect, the surface is considerably smoother at S9 than at S5 and S6 (Fig. 3).As a result the strongest seasonal cycle is found at S9 with monthly averaged summer wind speeds of 6 m s −1 and 11 m s −1 during February.Averaged over the K-transect, the modelled 10 m monthly wind speed deviates less than 1 m s −1 from the observations (Fig. 7b).Similar results are found for the different seasons.At S5 and S9 the averaged seasonal bias is uniform over the year and slightly negative (−0.4 and −0.3 m s −1 , respectively), but considerably smaller than the observed standard deviation (2.5 and 3.1 m s −1 ).At S6 the seasonal biases are close to zero, except for summer (bias = 0.7 m s −1 ), probably due to an inaccurate transition of snow to bare ice (see Sect. 4.4.1).At these lower elevations, the estimates of 10 m wind speed based on similarity theory may be more reliable, because enhanced turbulent mixing due to increasing wind speeds minimizes the stability effects.On a daily basis, the mean bias between the modelled and observed 10 m wind speed at S6 is 0.7 m s −1 for 2004 (Fig. 6b).The RMSE of daily means is 1.6 m s −1 for the 2003-2007 period.In summer, the daily 10 m wind speed is overestimated (bias = 1.1 m s −1 ) during both high and low www.the-cryosphere.net/4/511/2010/The Cryosphere, 4, 511-527, 2010 wind speed events, possibly due to a too low modelled surface roughness.A remarkable feature is the daily averaged wind speed, which is always above 1 m s −1 apart from a short period during which the sensor was frozen.This is because a continuous surface temperature inversion develops owing to negative net surface radiation in winter and a surface temperature restricted to the melting point in summer, causing a persistent katabatic wind throughout the year over the sloping surface of the ice sheet.
The wind regime on the ice sheet is dominated by semipermanent katabatic winds (Steffen and Box, 2001).Katabatic winds are characterised by (a) a maximum in wind speed close to the surface and (b) a constant wind direction.The directional constancy dc is a useful tool to detect local persistent circulations and is defined as the ratio of the vector-averaged wind speed to the mean wind speed usually taken at 10 m (Bromwich, 1989): where u and v are the horizontal components of the 10 m wind.A dc of zero implies that the near-surface wind direction is random.When dc approaches 1, the wind blows increasingly from the same direction.Close to the ice margin, the directional constancy and wind speed peak twice a year.In winter, the katabatic wind forcing is maintained by the radiation deficit at the surface, whereas in summer, the snow/ice at the surface melts and prevents the surface temperature from rising above melting point, so that katabatic winds persist.For S6, RACMO2/GR underestimates the persistence of the katabatic flow by ∼5% on average (Fig. 6c), but the double annual maximum is well (R = 0.9) represented.
The mean wind direction along the K-transect is southsoutheasterly (Fig. 8).This dominant wind direction is determined by storms and the persistent katabatic flow that is deflected to the right of the downslope direction due to the Coriolis force.A downslope (cross-isobar) component is maintained by friction.The wind direction is well simulated by RACMO2/GR, although it is too strongly (26 degrees on monthly basis) deflected at S9, possibly due to an underestimated surface roughness length.

Humidity at 2 m
The near-surface specific humidity is strongly controlled by air temperature.Along the K-transect, higher elevated sites have lower average specific humidity, modelled and observed.When specific humidity is high, temperatures are also high and visa versa, which follows the essential Clausius Clapeyron function.
Figure 9a, shows that at S6, the agreement between the daily RACMO2/GR values and observations of specific humidity is good (R = 0.98), both for the very low values during winter (<1 g kg −1 ) and for the maximum values during summer (≈4 g kg −1 ).The bias is rather constant throughout the year, also for the other years within the measurement period (bias = −0.05g kg −1 ; RMSE = 0.26 g kg −1 ).The seasonal variability is well captured as the daily modelled humidity follows the The Cryosphere, 4, 511-527, 2010 observations closely (Fig. 9a).The observed and modelled standard deviations are identical (1.42 g kg −1 ), and considerably larger than the above-mentioned bias and RMSE.At S9, the bias and RMSE are even smaller (bias = −0.005g kg −1 ; RMSE = 0.27 g kg −1 ).For S5, RACMO2/GR performs slightly worse (bias = −0.25 g kg −1 ; RMSE = 0.35 g kg −1 ).This bias is also persistent throughout the year.When analysing the 2 m relative humidity RH 2 m , it appeared that in the standard post-processing of RACMO2 data, the latent heat of vapourization is used for the computation of the saturated vapour pressure as prescribed by the WMO (World Meteorological Organization), whereas sublimation/deposition takes place at freezing winter temperatures.Since the observed RH 2 m is derived using the latent heat of sublimation, RACMO2/GR would significantly underestimate RH 2 m by −14.4%.Therefore, we recomputed modelled RH 2 m using the daily specific humidity model values and the latent heat of sublimation, which reduced the mean daily bias to −7.2%.The observed RH 2 m at S6 remains close to saturation throughout the year, while RACMO2/GR shows an unexpected decrease in wintertime (Fig. 9b).This discrepancy is also found for the observational years 2005 and 2006.In summer, both observed and modelled RH 2 m decrease towards the lower elevations (not shown).A possible explanation is that the katabatic wind transports colder, dry air downwards and that adiabatic compression and the associated heating results in a lower relative humidity downslope in summer.Measurement uncertainties at low temperatures are also a possible explanation.

Surface energy balance
The air temperature near the surface is strongly coupled to the surface temperature T s , which is determined by the surface energy balance (SEB).The SEB (Eq. 3) voor the GrIS is largely controlled by the radiative fluxes and the surface albedo, and to a lesser extent by the turbulent fluxes and the subsurface heat flux (Van den Broeke et al., 2008b,a).The performance of RACMO2/GR for different terms in the SEB will be discussed in this order.Few reliable measurements of SEB components on the ice sheet are available.We rely on SEB observations along the K-transect, where the AWSs are equipped with K&Z CNR1 radiation sensors that measure all four radiation components individually.

Net shortwave radiation and surface albedo
The SEB is strongly influenced by net shortwave radiation that is absorbed at the surface and which drives a clear seasonal and diurnal cycle unless the energy is used for melting.Along the K-transect, the model bias in SW ↓ is time-dependent.While RACMO2/GR estimates SW ↓ to be 126 W m −2 for all three sites, the observations are less uniform.A positive model bias of +14 W m −2 (11.2%) is found at S5 and a negative bias of −10 W m −2 (7.8%) at S9. Inaccuracies in modelled clear-sky transmissivity, clouds and/or cloud/radiation interactions in RACMO2/GR can cause these deviations from the observations.Quantification of a bias in www.the-cryosphere.net/4/511/2010/The Cryosphere, 4, 511-527, 2010 each of these processes separately cannot be clarified without detailed cloud-radiation observations and modelling.
The reflected shortwave radiation depends on the amount of incident shortwave radiation at the surface and the surface albedo.The latter is observed to be asymmetric through the year in the ablation zone (Van den Broeke et al., 2008a).Comparing daily model output with the K-transect observations reveals a too early decrease and a too late increase in modelled α, ranging from only a few days up to weeks (Fig. 10) for all evaluated years.In early summer, the winter snow pack melts, leading to a transition from a dry snow pack (modelled α of 0.825) to a wet snow pack with modelled α of ≈0.7, followed by the appearance of the underlying glacier ice with modelled α of ≈0.5.The rate of this transition process is hard for RACMO2/GR to capture, since the modelled surface albedo is determined based on the density of the upper 5 cm of dry snow, unaffected by the presence of water in the snow pack.Furthermore, in reality, some redistribution of falling snow by the wind occurs (Van den Broeke et al., 2008a).The radiation sensor is mounted on the AWSs that stands on top of an ice hummock (Fig. 3) and, thus, there is a likely sampling bias toward lower albedo, especially in the early melt season.
The observed daily variations in α associated with snowfall events are underestimated by the model (Fig. 10).In the observations, α rises more abruptly during a snowfall event, even if only a very thin layer of fresh snow covers the surface.In the model, α responds only to significant changes in the density of the upper 5 cm of the snow/firn/ice pack, which requires a more substantial snowfall event.The same discrepancy between model and observations is responsible for the late increase in model α during autumn, as fresh snow starts to cover the glacier ice.Similar systematic biases are found for the other years of the measurement period.The timing of the spring melt and of the fresh snowfall in autumn does change for the different years, but the time lag between the model and observations is similar (not shown).Overall, the surface albedo evolution through all four summers (2004)(2005)(2006)(2007) is captured reasonably well (quantified below) by RACMO2/GR (R = 0.73), taking into account that the ablation zone is characterised by a very inhomogeneous surface.
The underestimation of the albedo in early summer and autumn leads, on average, to a positive model bias in the reflected shortwave radiation of +9 W m −2 averaged over the K-transect (not shown).In the ablation zone, the positive biases in the reflected shortwave radiation lead to an overestimation in the net shortwave radiation, with the largest biases in the spring and summer months (Fig. 11b and Table 2).Figure 12a shows that RACMO2/GR significantly overestimates SW net at S6 by 31% in summer compared to the observations.As expected, the bias in SW net is smaller for most of the dry snow zone (GC-net stations in Fig. 11a), where the surface albedo remains relatively high and constant throughout the year.Only a significant deviation from the assumed  fresh snow α of 0.85 may result in an overestimation at the accumulation zone sites.

Net longwave radiation
At S6, the daily variation in net longwave radiation LW net is well captured by RACMO2/GR (Fig. 12b).The model tends to underestimate the lower range of values during the winter months (Table 2).This negative bias is caused by an underestimation of LW ↓ , with as largest bias −30 W m −2 .Van de Berg et al. ( 2007) encountered a similar problem over the Antarctic ice sheet using an earlier version of RACMO2, which they related to an underestimation of the clear-sky radiance, winter cloud cover and humidity.Similarly to biases in SW ↓ , detailed cloud observations are needed to quantify the effect of a potential bias in cloud properties on LW ↓ .At S6, the resulting winter negative bias in LW net is 16 W m −2 (Table 2), whereas the monthly average bias in LW ↑ is only ±5 W m −2 (Fig. 13a).In summer, the LW ↑ bias diminishes as the melting surface limits the surface temperature.For S9, the performance of RACMO2/GR is similar to S6.For S5 however, the cold bias (see Fig. 4b) results in an underestimation of LW ↑ in winter of 25 W m −2 , compensating for the bias in LW ↓ (Fig. 13a).

Net radiation
In Fig. 12c, the net result of the daily shortwave and longwave radiation fluxes is presented for S6.In wintertime, shortwave radiation is reduced to near zero and LW net drives the surface radiation budget.The negative bias in LW ↓ leads to an underestimation of net radiation and is thought to be the result of underestimated clear sky longwave radiance and/or of cloudiness (see Sect. 4.4.2).In summer, the positive bias in SW net is the dominant contribution to an overestimation of the net radiation absorbed at the surface.Figure 13b shows that for S6 the largest disagreement is found in spring, when the negative bias in albedo is largest.At S5 and S9, the bias in net radiation is smaller due to a better representation of the surface albedo variability in the summer half year.A similar bias is found for the GC-net sites JAR1 and Swiss Camp that are located in environments comparable to S9 (not shown).The correlation between net radiation observed at 20 ice sheet locations and modelled is 0.79 with climatological mean bias of 2.5 W m −2 and RMSE of 3.3 W m −2 .

Turbulent heat fluxes
Figure 14a shows that the daily sensible heat flux SHF at S6 is positive throughout the year, which indicates that the atmosphere continuously transfers heat to the surface.The double maxima (winter and summer) correspond to the maxima in wind shear and temperature gradient between the surface and atmosphere, which are coupled through the katabatic forcing.
During winter, RACMO2/GR simulates an excess SHF compared to observations of 20 W m −2 at S6 and S9 (Fig. 15a and Table 2).This balances most of the surplus in net LW cooling, explaining the realistic near-surface temperatures at these sites (Fig. 6a).It is known that the mixing scheme in RACMO2/GR is too active, especially under very stable atmospheric conditions (Van Meijgaard et al., 2008).The winter bias in SHF is smaller at S5 (10 W m −2 ), because this site is closer to the ice margin and affected by a deeper katabatic wind circulation, so the modelled and observed mixing layer depth are more similar.Here, the excess LW cooling during winter is only partly compensated by the overestimated SHF.During the summer, the largest positive bias is found at S6 (about +20 W m −2 ), while at S5 and S9 the biases (−4.9 and +4.1 W m −2 , respectively) are much smaller.The annual cycle of latent heat flux LHF is of importance to the SEB.Surface temperatures continuously below freezing lead to deposition (rime formation) in winter and sublimation in spring and summer (Fig. 14b).To obtain a realistic sublimation, it is important that at least the surface temperature is correctly represented.Differences in LHF between RACMO2/GR and observational sites along the K-transect are less than ±5 W m −2 in winter months and about 5 W m −2 during summer (Fig. 15b and Table 2).The annual bias is −2.0 W m −2 averaged over these 3 sites.The largest monthly biases are found at S5, coinciding with a large T 2 m bias.It should be noted here that "observed" turbulent fluxes are approximated by the bulk fluxes, which are also somewhat uncertain (Box and Steffen, 2001).

Summary and conclusions
An assessment of the performance of RACMO2/GR, a regional climate model with physical parameterizations optimized for use over the extensive ice sheets, is pre- sented using in situ observations on and around the Greenland ice sheet.This analysis has primarily focused on the near-surface atmospheric state (temperature, humidity, wind speed and direction), and the surface energy balance components including the radiative fluxes.We found a good correlation (bias = −0.8• C, R = 0.97, RMSE = 2.3 • C) between modelled and measured climatological value of T 2 m at 70 stations across the ice sheet.The temperature climatological bias seems correlated with land surface type, as a persistent warm/cold bias is found over the ice sheet/tundra of +0.9 and −1.5 • C, respectively.The largest monthly bias (−5 • C) occurs for winter near the ice margin, whereas in the higher ablation zone and in the percolation zone, the temperature is well captured.
The difference between modelled and measured wind speed appears to be substantial at several locations, caused by local topography, but generally the agreement is reasonable  (bias = 0.3 m s −1 , R = 0.74, RMSE = 1.9 m s −1 ).At about 60 out of the 70 stations, the difference in climatological mean 10 m wind speed is smaller than 2 m s −1 .Local topographical conditions at the stations and smoothing of steep terrain in the model make it difficult to directly compare the near-surface winds with model values, especially for the land sites.The force and persistency of the katabatic wind circulation is well captured by the model.The small deviations in wind direction in the ablation zone are probably caused by differences in surface roughness lengths in RACMO2/GR.The surface energy balance is evaluated using observations from three AWSs along the K-transect and AWSs of GC-net for which high-quality measurements were made available.The modelled net shortwave radiation flux matches the observations reasonably well (R = 0.79) in the dry snow zone, whereas it is overestimated in the ablation and percolation zone.The snow model has difficulties in simulating the instant decrease in surface albedo due to wetting and melting of snow and the sudden increase when a thin layer of fresh snow covers the bare glacier ice.Determining the albedo based on the microphysical properties of the upper snow/firn/ice layer would be preferable to the empirical correlation between snow density of the upper 5 cm and the albedo used here.Keeping in mind that the surface in the ablation zone is very inhomogeneous, which reduces the representation of the single point observations for the model grid box at 11 km resolution, the snow model captures the changing surface conditions under melting conditions reasonably well.
It is known that RACMO2 underestimates the downwelling longwave radiation at low atmospheric temperatures, which is related to an underestimation of the clear-sky component and/or of humidity and cloud cover.This is confirmed by the measurements at the higher elevated K-transect sites, where the model bias reaches 20 W m −2 in winter.Radiation budget errors suggest that the largest source of uncertainty next to the surface albedo is cloud-radiation interactions.During winter, an excess SHF of 15 W m −2 balances most of the excess LW cooling, except for the lower ablation zone, where S5 is located.Under very stable conditions, the vertical mixing scheme is too active, which introduces a compensating error.As a result, only a small bias is found in the surface and 2 m temperature.
The model evaluation described here demonstrates that RACMO2/GR is capable of realistically simulating presentday near-surface characteristics of the Greenland atmosphere on daily and monthly timescales, without post-calibration or reinitialization during the 51-year simulation.This makes RACMO2/GR a suitable and valid tool to study recent climate changes over the Greenland ice sheet.

Fig. 1 .
Fig. 1.Map of Greenland featuring the model domain, relaxation borders (the outer 16 grid points represented as dark gray dots), location of model grid points (light gray dots) and location of observational sites.The 51 DMI climate stations are indicated by triangles, the 15 GC-net automatic weather stations by squares and the three K-transect AWSs by circles.Thin dashed lines are 250 m elevation contours from Bamber et al. (2001).The thick black line represents the ice sheet contour as used in RACMO2/GR.

Fig. 2 .
Fig. 2. Schematic representation of modelled processes that determine the surface mass balance.Upper and lower blue surfaces denotes snow-air and snow-ice interfaces, respectively.

Fig. 3 .
Fig. 3. Images of the AWSs along the K-transect and their surroundings at S5, S6 and S9.Images taken at the end of the ablation season (end of August).Photos by Paul Smeets (UU/IMAU).

Fig. 6 .
Fig. 6.Comparison of simulated (gray lines) and observed (black lines) daily averaged (a) 2 m temperature [ • C], (b) 10 m wind speed [m s −1 ] at S6 for the year 2004, and (c) comparison of simulated (gray lines) and observed (black lines) monthly averaged directional constancy [−] of 10 m wind at S6 for the period January 2004-August 2007.

Fig. 8 .
Fig. 8.Comparison of simulated (open circles) and observed (solid circles) monthly averaged 10 m wind direction and speed at S5, S6 and S9 for the measurement period August 2004-August 2007.

Fig. 10 .
Fig. 10.Time evolution of the daily surface albedo [−] in the observations (black lines) and model output (gray lines) for the three AWSs (S5, S6 and S9) along the K-transect for the period April-October 2004.

Fig. 11 .
Fig. 11.Model performance for surface net shortwave radiation [W m −2 ] (a) model versus observations for GC-net (black) and K-transect (green) averaged over the available measuring period, (b) monthly model bias for S5 (black), S6 (red) and S9 (blue) over the period August 2003-August 2007.

Fig. 12 .
Fig. 12.Comparison of simulated (gray lines) and observed (black lines) daily averaged values of (a) net shortwave radiation flux SW net , (b) net longwave radiation flux LW net , and (c) net radiation flux in [W m −2 ] at S6 for the year 2004.Note the different vertical scales used in the panels.

Fig. 13 .
Fig. 13.Model performance for (a) the net longwave radiation and (b) the net radiation for S5 (black), S6 (red) and S9 (blue) along the K-transect in [W m −2 ] over the period August 2003-August 2007.

TheLHFFig. 14 .
Fig. 14.Comparison of simulated (gray lines) and observational based (black lines) daily averaged surface (a) sensible heat flux SHF, and (b) latent heat flux LHF in [W m −2 ] at S6 for the year 2004.

Table 2 .
Seasonal and annual bias between the modelled and observed surface energy fluxes [W m −2 ] for the stations S5, S6 and S9 over the period August 2003-August 2007.