Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 1: Greenland (1958–2016)

Abstract

We evaluate modelled Greenland ice sheet (GrIS) near-surface climate, surface energy balance (SEB) and surface mass balance (SMB) from the updated regional climate model RACMO2 (1958–2016). The new model version, referred to as RACMO2.3p2, incorporates updated glacier outlines, topography and ice albedo fields. Parameters in the cloud scheme governing the conversion of cloud condensate into precipitation have been tuned to correct inland snowfall underestimation: snow properties are modified to reduce drifting snow and melt production in the ice sheet percolation zone. The ice albedo prescribed in the updated model is lower at the ice sheet margins, increasing ice melt locally. RACMO2.3p2 shows good agreement compared to in situ meteorological data and point SEB/SMB measurements, and better resolves the spatial patterns and temporal variability of SMB compared with the previous model version, notably in the north-east, south-east and along the K-transect in south-western Greenland. This new model version provides updated, high-resolution gridded fields of the GrIS present-day climate and SMB, and will be used for projections of the GrIS climate and SMB in response to a future climate scenario in a forthcoming study.

Introduction

Predicting future mass changes of the Greenland ice sheet (GrIS) using regional climate models (RCMs) remains challenging . The reliability of projections depends on the ability of RCMs to reproduce the contemporary GrIS climate and surface mass balance (SMB), i.e. snowfall accumulation minus ablation from meltwater run-off, sublimation and drifting snow erosion . In addition, RCM simulations are affected by the quality of the re-analysis used as lateral forcing and by the accuracy of the ice sheet mask and topography prescribed in the models .

Besides direct RCM simulations, the contemporary SMB of the GrIS has been reconstructed using various other methods, e.g. positive degree day (PDD) models forced by statistically downscaled re-analyses , mass balance models forced by the climatological output of an RCM (HIRHAM4) and reconstruction of SMB obtained by combining RCM outputs with temperature and ice core accumulation measurements . In addition, and respectively used the Community Earth System Model (CESM) at 1${}^{\circ }$ resolution ($\sim$ 100 km) and the Goddard Earth Observing System model version 5 (GEOS-5) at 0.5${}^{\circ }$ resolution ($\sim$ 50 km) to estimate recent and future mass losses of the GrIS.

Polar RCMs have the advantage of explicitly resolving the relevant atmospheric and surface physical processes at high spatial (5 to 20 km) and temporal (subdaily) resolution. Nonetheless, good RCM performance often results from compensating errors between poorly parameterised processes, e.g. cloud physics and turbulent fluxes . Therefore, considerable efforts have been dedicated to evaluating and improving polar RCM output in Greenland , using in situ SMB observations , airborne radar measurements of snow accumulation and meteorological records , including radiative fluxes that are required to close the ice sheet surface energy balance (SEB) and hence quantify surface melt energy.

For more than two decades, the polar version of the Regional Atmospheric Climate Model (RACMO2) has been developed to simulate the climate and SMB of the Antarctic and Greenland ice sheets. In previous versions, snowfall accumulation was systematically underestimated in the GrIS interior, while melt was generally overestimated in the percolation zone . At the ice sheet margins, meltwater run-off is underestimated over narrow ablation zones and small outlet glaciers that are not accurately resolved in the model's ice mask at 11 km. Locally, this underestimation can exceed several m w.e. yr${}^{-\mathrm{1}}$, e.g. at automatic weather station (AWS) QAS$\mathit{\_}$L installed at the southern tip of Greenland . These biases can be significantly reduced by statistically downscaling SMB components to 1 km resolution . Computational limitations currently hamper direct near-kilometre-scale simulations of the contemporary GrIS climate, making it essential to further develop RACMO2 model physics at coarser spatial resolution. Important modelling challenges and limitations still need to be addressed in RACMO2 regardless of the spatial resolution used: e.g. cloud representation , surface albedo and turbulent heat fluxes (Sect. 6).

Here, we present updated simulations of the contemporary GrIS climate and SMB at 11 km resolution (1958–2016). The updated model incorporates multiple adjustments, notably in the cloud scheme and snow module. Model evaluation is performed using in situ meteorological data and point SEB/SMB measurements collected across the GrIS. We then compare the SMB of the updated model version (RACMO2.3p2) with its predecessor (RACMO2.3p1), discussed in for the overlapping period between the two simulations (1958–2015). Section 2 discusses the new model settings and initialisation together with observational data used for model evaluation. Modelled climate and SEB components are evaluated using in situ measurements in Sect. 3. Changes in SMB patterns between the new and old model versions are discussed in Sect. 4, as well as case studies in north-eastern, south-western and south-eastern Greenland. Section 5 introduces and evaluates the updated downscaled daily, 1 km SMB product. Section 6 discusses the remaining model uncertainties, followed by conclusions in Sect. 7. This paper is part of a tandem model evaluation over the Greenland (present study) and Antarctic ice sheets .

Model and observational data

The Regional Atmospheric Climate Model RACMO2

The polar (p) version of the Regional Atmospheric Climate Model (RACMO2) is specifically adapted to simulate the climate of polar ice sheets. The model incorporates the dynamical core of the High Resolution Limited Area Model (HIRLAM) and the physics package cycle CY33r1 of the European Centre for Medium-Range Weather Forecasts Integrated Forecast System . It also includes a multilayer snow module that simulates melt, liquid-water percolation and retention, refreezing and run-off , and accounts for dry snow densification following . RACMO2 implements an albedo scheme that calculates snow albedo based on prognostic snow grain size, cloud optical thickness, solar zenith angle and impurity concentration in snow . In RACMO2, impurity concentration, i.e. soot, is prescribed as constant in time and space. The model also simulates drifting snow erosion and sublimation following . Previously, RACMO2 has been used to reconstruct the contemporary SMB of the Greenland ice sheet and peripheral ice caps , the Canadian Arctic Archipelago , Patagonia and Antarctica .

Surface energy budget and surface mass balance

In RACMO2, the skin temperature ($T_{\text{skin}}$) of snow and ice is derived by closing the surface energy budget (SEB), using the linearised dependencies of all fluxes to $T_{\text{skin}}$ and further assuming, as a first approximate, that no melt occurs at the surface ($M=\mathrm{0}$). If the obtained $T_{\text{skin}}$ exceeds the melting point, $T_{\text{skin}}$ is set to 0 ${}^{\circ }$C; all fluxes are then recalculated and the melt energy flux ($M>\mathrm{0}$) is estimated by closing the SEB in Eq. (1), assuming that no solar radiation can directly penetrate the snow or ice interface.

$$\begin{array}{rl}M& ={\mathrm{SW}}_{\mathrm{d}}-{\mathrm{SW}}_{\mathrm{u}}+{\mathrm{LW}}_{\mathrm{d}}-{\mathrm{LW}}_{\mathrm{u}}+\mathrm{SHF}+\mathrm{LHF}+G_{\mathrm{s}}\\ & ={\mathrm{SW}}_{\mathrm{n}}+{\mathrm{LW}}_{\mathrm{n}}+\mathrm{SHF}+\mathrm{LHF}+G_{\mathrm{s}},\end{array}$$ where SW${}_{\mathrm{d}}$ and SW${}_{\mathrm{u}}$ are the shortwave down-/upward radiation fluxes, LW${}_{\mathrm{d}}$ and LW${}_{\mathrm{u}}$ are the longwave down-/upward radiation fluxes, SHF and LHF are the net sensible and latent turbulent heat fluxes, and $G_{\mathrm{s}}$ is the subsurface heat flux. SW${}_{\mathrm{n}}$ and LW${}_{\mathrm{n}}$ are the net short-/longwave radiation at the surface. All fluxes are expressed in W m${}^{-\mathrm{2}}$ and are defined positive.

In the percolation zone of the GrIS, liquid-water mass from melt (ME) and rainfall (RA) can percolate through the firn column, and is either retained by capillary forces as irreducible water (RT) or refreezes (RF). Combined with dry snow densification, this progressively depletes firn pore space until the entire column turns into ice (900 kg m${}^{-\mathrm{3}}$). The fraction not retained is assumed to immediately run-off (RU) to the ocean: $$\mathrm{RU}=\mathrm{ME}+\mathrm{RA}-\mathrm{RT}-\mathrm{RF}.$$ The climatic mass balance , hereafter referred to as SMB, is estimated as $$\mathrm{SMB}=P_{\mathrm{tot}}-\mathrm{RU}-{\mathrm{SU}}_{\mathrm{tot}}-{\mathrm{ER}}_{\mathrm{ds}},$$ where $P_{\mathrm{tot}}$ is the total amount of precipitation, i.e. solid and liquid, RU is meltwater run-off, SU${}_{\mathrm{tot}}$ is the total sublimation from drifting snow and surface processes, and ER${}_{\mathrm{ds}}$ is the erosion by the process of drifting snow. All SMB components are expressed in millimetre water equivalent (mm w.e.) for point-specific SMB values, or in Gt yr${}^{-\mathrm{1}}$ when integrated over the GrIS.

Model updates

In the cloud scheme, parameters controlling precipitation formation have been modified to reduce the negative snowfall bias in the GrIS interior ($\sim$ 40 mm w.e. yr${}^{-\mathrm{1}}$) . To correct for this, the critical cloud content ($l_{\mathrm{crit}}$) governing the onset of effective precipitation formation for liquid-mixed and ice clouds has been increased by factors of 2 (Eqs. 5.35 and 6.39 in ) and 5 (Eq. 6.42 in ). As a result, moisture transport is prolonged to higher elevations and precipitation is generated further inland. The values of $l_{\mathrm{crit}}$ adopted in RACMO2 were obtained after conducting a series of sensitivity experiments, i.e. 1-year simulations, to test the dependence of precipitation formation efficiency, spatial distribution and cloud moisture content on $l_{\mathrm{crit}}$ and other cloud tuning parameters. From these experiments, we found a linear relationship between $l_{\mathrm{crit}}$ for mixed and ice clouds, the vertical integrated cloud content, i.e. liquid and ice water paths that also affect the SEB through changes in cloud optical thickness, and the integrated precipitation over Greenland. These new settings were then tested for a longer period and proved to almost cancel the dry bias observed in RACMO2.3p1 (see Sect. 5.1). This led to larger but realistic vertical integrated cloud content and did not strongly affect the SEB and surface climate of the GrIS. For instance, the induced changes of surface downward shortwave and longwave radiation are only about $-\mathrm{4}$ and 7 W m${}^{-\mathrm{2}}$, peaking in central Greenland. While the obtained increase in $l_{\mathrm{crit}}$ is relatively large, especially for ice clouds, it is important to note that it is also strongly adjusted in the original ECMWF physics compared to commonly used values in the literature: e.g. set $l_{\mathrm{crit}}$ to $\mathrm{1}\times \mathrm{10}$${}^{-\mathrm{3}}$ kg kg${}^{-\mathrm{1}}$ for ice clouds, while the ECMWF physics, tuned for GCM sized grid cells, uses a value of $\mathrm{0.3}\times {\mathrm{10}}^{-\mathrm{4}}$ kg kg${}^{-\mathrm{1}}$ . As $l_{\mathrm{crit}}$ depends on model grid resolution, i.e. GCMs running at lower spatial resolution require lower values of $l_{\mathrm{crit}}$ , the use of a larger $l_{\mathrm{crit}}$, e.g. for ice clouds ($\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{4}}$ kg kg${}^{-\mathrm{1}}$) in RACMO2, is deemed reasonable. In addition, this value remains well within the range of values previously presented in the literature .

Furthermore, the previous model version overestimated snowmelt in the percolation zone of the GrIS . With the aim of minimising this bias, the following parameters have been tuned in the snow module:

• a.

  The model soot concentration, accounting for dust and black carbon impurities deposited on snow, has been reduced from 0.1 to 0.05 ppmv, more representative of observed values . A lower soot concentration yields a higher surface albedo; hence melt decreases .

• b.

  The size of refrozen snow grains has been reduced from 2 to 1 mm . Consequently, the surface albedo of refrozen snow increases, as smaller particles enhance scattering of solar radiation back to the atmosphere .

• c.

  In previous model versions, the albedo of superimposed ice, i.e. the frozen crust forming at the firn surface, was set equal to the albedo of bare ice ($\sim$ 0.55), underestimating surface albedo and hence overestimating melt. The snow albedo scheme now explicitly calculates the albedo of superimposed ice layers ($\sim$ 0.75), following .

• d.

  The saltation coefficient of drifting snow has been approximately halved from 0.385 to 0.190 . Saltation occurs when near-surface wind speed is sufficiently high to lift snow grains from the surface. In RACMO2, this coefficient determines the depth of the saltation layer, i.e. typically extending 0 to 10 cm above the surface, that directly controls the mass of drifting snow transported in the suspension layer aloft (above 10 cm). This revision does not affect the timing and frequency of drifting snow events, which are well modelled , but only reduces the horizontal drifting snow transport and sublimation, preventing a too-early exposure of bare ice during the melt season, especially in the dry and windy north-eastern GrIS (Sect. 4.2).

Initialisation and set-up

To enable a direct comparison with previous runs, RACMO2.3p2 is run at an 11 km horizontal resolution for the period 1958–2016, and is forced at its lateral boundaries by ERA-40 (1958–1978) and ERA-Interim (1979–2016) re-analyses on a 6-hourly basis over the model domain shown in Fig. . The forcing consists of temperature, specific humidity, pressure, wind speed and direction being prescribed at each of the 40 vertical atmosphere hybrid model levels. To better capture SMB inter-annual variability in this new model version, upper atmosphere relaxation (UAR or nudging) of temperature and wind fields is applied every 6 h for model atmospheric levels above 600 hPa, i.e. $\sim$ 4 km a.s.l. . UAR is not applied to atmospheric humidity fields in order not to alter clouds and precipitation formation in RACMO2. As the model does not incorporate a dedicated ocean module, sea surface temperature and sea ice cover are prescribed from the re-analyses . The model has about 40 active snow layers that are initialised in September 1957 using estimates of temperature and density profiles derived from the offline IMAU Firn Densification Model (IMAU-FDM) . These profiles are obtained by repeatedly running IMAU-FDM over 1960–1979 forced by the outputs of the previous RACMO2.3p1 climate simulation until the firn column reaches an equilibrium. The data spanning the winter season up to December 1957 serve as an additional spin-up for the snowpack and are therefore discarded in the present study.

SMB (mm w.e. yr${}^{-\mathrm{1}}$) modelled by RACMO2.3p2 at 11 km resolution for 2016. Black dots delineate the relaxation zone (24 grid cells) where the model is forced by ERA re-analyses. Ablation sites (213) are displayed as yellow dots, accumulation sites (182) as white dots, and AWS locations (23) are represented in green.

[Figure omitted. See PDF]

Relative to previous versions, the integration domain extends further to the west, north and east (Fig. ). This brings the northernmost sectors of the Canadian Arctic Archipelago and Svalbard well inside the domain interior and further away from the lateral boundary relaxation zone (24 grid cells, black dots in Fig. ). In addition, RACMO2.3p2 utilises the 90 m Greenland Ice Mapping Project (GIMP) digital elevation model (DEM) to better represent the glacier outlines and the surface topography of the GrIS. Compared to the previous model version, which used the 5 km DEM presented in , the GrIS area is reduced by 10 000 km${}^{\mathrm{2}}$ (Fig. a). This mainly results from an improved partitioning between the ice sheet and peripheral ice caps, for which the ice-covered area has, in equal amounts, decreased and increased, respectively. In RACMO2, a grid cell with an ice fraction $\ge \mathrm{0.5}$ is considered fully ice covered. The updated topography shows significant differences compared to the previous version, especially over marginal outlet glaciers where surface elevation has considerably decreased (Fig. b). Bare ice albedo is prescribed from the 500 m Moderate Resolution Imaging Spectroradiometer (MODIS) 16-day albedo version 5 product (MCD43A3v5) as the lowest 5 $\mathit{\%}$ surface albedo records for the period 2000–2015 (vs. 2001–2010 in older versions; Fig. c). In RACMO2, minimum ice albedo is set to 0.30 for dark ice in the low-lying ablation zone, and a maximum value of 0.55 for bright ice under perennial snow cover in the accumulation zone. In previous RACMO2 versions, bare ice albedos of glaciated grid cells without valid MODIS estimates were set to 0.47 .

Difference in (a) ice mask (b) surface elevation and (c) bare ice albedo between RACMO2.3p2 and RACMO2.3p1. In panel (a), the common ice mask for both model versions is displayed in grey, the ice sheet area is outlined in yellow; additional and removed ice-covered cells in RACMO2.3p2 are shown in red and blue.

[Figure omitted. See PDF]

Observational data

To evaluate the modelled contemporary climate and SMB of the GrIS, we use daily average meteorological records of near-surface temperature, wind speed, relative humidity, air pressure and down-/upward short-/longwave radiative fluxes, retrieved from 23 AWSs for the period 2004–2016 (green dots in Fig. ). Erroneous radiation measurements, e.g. caused by sensor riming, were discarded by removing daily records showing SW${}_{\mathrm{d}\mathrm{bias}}>\mathrm{6}{\mathit{\sigma }}_{\mathrm{bias}}$, where SW${}_{\mathrm{d}\mathrm{bias}}$ is the difference between daily modelled and observed SW${}_{\mathrm{d}}$, and ${\mathit{\sigma }}_{\mathrm{bias}}$ is the standard deviation of the daily SW${}_{\mathrm{d}}$ bias for all measurements. In addition, measurements affected by sensor heating in summer, i.e. showing LW${}_{\mathrm{u}}>\mathrm{318}$ W m${}^{-\mathrm{2}}$, were eliminated as these values represent $T_{\mathrm{s}}>\mathrm{0}$ ${}^{\circ }$C for $\mathit{ϵ}\approx \mathrm{0.99}$, where $T_{\mathrm{s}}$ is the surface temperature and $\mathit{ϵ}$ the selected emissivity of snow or ice. We only used daily records that were simultaneously available for each of the four radiative components. Eighteen of these AWS sites are operated as part of the Programme for Monitoring of the Greenland Ice Sheet (PROMICE, www.promice.dk) covering the period 2007–2016 . Four other AWS sites, namely S5, S6, S9 and S10 (2004–2016), are located along the K-transect in south-western Greenland (67${}^{\circ }$ N, 47–50${}^{\circ }$ W) . Another AWS (2014–2016) is situated in south-eastern Greenland (66${}^{\circ }$ N; 33${}^{\circ }$ W) at a firn aquifer site . The latter five sites are operated by the Institute for Marine and Atmospheric research at Utrecht University (IMAU).

We also use in situ SMB measurements collected at 213 stake sites in the GrIS ablation zone (yellow dots in Fig. ; ) and at 182 sites in the accumulation zone (white dots in Fig. ) including snow pits, firn cores , and airborne radar measurements . We exclusively selected measurements that temporally overlap with the model simulation (1958–2016). To match the observational period, daily modelled SMB is cumulated for the exact number of measuring days at each site.

For model evaluation, we select the grid cell nearest to the observation site in the accumulation zone. In the ablation zone, an additional altitude correction is applied by selecting the model grid cell with the smallest elevation bias among the nearest grid cell and its eight adjacent neighbours. One ablation site and seven PROMICE AWS sites presented an elevation bias in excess of $>$ 100 m compared to the model topography and were discarded from the comparison.

In addition, we compare modelled SMB with annual glacial ice discharge (D) retrieved from the combined Zachariae Isstrøm and Nioghalvfjerdsbrae glacier catchments in north-eastern Greenland (1975–2015; yellow line in Fig. a), presented in .

Comparison between modelled (RACMO2.3p2, 2004–2016) and observed (a) 2 m temperature ($T_{\mathrm{2}\mathrm{m}}$, ${}^{\circ }$C), (b) 2 m specific humidity ($q_{\mathrm{2}\mathrm{m}}$, g kg${}^{-\mathrm{1}}$), (c) 10 m wind speed ($w_{\mathrm{10}\mathrm{m}}$, m s${}^{-\mathrm{1}}$) and (d) surface pressure ($P_{\mathrm{surf}}$, hPa) collected at 23 AWSs (green dots in Fig. ). For each variable, the linear regression including all records is displayed as a red dashed line. Statistics including number of records ($N$), regression slope (b0) and intercept (b1), determination coefficient ($R^{\mathrm{2}}$), bias and RMSE are listed for each variable.

[Figure omitted. See PDF]

Results: near-surface climate and SEB

We evaluate the modelled present-day near-surface climate of the GrIS in RACMO2.3p2 using data from 23 AWS sites (see Sect. 2.5). Then, we discuss in more detail the model performance at four AWSs along the K-transect and compare RACMO2.3p2 outputs to those of RACMO2.3p1.

Near-surface meteorology

Figure  compares daily mean values of 2 m temperature, 2 m specific humidity, 10 m wind speed, and air pressure collected at 23 AWS sites with RACMO2.3p2 output. The modelled 2 m temperature is in good agreement with observations ($R^{\mathrm{2}}=\mathrm{0.95}$) and with a RMSE of $\sim$ 2.4${}^{\circ }$C and a small cold bias of $\sim$ 0.1${}^{\circ }$C (Fig. a). As specific humidity is not directly measured at AWS sites, it is calculated from measured temperature, pressure and relative humidity following . The obtained 2 m specific humidity is accurately reproduced in the model ($R^{\mathrm{2}}=\mathrm{0.95}$) with a RMSE $\sim$ 0.35 g kg${}^{-\mathrm{1}}$ and a negative bias of 0.13 g kg${}^{-\mathrm{1}}$ (Fig. b). The same holds for daily records of 10 m wind speed ($R^{\mathrm{2}}=\mathrm{0.68}$; Fig. c), with the model exhibiting a small negative bias and RMSE of $\sim$ 2 m s${}^{-\mathrm{1}}$. Surface pressure is also well represented ($R^{\mathrm{2}}=\mathrm{0.99}$) with a small negative bias of 0.8 hPa and RMSE $<\mathrm{6}$ hPa (Fig. d). A systematic pressure bias at some stations results from the (uncorrected) elevation difference with respect to the model, which can be as large as 100 m. To provide some regional insight into the model performance, Supplement Table S1 and Figs. S1–S4 compare modelled meteorological data from RACMO2.3p2 with AWS measurements (green dots in Fig. ) clustered in four sectors of the GrIS, i.e. NW, NE, SW and SE. These sectors correspond to the four quadrants delimited by 40${}^{\circ }$ W longitude and 70${}^{\circ }$ N latitude. These regional scatter plots unambiguously show that RACMO2.3p2 performs equally well in each of these four sectors of the GrIS.

Difference between daily modelled RACMO2.3p1 (2004–2015) or RACMO2.3p2 (2004–2016) and observed meteorological data and SEB components collected at 23 PROMICE AWSs (green dots in Fig. ). Statistics include model bias (RACMO2.3pX – observations), RMSE of the bias as well as the determination coefficient of daily mean data. All fluxes are set to positive.

PROMICE	23 AWS	RACMO2.3p1				RACMO2.3p2
variable	unit	bias	RMSE	$R^{\mathrm{2}}$		bias	RMSE	$R^{\mathrm{2}}$
$T_{\mathrm{2}\mathrm{m}}$	${}^{\circ }$C	$-$0.8	2.9	0.93		$-$0.1	2.4	0.95
$q_{\mathrm{2}\mathrm{m}}$	g kg ${}^{-\mathrm{1}}$	$-$0.2	0.4	0.93		0.1	0.4	0.95
$w_{\mathrm{10}\mathrm{m}}$	m s${}^{-\mathrm{1}}$	0.3	2.2	0.65		$-$0.02	2.0	0.68
Psurf	hPa	$-$0.2	7.7	0.98		$-$0.8	5.8	0.99
SW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	6.5	26.7	0.96		3.8	27.1	0.95
SW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	5.8	27.8	0.91		6.8	32.1	0.88
LW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	$-$13.8	25.7	0.79		$-$7.1	21.2	0.83
LW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	$-$7.4	14.6	0.91		$-$4.4	12.1	0.92

Table  and Fig. S5 compare the agreement of RACMO2.3p2 and version 2.3p1 with in situ measurements. We find an overall improvement in the updated model version, showing a smaller bias and RMSE as well as an increased variance explained. Notably, the remaining negative bias in 2 m temperature (Fig. S5a) and the systematic dry bias (Fig. S5b) in RACMO2.3p1 have almost vanished in the updated model version (Figs. a and b).

Comparison between daily average modelled (RACMO2.3p2, 2004–2016) and observed (a) shortwave downward, (b) shortwave upward, (c) longwave downward and (d) longwave upward radiation (W m${}^{-\mathrm{2}}$) collected at 23 AWSs (green dots in Fig. ). For each variable, a regression including all records is displayed as a red dashed line. Statistics including number of records ($N$), the linear regression slope (b0) and intercept (b1), determination coefficient ($R^{\mathrm{2}}$), bias and RMSE are listed for each variable.

[Figure omitted. See PDF]

Radiative fluxes

Figure  shows scatter plots of modelled and measured daily mean radiative fluxes, i.e. short-/longwave down-/upward radiation. Radiative fluxes are also well reproduced by the model with $R^{\mathrm{2}}$ ranging from 0.83 for LW${}_{\mathrm{d}}$ to 0.95 for SW${}_{\mathrm{d}}$ (Fig. ), showing relatively small biases of $-$7.1 and 3.8 W m${}^{-\mathrm{2}}$, and RMSEs of 21.2 and 27.1 W m${}^{-\mathrm{2}}$. The negative biases in LW${}_{\mathrm{d}}$ and 2 m temperature partly lead to an LW${}_{\mathrm{u}}$ underestimation of 4.4 W m${}^{-\mathrm{2}}$ with a small RMSE of 12.1 W m${}^{-\mathrm{2}}$. In combination with a positive bias in SW${}_{\mathrm{d}}$ this suggests an underestimation of cloud cover in the ice sheet marginal regions, where most stations are located. The larger biases and RMSEs in SW${}_{\mathrm{u}}$ of 6.8 and 32.1 W m${}^{-\mathrm{2}}$ can be ascribed to overestimated surface albedo, especially during summer snowfall episodes, when a bright fresh snow cover is deposited over bare ice. In RACMO2, precipitation falls vertically, i.e. no horizontal transport is allowed, and is assumed to be instantly deposited at the surface. Consequently, the spatial distribution of summer snow patches may be locally inaccurate, resulting in large albedo discrepancies when compared to point albedo measurements. Note that these AWS radiation measurements are also prone to potentially large uncertainties due to preferred location on ice hills, sensor tilt, riming and snow/rain deposition on the instruments, leading to spurious albedo and SW${}_{\mathrm{u}}$ data , e.g. the upper-left dots in Fig. b. By clustering AWS measurements within four sectors of the GrIS (Figs. S6–S9 and Table S1), RACMO2.3p2 shows good and equivalent agreement in NW, NE, SW and SE Greenland.

Compared to the previous model version (Table ), changes in the cloud scheme have significantly improved the representation of LW${}_{\mathrm{d}}$ (Figs. c and S10c), showing a reduced negative bias and RMSE. These modifications have also somewhat decreased the positive bias in SW${}_{\mathrm{d}}$ (Fig. a) relative to RACMO2.3p1 (Fig. S10a). In addition, LW${}_{\mathrm{u}}$ is notably improved in RACMO2.3p2: the remaining negative bias in LW${}_{\mathrm{u}}$ has almost vanished (Figs. d and S10d). This can be partly explained by the much better resolved 2 m temperature in RACMO2.3p2.

Seasonal SEB cycle along the K-transect

The K-transect comprises four AWS sites located in different regions of the GrIS: S5 and S6 are installed in the lower and upper ablation zones, respectively, S9 is situated close to the equilibrium line and S10 is in the accumulation zone. Figure  shows monthly mean modelled (continuous lines, RACMO2.3p2) and observed (dashed lines) SEB components, i.e. net short-/longwave radiation (SW${}_{\mathrm{n}}$/LW${}_{\mathrm{n}}$), latent and sensible heat fluxes (LHF and SHF, respectively), surface albedo and melt measured at these four AWS sites for the period 2004–2015. Tables – list statistics calculated at each individual AWS and for the two model versions.

Modelled and observed mean SEB components and statistics of the differences (2004–2015) at station S5 in the lower ablation zone (490 m a.s.l.). Statistics include means of measurements collected at S5, model bias (RACMO2.3pX – observations), RMSE of the bias as well as the determination coefficient of monthly mean data. All fluxes are set to positive.

AWS	S5	Observed	RACMO2.3p1			RACMO2.3p2
variable	unit	mean	bias	RMSE	$R^{\mathrm{2}}$	bias	RMSE	$R^{\mathrm{2}}$
SW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	109.5	26.2	33.1	0.99	20.7	27.2	0.98
SW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	70.9	15.8	25.0	0.93	4.5	34.3	0.74
LW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	241.4	$-$17.0	18.5	0.97	$-$11.8	13.4	0.97
LW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	278.3	$-$13.2	15.5	0.98	$-$12.1	14.2	0.98
SHF	W m${}^{-\mathrm{2}}$	41.1	$-$13.1	22.2	0.45	$-$15.3	22.4	0.66
LHF	W m${}^{-\mathrm{2}}$	5.3	2.6	5.6	0.72	3.4	6.5	0.64
M	W m${}^{-\mathrm{2}}$	42.6	$-$6.8	18.0	0.96	$-$0.4	11.9	0.97
ALB	(–)	0.74	0.03	0.09	0.75	$-$0.004	0.14	0.72
$T_{\mathrm{2}\mathrm{m}}$	${}^{\circ }$C	$-$6.4	$-$2.3	2.6	0.99	$-$2.0	2.2	0.992

Modelled and observed mean SEB components and statistics of the differences (2004–2015) at station S6 in the upper ablation zone (1010 m a.s.l.). Statistics include means of measurements collected at S6, model bias (RACMO2.3pX – observations), RMSE of the bias as well as the determination coefficient of monthly mean data. All fluxes are set to positive.

AWS	S6	Observed	RACMO2.3p1			RACMO2.3p2
variable	unit	mean	bias	RMSE	$R^{\mathrm{2}}$	bias	RMSE	$R^{\mathrm{2}}$
SW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	131.6	9.7	12.7	0.997	6.0	9.1	0.997
SW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	95.8	$-$2.9	16.3	0.97	$-$3.8	16.3	0.97
LW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	222.3	$-$8.7	11.4	0.96	$-$2.7	6.5	0.97
LW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	263.6	$-$1.6	4.0	0.991	$-$0.4	3.1	0.992
SHF	W m${}^{-\mathrm{2}}$	20.8	9.8	11.4	0.67	7.0	8.7	0.70
LHF	W m${}^{-\mathrm{2}}$	1.6	$-$3.9	5.2	0.42	$-$2.4	3.3	0.64
M	W m${}^{-\mathrm{2}}$	18.7	10.6	22.0	0.96	8.3	18.0	0.97
ALB	(–)	0.81	$-$0.02	0.06	0.89	$-$0.02	0.06	0.89
$T_{\mathrm{2}\mathrm{m}}$	${}^{\circ }$C	$-$10.9	0.4	0.8	0.994	0.7	1.0	0.995

Modelled and observed mean SEB components and statistics of the differences (2009–2015) at station S9 close to the equilibrium line (1520 m a.s.l.). Statistics include means of measurements collected at S9, model bias (RACMO2.3pX – observations), RMSE of the bias as well as the determination coefficient of monthly mean data. All fluxes are set to positive.

AWS	S9	Observed	RACMO2.3p1			RACMO2.3p2
variable	unit	mean	bias	RMSE	$R^{\mathrm{2}}$	bias	RMSE	$R^{\mathrm{2}}$
SW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	141.2	2.2	6.6	0.998	$-$1.5	7.8	0.997
SW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	106.5	3.5	9.4	0.991	3.5	7.6	0.995
LW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	217.8	$-$10.1	14.1	0.93	$-$3.1	8.9	0.94
LW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	255.2	$-$1.9	4.9	0.99	0.5	3.5	0.991
SHF	W m${}^{-\mathrm{2}}$	15.8	7.0	9.2	0.68	5.2	7.3	0.74
LHF	W m${}^{-\mathrm{2}}$	0.8	$-$3.8	5.4	0.20	$-$2.8	4.0	0.42
M	W m${}^{-\mathrm{2}}$	12.0	$-$0.7	7.8	0.89	$-$2.4	7.0	0.96
ALB	(–)	0.82	0.02	0.05	0.79	0.03	0.06	0.83
$T_{\mathrm{2}\mathrm{m}}$	${}^{\circ }$C	$-$13.3	$-$0.04	0.7	0.994	0.5	0.8	0.996

Modelled and observed mean SEB components and statistics of the differences (2010–2015) at station S10 in the accumulation zone (1850 m a.s.l.). Statistics include means of measurements collected at S10, model bias (RACMO2.3pX – observations), RMSE of the bias as well as the determination coefficient of monthly mean data. All fluxes are set to positive.

AWS	S10	Observed	RACMO2.3p1			RACMO2.3p2
variable	unit	mean	bias	RMSE	$R^{\mathrm{2}}$	bias	RMSE	$R^{\mathrm{2}}$
SW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	141.5	1.7	7.0	0.998	$-$2.1	8.3	0.998
SW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	113.8	$-$2.7	12.0	0.991	$-$1.7	7.3	0.997
LW${}_{\mathrm{d}}$	W m${}^{-\mathrm{2}}$	220.4	$-$14.4	17.2	0.93	$-$6.5	10.7	0.94
LW${}_{\mathrm{u}}$	W m${}^{-\mathrm{2}}$	252.5	$-$1.0	4.8	0.99	2.1	4.1	0.991
SHF	W m${}^{-\mathrm{2}}$	11.9	7.6	10.8	0.57	4.9	8.3	0.62
LHF	W m${}^{-\mathrm{2}}$	$-$2.7	$-$3.5	5.6	0.22	$-$2.1	3.5	0.62
M	W m${}^{-\mathrm{2}}$	8.9	2.5	6.6	0.89	$-$0.2	3.1	0.92
ALB	(–)	0.86	$-$0.01	0.04	0.69	0.01	0.03	0.73
$T_{\mathrm{2}\mathrm{m}}$	${}^{\circ }$C	$-$14.6	0.5	1.0	0.991	1.1	1.3	0.995

Observed and modelled (RACMO2.3p2) monthly mean (a) turbulent and net shortwave/longwave fluxes (W m${}^{-\mathrm{2}}$) and (b) surface albedo and surface melt energy (W m${}^{-\mathrm{2}}$) at site S5 for 2004–2015. Similar results are shown at S6 for 2004–2015 (c, d), S9 for 2009–2015 (e, f) and S10 for 2010–2015 (g, h).

[Figure omitted. See PDF]

Low ablation zone

At station S5 (490 m a.s.l.), surface melt is well reproduced in RACMO2.3p2, with a small negative bias of 0.4 W m${}^{-\mathrm{2}}$ (Table ; Fig. b). However, this good agreement results from significant error compensation between overestimated SW${}_{\mathrm{n}}$ (bias of 16.2 W m${}^{-\mathrm{2}}$) and underestimated SHF in summer (15.3 W m${}^{-\mathrm{2}}$; Fig. a). The bias in SW${}_{\mathrm{n}}$ is mostly driven by overestimated SW${}_{\mathrm{d}}$ (20.7 W m${}^{-\mathrm{2}}$; Table ) and to a lesser extent by underestimated SW${}_{\mathrm{u}}$ (4.5 W m${}^{-\mathrm{2}}$), resulting from underestimated cloud cover and ice albedo (Fig. b), respectively. AWSs are often installed on snow-covered promontories, i.e. hummocks, that maintain higher albedo in summer ($\sim$ 0.55) than their surroundings where impurities collect. Mixed reflectance from bright ice cover ($\sim$ 0.55) and neighbouring darker tundra, exposed nunataks or meltwater ponds ($<$ 0.30), located within the same MODIS grid cell, likely explains this underestimation. Another explanation stems from the deterioration of MODIS sensors in time, resulting in underestimated surface albedo records for the MCD43A3v5 product .

LW${}_{\mathrm{n}}$ is well reproduced in the model due to similar negative biases in LW${}_{\mathrm{d}}$ and LW${}_{\mathrm{u}}$ ($\sim$ 12 W m${}^{-\mathrm{2}}$), again indicating underestimated cloud cover. The large negative bias in SHF is attributed to an inaccurate representation of surface roughness in the lowest sectors of the ablation zone. show that observed surface roughness for momentum has a high temporal variability at site S5, with a minimum of 0.1 mm in winter, when a smooth snow layer covers the rugged ice sheet topography, and a peak in summer (up to 50 mm), when melting snow exposes hummocky ice at the surface. In RACMO2, surface aerodynamic roughness is prescribed at 1 mm for snow-covered grid cells and at 5 mm for bare ice, hence significantly underestimating values over ice in summer and thus causing too-low SHF . This bias in SHF at S5 is also partly ascribable to too-cold conditions (2 ${}^{\circ }$C). Although not negligible, LHF contributes little to the energy budget and shows a positive bias of 3.4 W m${}^{-\mathrm{2}}$, notably in winter.

Upper ablation zone

Station S6 is located at 1010 m a.s.l. in the GrIS upper ablation zone. There, summer melt is overestimated by $\sim$ 8 W m${}^{-\mathrm{2}}$ owing to both too-high SW${}_{\mathrm{n}}$ and SHF (9.8 and 7 W m${}^{-\mathrm{2}}$; Fig. c and Table ). As for S5, the bias in SW${}_{\mathrm{n}}$ results from overestimated SW${}_{\mathrm{d}}$ (6 W m${}^{-\mathrm{2}}$) and underestimated SW${}_{\mathrm{u}}$ (3.8 W m${}^{-\mathrm{2}}$). At the AWS location, surface albedo progressively declines from 0.60 to $\sim$ 0.40 when bare ice is exposed in late summer, whereas RACMO2.3p2 simulates bare ice at the surface throughout summer, with an albedo of 0.40. As a result, modelled surface albedo is systematically underestimated in summer, especially in July (Fig. d). Likewise, a small negative bias in LW${}_{\mathrm{n}}$ (2.3 W m${}^{-\mathrm{2}}$) is obtained as LW${}_{\mathrm{d}}$ and LW${}_{\mathrm{u}}$ are both slightly underestimated (Table ). Here, 2 m temperature is on average 0.7 ${}^{\circ }$C too high, causing SHF to be overestimated (7 W m${}^{-\mathrm{2}}$).

Equilibrium line

Close to the equilibrium line, RACMO2.3p2 slightly underestimates summer melt (2.4 W m${}^{-\mathrm{2}}$; Fig. f and Table ). At station S9 (1520 m a.s.l.), a perennial snow cover maintains a minimum albedo of 0.65 in summer, i.e. when melt wets the snow. A small positive bias in modelled snow albedo (0.03) combined with a slightly underestimated SW${}_{\mathrm{d}}$ (1.5 W m${}^{-\mathrm{2}}$) leads to an overestimated SW${}_{\mathrm{u}}$ (3.5 W m${}^{-\mathrm{2}}$), hence underestimating SW${}_{\mathrm{n}}$ (5 W m${}^{-\mathrm{2}}$). Although LW${}_{\mathrm{d}}$ is underestimated by 3.1 W m${}^{-\mathrm{2}}$ and LW${}_{\mathrm{u}}$ is overestimated by 0.5 W m${}^{-\mathrm{2}}$, especially in winter, LW${}_{\mathrm{n}}$ agrees well with the measurements. The 2 m surface temperature shows a 0.5 ${}^{\circ }$C positive bias, in turn causing a slightly too-large SHF (5.2 W m${}^{-\mathrm{2}}$; Fig. e and Table ).

Accumulation zone

All SEB components are well reproduced at site S10 (1850 m a.s.l.). Compensation of minor errors between underestimated SW${}_{\mathrm{d}}$ and SW${}_{\mathrm{u}}$ ($\sim$ 2 W m${}^{-\mathrm{2}}$) provides good agreement with observed SW${}_{\mathrm{n}}$ (Fig. g and Table ). Modelled surface albedo also compares well with the measurements, with only a small positive bias (0.01; Fig. h). LW${}_{\mathrm{n}}$ is underestimated by $\sim$ 9 W m${}^{-\mathrm{2}}$; this is mainly driven by a too-low LW${}_{\mathrm{d}}$ and a too-large LW${}_{\mathrm{u}}$ (Table ). The turbulent fluxes are well captured, although a significant positive bias in SHF persists ($\sim$ 5 W m${}^{-\mathrm{2}}$), especially in winter when LW${}_{\mathrm{d}}$ is underestimated. As biases in SHF and LW${}_{\mathrm{d}}$ are almost equal, modelled melt matches well with observations despite a small negative bias ($\sim$ 0.2 W m${}^{-\mathrm{2}}$).

Model comparison along the K-transect

Tables – compare statistics of SEB components between RACMO2.3p2 and 2.3p1. Although differences are relatively small, the new model formulation shows general improvements. The increased cloud cover over the GrIS reduced the biases in SW${}_{\mathrm{d}}$ and LW${}_{\mathrm{d}}$. Improvements in the representation of turbulent fluxes is partly attributed to the new topography prescribed in RACMO2.3p2 and the better resolved SW${}_{\mathrm{d}}$/LW${}_{\mathrm{d}}$, although significant biases remain at all stations.

At site S5 located in the low ablation zone (Table ), smaller SW${}_{\mathrm{d}}$ and lower ice albedo significantly reduce the SW${}_{\mathrm{u}}$ bias in RACMO2.3p2, and enhanced LW${}_{\mathrm{d}}$ decreases the negative bias in LW${}_{\mathrm{u}}$. As a result, melt increases substantially, reducing the negative bias compared to version 2.3p1. Note that SW${}_{\mathrm{d}}$ remains overestimated in RACMO2.3p2. This is compensated by underestimated SHF, i.e. partly caused by underestimated LW${}_{\mathrm{d}}$, providing realistic surface melt. In the upper ablation zone, similar improvements are obtained at site S6 (Table ). At site S6, all SEB components show smaller biases except for SW${}_{\mathrm{u}}$, as underestimated surface albedo increases the negative SW${}_{\mathrm{u}}$ bias.

Above the equilibrium line, enhanced cloud cover also reduces the SW and LW biases at sites S9 and S10 (Tables  and ). However, surface albedo overestimation in RACMO2.3p2 causes a small increase in melt underestimation.

Results: regional SMB

In Sect. 3, we discussed the overall good ability of RACMO2.3p2 to reproduce the contemporary climate of the GrIS, which is essential for estimating realistic SMB patterns. Here we compare SMB from RACMO2.3p2 and RACMO2.3p1 over the GrIS. For further evaluation, we focus on three regions where there are large differences in SMB between the two versions.

Changes in SMB patterns

Figure a shows SMB from RACMO2.3p2 for the overlapping model period 1958–2015. Differences with the previous version 2.3p1 are shown in Fig.  b and the changes in individual SMB components are depicted in Fig. . Owing to the modifications in the cloud scheme, clouds are sustained at higher elevations, enhancing precipitation further inland, while it decreases in low-lying regions. Changes are especially large in south-eastern Greenland where the decrease locally exceeds 300 mm w.e. yr${}^{-\mathrm{1}}$. Precipitation in the interior increases by up to 50 mm w.e. yr${}^{-\mathrm{1}}$ (Fig. a). This pattern of change is clearly recognisable in the SMB difference (Fig. b). In addition, the shallower saltation layer in the revised drifting snow scheme is responsible for reduced sublimation ($\sim$ 50 mm w.e. yr${}^{-\mathrm{1}}$; Fig. b) that reinforces the overall increase in SMB (Fig. b). Although drifting snow erosion changes locally, patterns are heterogeneous and the changes remain small when integrated over the GrIS (Fig. c). This process has only a limited contribution to SMB ($\sim$ 1 Gt yr${}^{-\mathrm{1}}$) resulting from drifting snow being transported away from the ice sheet towards the ice-free tundra and ocean.

(a) SMB (mm w.e. yr${}^{-\mathrm{1}}$) averaged for the period 1958–2015. The combined Zachariae Isstrøm and Nioghalvfjerdsbrae (79${}^{\circ }$ N) glacier basins are delineated by the yellow line. Yellow dots locate the K-transect measurement sites in western Greenland and the single AWS operated in south-eastern Greenland. (b) SMB difference (mm w.e. yr${}^{-\mathrm{1}}$) between RACMO2.3p2 and RACMO2.3p1 for the period 1958–2015. Areas showing significant difference are stippled in panel (b): difference exceeds 1 standard deviation of the difference between the two model versions.

[Figure omitted. See PDF]

Difference in SMB components (mm w.e. yr${}^{-\mathrm{1}}$) between RACMO2.3p2 and RACMO2.3p1 averaged for the period 1958–2015. Areas showing significant difference are stippled: the difference exceeds 1 standard deviation of the difference between the two model versions.

[Figure omitted. See PDF]

(a) Modelled basin-integrated SMB in RACMO2.3p2 (blue dots) and RACMO2.3p1 (red dots) and ice discharge estimates (black dots, ) from the glacier basins of Zachariae Isstrøm and Nioghalvfjerdsbrae (79${}^{\circ }$ N) in north-eastern Greenland (yellow line in b and c) for the period 1975–2015. Dashed lines represent average SMB for 1975–2001. Mean SMB as modelled by (b) RACMO2.3p2 and (c) RACMO2.3p1 in north-eastern Greenland for the period 1958–2015.

[Figure omitted. See PDF]

In the percolation zone, the decrease in run-off (Fig. d) is governed by reduced surface melt (Fig. e), mostly resulting from the smaller grain size of refrozen snow and the lower soot concentration in snow that have increased surface albedo (not shown), further increasing SMB (Fig. b). In western and north-eastern Greenland, this decrease in run-off even exceeds that of melt by 50 to 100 mm w.e. yr${}^{-\mathrm{1}}$, a result of combined enhanced precipitation and reduced summer melt (delaying the disappearance of the seasonal snow cover) that increased the snow refreezing capacity (Fig. f). At higher elevations, the decrease in refreezing is exclusively driven by melt reduction (Fig. e and f), while at the extreme margins of the GrIS, the lower ice albedo used in RACMO2.3p2 (Fig. c) locally increases run-off (Fig. d), in turn decreasing SMB (Fig. b).

North-eastern Greenland

For north-eastern Greenland's two main glaciers, Zachariae Isstrøm and Nioghalvfjerdsbrae (79${}^{\circ }$ N glacier; yellow line in Fig. a), solid ice discharge (D) estimates are available for the period 1975–2015 . Assuming that this glacier catchment draining of $\sim$ 12 $\mathit{\%}$ of the GrIS area remained in approximate balance until $\sim$ 2000 , i.e. D $\approx$ SMB, measurements of D at the grounding line of these marine-terminating glaciers can be used to evaluate modelled SMB.

In these two catchments, model updates significantly improve the representation of SMB, which was substantially underestimated in the previous version. Figure a compares ice discharge (black dots) with modelled SMB (RACMO2.3p2 as blue dots and 2.3p1 in red) integrated over the two glacier basins for 1958–2015. In a balanced system, i.e. before discharge accelerated in 2001, SMB equals ice discharge. Averaged over 1975–2001, modelled SMB in RACMO2.3p2 (20.5 Gt yr${}^{-\mathrm{1}}$) is similar to the estimated glacial discharge of 21.2 Gt yr${}^{-\mathrm{1}}$, significantly improving upon version 2.3p1 (15.8 Gt yr${}^{-\mathrm{1}}$). The negative bias in RACMO2.3p2 (0.7 Gt yr${}^{-\mathrm{1}}$; dashed blue line) is reduced by almost a factor of 8 relative to the previous version (5.4 Gt yr${}^{-\mathrm{1}}$) and SMB now equals discharge within the uncertainty. However, it is important to note that, while good agreement is obtained between averaged SMB and D before 2001, suggesting a glacier catchment in approximate balance as in , this does not necessarily confirm that spatial and temporal variability of north-eastern Greenland SMB is accurately resolved by the model. Averaging over 2001–2015 showed that basin mass loss accelerated due to enhanced surface run-off, decreasing SMB by 4.2 Gt yr${}^{-\mathrm{1}}$, and increasing ice discharge (2.8 Gt yr${}^{-\mathrm{1}}$).

Figures b and c show mean SMB for 1958–2015 as modelled by RACMO2.3p2 and 2.3p1. In the percolation zone, the difference between the two model versions primarily results from the smaller refrozen snow grain size that reduces melt and run-off through increased surface albedo in RACMO2.3p2. To a smaller extent, reduced soot concentration delays the onset of melt in summer. In the ablation zone, snow cover persists longer before bare ice is exposed in late summer, in turn reducing run-off (Fig. d). Superimposed onto this, precipitation has increased over the whole glacier basin (Fig. a), allowing for enhanced refreezing in snow (Fig. f) and increasing SMB by 4.7 Gt yr${}^{-\mathrm{1}}$ in RACMO2.3p2 (Fig. b). Note the large inter-annual variability in modelled SMB showing maximum and minimum values of approximately 30 and 8.5 Gt yr${}^{-\mathrm{1}}$ in RACMO2.3p2 vs. 25 and 0 Gt yr${}^{-\mathrm{1}}$ in the previous version and stressing the importance of accurately modelling individual SMB components. In this dry region, underestimation of snowfall accumulation in RACMO2.3p1 initiated a pronounced feedback decreasing SMB: active drifting snow processes erode the shallow snow cover, exposing bare ice prematurely and moving the equilibrium line too far inland (Fig. b and c).

(a) Observed and simulated SMB (m w.e. yr${}^{-\mathrm{1}}$) along the K-transect in western Greenland (67${}^{\circ }$ N), averaged for the period 1991–2015. The observed SMBs (grey dots) at S4, S5, SHR, S6, S7, S8, S9 and S10 are based on annual stake measurements; S10 observations cover 1994–2015. The coloured bars represent the standard deviation ($\mathrm{1}\mathit{\sigma }$) around the 1991–2015 modelled and observed mean value. Modelled SMB at stake sites are displayed for RACMO2.3p2 (blue dots) and RACMO2.3p1 (red dots). Panel (b) shows time series of modelled (continuous lines) and observed (dashed lines) annual SMB at stakes S4, SHR, S7 and S8 for the period 1991–2016. Similar time series are shown for S5, S6, S9 and S10 in panel (c). At S10, modelled SMB is estimated as the difference between total precipitation and melt.

[Figure omitted. See PDF]

K-transect

The K-transect in south-western Greenland consists of eight stake sites where SMB is measured annually (yellow dots in Fig. a) . Figure a compares modelled (RACMO2.3p2 as blue dots and RACMO2.3p1 in red) with observed SMB (black dots) along the transect, averaged for the period 1991–2015. Using mean annual SMB at each station, the updated model shows a decreased RMSE from 606 mm w.e. in RACMO2.3p1 to 424 mm w.e. in version 2.3p2, and reduced bias from $-$133 to $-$54 mm w.e., and an increased $R^{\mathrm{2}}$ from 0.92 to 0.97. In the low ablation zone ($<$ 600 m a.s.l.), the lower ice albedo increases run-off in summer, locally reducing SMB. Decreased run-off in the upper ablation zone, i.e. between 600 and 1500 m a.s.l., increases SMB, improving the agreement at all sites except SHR. A negative bias in SMB remains at site S6 where ice albedo in summer (0.45 in July) is underestimated by up to 0.1 (Fig. d). Above the equilibrium line ($>$ 1500 m a.s.l.), in situ stake SMB measurements systematically underestimate climatic SMB, as they do not or only partly account for internal accumulation, i.e. refreezing in the firn. For comparison at S10, we therefore use the difference between modelled total precipitation and melt instead of SMB, decreasing the bias by 260 to $-$40 mm w.e. yr${}^{-\mathrm{1}}$ and the RMSE by 200 to 210 mm w.e. yr${}^{-\mathrm{1}}$. Measured and modelled SMB-to-elevation gradients are estimated using a linear regression: 3.21 mm w.e. m${}^{-\mathrm{1}}$ from the observations, 2.62 mm w.e. m${}^{-\mathrm{1}}$ in RACMO2.3p1, and 3.16 mm w.e. m${}^{-\mathrm{1}}$ in RACMO2.3p2, indicating a notable improvement in model performance along the K-transect.

Figures b and c show time series of measured (dashed lines) and modelled SMB (continuous lines; RACMO2.3p2) at each site along the K-transect for the period 1991–2016. The model realistically captures inter-annual variability in the SMB signal, although substantial biases remain at stations SHR and S6 (Table ).

Modelled and observed mean annual SMB (m w.e. yr${}^{-\mathrm{1}}$) and statistics of the differences at S4, S5, SHR, S6, S7, S8 and S9 over 1991–2015; measurements at S10 are compared to modelled total precipitation minus melt for the period 1994–2015. Spatial coordinates of each site are listed.

Stakes	Observed	RACMO2.3p1				RACMO2.3p2				Coordinates
SMB	mean	bias	RMSE	$R^{\mathrm{2}}$		bias	RMSE	$R^{\mathrm{2}}$		long.	lat.	elevation
										(${}^{\circ }$ W)	(${}^{\circ }$ N)	(m a.s.l.)
S4	$-$4.2	0.64	0.84	0.40		$-$0.05	0.51	0.47		$-$50.20	67.10	383
S5	$-$3.7	0.64	0.79	0.45		$-$0.08	0.46	0.50		$-$50.09	67.10	490
SHR	$-$3.1	$-$0.32	0.57	0.53		0.41	0.62	0.51		$-$49.94	67.10	710
S6	$-$1.7	$-$0.68	0.87	0.30		$-$0.56	0.78	0.29		$-$49.40	67.08	1010
S7	$-$1.5	$-$0.65	0.75	0.64		$-$0.15	0.37	0.68		$-$49.15	66.99	1110
S8	$-$0.8	$-$0.31	0.49	0.62		$-$0.03	0.28	0.76		$-$48.88	67.01	1260
S9	$-$0.2	$-$0.13	0.21	0.83		0.07	0.16	0.88		$-$48.25	67.05	1520
S10	0.3	$-$0.25	0.33	0.44		$-$0.04	0.21	0.45		$-$47.02	67.00	1850

South-eastern Greenland

South-eastern Greenland experiences topographically forced precipitation maxima in winter, followed by high melt rates in summer, allowing for the formation of perennial firn aquifers . In April 2014, an AWS was installed in the aquifer zone of the south-eastern GrIS (yellow dot in Fig. a). In August 2015, the AWS was relocated from 1563 m a.s.l (66.18${}^{\circ }$ N and 39.04${}^{\circ }$ W) to 1663 m a.s.l (66.36${}^{\circ }$ N and 39.31${}^{\circ }$ W). Figure  shows time series of snow albedo and cumulative snowmelt energy (expressed in mm w.e.) modelled by RACMO2.3p2 (blue lines) and RACMO2.3p1 (red lines), and calculated from the AWS data (grey lines) for the summer of 2014. The comparison is limited to 2014 because of a 3-month data gap in summer 2015.

Time series of (a) daily snow albedo, and (b) cumulative surface melt (mm w.e. day${}^{-\mathrm{1}}$) modelled by RACMO2.3p2 (blue lines), RACMO2.3p1 (red lines) and measured (grey lines) at the south-eastern AWS (66${}^{\circ }$ N; 33${}^{\circ }$ W; 1563 m a.s.l.) during summer 2014.

[Figure omitted. See PDF]

Comparison between (a) modelled, i.e. RACMO2.3p2 (blue) and RACMO2.3p1 (red) at 11 km, and observed SMB (m w.e. yr${}^{-\mathrm{1}}$) collected in the GrIS accumulation zone (white dots in Fig. ). Regressions for RACMO2.3p2 (blue) and version 2.3p1 (red) are displayed as dashed lines. Comparison between SMB measurements from the GrIS ablation zone (yellow dots in Fig. ) and (b) original RACMO2.3p2 data at 11 km, (c) downscaled product at 1 km. Orange stars correspond to measurements collected at station QAS_L at the southern tip of Greenland. Regression including all records is displayed as orange dashed line in panels (b) and (c). Main statistics including number of records ($N$), regression slope (b0) and intercept (b1), determination coefficient ($R^{\mathrm{2}}$), bias and RMSE are listed for each graph.

[Figure omitted. See PDF]

As melt wets the snow in summer, surface albedo gradually decreases from values typical for dry fresh snow (0.85) to wet old snow ($\sim$ 0.75) in late summer, before sharply increasing again when a new fresh snow cover is deposited (grey line in Fig. a). In the previous model version, surface albedo could drop to values as low as $\sim$ 0.66 in summer (JJA), e.g. days 152 to 243, underestimating albedo by 0.04 on average. The bias is reduced to 0.01 in RACMO2.3p2 as combined lower soot concentration and decreased grain size of refrozen snow increase the surface albedo. The remaining small negative bias is mostly ascribable to a too-rapid snow metamorphism from fresh to old snow that leads to a premature drop in surface albedo, e.g. days 140 to 160. Sporadic fresh snow deposition over older snow, characterised by sharp peaks in surface albedo during summer, is well timed by the model. Consequently, the cumulative melt obtained at the end of summer (702 mm w.e.; blue line in Fig. b) is reduced by $\sim$ 100 mm w.e. relative to RACMO2.3p1 (red line), a significant improvement when compared to the observations (639 mm w.e.; grey line).

Results: SMB of the contiguous ice sheet

Modelled SMB at 11 km

In Fig. , we evaluate modelled SMB in RACMO2.3p2 using 182 measurements collected in the GrIS accumulation zone (white dots in Fig. ) and 1073 stake observations from 213 sites located in the ablation zone (yellow dots in Fig. ). The increased precipitation in the GrIS interior reduces the negative bias in the 11 km product (blue dots in Fig. a) compared to the previous model version (red dots in Fig. a). For the full data set, a significant bias of $-$22 mm w.e. yr${}^{-\mathrm{1}}$ and RMSE of 72 mm w.e. yr${}^{-\mathrm{1}}$ remain in RACMO2.3p2. Sites experiencing the highest precipitation rates on the steep slopes of south-eastern Greenland ($>$ 0.5 m w.e. yr${}^{-\mathrm{1}}$) primarily contribute to this bias. If only values $<$ 0.5 m w.e. yr${}^{-\mathrm{1}}$ are considered (156 measurements), the bias and RMSE decrease from $-$26 and 52 mm w.e. yr${}^{-\mathrm{1}}$ in RACMO2.3p1 to only $-$7 and 49 mm w.e. yr${}^{-\mathrm{1}}$ in RACMO2.3p2. In the ablation zone (Fig. b), the updated model performs as well as the previous version, i.e. with a bias of 1.20 m w.e. yr${}^{-\mathrm{1}}$ and RMSE of 0.47 m w.e. yr${}^{-\mathrm{1}}$ , although SMB remains overestimated in the lower sectors, caused by inaccurately resolved steep slopes, low ice albedo and relatively large turbulent fluxes at the GrIS margins, which require further downscaling (see Sect. 5.2).

Integrated over the GrIS, modelled SMB has increased by 66 Gt yr${}^{-\mathrm{1}}$ (415 Gt yr${}^{-\mathrm{1}}$; +19 $\mathit{\%}$) compared to the previous version. This difference is dominated by a significant increase in SMB in the percolation zone of the GrIS, driven by reduced meltwater run-off (61 Gt yr${}^{-\mathrm{1}}$ or $-$22 $\mathit{\%}$) and reduced sublimation (10 Gt yr${}^{-\mathrm{1}}$ or $-$24 $\mathit{\%}$), while precipitation decreased by less than 1 $\mathit{\%}$ (5 Gt yr${}^{-\mathrm{1}}$); the latter can be explained by the smaller GrIS area ($\sim$ 10 000 km${}^{\mathrm{2}}$ or 0.6 $\mathit{\%}$) in the new ice mask. We deem these changes in the 11 km fields to be realistic. For the poorly resolved marginal areas, the SMB product requires further statistical downscaling to reproduce the high melt rates in these rugged regions at the ice sheet margins. At 11 km resolution, run-off is locally underestimated by up to 6 m w.e. yr${}^{-\mathrm{1}}$, e.g. station QAS$\mathit{\_}$L in southern Greenland (orange stars in Fig. b).

Downscaled SMB to 1 km

To solve these issues at the margins, we apply the downscaling technique described in , which includes elevation and ice albedo corrections. As a result, modelled run-off increases by 82 Gt yr${}^{-\mathrm{1}}$ ($\sim$ 37 $\mathit{\%}$) to 305 Gt yr${}^{-\mathrm{1}}$ for the period 1958–2015, compared to the 11 km product, and the SMB biases and RMSEs in the GrIS ablation zone are reduced by 480 and 460 mm w.e. yr${}^{-\mathrm{1}}$. The error at QAS$\mathit{\_}$L is reduced to $\sim$ 2 m w.e. yr${}^{-\mathrm{1}}$ (orange stars in Fig. c), i.e. biases and RMSEs of 2.21 and 2.35 m w.e. yr${}^{-\mathrm{1}}$. A major improvement upon is that no additional precipitation correction is required here as the remaining negative bias in the GrIS interior has been almost eliminated in RACMO2.3p2 (Fig. a). At 1 km resolution, precipitation contributes 693 Gt yr${}^{-\mathrm{1}}$ to GrIS SMB. Relative to the 11 km product, GrIS-integrated SMB at 1 km decreases by 59 Gt yr${}^{-\mathrm{1}}$ ($-$14 $\mathit{\%}$) to 356 Gt yr${}^{-\mathrm{1}}$, in line with our previous estimate of 338 Gt yr${}^{-\mathrm{1}}$ (+5 $\mathit{\%}$) . This confirms once more that an 11 km resolution is insufficient to resolve run-off patterns over narrow ablation zones and small outlet glaciers, and that further downscaling is essential for obtaining realistic GrIS SMB.

Remaining limitations and challenges

Model resolution

Extensive model evaluation confirms that RACMO2.3p2 realistically reproduces the contemporary climate and SMB of Greenland, although significant biases remain. However, while a 11 km grid is sufficient to resolve large-scale inland SMB patterns, it does not well resolve irregular, low-lying regions at the GrIS margins where run-off peaks. There, the remaining issue is to accurately resolve total run-off of meltwater from the narrow ablation zone and small outlet glaciers. This demonstrates the need for higher-resolution (statistically or dynamically) downscaled products, e.g. the 1 km product as presented here, for regional mass balance studies.

An alternative approach is to carry out a dedicated Greenland simulation at higher spatial resolution, e.g. 5.5 km . This increase in resolution does lead to better resolved SMB gradients over marginal glaciers, without exceeding the physics constraints of a hydrostatic model like RACMO2. Subsequently applying the statistical downscaling technique to this 5.5 km product would likely result in further improvements.

Turbulent fluxes

Another model limitation stems from the turbulent fluxes scheme. While LHF remains generally small and contributes little to the energy budget, accurate SHF is crucial for capturing extreme melt events along the GrIS margins , such as those that occurred in summer 2012 . However, SHF shows significant biases in RACMO2.3p2 in low-lying regions at the GrIS margins. Improving the representation of the GrIS surface roughness and surface elevation using higher spatial resolution could reduce these biases.

Surface albedo

Snowmelt rate is highly sensitive to soot concentration in snow . Although assumed to be constant in time and space in RACMO2, show a heterogeneous distribution of impurities (soot, dust, microbiological material) over the GrIS, with a gradual increase towards lower elevations due to (a) the proximity of dust sources in the tundra region and (b) downslope transport of previously deposited soot by meltwater run-off.

Over bare ice, the accumulation of cryoconite and the growth of algae play a major role in reducing surface albedo . Therefore, explicitly modelling impurity concentration on ice, as described in , could substantially improve melt estimates. Future climate projections should include such a bio-darkening feedback .

Conclusions

We present a detailed evaluation of the regional climate model RACMO2.3p2 (1958–2016) over the Greenland ice sheet (GrIS). The updated model generates more inland precipitation at the expense of marginal regions, reducing the dry bias in the GrIS interior. Impurity concentration in snow, i.e. soot, has been decreased by a factor of 2, minimising the melt rate overestimation in the GrIS percolation zone. We demonstrate that the model successfully reproduces the contemporary climate of the GrIS compared to daily meteorological records and radiative energy flux measurements from 23 AWS sites. Apart from the ultimate margins, the model also accurately captures the seasonal cycle of radiative and turbulent heat fluxes as well as surface albedo along the K-transect in south-western Greenland. Compared to SMB observations, RACMO2.3p2 generally improves on the previous version, especially in the extensive GrIS interior. SMB improvements are also found along the K-transect as well as in north-eastern and south-eastern Greenland. This model version will be used for future climate scenario projections at 11 km resolution. Nonetheless, since run-off from narrow glaciers in the GrIS margins remains poorly resolved at this resolution, it is necessary to further statistically downscale present-day and future SMB fields to higher spatial resolutions for use in regional mass balance studies.

RACMO2.3p2 data at 11 km (1958–2016) and a daily downscaled product at 1 km resolution are available from the authors without conditions.

The Supplement related to this article is available online at https://doi.org/10.5194/tc-12-811-2018-supplement.

BN, WJB, JMW and MRB conceived this study, decided on the new model settings and performed the analysis and synthesis of the data sets. BN performed the model simulations and led the writing of the manuscript. JTML, EM, PKM and LHU contributed to the development of the model. DA, SL, CJPPS and RSWW processed and provided observational data sets. All authors contributed to discussions on the writing of this manuscript.

The authors declare that they have no conflict of interest.

Acknowledgements

Brice Noël, Willem Jan van de Berg, J. Melchior van Wessem, Roderik S. W. van de Wal and Michiel R. van den Broeke acknowledge support from the Polar Programme of the Netherlands Organization for Scientific Research (NWO/ALW) and the Netherlands Earth System Science Centre (NESSC), as well as the European Centre for Medium-Range Weather Forecasts (ECMWF) for hosting simulations and providing computation time. Edited by: Xavier Fettweis Reviewed by: two anonymous referees