A daily, 1 km resolution data set of downscaled Greenland ice sheet surface mass balance (1958–2015)

: This study presents a data set of daily, 1 km resolution Greenland ice sheet (GrIS) surface mass balance (SMB) covering the period 1958–2015. Applying corrections for elevation, bare ice albedo and accumulation bias, the high-resolution product is statistically downscaled from the native daily output of the polar regional climate model RACMO2.3 at 11 km. The data set includes all individual SMB components projected to a down-sampled version of the Greenland Ice Mapping Project (GIMP) digital elevation model and ice mask. The 1 km mask better resolves narrow ablation zones, valley glaciers, fjords and disconnected ice caps. Relative to the 11 km product, the more detailed representation of isolated glaciated areas leads to increased precipitation over the southeastern GrIS. In addition, the downscaled product shows a significant increase in runoff owing to better resolved low-lying marginal glaciated regions. The combined corrections for elevation and bare ice albedo markedly improve model agreement with a newly compiled data set of ablation measurements. Abstract. This study presents a data set of daily, 1 km resolution Greenland ice sheet (GrIS) surface mass balance (SMB) covering the period 1958–2015. Applying corrections for elevation, bare ice albedo and accumulation bias, the high-resolution product is statistically downscaled from the native daily output of the polar regional climate model RACMO2.3 at 11 km. The data set includes all individual SMB components projected to a down-sampled version of the Greenland Ice Mapping Project (GIMP) digital elevation model and ice mask. The 1 km mask better resolves narrow ablation zones, valley glaciers, fjords and disconnected ice caps. Relative to the 11 km product, the more detailed representation of isolated glaciated areas leads to increased precipitation over the southeastern GrIS. In in to better resolved low-lying marginal glaciated combined for and ice albedo improve model with a newly compiled data set of ablation


Introduction
During the last 2 decades, the Greenland ice sheet (GrIS) has experienced significant mass loss as a result of increased meltwater runoff and sustained high solid ice discharge from marine-terminating outlet glaciers ( Van den Broeke et al., 2009;Rignot et al., 2008Rignot et al., , 2011Sasgen et al., 2012;Shepherd et al., 2012;Enderlin et al., 2014). To fill spatial and temporal gaps in the scarce in situ observations, regional climate models (RCMs) are often used to produce maps of the GrIS surface mass balance (SMB; Burgess et al., 2010;Ettema et al., 2010a, b;Fettweis, 2007;Fettweis et al., 2005Fettweis et al., , 2011Noël et al., 2015;Lucas-Picher et al., 2012). RCMs explicitly calculate the individual SMB components , i.e. precipitation, runoff and sublimation, over the entire ice sheet (Fig. 1) at high spatial and temporal resolution and over extended periods. However, the current spatial resolution of RCMs, typically 5-20 km, remains too coarse to accurately resolve glaciated areas in topographically complex regions such as small isolated ice caps and marginal outlet glaciers flowing into narrow fjords. In these regions, the relatively coarse elevation and land ice masks used in RCMs might result in runoff underestimation (Franco et al., 2012;Noël et al., 2015), hampering realistic regional SMB estimates. Performing higher resolution simulations to address these issues would require a substantial computational effort and is thus restricted to case studies of small regions and relatively short time periods. This figure also depicts the location of 213 ablation measuring sites (yellow dots) and 182 accumulation sites (white dots) used for downscaled SMB evaluation as well as the four GrIS marginal regions (blue boxes), discussed in Sect. 5. Letters refer to the different transects shown in Fig. 9.
As an alternative, statistical downscaling can be applied to RCM output. Previously, this method has been applied to the GrIS using global reanalysis and climate data (Hanna et al., 2005(Hanna et al., , 2011. Machguth et al. (2013) downscaled nearsurface temperature and precipitation from three different RCMs (11-25 km spatial resolution) to force a glacier mass balance model on a 250 m grid derived from the Greenland Ice Mapping Project (GIMP) digital elevation model (DEM) , accurately resolving local glaciers and ice caps of Greenland. Vertical gradients of climate parameters were iteratively calibrated to enable the mass balance model to generate a realistic melt distribution for the period 1980-2010, but the very high resolution restricted the analysis to a few regions. Franco et al. (2012) statistically downscaled GrIS SMB by interpolating each component of the Modèle Atmosphérique Régional (MAR) from the original 25 km grid to a 15 km resolution. This method used local daily vertical gradients, except for precipitation, to correct for elevation differences between MAR and a down-sampled version of the 5 km DEM from Bamber et al. (2001). The elevation correction significantly reduced SMB biases. However, a resolution of 15 km remains insufficient to resolve the rugged topography at the ice sheet margins; to address this issue, near-kilometre resolution is necessary.
Here, we present a new data set of daily, 1 km resolution GrIS SMB components (precipitation, melt, runoff, refreezing, sublimation and snowdrift erosion) covering the period 1958-2015. The SMB product is statistically downscaled from data of the Regional Atmospheric Climate Model version 2.3 (RACMO2.3) at 11 km ( Fig. 1), using an elevationdependent technique based on the elevation and ice mask from the GIMP DEM , down-sampled to 1 km. The following section briefly describes RACMO2.3, the GIMP DEM, observational data sets and MODIS bare ice albedo product used to evaluate and correct the downscaled data set. The downscaling algorithm is explained in Sect. 3. Downscaled SMB is evaluated using ablation and accumulation measurements in Sect. 4. Section 5 discusses the downscaling results for four different regions and for the entire ice sheet. The added value, limitations and uncertainties of the downscaling method are argued in Sect. 6, followed by conclusions in Sect. 7.
2 Model and data 2.1 The regional climate model RACMO2 A detailed description of the Regional Atmospheric Climate Model (RACMO2) is presented in Van Meijgaard et al. (2008). RACMO2 incorporates the atmospheric dynamics and physics modules from the High Resolution Limited Area Model (HIRLAM) and the European Centre for Medium-range Weather Forecasts Integrated Forecast System (ECMWF-IFS, Unden et al., 2002). The polar version of RACMO2 is developed by the Institute for Marine and Atmospheric Research (IMAU), Utrecht University, and is especially adapted for use over ice sheets and other glaciated regions. Polar RACMO2 is interactively coupled to a multi-layer snow module, accounting for firn densification, meltwater percolation, refreezing and runoff (Ettema et al., 2010a); it also incorporates an albedo scheme with prognostic snow grain size (Kuipers Munneke et al., 2011) and a drifting snow module, simulating snow erosion and the drifting snow contribution to sublimation . Recently, RACMO2.1 has been updated to RACMO2.3 as discussed in Van Wessem et al. (2014) and Noël et al. (2015). Model evaluation against SMB measurements, collected in the accumulation and ablation zones of the GrIS, showed generally improved agreement (Noël et al., 2015). The native 11 km climate run is forced at the lateral boundaries by ERA-40 (1958-1978, Uppala et al., 2005 and ERA-Interim (1979, Stark et al., 2007Dee et al., 2011) reanalyses and uses the 5 km DEM and ice mask from Bamber et al. (2001).

GIMP DEM
To downscale RACMO2.3 output, we use the ice mask and topography from the GIMP DEM, described in Howat et al. (2014), and currently considered to be one of the most complete ice masks for Greenland (Rastner et al., 2012). A 1 km ice mask and DEM are obtained by averaging the original 90 m GIMP grid cells in each 1 km pixel covering Greenland. A 1 km resolution is deemed an acceptable trade-off between improved resolution, i.e. a 121-fold improvement compared to the 11 km grid, and manageable data handling given the daily time resolution, time span  and the number of SMB components. As an example, Fig. 2 shows the topography and ice mask from RACMO2.3 at 11 km in central east Greenland (left in Fig. 2; blue box 1 in Fig. 1) and the GIMP DEM at 1 km (right). The latter better resolves small scale landforms such as narrow fjords and calving glacier tongues. Integrated over the contiguous GrIS, the ice-covered area of 1.69 × 10 6 km 2 for the 1 km grid represents a 0.5 % decrease relative to the 11 km mask. For our SMB calculations, we only consider grounded ice, i.e. we discarded floating ice pixels using a 1 km version of the 90 m grounded ice mask used in Enderlin and Howat (2013).

Ablation and accumulation measurements
To evaluate the daily downscaled SMB product, we use 1155 SMB measurements collected in the GrIS ablation (1073) and accumulation (182) zones. The ablation data set  was compiled as part of the Programme for Monitoring of the Greenland Ice Sheet (PROMICE)  and includes stake and automatic weather station measurements retrieved from 213 sites (yellow dots in Fig. 1). Accumulation observations were derived from 182 sites including snow pits and firn cores (Bales et al., 2001(Bales et al., , 2009) as well as airborne radar measure-ments (Overly et al., 2016) (white dots in Fig. 1). We exclusively selected data that temporally overlap with RACMO2.3 simulations (1958-2015; 205 sites discarded). We rejected observations from sites with a > 100 m height bias relative to the representative elevation of the 1 km GIMP topography (one site discarded).
To compare modelled and downscaled SMB with observations, different selection approaches were applied in the ablation and accumulation zones, as described in Noël et al. (2015). In the accumulation zone, we select the closest grid cell on the 11 and 1 km grids to represent modelled and downscaled SMB, respectively. In the ablation zone, an altitude correction is applied by selecting the grid cell with the smallest elevation bias among the closest pixel and its eight adjacent neighbours.

MODIS bare ice albedo
A 1 km version of the 500 m MODerate-resolution Imaging Spectroradiometer (MODIS) 16-day Albedo product (MCD43A3) is used to retrieve estimates of bare ice albedo in the GrIS ablation zone. Bare ice albedo is estimated as the average of the 5 % lowest surface albedo measurements for the period 2000-2015. A similar ice albedo product is used in RACMO2.3 based on MODIS observations between 2001 and 2010 (Noël et al., 2015). In RACMO2.3, bare ice albedo ranges from 0.3, i.e. dark bare ice exposed in the low ablation zone, to 0.55 under persistent snow cover in the GrIS accumulation zone. Bare ice albedo of glaciated pixels with no valid MODIS estimate is set to 0.47.

Methods
The daily, 1 km SMB product consists of statistically downscaled output from a previously conducted RACMO2.3 simulation at 11 km, covering the period 1958-2015. RACMO2.3 settings and lateral forcing are described in Noël et al. (2015). The downscaling algorithm corrects the interpolated SMB components using their local regression to elevation. Figure 3 shows the spatial correlation of individual SMB components with elevation on the 11 km RACMO2.3 grid. The spatial correlation is calculated for each grid box using eight adjacent ice-covered pixels.
The elevation correction is exclusively applied to the SMB components, which show a significant and spatially homogeneous correlation with elevation, i.e. melt, runoff and sublimation (Fig. 3). These SMB components decrease with decreasing air temperature, represented by a negative correlation with elevation ( Fig. 3b, d, e). Although precipitation negatively correlates with elevation over most of the ice sheet, the correlation remains small and highly heterogeneous at the margins (Fig. 3a). Snowdrift erosion exhibits a noisy correlation pattern. Therefore, daily precipitation and snowdrift erosion are bilinearly interpolated to the 1 km ice mask without elevation corrections. Refreezing exhibits a marked bimodal correlation pattern (not shown), gradually increasing with height in the ablation zone, where pore space is more abundant, and decreasing towards the ice sheet interior due to limited meltwater supply. For this reason, and in order to have a consistent liquid water balance, daily refreezing is calculated as a residual: where RF is the residual refreezing, RA is rainfall, ME is surface melt, and RU is meltwater runoff.
Daily SMB values are obtained by summing the individually downscaled components: where P tot is total precipitation (liquid and solid), RU is meltwater runoff, SU is total sublimation (from surface and drifting snow) and ER is drifting snow erosion.

Elevation-dependent downscaling
The downscaling algorithm interpolates daily SMB components to the 1 km topography and ice mask in three successive steps (Fig. 4). First, the local dependence on elevation is calculated on the original RACMO2.3 11 km grid. Regression parameters are computed on a daily basis and are, therefore, only valid for that specific day. A local regression slope, b 11 km (mmWE m −1 , Fig. 4), is calculated for each ice-covered RACMO2.3 grid point using the maximum number of points available; i.e. we use a total of six to nine ice-covered grid cells, the current one and the minimum five to maximum eight adjacent pixels. This minimum number is chosen after testing the downscaling sensitivity to the number of regression cells used, as discussed in Sect. 3.2. An approximation of the SMB components at mean sea level, a 11 km (mmWE, Fig. 4), is then obtained using b 11 km and the current pixel. The regression is applied to the current grid cell to prevent local estimates of a 11 km to significantly differ from the original RACMO2.3 value. Local regression parameters for melt and runoff are only computed for pixels experiencing ablation. Moreover, erroneous positive regression slopes, i.e. increasing melt rates with altitude, are discarded until the following stage.
Next, valid estimates of b 11 km and a 11 km are extrapolated iteratively on the 11 km grid to fully cover the more extensive 1 km ice mask. To that end, daily regression parameters are extrapolated outwards of the 11 km ice mask by averaging b 11 km from at least three ice-covered pixels from the eight cells surrounding the current one.
Finally, the extrapolated fields of b 11 km and a 11 km are bilinearly interpolated to the 1 km ice mask, providing estimates of b 1 km and a 1 km . The downscaled SMB components Current grid cell Adjacent grid cells a 1 1 k m = yb 1 1 k m * x Figure 4. Elevation-dependent downscaling procedure: b 11 km and a 11 km are respectively the daily local estimates of the SMB components regression to elevation and the SMB components value at mean sea level obtained on the RACMO2.3 grid at 11 km. The red line corresponds to the regression (b 11 km ) calculated using the current grid cell (blue dot) and the adjacent ones (red dots). The dashed green line applies the regression slope to the current grid cell to estimate a 11 km .
(X v0.2 ), i.e. runoff, melt and sublimation, are then computed as a linear function of the high-resolution topography as The downscaled data set that is based on the above elevation-dependent technique is hereafter referred to as version v0.2. Figure 5 shows the difference between 11 km and downscaled GrIS integrated daily runoff in summer 2011. Each line represents a different minimum number of grid cells, ranging from three to nine, used to estimate the local regression of runoff with elevation (b 11 km ; Fig. 4). The results are moderately sensitive to the number of regression points used except for the nine cells setting, which systematically underestimates runoff at the beginning and the end of the melt season as it discards all low-lying glaciated pixels at the edge of the GrIS that experience early melt and the largest values of runoff. The standard deviation between the different settings (∼ 0.2 Gt day −1 ) is smaller than the difference between 11 and 1 km runoff (∼ 0.6 Gt day −1 ). The more regression points used, the smoother the runoff to elevation gradient field becomes, lowering the downscaled runoff and bringing it closer to the 11 km model output. Conversely, a small number of regression points can lead to spuriously large local gradients. To prevent the downscaling algorithm from substantially converging to, or diverging away from, 11 km RACMO2.3 output, we adopted a setting of a minimum of 2011 -1 Figure 5. Summer 2011 time series of daily, ice sheet-integrated runoff difference (Gt day −1 ) between the downscaled product at 1 km, using a minimum threshold of three-nine regression points (legend), and the RACMO2.3 model at 11 km. six regression points, which is closest to the average value of the different experiments (±0.1 Gt day −1 ).

Melt and runoff adjustments
RACMO2.3 uses a prescribed bare ice albedo field, typically ranging from 0.30 in the low ablation zone to 0.55 under persistent snow cover. It is based on the lowest 5 % of MODIS surface albedo values averaged for the period 2001-2010 (Noël et al., 2015). A comparison with a similar 1 km MODIS product averaged for 2000-2015, ranging from 0.15 to 0.55, shows a systematic overestimation of ice albedo at 11 km, especially for low-lying marginal glacier tongues (Fig. 13i). This causes melt energy to be underestimated during the melt season. To correct for this, downscaled melt and runoff are adjusted by estimating the missing amount of ice melt (ME add ) resulting from underestimated absorption of downward shortwave radiation (SW d ). In addition, as RACMO2.3 calculates radiative fluxes on a horizontal plane, the direct fraction of SW d is corrected for the slope and orientation of each 1 km glaciated grid cell, as described in Weiser et al. (2016). For simplicity, we assume SW d to be equally partitioned between diffuse and direct radiation, and that the sun is exactly in the south at noon. This assumption is purely pragmatic; on the basis of data availability, it could be further refined in future versions of the downscaling procedure. Figures 6b and 7 show that ablation underestimation in v0.2 is restricted to the low ablation zone (SMB < −4 mWE), where bare ice is exposed for long episodes in summer. Therefore, the following corrections are only applied to the ablation zone on days of melting bare ice when both surface runoff and melt are non-zero in the downscaled product v0.2: where ME add (mmWE day −1 ) is the additional amount of ice melt calculated at 1 km; α (-) is the difference be-  tween the averaged bare ice albedo retrieved from the set of regression cells used to downscale runoff at 11 km and the MODIS albedo product at 1 km; SW d, 1 km is the modelled daily cumulated downward shortwave radiation bilinearly interpolated to 1 km; L f is the latent heat of fusion (3.337 × 10 5 J kg −1 ); and ξ (-) is the correction factor for a tilted plane (Fig. 8), applied to the direct component of downward shortwave radiation: where ζ * is the solar angle of incidence for a tilted plane, ζ is the solar zenith angle, a is the azimuth of the tilted plane, σ is the local surface slope, is the orientation, φ is the latitude, δ is the solar declination and H is the hour angle set to 0 at noon (Fig. 8). All angles are expressed in radians. The bare ice albedo bias correction aims at minimizing the misfit between downscaled SMB v0.2 and in situ measurements (Fig. 6b) by estimating the missing runoff in the low ablation zone. Additional runoff RU add is calculated by applying a daily specific fraction (-) to ME add , estimating the melt contribution to surface runoff. is defined as the ratio between daily downscaled runoff and melt in v0.2 estimated using elevation dependence only: Figure 8. Scheme of a tilted plane as described in the GIMP DEM at 1 km. SW d is the downward shortwave radiation, ζ * is the solar angle of incidence for a tilted plane, ζ is the solar zenith angle, a is the azimuth of the tilted plane, σ is the local surface slope and is the orientation. n and z are respectively the vector normal to the tilted plane and the local vertical axis.
Assuming that the residual misfit between reconstructed and observed SMB ( SMB, Fig. 6b) for the different ablation sites can be ascribed to underestimated runoff in the low ablation zone of the GrIS, RU add is then scaled by a factor f scale (-), obtained by computing a least-squares fit, minimizing the difference between SMB and RU add using all ablation measurements: The least-squares fit yields a value of f scale = 1.176 for the GrIS. This means that RU add , i.e. accounting for elevation and bare ice albedo corrections, has yet to be increased by ∼ 18 % to optimize the agreement between downscaled and in situ SMB (Fig. 6c). The fact that f scale > 1 strongly suggests that additional processes might play a role in enhancing surface ablation, e.g. underestimation of modelled sensible heat flux from warm air advection along the GrIS periphery (Noël et al., 2015;Fausto et al., 2016) and uncertainties in cloud representation (Van Tricht et al., 2016). However, as the statistical downscaling approach is not designed to correct for these physical processes, we adopted the empirical approach presented above. The adjusted amount of runoff (RU v1.0 ) is obtained by adding the missing runoff to the downscaled runoff (RU v0.2 ): The corrected melt (ME v1.0 ) is obtained in a similar fashion, and refreezing (RF v1.0 ) is estimated as a residual between adjusted melt, runoff and rainfall: ME v1.0 = ME v0.2 + ME add , The downscaled SMB data set resulting from the combined elevation correction and runoff adjustment is referred to as version v1.0 in the following sections.
4 Evaluation of daily downscaled SMB Figure 6 evaluates the original RACMO2.3 SMB at 11 km (a), the 1 km raw downscaled SMB version v0.2 (b) and the 1 km corrected downscaled SMB version v1.0 (c) (mWE yr −1 ) with 1073 observations from 213 ablation sites (yellow dots in Fig. 1). The observational period was matched with the modelled and downscaled SMB using the exact number of days. Each blue star corresponds to the cumulative SMB for a duration ranging from 10 days to a full hydrological year. The downscaled SMB v0.2 agrees better with observations compared to the RACMO2.3 output at 11 km ( Fig. 6a, b): we find a significant decrease of the root mean square error (RMSE) (190 mmWE or −16 %) and a smaller bias (100 mmWE or −21 %). The deviation from unity of the regression slope decreases from 0.28 to 0.21 (−25 %), and the variance explained increases from 47 to 61 %. When applying the bare ice albedo and local orientation corrections, we find further significant improvements relative to version v0.2 (Fig. 6c) Fig. 12a), situated in an extremely narrow ablation zone (∼ 10 km) at the southwestern tip of Greenland. Here, modelled ablation gradients at 11 km are strongly underestimated in RACMO2.3 and are only marginally better resolved at 1 km. At this site, the additional corrections are especially important to obtain agreement with observations. Figure 9 compares annual mean observed and downscaled SMB (v1.0) along eight different SMB transects. There is good agreement for most transects, except for Helheim Glacier (66.41 • N, −38.34 • W). The downscaled product fails at reproducing the quasi-constant ablation rate (∼ −1 mWe) characterizing the Helheim transect. The reason for this low SMB gradient is not clear at present; it may be due to uncertainties in individual observation covering relatively short periods, i.e. 1 or 2 months, which are only limited to the melt season (July-August). Another possible explanation is that Helheim Glacier experiences large and variable winter accumulation at low elevations, potentially caused by drifting snow transport, limiting summer ablation. In addition, Nioghalvjerds-fjorden and Storstrømmen transects (Fig. 9a, b) also show significant remaining biases between in situ and downscaled SMB at elevations lower than 200 m. We hypothesize that these SMB measurements are located on floating glacier tongues with melt ponds, resulting in very low satellite albedo, while stake measurements are -1 -1 Figure 9. Annual mean observed (red dots) and downscaled (blue dots, v1.0) SMB for eight selected transects in the GrIS ablation zone (mWE yr −1 ). Name and locations of these transects (Fig. 1) are listed at the bottom of each graph. Graphs also list the number of sites used for each transect, linear SMB-to-elevation regression retrieved from observations and downscaled (v1.0) data in mmWE yr −1 m −1 , the RMSE and the mean bias.
performed between ponds on brighter surfaces. As a result, the bare ice albedo correction could be overestimated. In the accumulation zone, a small improvement is also found in v0.2 compared to RACMO2.3 (Fig. 7), but accumulation remains underestimated. The SMB bias and RMSE are reduced by 0.7 mmWE (−2 %) and 1.8 mmWE (−3 %), whereas the regression slope and variance explained remain unchanged. In the accumulation zone, SMB is mostly driven by precipitation which is bilinearly interpolated to 1 km without elevation correction. In addition, changes in sublimation are small due to the relatively homogeneous topography of the ice sheet interior, limiting SMB changes through downscaling. Despite significant improvements in the cloud scheme of RACMO2.3 (Noël et al., 2015), clouds become saturated and start to produce precipitation at elevations that are too low, resulting in overestimated precipitation at the margins, e.g. southeastern Greenland, while the ice sheet interior experiences too dry conditions. This precipitation bias is currently being investigated, and we aim to resolve it in the upcoming version RACMO2.4. To overcome the systematic negative SMB bias of RACMO2.3 in the GrIS accumulation zone (−37.5 mmWE yr −1 , Fig. 7), the daily total precipitation v0.2 is adjusted to correct for underestimation in the ice sheet accumulation zone (SMB > 0 mmWE yr −1 ): where PR v1.0 is the daily adjusted total precipitation v1.0, PR v0.2 is the daily bilinearly interpolated total precipitation v0.2, PR a v0.2 is the annual cumulative bilinearly interpolated total precipitation v0.2 and σ SMB is the accumulation zone SMB bias in the downscaled product v1.0.
The final SMB v1.0 product is reconstructed as 5 High-resolution SMB patterns: case studies Table 1 lists annual mean modelled and downscaled SMB components (Gt yr −1 ) integrated over four different regions (blue boxes in Fig. 1) as well as over the entire GrIS. These regions were selected for their specific climates, rough topography and narrow glaciated features, which were not well resolved at 11 km. Figures 10, 11, 12 and 13 show the ice sheet mask for the selected regions at 11 km (red cells) and 1 km (orange cells) as well as peripheral glaciers and ice caps at 1 km (blue cells), the elevation bias between the 1 and 11 km DEMs and the bare ice albedo bias between the 1 km MODIS product and RACMO2.3 at 11 km; the latter figures moreover show the main SMB components at both resolutions for the two downscaled products (v0.2 and v1.0). In the following sections, we discuss the impact of downscaling on regional SMB. Here, SMB components are exclusively integrated over the contiguous GrIS; the SMB of detached ice caps will be discussed in a forthcoming paper.

Central east Greenland
Central east Greenland (blue box 1 in Fig. 1) is characterized by a large body of interconnected valley glaciers, mostly terminating in narrow glacial fjords. Figure 10a, e and i underline the inability of the 11 km mask to properly represent many glaciated areas, local topography or bare ice albedo. In the 1 km mask, the ice-covered area increases (∼ 2 %), while the elevation bias can locally exceed 500 m over glacial valleys and small-scale promontories (Table 1, Fig. 10e); the av-erage elevation bias is 80 m. These differences affect SMB in two ways. First, precipitation increases by 2.6 Gt yr −1 (12 %) in v0.2 (Table 1, Fig. 10b, f), exclusively caused by the expansion of glaciated area (no elevation correction is applied). Another 1.6 Gt yr −1 (6 %) of precipitation is added in v1.0 (Fig. 10j) to compensate for the systematic negative SMB bias in the GrIS accumulation zone, as discussed in Sect. 4. For both downscaling versions, changes in runoff mirror the elevation change between the two resolutions (Fig. 10e), -1 -1 -1 Figure 11. Centre west: same as Fig. 10 but for central west Greenland (blue box 2 in Fig. 1).
highlighting the high sensitivity of runoff to elevation. In version v0.2, integrated runoff increases by 7.7 Gt yr −1 (48 %) (Fig. 10c, g). Furthermore, Fig. 10i reveals a systematic overestimation of bare ice albedo at 11 km. Correcting for this further increases runoff over the glaciers tongues (Fig. 10k), accounting for ∼ 13 Gt yr −1 (55 %) of additional runoff with respect to v0.2 (Table 1). Negligible changes in sublimation and drifting snow are found (Table 1). As a consequence, integrated SMB on the 1 km mask decreases by 5.3 Gt yr −1 in version v0.2 (Fig. 10d, h) and by 16.6 Gt yr −1 in version v1.0 (Fig. 10l). This analysis for central east Greenland demonstrates the importance of accurately reproducing small-scale topography and ice albedo to realistically capture local SMB variations.

Central west Greenland
The 11 km resolution DEM provides a reasonable representation of the wide, gently sloping western ablation zone of the GrIS, where most glaciers are land-terminating. The northern part of the selected area includes several marine-terminating glaciers which are better represented at 1 km (Fig. 11d, h). Owing to negligible difference in glaciated area, precipitation remains almost unchanged for the two resolutions and versions (∼ 15 Gt yr −1 ). In both downscaled versions, enhanced runoff is mostly obtained over narrow, low-lying glaciers tongues and detached ice caps (Fig. 11c, g, k) where most of the elevation and ice albedo biases are found (Fig. 11e, i). On the ice sheet, the elevation correction increased runoff by about 1 Gt yr −1 (5 %) (Fig. 11g), while an additional ∼ 2 Gt yr −1 (10 %) (Fig. 11k) can be ascribed to the ice albedo correction (Table 1).

South Greenland
Southeast Greenland (blue box 3 in Fig. 1) is a rugged region (Fig. 12e), characterized by multiple topographically forced precipitation maxima (Fig. 12b, f) and narrow marginal ablation zones (Fig. 12c, g, k). Similar to central east Greenland, the larger glaciated area (+6.5 %, Fig. 12a) at 1 km enhances integrated precipitation by ∼ 6 Gt yr −1 (+7 %) in v0.2 and 8.4 Gt yr −1 (+9 %) in v1.0. Increased runoff (2.2 Gt yr −1 or 5 % in version v0.2) at the southern margins can be ascribed to additional melt production over the better resolved narrow ablation zones (Fig. 12g, k) combined with a moderate mean elevation difference (∼ 17 m) between both resolutions. In v0.2, the ice mask expansion explains most of the integrated SMB changes, leading to an overall mass gain of 3.3 Gt yr −1 . Figure 6b reveals considerable ablation underestimation in south Greenland, expressed as a systematic SMB bias of 2-4 mWe relative to measurements collected at PROMICE station QAS_L (red dots in Fig. 6a). The main reason for this underestimation is that SMB at this location is characterized by a rare combination of high snowfall and strong summer melt.
The remaining ablation underestimation in v0.2 can be partly ascribed to an overestimated bare ice albedo (0.47) prescribed in RACMO2.3 (Noël et al., 2015); observed albedo at QAS_L frequently falls to 0.2 during the melt season . As a result, the additional bare ice albedo correction significantly improves runoff at station QAS_L (Fig. 6c). Integrated over region 3, runoff increases by another ∼ 13 Gt yr −1 (29 %) relative to v0.2 (Fig. 12k). The increased marginal mass loss leads to the expansion of the southern ablation zone towards higher elevations (Fig. 12k, l), in line with local observations (Fig. 6c).

North Greenland
In north Greenland (blue box 4 in Fig. 1), the climate is dry, and most glaciers are marine-terminating. The ice sheet surface is relatively smooth and homogeneous. The wide ablation zone is reasonably well captured at 11 km, leading to a modest deviation in elevation (∼ 43 m) (Fig. 13e). However, the ice-covered area decreases by ∼ 11 % between both resolutions as the 11 km grid contained erroneous floating glacier tongues (Fig. 13a). The ice area reduction at 1 km affects precipitation (−0.8 Gt yr −1 or −12 %) (Fig. 13b, f) and runoff (−3.1 Gt yr −1 or −35 %) (Fig. 13c, g), resulting in a small SMB increase (2.3 Gt yr −1 ) in version v0.2 Figure 13. North: same as Fig. 10 but for north Greenland (blue box 4 in Fig. 1). The green line in (a) shows the grounded ice mask at 1 km. The yellow dot in (a) locates the Petermann glacier site settled on a floating ice tongue.

Greenland ice sheet
Although similar in area, the 1 km ice sheet mask better resolves peripheral glaciers at the GrIS margins than RACMO2.3 at 11 km. GrIS-integrated precipitation increases by 16.6 Gt yr −1 (+2 %) in v0.2, most of which can be ascribed to ice area expansion in the east (2.6 Gt yr −1 ) and south of Greenland (5.8 Gt yr −1 ), where precipitation is large. An additional 56.2 Gt yr −1 (+8 %) is obtained in v1.0 when correcting for the accumulation zone SMB bias. The smooth topography of the ice sheet interior results in a small elevation difference of 4 m between both resolutions. Significant elevation biases are mostly restricted to peripheral glaciers and narrow ablation zones at the GrIS margins. As a result, runoff increases by 13.6 Gt yr −1 (+5 %) in version v0.2. Accounting for the bare ice albedo bias in RACMO2.3 further increases runoff by 69.3 Gt yr −1 in version v1.0, leading to a much improved agreement with ablation measurements. Of our selected areas, central east and south Greenland contribute 25 and 18 % to the total runoff increase in the downscaled product v1.0 owing to the many low-lying glaciers tongues that can only be resolved at 1 km. Due to their smoother topography, north and centre west Greenland contribute much less to the runoff change (∼ 3 and 1 %, respectively). Integrated over the contiguous ice sheet, SMB is not significantly affected by the elevation dependence for which enhanced precipitation (16.6 Gt yr −1 ) yearly balances the moderate increase in runoff (13.6 Gt yr −1 ). In con-trast, the bare ice albedo and precipitation corrections substantially increase marginal runoff (82.9 Gt yr −1 ) and accumulation (72.8 Gt yr −1 ), resulting in a decrease of SMB of −11.1 Gt yr −1 (−3 %) relative to the 11 km product.

Added value, limitations and uncertainties
The downscaled SMB v1.0 is the first data set to provide daily SMB estimates for all outlet glaciers of the GrIS at a 1 km resolution for 58 years . Relative to the original RACMO2.3 output, this data set improves local SMB values (Fig. 9) and produces more realistic SMB patterns over rugged glaciated areas along the GrIS margins (Figs. 10,11,12,13). Figures 6 and 7 show that SMB v1.0 is an overall improvement on the original RACMO2.3. To further investigate this, Fig. 14 shows the annual mean SMB RMSE (model vs. observations) of the 11 km SMB field in RACMO2.3 (red), the downscaled product v0.2 (green) and v1.0 (blue) as a function of observed SMB, binned in 0.5 mWE intervals. In the ablation zone (SMB < 0), the SMB RMSE is reduced by 29-65 % in v1.0 relative to the 11 km product, owing to the elevation correction in v0.2 (9-23 %) and the additional albedo correction (20-42 %). In the accumulation zone, the elevation dependence (9 %) and the precipitation adjustment (19 %) also contribute to reduce the SMB RMSE by 28 % in v1.0. The largest RMSE reduction occurs in the lower GrIS ablation zone, where improvements in topography and bare ice albedo in v1.0 are greatest. Although significantly improved, the downscaled SMB v1.0 is likely to be locally underestimated for four reasons: (a) the bare ice albedo correction is evenly applied to both snow-covered and bare ice regions experiencing surface melt and runoff, as no relevant proxy, reflecting day-to-day snow coverage, could be derived from RACMO2.3. However, this issue should have a limited effect on the magnitude of downscaled melt and runoff since the albedo correction is most efficient in summer, when the snow cover of low-lying glaciers has likely melted. (b) The MODIS ice albedo product at 1 km becomes less accurate at high latitudes, likely suffering from bare soil contamination resulting from mixed reflectance signals recorded in both the tundra-and ice-covered regions. Note that floating glacier tongues also show surface albedo that is too low, e.g. Petermann glacier (yellow dot in Fig. 13a), resulting from mixed signals from adjacent dark melt pond and brighter dry ice. The resulting albedo underestimation over low-lying floating tongues below 200 m leads to overestimated ablation (∼ 0.2 mWE yr −1 ; Fig. 9a, b).  Figure 14. Annual mean modelled SMB RMSE (model vs. observations) of the 11 km SMB field in RACMO2.3 (red dots), the downscaled SMB data set v0.2 (green dots) and v1.0 (blue dots) as a function of observed SMB (395 observations). Modelled SMB is grouped in 0.5 mWE yr −1 bins except for the first bin, which ranges from −6.00 to −3.75 mWE yr −1 . Numbers indicate the number of observations used in each bin. initialize the downscaling procedure, i.e. RCM version used, the resulting modelled SMB components, bare ice albedo records, ablation measurements, topography and ice mask. The downscaled SMB v1.0 presents an estimated uncertainty of ∼ 12 Gt yr −1 in the GrIS ablation zone, which was estimated by integrating the SMB bias in v1.0 (60 mmWE, Fig. 6c) over the ablation zone of the contiguous ice sheet (∼ 202 000 km 2 ).
We anticipate that the new, 1 km Greenland SMB product is especially useful for studies that address the mass balance of Greenland outlet glaciers that are too steep and/or narrow to be properly resolved at the typical horizontal resolution of regional climate models (∼ 5-20 km). Future downscaled products can have even higher resolution (100 m) and will be based on further improved RCM output fields of precipitation and melt.

Conclusions
The relatively coarse spatial resolution currently used in RCMs remains insufficient to properly resolve small-scale variations in elevation and ice cover at the ice sheet margins, significantly affecting the calculation of melt and runoff. In the present study, we statistically downscale individual SMB components from RACMO2.3 at 11 km to a 1 km ice mask and topography derived from the GIMP DEM, using a daily specific elevation dependence. Moreover, runoff and melt are corrected for biases in bare ice albedo in RACMO2.3. Precipitation and snowdrift erosion are bilinearly interpolated without applying an elevation correction. Total precipitation is also adjusted to compensate for the dry accumulation bias of RACMO2.3 in the ice sheet interior. Downscaled daily SMB is then retrieved for the period 1958-2015 by summing daily downscaled precipitation, runoff, sublimation and drifting snow erosion. An evaluation of the downscaled SMB product against observations, collected both in the ablation and accumulation zones of the GrIS, shows improved agreement. In the ablation zone, the variance explained by the downscaled product v1.0 increased by 31 % relative to the original RACMO2.3 11 km output, mainly through better resolved narrow outlet glaciers at the GrIS margins.
Integrated over the GrIS, precipitation increased by 16.6 Gt yr −1 due to the larger glaciated area in south and east Greenland at 1 km; an additional correction of 56.2 Gt yr −1 must account for the accumulation bias in the ice sheet interior in RACMO2.3. Likewise, a 13.6 Gt yr −1 increase in runoff is attributed to elevation corrections on the 1 km topography, and another 69.3 Gt yr −1 extra runoff can be ascribed to underestimated bare ice albedo over narrow outlet glaciers at the GrIS margins. A small area in central east Greenland alone, characterized by multiple narrow glacier tongues poorly resolved at 11 km, accounts for ∼ 25 % of the total additional runoff.

Data availability
The daily, 1 km SMB data set v1.0 presented in this study is freely available from the authors without conditions.