Regional climate model

Output of the Regional Atmospheric Climate Model (RACMO2.3) is used as input for the downscaling procedure14,16. RACMO2.3 combines the atmospheric dynamics from the High-Resolution Limited Area Model (HIRLAM) and the physics from the European Centre for Medium-range Weather Forecasts Integrated Forecast System (ECMWF-IFS)27. The polar version of RACMO2.3 is developed by the Institute for Marine and Atmospheric Research Utrecht University (IMAU), to simulate the evolution of SMB over ice sheets and surrounding smaller glacierized regions. Polar RACMO2.3 incorporates a multi-layer snow module to simulate firn compaction, meltwater retention and percolation, refreezing and runoff9. In RACMO2.3, the excess energy available at the surface, resulting from closure of the surface energy budget, is used to melt snow and ice. Liquid water from melt and rain percolates through the firn column, and is either held as irreducible water or refreezes, progressively reducing pore space from bottom to top layers until the entire firn column turns into ice (900 kg m−3) and no additional water can be stored. At this point, any additional water is assumed to run off. The model also includes a snow albedo scheme using prognostic snow grain size28; a drifting snow routine accounting for sublimation and snow erosion29. For the contemporary Arctic simulation, RACMO2.3 was run at 11 km and forced on a six-hourly basis by ERA-40 (ref. 30; 1958–1978) and ERA-Interim31 (1979–2015) re-analyses. The ice mask and topography at 11 km are based on a 5 km Greenland DEM32. For more information about RACMO2.3, recent updates and evaluation see refs 14, 33, 34.

Topography and ice mask

We used a downsampled version of the GIMP DEM17 to downscale the output of RACMO2.3 to a 1 km topography and ice mask (Supplementary Fig. 4a,b), obtained by averaging the original GIMP DEM at 90-m resolution. To distinguish between the ice sheet, including glaciers strongly connected to the ice sheet (corresponding to connectivity level CL2 in ref. 21), and GICs that are physically or dynamically separated from the ice sheet (respectively CL0 and CL1 in ref. 21), we used the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) ice classes projected on the 1 km GIMP ice mask (Supplementary Fig. 4b). In addition, floating glacier tongues were eliminated using a 1 km ice grounding line35. This results in a GICs area of ∼81,400 km2, ∼8% less than previous estimates22, owing to unresolved small glaciers in the original GIMP DEM at 90-m (ref. 17).

Bare ice albedo

To correct for the systematic bare ice albedo overestimation of RACMO2.3 in low-lying glaciated regions at 11 km, we used a 1 km version of the 500 m MODerate-resolution Imaging Spectroradiometer (MODIS) 16-day Albedo product (MCD43A3). Bare ice albedo is defined as the average of the 5% lowest surface albedos recorded for 2000–2015. At 1 km, bare ice albedo values range from 0.15 for dark ice surfaces at the ablation zone edges and local glacier tongues, to 0.55 under persistent snow cover in the accumulation zone16 (Supplementary Fig. 4c). In RACMO2.3, bare ice albedo is prescribed from a similar 11 km product (2001–2010) with a lower threshold of 0.30 (ref. 14). A value of 0.47 is assigned to glaciated pixels showing no valid MODIS estimate.

In-situ measurements

A total of 965 SMB measurements collected at 101 stake sites4 was used to evaluate the downscaled GICs SMB product at 1 km. These records were combined with another 1,073 observations from 230 sites on the GrIS (yellow dots in Fig. 1a) to adjust runoff and melt in the downscaling procedure. The ablation data set4 was compiled as part of the PROMICE36. For consistency, we only selected data with a temporal overlap with RACMO2.3 (1958–2015), and rejected sites with a >100 m height difference relative to the GIMP DEM at 1 km. To compare downscaled and observed SMB (Fig. 1b), we selected the grid cell with the smallest elevation bias among the closest pixel and its eight neighbours.

Remote sensing

Elevation changes for 2003–2009 and 2010–2014 were derived from ICESat and CryoSat-2 measurements, following the methods described in refs 6, 37. For ICESat, observations were grouped every 700 m along repeated ground tracks, whereas for CryoSat-2, neighbouring observations are collected within 1 km of each individual echo location. To these clusters of elevation observations, a model is fitted to estimate the local surface topography and elevation rate at the central point, where outliers are removed in an iterative procedure. For full details, see ref 37. After estimating the local topography and elevation rate for the ICESat and CryoSat-2 periods, local elevation anomalies at the echo locations can be estimated by adding back the constant elevation rate of the fitted model to the residuals. These anomalies are subsequently used to compute monthly volume anomalies for selected regions. We do so by parameterizing the elevation anomalies as a function of absolute elevation using a third-order polynomial. The resulting fit is then used to derive regional volume anomalies within 100 m elevation intervals, by multiplying the polynomial value at each interval’s midpoint with the total glaciated area within this elevation bin38. We perform the polynomial fit for nine regions individually6 and sum the results to obtain the GICs volume anomaly. Finally, volume anomalies are converted to mass anomalies by assuming a constant density profile, using the density of ice below the equilibrium line altitude (ELA), and a density of 600±250 kg m−3 above the ELA. Figure 2 shows the cumulative GICs mass change; using cumulative values suppresses noise in the ICESat and CryoSat-2 time series.

Downscaling procedure

The daily, 1 km SMB product is statistically downscaled from the output of RACMO2.3 at 11 km resolution (1958–2015) to the 1 km GIMP topography and ice mask (Supplementary Fig. 4a,b), using elevation dependence. Elevation correction is only applied to SMB components showing a significant correlation with height: runoff (RU), melt (ME) and sublimation (SU)16; total precipitation (PR), that is, rainfall (RA) and snowfall (SF), and snowdrift erosion (ER) are bi-linearly interpolated to the 1 km ice mask, without elevation correction. Daily SMB is then reconstructed as:

Elevation dependence

Daily regression parameters are calculated for the dependence of modelled SMB components on the 11 km RACMO2.3 topography. A local regression slope, b 11 km (mm w.e. per m, Supplementary Fig. 5), is estimated for each glacierized 11 km pixel using at least five adjacent ice-covered pixels. By applying the obtained b 11 km to the current pixel, the SMB component is approximated at mean sea level, a 11 km (mm w.e., Supplementary Fig. 5). To fully cover the GrIS and detached GICs, the latter regression parameters are extrapolated outwards by averaging b 11 km from at least three glaciated pixels. An estimate of a 1 km and b 1 km is then obtained by interpolating bi-linearly the 11 km regression parameters to the 1 km ice mask. Runoff, melt and sublimation (X v0.2 ) are calculated using a linear function of the high-resolution topography (h 1 km ) as:

The downscaled product based on elevation dependence only is hereafter referred to as version v0.2.

Runoff and melt adjustments

To correct for bare ice albedo overestimation in RACMO2.3 at 11 km, melt and runoff (v0.2) are adjusted to account for additional ice melt (ME add ) observed in low-lying regions compared to the downscaled product v0.2.

where ME add (mm w.e. per day) is the additional amount of ice melt calculated at 1 km (Supplementary Fig. 4c); Δα (dimensionless) is the difference between 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 direct 1 km is the modelled daily cumulated downward shortwave radiation bi-linearly interpolated to 1 km; L f is the latent heat of fusion (3.337 × 105 J kg−1). To account for the slope of the glacier surface, the dimensionless correction factor for a tilted plane ξ is applied to the direct component of modelled downward shortwave radiation. This correction is required as RACMO2.3 models radiation assuming a horizontal surface. This factor ranges from 0 for north sloping glaciers to values larger than 1 for south-oriented slopes.

A daily ratio Γ, between downscaled runoff and melt in v0.2, is applied to ME add to calculate the additional runoff (RU add ).

Assuming that the residual misfit between reconstructed and observed SMB (ΔSMB) for the different GICs ablation sites is ascribable to underestimated runoff in narrow ablation zones, RU add is then scaled by a factor f scale (0.682), obtained by computing a least-squares fit minimizing ΔSMB using both GICs and GrIS ablation measurements. The justification for including GrIS observations is that GICs SMB is not independent of the GrIS as the regression parameters were extrapolated from the ice sheet margins onto the GICs. The fact that f scale <1 indicates a slight overestimation of the melt adjustment calculated in equation (3), which could be due to the uncertainties in clouds39 and/or ice albedo underestimation at the ice caps margins.

The adjusted amount of runoff (RU v1.0 ) is obtained by adding RU add to the downscaled runoff (RU v0.2 ):

The corrected melt (ME v1.0 ) is obtained in a similar manner:

Refreezing (RF v1.0 ) is estimated as a residual between adjusted melt, runoff and rainfall:

The downscaled SMB data set resulting from the combined elevation correction and runoff adjustment is referred to as version v1.0.

Precipitation correction

To eliminate the systematic negative SMB bias of RACMO2.3 in Greenland’s accumulation zones16 (−37.5 mm w.e. yr−1), daily precipitation totals from v0.2 were adjusted in regions presenting a positive annual cumulative SMB in v1.0:

where PR v1.0 is the daily adjusted total precipitation v1.0, PR v0.2 is the daily bi-linearly interpolated total precipitation v0.2, PR v0.2 a is the annual cumulative bi-linearly interpolated total precipitation v0.2 and σ SMB is the accumulation zone SMB bias in the downscaled product v1.0.

Product uncertainty

Assuming that the remaining discrepancies in Fig. 1b consist of a systematic bias due to model uncertainty, combined with additional random scatter, attributed to observations, a product uncertainty can be obtained by integrating the average accumulation and ablation zone biases and adding them as if they were independent. To that end, SMB measurements (Fig. 1b) were binned in 500 mm w.e. bins for which the mean bias was estimated, that is, modelled minus measured SMB. The average uncertainty that results is 247 mm w.e. for the ablation zones (8 bins) and 135 mm w.e. for the accumulation zones (6 bins). The product SMB uncertainty of 15.7 Gt yr−1 (∼40%) is obtained by summing the mean ablation bias integrated over the GICs ablation zones (45,600 km2) and the mean accumulation bias over the whole GICs area (81,400 km2), to account for a potential precipitation bias in the ablation zones.

Average biases for the GrIS accumulation (22 mm w.e.) and ablation zones (170 mm w.e.) have been calculated in a similar fashion and provided a product SMB uncertainty of 52.5 Gt yr−1 (∼20%). These calculations were repeated for 250 mm w.e. SMB bins, which yielded similar results.

Breakpoint analysis

We applied the segmented regression method described in ref. 40 to retrieve breakpoints, that is, structural changes, in the GICs-integrated refreezing time series. To detect a break point, the technique fits a regression function, consisting of segments with different slopes, to the studied data set. The algorithm calculates multiple continuous regressions before and after the break points and determines an optimized solution. A confidence interval can be estimated for the resulting break point by using the 95% Wald-based statistics. In this study, we used the segmented regression method to identify a tipping point in the GICs refreezing capacity at 1997±5 years.

Data availability

The daily, 1 km SMB data set v1.0 presented in this study is freely available from the authors without conditions. The background in Supplementary Fig. 1 stems from Landsat satellite imagery freely available at https://landsat.usgs.gov/.