Interactive comment on “ On the recent contribution of the Greenland ice sheet to sea level change ”

This is an interesting and useful up-to-date analysis of recent changes in the Greenland Ice Sheet mass balance, considering independent estimates of surface mass balance, solid ice discharge and total mass change from GRACE, that is set in a longer-term context of the last 58 years. The paper is generally well-written and -presented, and the underlying analysis seems proficient, although a few aspects (detailed below) need clarification. I’d be happy to recommend full publication in TC once the authors have addressed the following points.


Introduction
Three methods are commonly used to assess the mass balance (MB) of the Greenland ice sheet (GrIS): gravimetry, radar/laser altimetry and the mass budget method (MBM).Each method has distinct advantages and disadvantages.The main advantage of gravity observations is that they provide, once corrected and deconvolved, a relatively direct measure of glacial mass change; as a result, several years after its launch, the dual-satellite Gravity Recovery and Climate Experiment (GRACE) mission was able to confirm that the GrIS is losing mass (Velicogna and Wahr, 2006).The main drawbacks of this technique are the relative brevity of the time series (since late 2002), the large footprint of ∼ 200 km and the fact that corrections must be made for mass movements in atmosphere, ocean, soil and solid earth (e.g.glacial isostatic adjustment, or GIA) in order to isolate ice mass changes.
Like GRACE, satellite altimetry provides full coverage of the GrIS.The time series are longer, starting in the mid-1990s, with aerial photography extending some of the records back to the mid-1980s and earlier (Kjaer et al., 2012;Kjeldsen et al., 2015).Difficulties in interpreting radar altimeter data arise from the variable penetration depth of the radar signal in firn (Thomas et al., 2008) and, especially for the earlier instruments, signal loss along the steep coastal margins.The radar altimeter on board CryoSat-2, launched in 2010, partly mitigates these issues and shows GrIS elevation changes in unprecedented detail (Helm et al., 2014).ICE-Sat's laser altimeter measured the surface elevation change Published by Copernicus Publications on behalf of the European Geosciences Union.
accurately but was sensitive to cloud cover and had a relatively large ground track separation.The various altimeter records must be intercalibrated and spatially interpolated between ground tracks before a continuous time series is obtained.Both radar and laser altimeters measure ice sheet volume changes, which must be converted to mass changes using a model that accounts for vertical bedrock motion and variability in the depth and mass of the firn layer, which introduces additional uncertainties (Kuipers Munneke et al., 2015).
The MBM estimates the difference between individual mass sources (mainly snowfall) and sinks (mainly meltwater runoff and solid ice discharge).Because it resolves the individual components of the MB, this method has the advantage that it identifies the physical processes that are responsible for the mass change.However, because the mass change represents the relatively small difference between three large source and sink terms, it is very sensitive to uncertainties in any of these.This is especially true for surface mass fluxes such as snowfall and meltwater runoff; because these cannot be measured from space, they must be interpolated from scarce in situ measurements and/or simulated using dedicated regional climate models, which introduces potentially large uncertainties (Vernon et al., 2013).Shepherd et al. (2012) reconciled results of the three methods for the GrIS to obtain an average GrIS 1992-2011 mass loss of 142 ± 49 Gt yr −1 .The fifth assessment report of Working Group I of the Intergovernmental Panel on Climate Change (IPCC), while weighing the results of the various studies somewhat differently, arrived at a similar conclusion and shows that in 2012 the GrIS had become the largest single contributor to sea level rise (Vaughan et al., 2013).A compilation of GrIS MB studies covering sub-periods of the satellite era confirms that the mass loss of the GrIS is accelerating (Hanna et al., 2013b).
Combining two or more methods may reduce uncertainties and may provide additional insights in the physical processes causing the mass loss.Sasgen et al. (2012) used GRACE, the MBM and altimetry to assess mass changes in seven GrIS regions; for the GrIS as a whole, between August 2002 and August 2010 they found a mass trend of −228 ± 22 Gt yr −1 , with an acceleration of −15 ± 7 Gt yr −2 .Rignot et al. (2011) used a combination of GRACE and the MBM to show that GrIS mass loss between 1992 and 2010 accelerated by 17-22 Gt yr −2 .When combined with surface mass balance (SMB) fields, satellite altimetry can be used to spatially separate mass losses arising from ice dynamical (i.e.driven by ice flow) and surface (i.e.driven by the atmosphere) processes (Csatho et al., 2014;Kjeldsen et al., 2015).Using anomalies relative to a reference period, Van den Broeke et al. (2009) showed that 1996-2008 GrIS mass loss was approximately equally partitioned between increased surface meltwater runoff and ice discharge.Enderlin et al. (2014) showed that the relative contribution of ice discharge to total GrIS mass loss decreased from 58 % before 2005 to 32 % between 2009 and 2012, indicating an increasingly important role for surface processes.
In this paper we combine GRACE data and the MBM, using new SMB and D data that allow updating the time series to 2015, to identify the causes and temporal evolution of recent (1958-2015) GrIS mass loss and its contribution to sea level rise.In Sect. 2 we discuss methods, in Sect. 3 we discuss results and in Sect. 4 we identify outstanding problems and avenues for future research.

Definitions
The MB of an ice sheet, usually expressed in Gt yr −1 , represents the change of its mass in time dM/dt.Neglecting basal melting of grounded ice, which typically does not exceed several mm per year, and assuming the grounding line position to remain unchanged, the MB of the grounded ice sheet is governed by the difference between SMB and solid ice discharge across the grounding line (D): SMB represents the sum of mass fluxes towards and away from the ice sheet surface: where P tot is total precipitative flux (sum of snowfall, SN, and rainfall, RA), SU tot is total sublimation (from the surface and from drifting snow particles), ER ds is erosion of snow by divergence of the drifting snow transport and RU is meltwater runoff. 1 The accumulation/ablation zones of an ice sheet are defined as the areas where SMB > 0 and SMB < 0, respectively.These two zones are separated by the equilibrium line, where SMB = 0.The amount of RU from the ice sheet is determined by the liquid water balance (LWB), the sum of all sources (mainly rainfall and melt) and sinks (mainly refreezing and capillary retention) of liquid water in the column of firn and/or ice: where CO is condensation of water vapour at the ice sheet surface, ME is surface meltwater production, RT is retention of liquid water in the snow/firn by capillary forces and RF is refreezing of liquid water at or below the surface.Equation (3) in first order includes processes associated with the formation of perennial firn aquifers (Forster et al., 2013;Kuipers Munneke et al., 2014) but neglects the delay in runoff by storage of meltwater in semipermanent supra-, suband englacial lakes and channels, which is potentially significant (Rennermalm et al., 2013).In summary, to quantify the MB of the GrIS, the MBM relies on the quantification of mass sinks and sources in Eqs.(1), ( 2) and (3).

Data sources
To calculate GrIS MB using the MBM we combine SMB components calculated with the regional atmospheric climate model RACMO2.3 with annual estimates of D, updated from Enderlin et al. (2014).The latter data represent discharge summed over all marine-terminating glaciers wider than 1 km, and cover the 16-year period 2000-2015.Between 1996 and 2000 we assume a linear increase in D of in total 38 Gt yr −1 following Rignot et al. (2008).In the absence of data, no changes in D are assumed between 1958 and 1996.A seasonal cycle in ice-sheet-wide D is not considered, although it is well known that marine-terminating outlet glaciers can show pronounced (sub-)seasonal velocity oscillations (Moon et al., 2014).Because seasonal acceleration leads to thinning, the effect on D is smaller than based on velocity changes alone.Moreover, by using the median annual velocity, seasonal variability is accounted for to a large extent (Enderlin et al., 2014).Land-terminating glaciers also exhibit a seasonal cycle in their velocity (Van de Wal et al., 2008;Joughin et al., 2008;Bartholomew et al., 2011;Sole et al., 2013), but this does not influence D.
SMB components are derived from a run with RACMO2.3 for the period January 1958-December 2015, using 40 vertical layers and a horizontal resolution of ∼ 11 × 11 km 2 (Noël et al., 2015).Figure 1 shows the Greenland domain of RACMO2.3, which consists of 312 (latitude) × 306 (longitude) grid cells and includes Iceland, the Svalbard archipelago and the Canadian Arctic.The model is forced at the lateral boundaries by the 40-year European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis (ERA-40) for the period January 1958-December 1979 and the ECMWF Interim Reanalysis (ERA-Interim) afterwards.In previous versions of RACMO2, the impact of using inhomogeneous forcing before and after 1979 (ERA-40 vs. ERA-Int) was found to be small (Ettema et al., 2009).According to Fettweis et al. (2013a), the regional climate model MAR shows precipitation to be 5 % greater when forced by ERA-40 compared to ERA-Interim for the period 1980-1999, probably a result of biases in the ERA-40 humidity scheme that were corrected in ERA-Interim.This uncertainty falls within the error bar used here.
The polar version of RACMO2.3 has been especially developed to simulate SMB of glaciated regions and is an 3 Greenland domain with ice sheet mask (white), land mask (green) and elevation contours (dashed lines, every 500 m a.s.l.).Boundary relaxation zone indicated by dots every other grid point.update of RACMO2.1 (Ettema et al., 2009;Van Angelen et al., 2014).It is interactively coupled to a multilayer onedimensional (column) firn model that simulates (sub-)surface processes like vertical heat transport, grain growth, firn densification, meltwater percolation, retention and refreezing.RACMO2.3uses a prognostic calculation of snow grain size from which broadband snow albedo is derived (Kuipers Munneke et al., 2011).Ice albedo is not explicitly calculated, as it is influenced by poorly known processes such as dust accumulation and biological activity; to account for its considerable spatial heterogeneity (Bøggild et al., 2010), ice albedo is prescribed from the Moderate Resolution Imaging Spectroradiometer (MODIS) on board the Terra and Aqua satellites (Stroeve et al., 2005).
RACMO2.3 includes a routine for drifting snow sublimation, which removes on average 25 Gt yr −1 from the ice sheet (Lenaerts et al., 2012), with little interannual variability; when compared to scarce observations, this scheme was found to accurately predict the occurrence of drifting snow, but overestimate drifting snow transport, and thus likely also overestimate drifting snow sublimation, although no direct observations are available to verify this.For a more detailed description of recent changes in RACMO2.3model physics and how they impact the modelled SMB of the GrIS, the reader is referred to Noël et al. (2015) and references therein.
Other RACMO2.3 modelled SMB components and atmospheric parameters have been extensively evaluated against in situ observations both in the accumulation and ablation zone of the GrIS (Ettema et al., 2009;Noël et al., 2015).From these comparisons, typical uncertainties of 9 and 15 % were found for ice-sheet-wide integrated accumulation and ablation, respectively.These are combined into an uncertainty for ice-sheet-wide SMB assuming accumulation and ablation to be independent.We note that this assumption is debatable, as ablation and accumulation tend to be connected via surface albedo, especially in summer (Fyke et al., 2014;Noël et al., 2015).The scarcity of accumulation and ablation measurements does not allow for a further regionalization of the uncertainty, but obviously uncertainties can be significantly larger for smaller areas and sub-climatological time periods.For the trend in cumulative SMB-D, an uncertainty is derived following Van den Broeke et al. (2009).
We compare ice sheet MB obtained from the MBM, i.e.MB = SMB-D, with monthly GRACE gravity field solutions.Mass variations for the GrIS were derived from the CSR RL05 data release (April 2002-September 2015), following the method described in Wouters et al. (2008).In brief, regional mass anomalies are adjusted in a model consisting of eight pre-defined GrIS basins and the resulting gravity disturbance is computed.The modelled mass anomalies are then adjusted until convergence with the actual GRACE observations is reached.All standard corrections are applied to the GRACE data, including a correction for GIA, based on the model of A et al. (2013).The uncertain- ties in monthly GRACE values and mass trend are estimated at 40 Gt and 20 Gt yr −1 , respectively.For the monthly values, the uncertainty is computed by conservatively assuming that all observed signals with a periodicity smaller than 6 months are due to noise (Wahr et al., 2006).For the trends, the quoted uncertainty takes into account methodological differences in the processing of the level-1 GRACE data, the uncertainty in the GIA correction, the formal error of the least-square fit and aliasing of ocean tides (Wouters et al., 2013).To assess the methodological uncertainty in the GRACE time series, results were compared to mass anomalies from the Jet Propulsion Laboratory (JPL) mascons using a fully independent approach (Watkins et al., 2015).No significant differences in trend and interannual variability were found.GRACE data are typically provided for mid-month, while cumulated SMB values from RACMO2.3 represent the end of the month, so we linearly interpolated the GRACE data to the end of the month.Missing monthly values were linearly interpolated.Because GRACE provides mass anomalies in time δMB(t) rather than mass balance dM/dt, SMB and D must be integrated in time before being subtracted and compared with GRACE.Alternatively, we could also use the temporal derivative of the GRACE time series to obtain dM/dt, but given the inherent noise in the GRACE data this would introduce large uncertainties.

Comparing MBM and GRACE
In this paper we are principally interested in the MB of the contiguous ice sheet, the main reason being that Greenland's peripheral glaciers and ice caps (GIC) are usually as-sumed to be part of the global population of glaciers when a MB assessment is made, e.g.Radic and Hock (2011).However, because GRACE has a footprint similar to the maximum width of the ice-free tundra in Greenland (∼ 200 km), it cannot readily separate mass changes of the contiguous ice sheet from those of GIC.Moreover, GRACE mass signals also include the waxing and waning of the tundra seasonal snow cover and hydrological signals.Although the latter two processes in principle do not contribute to long-term mass changes, because of their seasonal character they do represent a significant seasonal cycle in mass loading over Greenland that varies from year to year, modifying the amplitude of the GRACE signal (Bevis et al., 2012).To enable a meaningful comparison with GRACE therefore requires integrating surface mass fluxes over the entire island, including the GIC and ice-free tundra.This is rather straightforward, because both are explicitly modelled in RACMO2.3,albeit the seasonal snow cover with a simpler (single-layer, no refreezing) snow model (Dutra et al., 2010).This islandintegrated SMB, only used in Sect.(3.1), is indicated by SMB Greenland .No marine-terminating outlet glaciers wider than 1 km wide originate from detached ice caps, so we neglect their contribution to D and assume mass changes in GIC to be solely caused by SMB, i.e.D Greenland = D and MB Greenland = SMB Greenland -D.
The most direct way to compare MBM and GRACE is to cumulate SMB Greenland and D in time and calculate the difference to get cumulative MB Greenland , i.e. δMB relative to 1 January 1958.Figure 2 shows these cumulative values (expressed in Gt), starting at 0 on 1 January 1958 when our time series of SMB Greenland starts.For reference, the equivalent mass required for 5 cm global mean sea level change is indicated on the left.Until mid-1996, cumulative D rep-resents a straight line, because its annual value is assumed constant at the 1996 value of 411 Gt yr −1 .The average value of SMB Greenland over the period 1958-1995 is within 2 % of this value (418 Gt yr −1 ), resulting in an estimated pre-1996 cumulative MB (red line) that remains close to 0, in line with previous results of Howat and Eddy (2011).The fact that the estimated 1995 value of D and the 1958-1995 average value of SMB are similar suggests that ice flow in the mid-1990s was well adjusted to the average annual mass input of the previous decades, reminiscent of an ice sheet in approximate balance (Hanna et al., 2013b).Because we do not include a seasonal cycle in D, the mass curve shows a gradual wintertime increase, when SMB Greenland exceeds D in magnitude and Greenland gains mass, and a steep summer drop, when SMB Greenland becomes strongly negative and acts together with D to remove mass from Greenland.
After 1995, following the acceleration of several large outlet glaciers in the southeast and northwest (Rignot and Kanagarathnam, 2006;Moon et al., 2012), D increased while simultaneously SMB Greenland decreased.As a result, their cumulated values in Fig. 2 curve upward and downward, respectively, and cumulative MB Greenland becomes persistently negative as a result.According to Fig. 2, the most significant mass loss derives from the last 1-2 decades, and it is a fortunate coincidence that the GRACE mission covered most of this period.The recent Greenland mass evolution from the MBM agrees qualitatively well with the GRACE observations, represented by the dark grey line in Fig. 2. Note that, because GRACE measures mass changes rather than absolute mass, the time series have been vertically offset by 1000 Gt for clarity without losing information.
Figure 3a zooms in on the time series of cumulative MB Greenland and GRACE during the overlapping period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).Although the seasonal oscillations in the GRACE time series have lost some of their amplitude because of interpolation, Fig. 3b confirms the good agreement (R 2 > 0.995) between the fully independent time series.The seasonal and interannual variations in the GRACE time series are qualitatively well reproduced, with the largest summertime mass losses in 2007, 2010 and 2012, limited mass loss in 2013 and large interannual variability in wintertime accumulation.The magnitude of a fitted linear trend over the period 2003-2014 is also similar, −294 ± 5 in MB Greenland and −270 ± 4 Gt yr −1 in GRACE.These errors represent fitting uncertainties; the real uncertainties in the trends are estimated at 20 Gt yr −1 in both time series (Van den Broeke et al., 2009).This good agreement between both methods is partly fortuitous due to compensating biases in distinct SMB components (see Sect. 4).Nonetheless, Figs. 2 and 3 inspire sufficient confidence to support a more quantitative analysis of the components of year-to-year contributions of the GrIS to global average sea level rise.

Temporal SMB variability
In order to combine SMB with annual D in Sect.3.4 we integrated SMB components over the contiguous GrIS ice mask of RACMO2.3 and over calendar years.Figure 4 shows the resulting time series of the main SMB components, with 1991-2015 trends included as dashed lines.Table 1 lists the averages and trends over the climate period 1961-1990, for which the GrIS was assumed to be in approximate balance, and the recent melting period 1991-2015.
Total annual precipitation (P tot ) typically varies between 600 and 800 Gt yr −1 , without significant trend over the full period.A small positive trend in P tot of 0.3 % per year for the period 1961-1990 is followed by an insignificant negative trend in the subsequent decades.Although total precipitation on the GrIS has not significantly changed over the last 6 decades, the RA fraction has increased in response to a warmer atmosphere.During 1961-1990, 3.3 % or 23 Gt yr −1 of the modelled precipitation on the ice sheet fell as rain, increasing to 3.9 % or 28 Gt yr −1 in 1991-2015.For the full period, the annual rain fraction varied between 2.0 and 7.0 %, the latter value occurring in 2012.Total sublimation and drifting snow erosion are relatively constant from year to year and do not show significant trends.
During 1961-1990, the small positive trends in ME and RU were offset by a similarly small precipitation increase, resulting in an insignificant trend in SMB.This changed dramatically in the ensuing period , during which ME and RU trends increased 6-to 7-fold.In combination with a small decrease in precipitation, this led to a sharp decrease in SMB of 3.3 % per year.A synchronous increase in RF has limited the mass loss: based on annual values, the refrozen fraction RF/(RA+ME) varied between 32 % and 57 %, underlining the great importance of firn processes for contemporary GrIS MB.During 1991-2015, the refrozen fraction averaged 41 %, down from 44 % in 1961-1990, signalling a decrease in the retention efficiency of the GrIS firn layer.In RACMO2.3 this is mainly caused by a decrease in firn pore space (Van Angelen et al., 2014), but in reality this effect is exacerbated by the formation of impenetrable ice lenses during warm summers, preventing the surface meltwater from reaching the deeper firn layers and using the full retention potential (De la Peña et al., 2015;Machguth et al., 2016).
Figure 4 clearly demonstrates the large interannual variability of the major GrIS SMB components P tot , ME and RU.For ME and RU, the peak of 2012 stands out, with a modelled melt flux in excess of 1000 Gt, i.e. 1 Tt, exceeding the previous record of 2010 (∼ 800 Gt) by a wide margin.Remarkably, the summer of 2013 saw a return to near-normal melt conditions, with melt close to the 1961-1990 average, while summer 2015 saw record melting in the northern reaches of the ice sheet (Tedesco et al., 2016).This exceptional interannual variability in the melt climate of the GrIS points towards important roles Table 1.Contiguous ice sheet (GrIS) averages , Gt yr −1 with standard deviation) and trends (1961-1990 and 1991-2015, Gt yr −2 , with standard error) of SMB components, discharge (D) and mass balance (MB).

SMB
Average Trend Average Trend component (1961-1990) (1961-1990) (1991-2015)   for large-scale atmospheric drivers (Fettweis et al., 2013b;Hanna et al., 2013aHanna et al., , 2014Hanna et al., , 2016;;McLeod and Mote, 2016) and local feedback processes.Especially important is the albedo-melt feedback (Box et al., 2012), which constitutes the darkening of snow once it has melted, as well as the lengthening of the exposure of dark, bare ice in the ablation zone once the winter snow has melted away (Tedesco et al., 2011).However, precipitation, where feedbacks play a lesser role, is also highly variable from year to year; for instance, P tot increased by ∼ 300 Gt yr −1 between 1971 and 1972, a year-to-year change equivalent to 40 % of the longterm average.Fitting a linear trend to the standard deviation of running decadal values reveals that precipitation variability decreased by ∼ 30 Gt yr −1 , while that of runoff increased by approximately the same amount.The reasons for this are presently not clear.

Spatial SMB variability
In this section we discuss the spatial distribution of changes in LWB and SMB components between the climatic period 1961-1990 and the recent period of GrIS mass loss.For a description of spatial differences in D the reader is referred to e.g.Enderlin et al. (2014) and Csatho et al. (2014).To maximise the length and to avoid spurious trends, we let the recent period start in 1991 rather than in 1995, when the first changes became noticeable.Figures 5 to 8 show the average for 1961-1990 (a) and the difference between 1991-2015 and 1961-1990 (b)  spectively.In these maps, mass flux (difference) is expressed as kg m −2 yr −1 , equivalent to mm w.e.yr −1 .Note that over non-glaciated areas, RACMO2.3uses a simpler snow model that does not calculate refreezing, only melt and runoff; SMB is therefore only calculated and physically meaningful over glaciated areas.
Figure 5a shows that over tundra, the average annual ME rate is limited by the annual snowfall, which explains the generally lower values when compared to the adjacent glaciated areas, where ablated ice is continuously replenished by glacier flow.Over glaciated areas, ME increases strongly with decreasing elevation and latitude, reaching 3500 mm w.e.yr −1 in the low-lying parts of the southwestern GrIS.Although higher observed melt rates have been reported for the GrIS ablation zone, this concerns mainly locally very dark ice surfaces and/or isolated glacier tongues surrounded by ice-free land that are not well resolved at the model resolution of 11 km (Fausto and van As, 2012).
Most of the meltwater produced at higher elevations refreezes in the cold firn.Figure 6a shows that RF peaks in the lower accumulation zone between 1000 and 2000 m a.s.l., where significant summer melt occurs yet the firn layer still has sufficient pore space to store the meltwater.As a result, all RU from the GrIS occurs from ice sheet regions roughly below 2000 m a.s.l. in the south and 1500 m a.s.l. in the north (Fig. 7a).The resulting 1961-1990 SMB distribution (Fig. 8a) shows relatively wide (50-150 km) ablation zones in the dry southwest, north and northeast of the GrIS and narrow (10-50 km) ablation zones in the wetter and therefore steeper-sloped northwest and southeast of the GrIS.
Relative to 1961-1990, melt in 1991-2015 has increased everywhere over the GrIS (Fig. 5b).The change is not statistically significant on the higher ice sheet domes, where melt occurs intermittently.Integrated over the contiguous GrIS, melt has increased from 433 Gt yr −1 in 1961-1990 to 581 Gt yr −1 in 1991-2015, an increase of 34 % (see Table 1).Not all additional liquid water reaches the ocean: part of the mass loss is buffered by increased RF.Table 1 and Fig. 4 show that 29 % of the combined increase in ME and RA is buffered by increased RF. Figure 6b clearly shows that this increase in RF is confined to areas where the firn layer has sufficient storage capacity, i.e. well above the equilibrium line.As a result, the increase in RU is confined to the ablation zone and the lower accumulation zone (Fig. 7b).Resulting changes in the SMB (Fig. 8b) have two components: a significant decrease in the ablation zone and lower accumulation zone that mirrors the change in RU, and a partially significant increase in the interior, owing to an increase in snowfall.The result is that SMB gradients have steepened along the perimeter of the ice sheet, which is also visible in high-resolution altimetry (Helm et al., 2014).
In line with observations, the neighbouring ice caps in the Canadian Arctic, Iceland and Svalbard have also experienced strongly increased melt rates.ME decreased only over some non-glaciated areas, mainly as a result of decreased (winter) snow accumulation.Only the interior parts of the highest and coldest ice caps in the northern Canadian Arctic remain free of runoff, while ice caps in the southern Canadian Arctic, Iceland and Svalbard all produce runoff from the entire ice surface.

Temporal MB variability
Figure 9 combines GrIS integrated values of SMB and D into ice sheet MB with uncertainties as defined in Sect.2.2.Linear trends for the period 1991-2015 are indicated by dashed lines.The equivalent sea level rise (eq.SLR) for negative MB is provided on the lower right axis.The MB values be-fore 1996 are uncertain because reliable estimates of D are missing, although previous work reported little difference between discharge estimates from the early 1960s and the mid-1990s (Rignot et al., 2008).Before 1995, under the assumption of constant ice discharge, we see that MB typically varied between +200 and −200 Gt yr The period-average mass loss we obtain here can for instance be compared to the Ice sheet Mass Balance Intercomparison Experiment (IMBIE) results in Shepherd et al. (2012), which used data of gravimetry, altimetry and MBM to reconcile differences in ice sheet MB.For the periods 1992-2000 and 2000-2011, average values for GrIS MB in that (this) study are −51 ± 65 (−35 ± 79) and −211 ± 37 Gt yr −1 (−236 ± 86 Gt yr −1 ).Although not fully independent, this good agreement suggests that the uncertainties in these numbers may actually be smaller than stated.
The decrease in MB since 1991 is significant and indicates an acceleration of the mass loss of 16.8 ± 2.8 Gt yr −2 .This is somewhat less than the 21.9 ± 1 Gt yr −2 reported by Rignot et al. (2011) for the period 1992-2010, which is obviously caused by the inclusion of years after 2012 with higher MB.For the same period 1992-2010 we obtain 21.8 ± 3.7 Gt yr −2 , i.e. a number very close to Rignot et al. (2011).Using trends in SMB and D as reliable indicators for mass loss partitioning, the acceleration is caused for 61 % by a decrease in SMB and for the remainder by an increase in D. Again owing to the inclusion of the years 2013-2015, this partitioning between the two components is somewhat closer to equality than reported in Enderlin et al. (2014).

Discussion and conclusions
These results show that GrIS MB has been persistently negative since 1998 and continues to be negative in spite of a temporary rebound in 2013.The mass loss that occurred each year between 2006 and 2012 was unprecedented since 1958.How significant are these recent mass losses and how do they impact the future resilience of the ice sheet?
First we note that the capacity of the firn layer to buffer runoff (Harper et al., 2012) is compromised by increased meltwater refreezing.The associated release of latent heat warms the firn, reducing its cold content and enhances dry firn compaction, which together with the refrozen mass reduces the pore space where meltwater can be stored and refrozen.Figure 10 shows results of calculations with an offline firn densification model forced with RACMO2.3 output.It shows that firn temperatures in response to enhanced melt and refreezing have increased by up to 5 K in the lower accumulation area (Fig. 10a).Firn air content, defined as the vertically integrated depth of the air column (expressed in m), decreased by up to 6 m in the same areas (Fig. 10b).These changes are highly significant considering that typical values for total firn air content in the dry snow zone range between 20 and 25 m (Kuipers Munneke et al., 2015).Recent research has shown that warm summers can generate thick ice layers that prevent meltwater from reaching the deeper pores, further reducing the meltwater buffering capacity (De la Peña et al., 2015;Machguth et al., 2016).
The mass overturning rate of an ice sheet is defined as its total mass divided by the annual mass gain.We approximate the latter by the 1961-1990 average of SN-SU tot = 632 Gt yr −1 , assuming most rainfall to fall on the lower ice sheet margins and to runoff quickly.Combined with the estimated volume of the GrIS of 2.67 × 10 6 km 3 (Vaughan et al., 2013) and an ice density of 900 kg m −3 , we obtain a mass overturning rate of 4200 years.If we interpret this as the reaction timescale of an ice sheet to adjust its dynamics to changes in accumulation, it is clear that the recent changes owing to melt, with a typical timescale of decades, by far outpace a potential ice sheet thickening owing to increased snowfall in a warmer climate.In spite of that, at the current rate of mass loss, it would still take over 10 000 years to melt the entire ice sheet.
Alternatively, we may state that the GrIS is significantly out of balance, noting that the average 1991-2015 mass loss (171 Gt yr −1 ) represents a sizeable fraction (27 %) of the annual mass gain.Because D is definite positive, the situation in which SMB becomes persistently negative leads to a definite negative MB, even when the ice sheet has lost contact with the ocean, i.e.D = 0.This is therefore sometimes labelled a "tipping point" for GrIS MB, beyond which the ice sheet will not be able to recover.At the current rate of SMB decrease (10.2 ± 2.3 Gt yr −2 ), this tipping point would be reached between 2024 and 2043.The limited length of the time series and our incomplete knowledge of the main drivers of changes in SMB and D for now preclude firm statements about how realistic such a scenario is.It is therefore desirable that both MB components be reconstructed as far back in time as possible; this will require the smart use of climate archives, such as firn cores from the accumulation zone of the ice sheet (Box and Colgan, 2013), robust reanalysis products that cover the full 20th century (Hanna et al., 2011;Lee and Biasutti, 2014) and early satellite products and photogrammetry (Kjeldsen et al., 2015).
Because pre-1995 values of SMB and D are similar, in this study it was not necessary to define a reference period from which cumulative anomalies are defined (Van den Broeke et al., 2009).Instead, absolute mass fluxes could simply be integrated in time and subtracted to obtain ice sheet MB for comparison with GRACE (Fig. 2).However, this pre-1995 agreement is almost certainly in part fortuitous, because uncertainties in especially SMB, which is a modelled quantity, and to a lesser extent in D, which is largely observed, remain significant.For instance, current model horizontal resolution of RACMO2.3 (11 km) is insufficient to resolve the individual, low-lying outlet glaciers of the GrIS where runoff is especially large; as a result RU increases  when the 11 km field is statistically downscaled to 1 km resolution (Nöel et al., 2016).This unresolved mass loss is likely in part error-compensated by snowfall in RACMO2.3being underestimated in some regions of the ice sheet (Overly et al., 2016).In a recent study, it was moreover demonstrated that while RACMO2.3tends to time drifting snow events well, the model likely overestimates drifting snow transport and therewith drifting snow sublimation (Lenaerts et al., 2012).
This leads to uncertainties in SMB of 60-100 Gt yr −1 , clearly dominating the uncertainty in MB (Fig. 9).
To reduce these biases and increase our diagnostic and prediction skills of GrIS MB, it is imperative that SMB and firn models are further improved and their horizontal resolution enhanced.This can be achieved through statistical/dynamical downscaling in combination with targeted in situ observations.Examples of important processes that are poorly or not at all represented in current models are interactive snow/ice darkening by future enhanced dust/black carbon deposition or microbiological processes (Stibal et al., 2012), and sub-, supra-and englacial hydrology, including vertical and horizontal flow of meltwater in firn or over ice lenses (De la Peña et al., 2015;Machguth et al., 2016).Other emerging research topics of GrIS melt climate are the impact of atmospheric circulation changes on Greenland melt (Hanna et al., 2013a(Hanna et al., , 2014(Hanna et al., , 2016;;McLeod and Mote, 2016;Tedesco et al., 2013), the impact of rain on ice sheet motion (Doyle et al., 2015), the effect of liquid water clouds on the surface energy balance and melt (Bennartz et al., 2013;Van Tricht et al., 2016) and the increased role of turbulent heat exchange during strong melting episodes over the margins of the GrIS (Fausto et al., 2016).Finally, it is desirable that, once developed and tested, a single, sophisticated snow model is used to simulate both the deep firn layer over the ice sheet and the seasonal snow cover over the tundra.

Data availability
All data presented in this study are available without conditions from the authors.

Figure 1 .
Figure 1.RACMO2.3 Greenland domain with ice sheet mask (white), land mask (green) and elevation contours (dashed lines, every 500 m a.s.l.).Boundary relaxation zone indicated by dots every other grid point.

Figure 2 .Figure 3 .
Figure2.Cumulative surface mass balance (SMB Greenland ) for the full island (indicated by subscript "Greenland", amber line), i.e. including marginal glaciers and ice caps (GIC) and seasonal snow on the tundra; cumulative ice discharge D (blue line) and resulting cumulative mass balance MB Greenland (red line).GRACE time series included (grey line) has been offset by 1000 Gt for clarity.

Figure 9 .
Figure 9. Annual values of D, SMB and MB integrated over the contiguous GrIS.Dashed lines indicate 1991-2015 trends.The equivalent sea level rise (eq.SLR) for negative MB is provided on the right axis.

Figure 10 .
Figure 10.Modelled 1991-2015 minus 1961-1990 difference in 2-10 m average firn temperature (a) and change in firn air content (b) for the contiguous GrIS.Dashed contours are 500 m elevation intervals; thick solid contour represents glacier mask. 
Modelled 1961Modelled  -1990Modelled   average refreezing (a) and 1991Modelled  -2015Modelled   minus 1961Modelled  -1990 difference (b) difference (b).Stippled areas indicate differences that are not significant at the 95 % level.Dashed contours are 500 m elevation intervals; thick solid contour represents glacier mask.Note that refreezing is not considered outside of the glacier mask.
of ME, RF, RU and SMB, rewww.the-cryosphere.net/10/1933/2016/The Cryosphere, 10, 1933-1946, 2016 M. R. van den Broeke et al.: Greenland ice sheet and sea level rise −1 , with an average close to 0. After 1995, MB becomes persistently negative, with a www.the-cryosphere.net/10/1933/2016/The Cryosphere, 10, 1933-1946, 2016 minimum in 2012 of −446 ± 114 Gt yr −1 , equivalent to an SLR of 1.2 ± 0.3 mm yr −1 .In 2013, MB sharply increased in response to a summer with near-normal surface melt conditions; this temporarily limited GrIS mass loss but did not eliminate it, because D remained elevated.After 2013, ME and RU increased once more, reducing SMB and increasing the mass loss to values again approaching 1 mm eq.SLR in 2015.