Timing and origin of recent regional ice-mass loss in Greenland

Gravity Recovery and Climate Experiment (GRACE; Tapley et al., 2004), surface-ice velocities from Inteferometric Synthetic Aperture Radar (InSAR; Rignot and Kanagaratnam, 2006) together with output of the regional atmospheric climate modelling (RACMO2/ GR; Ettema et al., 2009), and surface-elevation changes from the Ice, cloud and land elevation satellite (ICESat; Sørensen et al., 2011). We show that changing ice discharge (D), surface melting and subsequent run-off (M/R) and precipitation (P) all contribute, in a complex and regionally variable interplay, to the increasingly negative mass balance of the GrIS observed within the last decade. Interannual variability in P along the northwest and west coasts of the GrIS largely explains the apparent regional mass loss increase during 2002–2010, and obscures increasing M/R and D since the 1990s. In winter 2002/2003 and 2008/2009, accumulation anomalies in the east and southeast temporarily outweighed the losses by M/R and D that prevailed during 2003–2008, and after summer 2010. Overall, for all basins of the GrIS, the decadal variability of anomalies in P, M/R and D between 1958 and 2010 (w.r.t. 1961–1990) was signiﬁcantly exceeded by the regional trends observed during the GRACE period (2002–2011).


Introduction
The GrIS is an important part of the climate system that interacts with other components in many ways, including the discharge of freshwater by surface melting and iceberg calving, changing sealevel and potentially affects the ocean's thermohaline circulation. Satellite observations have led to a general consensus about the increasingly negative mass balance for the GrIS since the late 1990s (e.g. Krabill et al., 2004;Cazenave, 2006;Hanna et al., 2005;Thomas et al., 2006;Zwally and Giovinetto, 2011), with a significant contribution of $ 0:77 0:1 mm=yr (2003van den Broeke et al., 2009) to the observed global sea level rise of 3.170.6 mm/yr (1993Ablain et al., 2009). A negative ice sheet mass balance means that the amount of mass gained by surface accumulation through rain and snow (P) is smaller than that lost by ablation, i.e. the combined effect of melting and subsequent meltwater run-off (M/R) and sublimation (SU), as well as ice discharge across the grounding line (D). The difference of P-SU-R approximates the surface-mass balance (SMB) neglecting snowdrift sublimation (o 6% of SMB; Ettema et al., 2009), SMB-D represents the mass balance of the ice sheet, whereas the time integrated (i.e. cumulative) SMB and D anomalies (here, w.r.t. 1961-1990) represent change in mass stored by the ice sheet, that are measured by GRACE and ICESat. For the GrIS, abrupt changes in ice motion have locally been linked to penetration of surface meltwater to the bed (Zwally et al., 2002;Bartholomew et al., 2010). For individual marineterminating glaciers, intruding warm ocean waters have been identified as a cause for enhanced submarine melting that leads to the thinning or disintegration of floating ice tongues and grounding line retreat, with a resultant increase in discharge (Rignot et al., 2010;Holland et al., 2010). Regional atmospheric climate models supported by passive microwave data indicate a changing ratio between surface melt and precipitation, leading to an increasing imbalance between accumulation and ablation over the last two decades (Fettweis et al., 2011). Enhanced meltwater production is consistent with a general warming of the Arctic atmosphere in the past 30 yr (Mote, 2007;Hanna et al., 2008).
Here, we derive regional mass variations of the GrIS using three methods: (i) measuring variations in the Earth's gravitational potential with GRACE, (ii) contrasting SMB from the regional atmospheric climate model RACMO2/GR  with D over the grounding line from InSAR (hereafter SMB-D), and (iii) determining surface elevation changes with the Geoscience Laser Altimeter System on ICESat. Each method has its strengths and weaknesses; e.g. ICESat altimetry measurements are accurate, but yield geometric volume changes that need to be converted to mass changes, relying on the accuracy of a snow/ice density model. SMB-D gives the best physical insight into the processes driving mass change, but it relies on subtracting two large numbers, where uncertainties of the difference can be significant w.r.t. SMB or D. In this respect, GRACE is advantageous, because it provides at regular (sub-)monthly intervals the gravity field disturbances induced by net ice-mass changes. However, the GRACE resolution of $ 400 km is much coarser than that of ICESat (track-spacing o 30 km for GrIS) or the regional climate model RACMO2/GR (11 km). Also, GRACE is influenced by all mass changes in the Earth system, ranging from short-term atmospheric pressure variations to the long-term redistribution of mantle material caused by the glacial-isostatic adjustment (GIA). Despite these limitation, van den Broeke et al. (2009) and Rignot et al. (2011) have demonstrated a promising agreement of SMB-D with GRACE for mass trends of the entire GrIS, which are also consistent with ICESat (Sørensen et al., 2011). As shown below, we find excellent basin-scale agreement of the temporal mass variations of GrIS from GRACE and SMB-D, and trends from ICESat, allowing us, for the first time, to provide a robust regional separation of the relative importance of P, M/R and D in causing recent GrIS mass loss.

GRACE
In this study, 103 unconstrained monthly mean GRACE gravity field solutions of the Centre for Space Research at University of Texas, Austin, USA (CSR RL04; Bettadpur, 2007) were used, covering the time period from August 2002 to September 2011. The solutions of CSR RL04, available from http://isdc.gfz-potsdam.de/, are represented as Stokes potential coefficients complete to spherical-harmonic degree and order 60. During post-processing of the GRACE data, each monthly set of Stokes coefficients is smoothed with an isotropic spatial averaging function (Sasgen et al., 2006), which optimizes the trade-off between spatial resolution -typically between 400 and 500 km -and noise with respect to the expected signals. For further regionalization, the weight of low degree and order coefficients are reduced according to the correlation of predicted GrIS mass signal (from ICESat) and its superposition with GIA in response to the history of the Laurentide ice sheet. The band-pass filtered gravity fields are corrected for GIA and simultaneously inverted for monthly mean mass anomalies within seven major GrIS drainage basins and Ellesmere Island with a forward modeling approach (Appendix A1). The uncertainty associated with the filtering and inversion procedure for the entire GrIS was estimated to be 75 mm equivalent water height for the monthly solutions and about 77% for the linear trends (Appendix A2). Following Paulson et al. (2007), the GIA correction is a version of ICE-5G (Peltier, 2004), adjusted to reconcile with the GRACE estimates, assuming upper-and lower-mantle viscosities of 6 Â 10 20 Pa s and 2 Â 10 22 Pa s (Wolf et al., 2004), and an elastic lithosphere with a thickness of 100 km. The uncertainty of the GIA corrections is assessed with the alternative retreat history of Lambeck and Chappell (2001), and by varying upper-and lowermantle viscosities between 2 Â 10 20 and 8 Â 10 20 Pa s and between 5 Â 10 21 and 4 Â 10 22 Pa s, respectively, and the thickness of the lithosphere by 720 km. Depending on the applied GIA scenario (i.e. load history and viscosity distributions), GrIS mass loss presented here could be stronger by À 3 Gt/yr or weaker by 20 Gt/yr (Table 1 and Appendix A3); this uncertainty agrees with the range (difference between minimum and maximum correction) of 2976 Gt/yr provided in the more complete GIA assessment of Barletta et al. (2008).

SMB-D
The surface mass balance (SMB) estimate is obtained from the regional climate model RACMO2/GR  providing monthly means of the individual SMB components, i.e. precipitation, evaporation/sublimation, run-off, melting and re-freezing. The model is forced at its lateral boundaries and at the sea surface with atmospheric reanalysis data of the European Centre for Medium-Range Weather Forecasts (ECMWF) and run at a spatial resolution of $ 11 km. Due to unavailability of data, our SMB time series is 1 yr shorter than that of GRACE, ending August 2010. The uncertainty of the SMB provided by RACMO2/GR is estimated to be 718% . It should be noted that RACMO2/GR and other regional climate models, e.g. MAR (Fettweis et al., 2011), Polar MM5 (Box et al., 2004), or ECMWF reanalysis (Hanna et al., 2011), may differ considerably in the absolute values provided for P and M/R; the temporal behaviour, however, is captured by all models, and the uncertainty is considerably reduced by building anomalies w.r.t. the reference period 1961-1990, as we will do in the following. The discharge (D) estimate is based on InSAR satellite observations of ice flow within 38 drainage basins (Rignot and Kanagaratnam, 2006) using the speckle tracking method, updated until the year 2009; due to unavailability of data, D is assumed to remain unchanged in 2010, compared to 2009. Ice-velocity is measured at flux gates located upstream of the grounding line and therefore a SMB estimate based on RACMO2/GR for the region located in between is added to the measurement at the flux gates. Together with ice thickness estimates, and assuming a volume-mass density of ice of 910 kg/m 3 , a fixed grounding line position and a constant velocity within the ice column, the ice exported by outlet glaciers is quantified with an accuracy of about 714% (Rignot et al., 2011). The discharge data of the 38 drainage basins are aggregated into seven drainage basins optimized for the GRACE analysis by interpolating discharge values by area, if necessary. Details on the extrapolation of D to unmeasured in regions is provided in Appendix A to this letter. The annually resolved discharge for each basin is then re-sampled at monthly time intervals by linear interpolation. Then, anomalies of SMB, as well as of D are calculated w.r.t. their respective 1961-1990 mean fields, and integrated in time to obtain the cumulative storage changes of the ice sheet, comparable with the GRACE observations. From the deviation of the SMB and D reference fields  an additional uncertainty of the trends in cumulative difference SMB-D is estimated to be below 78 Gt/yr for all basins, except basin B (mean D is 17 Gt/yr larger) and 730 Gt/yr for the entire GrIS (Appendix A4).

ICESat
Surface elevation estimates from ICESat in the period October 2003-October 2009 (Zwally et al., 2010) are analysed to derive the mean rate of volume change of the GrIS. Data culling based on the shape of the return signal, and various quality indicators and warnings are applied to remove less reliable surface elevation estimates. The rate of elevation change, dH/dt, is derived at a resolution of 500 m along-track, by applying a least squares linear regression on all data in each 500 m segment, with the assumption that the surface elevation varies linearly with position, time and a cosine and sine terms. The total volume change is found by fitting a smooth surface, which covers the entire GrIS, through the dH/dt estimates and using the variances from the least squares regression as weights (Sørensen et al., 2011). The error on the volume change is estimated using a bootstrap approach. The observed geometric rates are corrected for elevation changes not related to ice mass changes, including firn compaction, vertical bedrock movement, and an ICESat intercampaign elevation bias. The corrected volume rate is converted into a mass balance estimate by snow/ice density modelling (Sørensen et al., 2011). The firn compaction and density models are forced by climate parameters from a high resolution regional climate model (HIR-HAM5 RCM; Christensen et al., 2006) that is independent from RACMO2/GR. The uncertainty of the total mass balance is estimated to be 728 Gt/yr (Appendix A4).

GrIS mass balances from GRACE, SMB-D and ICESat
For the entire GrIS, our GRACE analysis results in a mass balance of À 238729 Gt/yr for October 2003-October 2009, the period with ICESat data available ( Fig. 1 and Table 1). This value is consistent with our SMB-D estimates of À260753 Gt/yr, and our ICESat estimates of À 245728 Gt/yr, as well as with another recent GRACE estimate of À 230733 Gt/yr obtained with a different GRACE inversion strategy ( loss of À 240718 Gt/yr, which is explained by significant acceleration of mass loss, discussed later. Earlier GRACE estimates are to some extent discordant with our inferred rate of GrIS mass loss; e.g. Barletta et al. (2008) present a preferred ice mass loss of À101722 Gt/yr (from 2003 to early 2006), which is in good agreement with analyses from Velicogna and Wahr (2005) (from April 2002 to July 2004), Ramillien et al. (2006) and Luthcke et al. (2006), but by a factor of about two lower than the estimates of Velicogna and Wahr (2006) (from April 2002to April 2006 and Chen et al. (2006) (from April 2002to November 2005. Part of this discrepancy may be related to the observation period, owing to the strong interannual mass variability of the GrIS. However, truncating our time series from January 2003 to July 2005, results in GrIS mass loss of À171731 Gt/yr, which still does not lead to an agreement with Barletta et al. (2008) or Luthcke et al. (2006), despite similar observation periods. An updated mascon estimate shows a stronger mass loss of À 18076 Gt/yr (August 2003(August -2009Pritchard et al., 2010), yet it remains considerably below our GRACE/SMB-D/ICESat estimate for roughly the same period. The earlier ICESat estimate of À18174 Gt/yr (October Shown are GRACE analysis (grey circles), as well as the relative contribution of cumulative SMB anomalies (red), as well as positive (blue) and negative (magenta; À 1 to À 2 Gt/yr for regions A and B) cumulative D anomalies, and the trends derived from ICESat (yellow circles). Light colouring indicates uncertainties provided in Table 1. Underlain is the spatial rate of geoid-height change over Greenland calculated from the forward model adjusted to the GRACE observations. An elastic-compressible earth model is assumed and the lower and upper spherical-harmonic cut-off degrees are 2 and 340, respectively. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.) (2011), is broadly consistent, yet somewhat lower than our GRACE estimate for the same period (À 218728 Gt/yr). Compared at the regional scale, our GRACE analysis (2003)(2004)(2005)(2006)(2007) yields stronger mass loss in the west and north of the GrIS with regard to ICESat (Zwally and Giovinetto, 2011); but except for these two regions, deviations between both data sets are o 10 Gt=yr.

2003-2007) from Zwally and Giovinetto
Although a dedicated inter-comparison of GRACE analyses is yet to be undertaken, Schrama and Wouters (2011) suggest that differences can be related to treating the leakage of the gravity field signal from and to regions outside Greenland, as well as to compensating for signal attenuation due to filtering. Improved strategies for the gravity field determination of the GRACE processing centres are also relevant; the baseline accuracy of the GRACE gravity fields improved by a factor of four from RL03 to RL04. Earlier GRACE estimates based on the gravity field inversion with optimal basin averaging functions may invoke uncertainties associated with the attenuation by filtering and leakage compensating factor at least 716% (Schrama and Wouters, 2011). Here, however, no calibration factor is applied, and the error on the trends associated with the filtering and inversion procedure is estimated to be o 5% for the entire GrIS (Appendix A3).

Regional ice-mass balances from GRACE, SMB-D and ICESat
At the basin scale, we find agreement of GRACE, SMB-D and ICESat within their respective ranges of uncertainties ( Fig. 1), for the common time interval October 2003-October 2009. In the east (C), southeast (D), southwest (F) and northwest (G) of the GrIS, GRACE and ICESat consistently exhibit mass loss stronger than $ À40 Gt=yr, while mass loss in the northern basins A and B are much weaker ( $ À15 Gt=yr). In the northwest (G), SMB-D results in smaller mass loss compared to GRACE and ICESat, which is most likely due to an underestimated D anomaly, to be addressed later. For all regions, except the east (C) and southeast (D), anomalous SMB markedly exceeds anomalous D (w.r.t. 1961-1990 mean). In the north (A) and northeast (B), the anomaly in D (w.r.t. 1961-1990 mean) is close to zero ( $ À1 Gt=yr) and the mass loss is entirely explained by a distinctly lower SMB compared to the 1961-1990 mean caused by enhanced meltwater production.
3.3. Acceleration of regional ice-mass change observed by GRACE and SMB-D The GRACE time series for the entire GrIS exhibits a statistically significant mass loss acceleration of À15 77 Gt/yr 2 , which is in good agreement with that of SMB-D ( À 2177 Gt/yr 2 ; common time interval August 2002-August 2010). Our estimate confirms another recent GRACE estimate of À17 78 Gt/yr 2 (from April 2002 to June 2010; Rignot et al., 2011) derived with a different gravimetric inversion approach, even though, due to the large interannual mass variability of the GrIS, the exact coverage of the short time interval (9 yr for GRACE), may influence the acceleration estimates. Extending our GRACE time series until September 2011, we obtain a mass loss acceleration of À 1877 Gt/yr 2 , which is in good agreement with the long-term SMB-D estimate of À22 71 Gt/yr 2 for the years 1992-2009 (Rignot et al., 2011). Our regional analysis locates the origin of increasingly strong mass loss primarily in the west (basin F, À16 76 Gt/yr 2 ) and northwest (basin G, À 1074 Gt/yr 2 ), for 2002-2011. Comparison with the accelerations from SMB-D (Table 2) reveals that on one hand, a SMB anomaly -positive from 2002 until 2005 and negative afterwards (w.r.t. 1960-1961 mean), constitutes a large part ( À 2473 Gt/yr 2 ; 2002-2010) of the apparent acceleration of for basins F and G observed with GRACE. But also increasing D by glacier speed-up accelerated mass loss by another 2 71 Gt/yr 2 for basin F following a long-term trend starting in the late 1990s, and 371 Gt/yr 2 for basin G (Fig. A2 of Appendix).

Separation of GRACE-observed mass loss into P, M/R and D
In addition to the trend and acceleration terms, Fig. 2 reveals a very high level of agreement between the time series of regionalscale GrIS mass change from GRACE and SMB-D for the common time interval August 2002-August 2010 (Fig. 2). The de-trended monthly time series (i.e. seasonal harmonic and linear trend removed) exhibit zero-lag correlations of 4 0:6, except for basin A (0.5) and E (0.3), which are statistically significant at the 95% confidence level, allowing us to identify the relative importance of P, M/R and D in the regional ice-mass budget. Pronounced steps in the GRACE mass loss observed in the west (basin F) are clearly related to episodes of enhanced M/R in 2002 and 2003, and from 2006 to 2008 (Fig. 2)-this is consistent with elevated seasonal maximum daily temperature anomalies (June-August) at Aassiaat, West Greenland (68.71N, À 52.81) (Mote, 2007(Mote, ) in 2003(Mote, , 2005(Mote, and 2007, as well as record melt extent for the entire GrIS in these years (Fettweis et al., 2011).

Changes in the northwest of the GrIS
In the northwest (basin G), this and other GRACE analysis (Khan et al., 2010;Schrama and Wouters, 2011;Chen et al., 2011) indicate a rapid transition to stronger negative mass balances after summer 2005, an area which has been dynamically thinning since the 1990s (Krabill et al., 2000). Regarding the main SMB constituents our analysis indicates that P is important in explaining the abrupt change in trend as a consequence of higher accumulation (i.e. higher SMB) before the year 2005, and lower afterwards w.r.t. the 1961-1990 mean (Fig. 2). Starting in winter 2009/2010, enhanced accumulation again slightly weakens the negative ice-mass balance observed with GRACE; the mass loss in the west (basin F), however, remains strongly negative for the year 2011.
It should be noted, however, that the agreement between SMB-D and GRACE is somewhat poorer in the northwest (basin G) than in other basins. We consider the GRACE trend to be more reliable due to the agreement with our ICESat estimate in the period 2003-2009 (Fig. 1). In addition, the distinct shift in the GRACE time series to stronger negative mass balances is recorded with the GPS station Thule (76.51N, À 68.81) showing a similar temporal behaviour in the deformation of the elastic lithosphere, induced by changes in the ice mass balance (Khan et al., 2010). Moreover, localized surface-lowering in this area detected with ICESat (Fig. A3 of Appendix) between 2003 and 2008 is also an indication of considerable glacier thinning. This may have been missed in the estimates of D due to incomplete ice thickness mapping in this region (Rignot and Kanagaratnam, 2006). Considering that glacier speed-up and ice-front retreat between 2000/2001 and 2005/2006 have been documented for the northwest (Joughin et al., 2010), in addition to P, part of the change in trend observed with GRACE could be related to, potentially rapid, changes in D. If we consider GRACE minus SMB-D as an estimate for the missing portion of D, we deduce that non-surveyed dynamic thinning could be responsible for an additional 24713 Gt/yr of mass loss, implying an imbalance for the entire region of about 40% (w.r.t. 1961-1990).

Changes in the north of the GrIS
The moderate negative mass balances in the north (A) and northeast (B) of the GrIS (À 1675 Gt/yr and À 1275 Gt/yr, respectively; GRACE for 2002-2011) are mainly caused by reduced SMB due to enhanced M/R, modulated by a 1-(A) and 2-yr (B) period of mass gain starting in winter 2005/2006. For Ellesmere Island, located $ 40 km northwest of the GrIS at its closest point, mass loss is À3575 Gt/yr and significantly accelerating at À 875 Gt/yr 2 (GRACE for 2002-2011), which is consistent with another recent estimate from GRACE, SMB-D and ICESat (Gardner et al., 2011). On Ellesmere Island, however, a gradual increase in M/R in direct response to increasing summer temperatures (Gardner et al., 2011) induces most of this acceleration (þ671 Gt/yr 2 ; 2002-2010), indicating a sensitive response of this region to the climatic forcing.

Conclusion
For all regions and the GrIS as a whole, the trends imposed by anomalies in M/R and D after 9-yr (2002-2011) significantly exceed decadal variability of trends in P for 1958-2010, meaning that GRACE and ICESat record long-term changes of the GrIS. For the west and northwest (basins F and G, respectively), the joint acceleration of M/R and D, significantly exceed the interannual variability in P, suggesting that, despite a strong contribution of P, 9-yr of GRACE data contain a long-term climate signal of mass loss acceleration in these regions ( À 2677 Gt/yr 2 ; basin F þG; 2002-2011). Part of the recent apparent acceleration is, however, due to gradual variations in P. For the GrIS as a whole, interannual mass variability during the 9-yr period mainly arises from P variations in the east (basin C), southeast (basin D) and west (basin F), and obscures the acceleration of M/R and D persisting since the 1990s (Fig. A2 of Appendix). Weaker mass loss rates observed with GRACE in the east (basin C) and to a lesser extent in the southeast (basin D) after winter 2008/2009, which have previously been attributed to decreasing ice-dynamic flow (Chen et al., 2011), are mainly a consequence of enhanced accumulation.
In summary, we show that GRACE, SMB-D and ICESat data (October 2003-October 2009) are sufficiently mature to provide consistent regional scale mass balances of the GrIS. This yields new insight into the complex mass loss behaviour of the GrIS showing that all mass balance components (P, M/R and D) are relevant in the regional mass budget. In particular, we demonstrate that variations in P substantially influence the year-to-year ice-mass balance, and care has to be taken in interpreting GRACEderived net mass changes in terms of changing ice dynamics. We find that currently (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) that the M/R production is enhanced by at least 39-59% compared to 1961-1990 for all basins, but most significantly in the northwest (basin G; 90-134%), whereas D is up to 40-60% larger in the east (basin C). For all basins, M/R has more strongly increased than D for the entire GrIS (50-74% for M/R and 18-27% for D, w.r.t. 1961-1990), owing to the sensitivity of M/R compared to D to increasing surface-air temperature, observed in the Arctic during last decade (Hanna et al., 2008). It remains to be seen, whether dynamic thinning will follow this melting trend, particularly in the north (basin A), where M/R increased by 75-111% (compared to 1961-1990 mean), and so far only minor changes in ice flow are observed (Moon and Joughin, 2008). Continuous satellite and in situ observations will be important to place increasing confidence on the long-term trends and follow the ongoing evolution of the GrIS.

Acknowledgments
Ingo Sasgen would like to acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through Grant SA 1734/2-2 (SPP1257). We would like to thank the German Space Operations Center (GSOC) of the German Aerospace Center (DLR) for providing continuously and nearly 100% of the raw telemetry data of the twin GRACE satellites. This Table 2 Regional mass trends and accelerations of the GrIS from GRACE and SMB-D for the common time interval August 2002-August 2010, as well as the trends and acceleration for the individual components of SMB-D; cumulative anomalies (w.r.t. 1961-1990 mean) in total precipitation (P), melting/run-off (M/R) and discharge (D). The GRACE results are provided for CSR RL04 and the GIA correction ICE-5G n (Appendix A3).

Drainage basin
Area (10 3 km 2 ) GRACE CSR RL04 RACMO2/GR þInSAR  methodology and provided his GIA modelling code; Sebastian B. Simonsen converted ICESat volume changes to mass changes by firn compaction and density modelling; all authors discussed the material presented and contributed to the writing of the manuscript.

A.1. GRACE post-processing
In this study, 103 unconstrained monthly mean GRACE gravity field solutions of the Centre for Space Research University of Texas, Austin, USA (CSR RL04; August 2002 to September 2011 Bettadpur, 2007) are filtered and inverted for mass changes using a forward modeling approach (Sasgen et al., 2010). This involves calculating the gravity field signal for each drainage basin and Ellesmere Island for a predefined mass change distributionhere, is based on the standard deviation of the mass variability of RACMO2/GR for the years January 2002 to December 2008 (74 months), assuming an elastic-compressible earth (Farrell, 1972;Wahr et al., 1998). Our forward model localizes the mass changes more in the coastal areas of GrIS and for Ellesmere Island in the ice covered areas compared to a uniform mass distribution, which leads to a less biased inverse solution. Changes in D are not prescribed in the forward model, and the uncertainty caused by this simplification is assessed by inverting the GRACE gravity fields with an alternative forward model based on ICESat data, which reflects both, changes in SMB and D. The errors associated with our inversion method are detailed in Appendix A.2.
The prescribed mass distribution for all drainage basins is simultaneously adjusted in magnitude such that the misfit between the predicted and GRACE-observed spatial gravity field changes, both identically filtered, is minimized in a least-squares sense. The seven drainage basins (henceforth, drainage basins A to G, Fig. 1) are defined according to the principal ice divides inferred from balance ice-velocities (Hardy et al., 2000). For Ellesmere Island, we adopt all areas that are ice covered in RACMO2/GR. The number of GrIS basins and their spatial separation is chosen to allow for recovering nearly independent estimates of mass change (Schrama and Wouters, 2011). The monthly GRACE coefficients are smoothed with an isotropic (i.e. degree-dependent) spatial averaging function (Sasgen et al., 2006), which optimizes the trade-off between spatial resolution -typically between 400 and 500 km -and noise with respect to the expected signals. The monthly coefficients of degree 2 and order 0 from GRACE are not reliable and are substituted by an estimate from Satellite Laser Ranging (Cheng and Tapley, 2004). For further regionalization, the weight of low degree and order coefficients are reduced according to the correlation of the predicted GrIS mass signal and its superposition with GIA in response to the history Laurentide ice sheet (Fig. A.1). Our inversion scheme for the GRACE solutions is validated with a simulated time series of monthly gravity-field variations based on cumulative RACMO2/GR mass anomalies for the years 2002 to 2008, w.r.t 1961 to 1990. To derive the long-term trends, the annual-oscillating component is removed from the monthly time series of GRACE mass anomalies, then, the time series are expanded into the first three Legendre polynomials. The advantage of a fit with Legendre polynomials is that the resulting trend and acceleration time series are orthogonal, providing independent estimates that are insensitive to a small change in the observation interval. The statistic significance of each term is determined using a Student ttest, estimating errors from the regression residuals and assuming a confidence limit of 5%.

A.2. GRACE error budget
Formal, as well as calibrated errors provided with the GRACE solutions of Stokes potential coefficients underestimate the true uncertainties, as they do not account for uncertainties associated with the de-aliasing models of atmospheric and oceanic shortterm variability. Therefore, we derive an empirical error estimate from the difference between two largely independent sets of GRACE coefficients, the CSR RL04 (Bettadpur, 2007) and GFZ RL04 (Flechtner, 2007). Although both GRACE releases rely on similar processing standards and background models, differences in the details of the gravity field determination cause deviations of the Table A1 Uncertainties of GRACE-inferred mass change for monthly solutions and linear trends for the seven major drainage basins (Fig. 1). Total uncertainties (Total) constitute the empirical errors in the GRACE gravity field solutions estimated from the difference between CSR RL04 and GFZ RL04 (Rel.), the filtering and inversion procedure (Filt./Inv.), the a priori assumption on the mass distribution (Forw. Model), and, in the case of the linear trend and acceleration terms, of formal regression errors (Regr.). Also provided is the range of plausible GIA corrections (GIA).  Filter response function applied to the monthly GRACE gravity field coefficients. Damping in the lower spectral part (medium shaded) accounts for the degree-correlation between the predicted GrIS mass change signal and its superpositions with the Laurentide GIA signal, whereas damping in the medium spectral part (light shaded) follows the signal-to-noise ratio in the GRACE coefficients. The upper limit of cut-off degree and order 60 (dark shaded) is governed by availability of the GRACE coefficients in CSR RL04.  GrIS (w.r.t. 1961GrIS (w.r.t. -1990, for the seven major drainage basins (Fig. 1). The starting date for the availability of GRACE monthly solutions (April 2002) is indicated by the vertical grey line. gravity field solution, and the inferred mass-change estimate (Table A.1). In addition, we estimate the uncertainty associated with our filtering and inversion scheme by performing our calculation on simulated (error-free) GRACE solutions. We also test the sensitivity of our estimates on the predefined mass change distribution (here, based on RACMO2/GR), by inverting the actual GRACE solutions with an alternative forward model based on mass trends inferred from ICE-Sat elevation changes (Fig. A.3). The uncertainty introduced by GIA is discussed in the following section. It should be stated that the applied inversion scheme and its accuracy stated in Table A.1 reflects our current understanding of GRACE processing, although alternative methods using e.g. Mascon (Luthcke et al., 2008) or optimal basin averaging functions (Velicogna, 2009) may perform similarly well.

A.3. Correcting for glacial-isostatic adjustment
In the GRACE trends over Greenland, GIA signature require removal, that are associated with the history of GrIS itself, as well as with the recession of the Laurentide ice sheet in North America, following the last glacial maximum. For the Laurentide GIA signal the GRACE trends can serve as a new constraint (e.g Tamisiea et al., 2007), augmenting more traditional measurements like GPS and tide-gauges (Davis and Mitrovica, 1996;Braun et al., 2008), terrestrial gravity, satellite radar altimetry (Lee et al., 2008), traditional leveling or sea-level indicators (e.g. Tushingham and Peltier, 1992;Lambeck et al., 1998).
To minimize the influence of GIA on the GrIS mass trends from GRACE, and to assess their uncertainty on the GrIS mass change derived from GRACE, we proceed as follows. First, GIA induced by the GrIS and Laurentide ice sheet are modeled independently; for the GrIS the glacial reconstruction of Fleming and Lambeck (2004) is adopted, for North America, the Lauren-tide ice sheet component of ICE-5G is isolated. Following the approach of Paulson et al. (2007), the total glacial load of the Laurentide ice sheet component of ICE-5G is adjusted in order to minimize the misfit between the predicted and GRACE-inferred trends in gravity fields over the Hudson Bay, North America (adjusted version is denoted as ICE-5G n ). In advance to the adjustment, the GIA prediction is filtered identically to the GRACE trends, and the GRACE trends are corrected for variations in total water storage estimated with the hydrology model WGHM (Döll et al., 2003). It should be mentioned, that the adjustment of GIA over North America is optimal only with respect to the band-pass filtered GRACE coefficients.
To estimate the uncertainty of our GIA correction, the scaling procedure is repeated for a variety of viscosity models ranging between 2 Â 10 20 Pa s and 8 Â 10 20 Pa s for the upper mantle and 5 Â 10 21 Pa s and 4 Â 10 22 Pa s for the lower mantle, for thickness of the lithosphere of 80, 100 and 120 km, as well as for the glacial history ANU (Lambeck and Chappell, 2001) developed at the Australian National University (adjusted version is denoted as ANU n ). In terms of the GIA related to the glacial history of the GrIS, GRACE provides no useful additional constraint, because the GIA signal directly underlays that of ongoing ice-mass changes. We estimate the related uncertainty by varying upper-and lower-mantle viscosities in the range stated above. For ICE-5G n and the glacial history of GrIS, upper-and lower mantle viscosities of 6 Â 10 20 Pa s and 2 Â 10 22 Pa s are adopted as reference values; for the Laurentide component of ANU* we employ 5 Â 10 20 Pa s and 1 Â 10 22 .
A.4. Regional climate modelling and InSAR satellite observations The attribution of GRACE-inferred mass losses to anomalies in SMB or D requires the definition of a reference period during which the Greenland ice sheet was approximately in balance, i.e. SMB E D. The 30-year time interval from 1961 to 1990 may serve as such a period, as it provides a meaningful multi-decadal climatology that excludes large SMB variations and imbalance starting in the year 1996. Therefore, for all basins, the anomalies in SMB and D are calculated separately with respect to the means for the years . The deviation of the SMB and D mean fields  provides an indication of the systematic uncertainties in the SMB and D data sets, and the approach of combining their anomalies in SMB-D; the difference of mean D and mean SMB is o8 Gt/yr for all basins, except basin B (mean D is 17 Gt/yr larger), for the GrIS the uncertainty is estimated to 7 14%. In the case of D, estimates are not continuous prior to 1996 (Rignot et al., 2008;Rignot and Kanagaratnam, 2006). To interpolate D between years with observations we correlated D with averaged runoff for the preceding four years and current year. This is similar to a previous approach (Rignot et al., 2008) except here we use runoff rather than SMB and we use estimates from RACMO2/GR. For the ice covered area in RACMO2/GR, that is unsurveyed by ( $ 12%; Rignot and Kanagaratnam, 2006), the D anomalies are assumed to be 0. In addition, D data are extrapolated to areas of basin F and basin G, which were problematic to survey with InSAR (basin 38; Rignot and Kanagaratnam, 2006), as follows: first, the mean D (years 1961-1990) is scaled by the factor ''total area''/''surveyed area'', resulting in 1.1 for basin F and 1.5 for basin G. Here, ''total area'' refers to the spatial coverage of the InSAR survey provided in (Rignot and Kanagaratnam, 2006), covering 99% (F) and and 92% (G) of the ice covered area in RACMO2/GR. Then, a regional ''imbalance'' factor for each year between 2002 and 2008 is applied to the mean D of the unsurveyed area; for basin G, the imbalance factor for drainage basin number 25 in the InSAR data set of Rignot and Kanagaratnam (2006) is adopted, whereas for basin F the factor is based on the average of the neighbouring basins number 21, 22 and 23 (Rignot and Kanagaratnam, 2006). Then, the extrapolated D minus the mean D for 1961-1990 of the unsurveyed region is calculated, and added to the surveyed D anomalies. For all basins, due to the unavailability of data, D is assumed to remain unchanged in 2010, compared to 2009. The time series of the SMB components and D evaluated for each basin from 1958 to 2010 are shown in Fig. A.2. For Ellesmere Island, D anomalies are assumed to be zero due to lack of recent data. Since the 1999-2003 mean D is estimated to be 5 72 Gt/yr (Gardner et al., 2011, Suppl. Information), the error introduced by our assumptions is small compared to the GRACE-estimated trend.

A.5. ICESat observations
NASA's ICESat carried the Geoscience Laser Altimeter System (GLAS) instrument (Zwally et al., 2010), which performed ice sheet elevation measurements in the period 2002-2009. Here, ICESat altimetry measurements from the GrIS and the surrounding glaciers and ice caps are considered in the elevation change analysis. The rejection of problematic measurements reduced the number of measurements by approximately 13%. The total volume change of the GrIS from October 2003 to October 2009 is found by fitting a smooth surface covering the entire ice sheet, through the ICESat derived rate of surface-elevation change. This geometric rate of ice sheet loss is related to mass, by correcting for volume losses not contributing to a mass loss and applying the appropriate conversion density to the corrected volume (Sørensen et al., 2011). The correction terms are firn compaction, vertical bedrock movement, and ICESat intercampaign elevation bias. The derived rate of elevation changes, corrected for these terms are shown in Fig. A.3. The firn compaction term is by far the largest of the corrections applied. Both the surface density model, and the firn compaction model are forced by surface temperature, accumulation and run-off from a high resolution regional climate model (HIRHAM5 RCM). The HIRHAM5 RCM simulation covering the period from 1989 to 2009, and are driven by the ERA-Interim re-analysis product, at the lateral boundaries.
An independent error analysis has been conducted for each of the drainage basins to estimate the regional error associated with the derived mass balance of the GrIS. This error analysis consists of six error sources, which are assumed to be independent; ICESat volume estimates, ICESat intercampaign bias, vertical bedrock movement, firn compaction, neglecting basal melt and using simplified assumptions on ice dynamics (Sørensen et al., 2011). A bootstrap procedure has been applied to estimate the regional error from the sampling and interpolation of the measured surface elevation change in the ICESat tracks. Both the firn compaction modelling and the vertical bedrock correction are associated with a volume error estimate. The error on the corresponding mass change has been estimated by regionally adding a perturbation of the volume error and applying the density model. By neglecting basal melt and using simplified assumptions on ice dynamics-possible build up of ice above the equilibrium line altitude (ELA), an error may be present in areas above the ELA, which are associated with an elevation increase. This error contribution is determined by the difference between the applied density of firn and the density of ice, and the magnitude of the ice dynamics are determined by a simple flow model (Sørensen et al., 2011).