Improved GRACE regional mass balance estimates of the Greenland ice sheet cross-validated with the input–output method

. In this study, we use satellite gravimetry data from the Gravity Recovery and Climate Experiment (GRACE) to estimate regional mass change of the Greenland ice sheet (GrIS) and neighboring glaciated regions using a least squares inversion approach. We also consider results from the input–output method (IOM). The IOM quantiﬁes the difference between the mass input and output of the GrIS by studying the surface mass balance (SMB) and the ice discharge ( D ). We use the Regional Atmospheric Climate Model version 2.3 (RACMO2.3) to model the SMB and derive the ice discharge from 12 years of high-precision ice velocity and thickness surveys. We use a simulation model to quantify and correct for GRACE approximation errors in mass change between different subregions of the GrIS, and investigate the reliability of pre-1990s ice discharge estimates, which are based on the modeled runoff. We ﬁnd that the difference between the IOM and our improved GRACE mass change estimates is reduced in terms of the long-term mass change when using a reference discharge derived from runoff estimates in several sub-areas. In most regions our GRACE and IOM solutions are consistent with other studies, but differences remain in the northwestern GrIS. We validate the GRACE mass balance in that region by considering several different GIA models and mass change estimates derived from data obtained by the Ice, Cloud and land Elevation Satellite (ICESat). We con-clude that the approximated mass balance between GRACE and IOM is consistent in most GrIS regions. The difference in the northwest is likely due to underestimated uncertainties in the IOM solutions.


Introduction
During the last decade, the ice mass loss from the Greenland ice sheet (GrIS) has become one of the most significant mass change events on Earth. Because of its ongoing and potentially large future contribution to sea level rise, it is critical to understand the mass balance of the GrIS in detail. As a result of increasing runoff and solid ice discharge (D), the GrIS has been experiencing a considerable and increasing mass loss since the mid-1990s (Hanna et al., 2005;Rignot and Kanagaratnam, 2006;. The changes in mass loss rates are due to different processes. For instance, mass loss acceleration in the northwestern GrIS is linked to the rapidly increasing discharge in this region (Enderlin et al., 2014;Andersen et al., 2015), while in the southeast the increase in mass loss rate after 2003 is mainly due to enhanced melting and less snowfall (Noël et al., 2015).
To quantify the recent changes in the GrIS mass balance, three methods are used: satellite altimetry, satellite gravimetry and the input-output method (IOM) (Andersen et al., 2015;Colgan et al., 2015b;Sasgen et al., 2012;Shepherd et al., 2012;Velicogna et al., 2014;Wouters et al., 2013). We will concentrate on the latter two methods in this study, using results from satellite altimetry for validation purposes.
The IOM is used to evaluate the difference between mass input and output for a certain region. It considers two major mass change entities: surface mass balance (SMB) and solid ice discharge. SMB is commonly estimated using climate models Fettweis, 2007;Tedesco et al., 2013;van Angelen et al., 2012), whereas D can be estimated from combined measurements of ice velocities and ice Published by Copernicus Publications on behalf of the European Geosciences Union. 896 Z. Xu et al.: Improved GRACE regional mass balance estimates of the Greenland ice sheet thicknesses, e.g., Rignot and Kanagaratnam (2006), Enderlin et al. (2014) and Andersen et al. (2015). To reduce the uncertainties in the mass change of SMB and D, the SMB and D from 1961 to 1990 are sometimes used as a reference when applying the IOM method Sasgen et al., 2012). However, introducing SMB and D as reference may introduce new uncertainties in the IOM. We will discuss the details of the IOM, as well as the uncertainties of the reference SMB and D, in Sect. 2.
The satellite gravity observations from GRACE (Gravity Recovery and Climate Experiment) provide snapshots of the global gravity field at monthly time intervals, which can be converted to mass variations. Mass variation solutions of a given area that are obtained from GRACE observations are, however, influenced by measurement noise and leakage of signals caused by mass change in neighboring areas. Furthermore, the GRACE monthly gravity fields contain northsouth oriented stripes as a result of measurement noise and mis-modeled high-frequency signal aliasing. Therefore, to estimate the mass balance for GrIS subregions from GRACE data, we apply the least squares inversion method (Schrama and Wouters, 2011) in this study with an improved approach (Xu et al., 2015). However, as shown by Bonin and Chambers (2013) in a simulation study, the least squares inversion method introduces additional approximation errors.
Previous studies have compared regional GrIS mass change from different independent methods. In Sasgen et al. (2012), the mass balance in seven major GrIS drainage areas was derived from the IOM and GRACE data using a forward modeling approach developed by Sasgen et al. (2010). When separating out the IOM components and comparing them with the seasonal variability in the derived GRACE solution, the relative contributions of SMB and D to the annual mass balances were revealed. In the northwestern GrIS, important differences between IOM and GRACE were noted, which were attributed to the uncertainty in the regional discharge component in this area, where detailed surveys of ice thickness are lacking. The comparison between two approaches shows a mass loss difference of 24 ± 13 Gt yr −1 in this region, and as a result the uncertainty in the regional mass balance estimate is estimated at ∼ 46 %. However, using new discharge estimates and the corresponding IOM regional mass change in the northwestern GrIS, Andersen et al. (2015) found that the difference between GRACE and IOM mass loss estimates fell within the combined uncertainty range. However, using the least squares based inversion approach of Schrama and Wouters (2011), we find that the mass change differences between results obtained from GRACE and the IOM in the southern GrIS increase and cannot be explained by the assumed uncertainties. An example of the regional differences between the GRACE data and the IOM solution can be seen in Fig. A1. The details of this difference will be discussed in Sect. 4. In this study, we aim to investigate the two aforementioned sources of uncertainties in GRACE and IOM mass balance  Zwally (2012). A mascon with the same digits refers to a region belonging to the same drainage system. The "a" and "b" terms indicate the GrIS margin (< 2000 m) and GrIS interior (≥ 2000 m), respectively. There are 16 GrIS mascons and 4 neighbouring Arctic mascons. The locations of the three largest discharge outlets are marked with a star, i.e., Jakobshavn (green), Kangerdlugssuaq (red) and Koge Bugt (blue) glaciers. The glacier surface area is defined in the RACMO2.3 model. estimations: (i) we present a way to reduce the error from the inversion approach and (ii) we investigate different discharge estimates. We then evaluate our results by comparing the GRACE and IOM estimates both with each other and with published estimates from satellite altimetry. The GrIS drainage system (DS) definition of Zwally (2012) is employed here to investigate the mass balance in GrIS subregions. This definition divides the whole GrIS into eight major drainage areas, and each drainage area is further separated by the 2000 m elevation contour line, creating interior and coastal regions for each drainage area. This GrIS DS definition is employed by several other studies (Andersen et al., 2015;Barletta et al., 2013;Colgan et al., 2015a;Luthcke et al., 2013;Sasgen et al., 2012). The regional GrIS mass change estimated with GRACE are influenced by mass change from areas outside the ice sheet, i.e., from Ellesmere Island, Baffin Island, Iceland and Svalbard (EBIS) (Wouters et al., 2008). Therefore, we include four additional DSs to reduce the leakage from these regions. The overall mascon definition used in this study is shown in Fig. 1.
The main topic of this study is to provide improved GrIS regional mass balance estimates from GRACE and the IOM. We show that the improved GRACE solution reduces the regional differences between two mass change estimates, especially in the southeastern GrIS region. Furthermore, we compare the GRACE solution with the IOM, which employs different reference discharge estimates, showing that the uncertainties in the reference discharge can result in an under-The Cryosphere, 10, 895-912, 2016 www.the-cryosphere.net/10/895/2016/ estimated mass loss rate in the IOM regional solution, in particular in the northwestern GrIS region. In Sect. 2, we present SMB mass change from a recently improved regional atmospheric climate model (RACMO2.3) (Noël et al., 2015) and discharge estimates of Enderlin et al. (2014), which are based on a near-complete survey of the ice thickness and velocity of Greenland marine-terminating glaciers. In Sect. 3, we introduce the least squares inversion approach. In Sect. 4, we start by investigating different methods to calculate mass change in GrIS drainage areas using the modeled SMB and D estimates. Subsequently, we identify the approximation errors in regional mass change estimates from GRACE data, followed by a comparison of mass change estimates from GRACE and the IOM and a discussion of the remaining differences. The conclusions and recommendations are given in Sect. 5.

SMB and D
For the GrIS, precipitation (P ) in the form of snowfall is the main contribution to the mass input, while mass loss is a combination of sublimation (S), meltwater runoff (R) and D. The SMB equals to P − S − R, and subtracting D from SMB yields the total mass balance (TMB). In this study, we use the Regional Atmospheric Climate Model, version 3 (RACMO2.3), to model the SMB of the GrIS. RACMO2.3 Van Angelen et al., 2012;) is developed and maintained at the Royal Netherlands Meteorological Institute (KNMI) and has been adapted for the polar regions at the Institute for Marine and Atmospheric Research, Utrecht University (UU/IMAU). RACMO2.3 model output is currently available at ∼ 0.1 • spatial resolution for January 1958 to December 2014. The differences between a previous model version (RACMO2.1) and other SMB models are discussed by Vernon et al. (2013). For RACMO2.3, we assume 18 % uncertainties for the P and R components in each grid cell. Assuming both components to be independent, the uncertainty of the SMB is the quadratic sum of uncertainties of P and R. The magnitude of S is small and its absolute uncertainty is negligible compared to those in P and R. The RACMO2.3 model also provides estimates of the SMB in the peripheral glacier areas, which we have included in this study.
D estimates from Enderlin et al. (2014) (hereafter Enderlin-14, with the associated discharge estimates D-14) are used in this study. In Enderlin-14, the ice thickness of 178 glaciers is estimated from the difference in ice surface elevations from repeat digital elevation models and bed elevations from NASA's Operation IceBridge airborne ice-penetrating radar data. The ice surface velocity is obtained from tracking the movement of surface features that are visible in repeat Landsat 7 Enhanced Thematic Map-per Plus (ETM+) and Advanced Spaceborne Thermal and Reflectance Radiometer (ASTER) images. For glaciers with thickness transects perpendicular to ice flow (i.e., flux gates), the ice flux is estimated by summing the product of the ice thickness and surface speed across the glacier width. For glaciers with no thickness estimates, we use empirical scaling factors for the ice flux, as derived in Enderlin et al. (2014). Because the ice fluxes are calculated within 5 km of the estimated grounding line locations, the SMB gain or loss between the observations and the grounding lines will be small and the ice discharge is estimated directly from the fluxes (Enderlin et al., 2014). The resulting estimation of discharge uncertainty of 1-5 % for each glacier is smaller than in previous studies, e.g., Rignot et al. (2008) (hereafter Rignot-08, and the associated estimates are denoted by D-08), which relied on interior ice thickness estimates that were assumed to be constant in time.

Cumulative TMB anomaly
For both the whole GrIS and a complete drainage area from ice sheet maximum height to the coast, the total mass balance is (1) In this study, we further separate each GrIS drainage area in a downstream (I) and upstream (II) region, separated by the 2000 m surface elevation contour line. Thus, for the subdivided regions Eq. (1) becomes where and in which F II refers to the ice flux across the 2000 m elevation contour, and F I refers to the ice flow across the flux gate. Note that F II is canceled if the study area includes both the regions below and above the 2000 m contour but F II has to be considered when the upstream and downstream regions are considered separately. As described above, we assume that SMB changes downstream of the Enderlin-14 flux gates are negligible and that F I = D. In order to fit the temporal resolution of the modeled SMB data, we interpolate the yearly D on a monthly basis. Significant seasonal variations in ice velocity have been observed along Greenland's marine-terminating outlet glaciers (Moon et al., 2014). However, since we focus mostly on long-term mass change in this study, monthly variations in D should have a negligible influence on our analysis, and we assume that D is approximately constant throughout the year. Contrary to the total mass represented in the GRACE data, the www.the-cryosphere.net/10/895/2016/ The Cryosphere, 10, 895-912, 2016 SMB, D and TMB are estimates of rates of mass change (i.e., mass flux) in "Gt per month" or "Gt per year". Hence, in order to compare with GRACE, one has to integrate the SMB and D from a certain month (or year) onwards, which yields where TMB i is the cumulative mass change at month i in the IOM (unit is Gt) and the integration time period is from a certain initial month i 0 to month i. In a previous study of mass balance with the IOM, for which estimates of D were not available for certain regions (Rignot et al., 2008), the 1961 to 1990 reference SMB is used to approximate the missing regional D (Sasgen et al., 2012). Also, due to the uncertainties in the SMB model, accumulating the TMB over decades may also lead to unrealistic mass gains or losses . By removing the reference, the influence of large uncertainties and interannual variability in SMB and D can be reduced . For instance, the uncertainties due to model configurations could be similar in every monthly SMB estimate, and accumulating those over a long period may result in a large uncertainty. The choice of reference period is based on the assumption that the mass gain from the surface mass balance during that period is compensated by ice discharge, i.e., the GrIS was in balance during that period (no mass change).
For the reference period we defined the period for the integral in Eq. (5) to be from years 1961 to 1990. For the subsequent period the lower and upper bounds of the integral are 1991 and i, respectively. Since we assume the GrIS was in balance during the reference period, 1990 1961 (SMB t − D t ) dt = 0.
By removing the reference SMB and D (i.e., SMB 0 and D 0 ) Eq. (5) becomes where δSMB t = SMB t − SMB 0 and δD t = D t − D 0 . Note that SMB 0 and δSMB t are both rates of mass change (and similarly for the discharge). As explained before, when Eq. (6) is used to compute the mass balance for the regions below and above 2000 m separately, the ice flux across the 2000 m contour (F II ) has to be considered. Therefore we introduce two assumptions: (1) F II is constant over time, which means F II = F II 0 (F II 0 is the F II during the reference period), so that ( (1) is necessary since there is a lack of yearly measurements of ice velocity across the 2000 m contour. An estimate of decadal change by  suggests it is reasonable to assume a constant F II for the entire GrIS, except for a few glaciers, such as the Jakobshavn glacier where, after 2000, the F II may be higher than F II 0 . In Andersen et al. (2015), the mass balance of the interior GrIS (in their study defined as the ice sheet above the 1700 m elevation contour) was 41 ± 61 Gt yr −1 during the 1961-1990 reference period and in Colgan et al. (2015a) the ice flux across the 1700 m contour was estimated to be 54 ± 46 Gt yr −1 for the same time period, indicating the assumption of balance approximately holds within the given uncertainties.
Based on these two assumptions, we apply Eq. (6) for the interior and coastal GrIS regions, yielding and We quantify the combined uncertainties of assumptions (1) and (2) by comparing the results from Eq. (7) to the regional mass balance derived from GRACE by Wouters et al. (2008), as well as those derived from Ice, Cloud and land Elevation Satellite (ICESat) by Zwally et al. (2011), resulting in an uncertainty of ±15 Gt yr −1 for the entire interior GrIS. The regional uncertainties are summarized in Table A2. Note that for each region, the same uncertainty is applied to both the interior and coastal areas. For the whole drainage area, the uncertainties associated with assumption (1) and (2) will vanish, because these two assumptions are needed only when we separate the coastal and interior regions.

Post-processing GRACE data
In this study, we use the GRACE release 5 level 2 monthly spherical harmonics coefficients C lm and S lm (GSM) produced by the University of Texas Center for Space Research (CSR). The monthly coefficients are available from January 2003 to January 2014 and have a maximum spherical harmonic degree l = 60. We add C 10 , C 11 and S 11 coefficients to these GRACE solutions (related to the motion of the Earth's geocenter) obtained from GRACE data and independent oceanic and atmospheric models (Swenson et al., 2008). The geopotential flattening coefficients (C 20 ) from The Cryosphere, 10, 895-912, 2016 www.the-cryosphere.net/10/895/2016/ GRACE data are less accurate than those from the satellite laser ranging (SLR) (Chen et al., 2004). Consequently, we replace these coefficients with the ones obtained from SLR data by Cheng et al. (2013). The GRACE potential coefficients are averaged between January 2003 and January 2014 to obtain the reference field with respect to which monthly anomalies C lm and S lm are defined. GRACE observations of mass change within a subregion of the GrIS are affected by the mass change in neighboring areas, a phenomenon known as leakage (Wahr et al., 1998). Furthermore, GRACE data should be corrected for known oceanic and atmospheric mass motions, continental hydrology and glacial isostatic adjustment (GIA). The oceanic and atmospheric mass change is already removed from the coefficients provided by CSR. The Global Land Data Assimilation System (GLDAS) model (Rodell et al., 2004), specifically the 1 • spatial resolution version 2 monthly data obtained from Goddard Earth Sciences Data and Information Services Center, is employed here to simulate the continental hydrology. Before removing the output of this hydrology model from the GRACE monthly coefficients, permafrost and glaciated regions are excluded.
We correct for the GIA effect in the GRACE data of the GrIS by using the output of the model by Paulson et al. (2007), which is based on the ICE-5G ice loading history and the VM2 Earth model (Peltier, 2004). Hereafter, we refer to this model by . In addition to this approach, three GIA models with a total of 11 variations are employed based on different ice histories and viscosity models to determine the uncertainty in the GIA correction. For instance, the models of van der Wal et al. (2013) include 3-D changes in viscosity and the model of Simpson et al. (2009) uses a different ice loading history. A summary of the GIA models used in this study is given in Table A3. An isotropic Gaussian filter is employed to reduce the noise in GRACE data (Wahr et al., 1998), with a half width of r 1/2 = 300 km.

Inversion of the regional mass balance
To estimate the regional mass balance in separate GrIS drainage area, we use a constrained least squares inversion approach (Bonin and Chambers, 2013;Schrama and Wouters, 2011): The vector y contains the monthly GRACE data. To compute the influence functions in the design matrix H we assume a layer of water with unit height uniformly distributed over the mascon and subsequently express the mass change in spherical harmonic coefficients up to d/o 60, similar to the GRACE data. The vectorx represents the scale factors for the unit mass change in each drainage area that we aim to find. P is the a priori covariance matrix of the mass change in each mascon. When assuming that the mass change in each equally weighted mascon is independent, P = λI, with λ the a priori variance of the regional mass change. In our previous study, we demonstrated that three different prior variances for the GrIS regions below and above 2000 m, as well as for the surrounding Arctic regions, improved the recovery of regional mass change (Xu et al., 2015). Using a simulation model based on the IOM (see Sect. A3), optimal regional constraints were determined. We find λ a = 13 m 2 for coastal mascons, λ b = 0.1 m 2 for inland mascons and λ EBIS = 11 m 2 for the nearby surrounding EBIS regions. It has to be highlighted that the GRACE method used in this study is not directly constrained by the mass change derived from the RACMO2.3 and the ice discharge estimates. Instead, from the RACMO2.3 and the ice discharge estimates, we derived the variance of the mass change in each region as the constraint. Moreover, in our previous study (Xu et al., 2015), we showed that the derived mass change in the coastal regions is largely obtained from the GRACE data, with only a small influence from the constraints. Conversely, for the interior regions where the mass change rates are small, they are mainly determined from the constraints.

Reference SMB and D
In this study, the total error in SMB 0 , hereafter σ SMB 0 , is a result of the systematic error caused by the assumption of a reference period and the averaging within the chosen reference period. Both components will be explained hereafter. The systematic error is the uncertainty in the SMB derived from model output, whereas the averaging error is related to the variability of the reference SMB 0 during 1961-1990. To quantify the latter, we perform a Monte Carlo simulation to evaluate the standard deviations of the SMB 0 that result from using different combinations of 20-year averages of the SMB, following . The sampled combinations are randomly chosen from the months between 1961 and 1990. For RACMO2.3, we find 20 Gt yr −1 averaging errors in σ SMB 0 . The SMB 0 from RACMO2.3 yields 403 Gt yr −1 ; hence the systematic error is approximately 73 Gt yr −1 (considering 18 % uncertainty in RACMO2.3). If we assume both errors are independent, then σ SMB 0 = 75 Gt yr −1 .
We also investigate the uncertainties of the 1961-1990 reference discharge. In this study we employ D-14 as the D estimate in the IOM. However, the D-14 time series starts from the year 2000, when the GrIS already was significantly out of balance. To retrieve D 0 for D-14, we use D 0 = 413 Gt yr −1 in 1996 from D-08 (D 0 -08) for the entire GrIS and assume that the regional D changes from 1990 to 2000 in D-08 are proportional to the changes in D-14 in each region, i.e., that D-14 and D-08 are linearly correlated. The details of the interpolation of regional D 0 are given in Sect. A1. Note that the averaging error in D 0 is minimized via an iteration process, the details of which can be found in Rignot et al. (2008). Due to the lack of ice thickness information before 2000, the reference D 0 in Rignot-08 has high uncertainty, especially in the northwest of the GrIS. Another way to obtain historic discharge estimates is by using the presumed correlation between discharge and SMB or runoff (Rignot et al., 2008;Sasgen et al., 2012). This approach assumes that the anomaly of the discharge with respect to a reference SMB (δD = SMB 0 − D) is correlated with a reference runoff (δR = R − R 0 ), which is based on the anomaly of the 5-year averaged runoff. Note that the lagging correlation is discussed in Bamber et al. (2012) and Box and Colgan (2013). In this study we choose to use the runoff output from the RACMO2.3 model. We consider three estimates of D, i.e., by Rignot-08, Enderlin-14 and Andersen et al. (2015), which use different types of measurements of the ice thickness and flux velocity changes, integration areas (areas between the flux gate and the grounding line), SMB and ice storage corrections. Additionally, they differ in whether the peripheral glaciers are included or not. In this study we provide runoff-based estimates for D 0 only for those GrIS drainage areas where the correlation coefficients between δD and δR are equal or higher than 0.7 (Fig. 2), i.e., for DS1, DS3, DS7 and DS8. For the entire GrIS, we obtain a high correlation (R 2 = ∼ 0.86), similar to the correlation found by Rignot et al. (2008). However, the high correlation cannot be obtained in all the GrIS drainage areas; the regional correlation coefficients vary from 0.19 to 0.94. In DS7 and DS8, the discharge anomaly is obviously correlated with the runoff anomaly (R 2 > 0.9), while in other regions (i.e., in DS2, DS4, DS5 and DS6), the correlation is low (R 2 < 0.5). In DS3a, when we consider only the D estimates from Enderlin et al. (2014) and Andersen et al. (2015), the correlation increases to R 2 = 0.72. Note that, the regions with high correlation are also those that have a large fraction of marineterminating glaciers. We derive the linear relation between δD and δR for eight major GrIS DSs and calculate the regional annual δD from 1961 to 2013 using this linear relation.
Hereafter, the regional cumulative discharge anomaly (δD), which is derived from the RACMO2.3 runoff, is denoted as δD R , while δD-08 and δD-14 refer to δD based on Rignot-08 and Enderlin-14, respectively. We compare δD R , δD-08 and δD-14 in Fig. 3 for the time interval 2000 to 2007, which is common to both δD-08 and δD-14. In DS7 we find a correlation of R 2 = 0.94. In that region, δD-08 and δD-14 are similar, i.e., 20.1 ± 1.9 and 17.6 ± 2.2 Gt yr −1 , respectively, but δD R is 8.9 ± 4.7 Gt yr −1 . The difference between the runoff-derived and flux gate D estimates may indicate that the reference D 0 for this region should be ∼ 9 Gt yr −1 lower than D 0 estimated by Rignot-08. A similar difference can be seen in DS4 where we obtain 36.2 ± 2.5 Gt yr −1 for δD-14 and 37.9 ± 2.8 Gt yr −1 for δD-08, but δD R is 8.4 ± 3.3 Gt yr −1 . However, in DS4, δD R is probably not reliable, since the runoff-to-discharge correlation is weak in this region (R 2 = 0.38). For the entire GrIS, the reference discharge D 0 is 427 ± 30 Gt yr −1 for δD-08 and 414 ± 44 Gt yr −1 for δD-14. becomes 410 ± 37 Gt yr −1 and all three versions of reference discharge agree within their respective uncertainties.
To evaluate the SMB 0 and D 0 used in this study, we compare the IOM regional mass balance in eight major drainage areas (interior and coastal regions are combined) and apply both Eqs. (5) and (6). The latter equation relies on the determination of the SMB 0 and D 0 while Eq. (5) does not, so the comparison can provide an indication about the reliability of the SMB 0 and D 0 for some drainage area. For the application of Eq. (6) we use two methods. Method 2 uses δD-14 while method 3 uses δD R in DS1, DS3, DS7 and DS8. As can be seen in Fig. 3, the three methods agree for the whole GrIS, as well as for most of the drainage areas, within their uncertainties. In DS4, DS7 and DS8, however, methods 1 and 2 are significantly different, which may be caused by underestimated cumulative errors in Eq. (5) or a less accurate reference surface mass balance SMB 0 and reference discharge D 0 . This is further discussed in Sect. 4.3.

Approximation errors
In the solution ofx, two types of errors occur: (a) systematic errors caused by measurement errors propagated through the least squares approach and (b) an additional error that is introduced when applying Eq. (9). For the type (b) error, Bonin and Chambers (2013) show that the use of Eq. (9) results in a noticeable difference between the approximationx and the "truth" (a GrIS mass change simulation), especially in GrIS subregions, which we recognize as an error source; see also the discussion in Tiwari et al. (2009). Hereafter the type (b) error is denoted as "approximation error" or ε. We estimate ε by using simulations of GrIS for x, following Bonin and Chambers (2013), so that the approximation error becomes ε = x −x. The simulated regional mass change on the mascons is x = [x 1 , x 2 , x 3 , . . . , x n ], where n is the total number of mascons. We will show that there is a relation betweenx and x, which can be used to correct for the approximation errors.
The simulation model y = f (θ λ) is based on the 10-year linear trend (2003 to 2012) of mass change of SMB and D estimates (see Sect. A3), with uncertainties of the simulation model written as σ (θ λ). We employ a Monte Carlo approach to simulate a sample of 1000 randomly distributed observations, according to y l = y + k l σ with k l = k l (θ λ) a vector of random scaling factors that vary from −1 to 1 and index l www.the-cryosphere.net/10/895/2016/ The Cryosphere, 10, 895-912, 2016 running from 1 to 1000. It is important to note here that we assume that measurement errors do not exist (i.e., the simulation model is assumed to be reality). In addition, we assume that the generated samples in the simulation (σ ) are normally distributed.
Next we apply Eq. (9) to yield approximated regional mass changex = [x m ], in which m is the index of the mascons (see Fig. 1). From the simulation we can derive the real regional mass change rate x = [x m ]. As mentioned above, the difference betweenx and x equals the approximation error. In Fig. 4 we show that the x i are linearly correlated withx m . By applying this correlation to the approximations derived from GRACE data, one can reduce the approximation errors in the GRACE-based regional mass balance approximations.
The simulated trends in regional mass changes x and the corresponding approximationx are shown in Fig. 4. It can be noticed that the approximations are strongly correlated over time with the simulation in the coastal regions, having an average correlation coefficient R 2 of ∼ 0.9. This means that the approximated regional solutions are close to the simulation. The correlation in region DS1a is weaker (∼ 0.6), which suggests that the approximation for region DS1a is influenced more by mass change in neighboring regions such as region DS8a. In the simulation, the inter-region correlation between DS1a and DS8a is ∼ −0.1, while in the approximations, the correlation rises to ∼ −0.5. By compari-son, another neighbor of DS8a, DS7a, has a very weak interregion correlation with DS8a (∼ 0.04), both in the simulation and in the approximation. The inter-region correlation errors are a systematic result of the least squares inversion (Bonin and Chambers, 2013;Schrama and Wouters, 2011). Previous work shows that the regional approximation errors can be reduced by specifying constraints for the GrIS coastal and inland regions. However, in this study all the sub-drainage areas within the coastal region are constrained by the same prior variance, resulting in the large remaining correlation between DS1a and DS8a.
For the coastal regions, there is a linear relationship between the simulations x and the approximationx, as can be seen in Fig. 4. We fit this relationship by x = α 1x + α 0 , with α 1 and α 0 given in Table A1. The linear relationship between the simulated and the approximated regional mass change rates is found to be stable; even when the simulation uncertainties are multiplied with a factor of 5 (light green marks in Fig. 4), the regression parameters (α 1 and α 0 ) vary by less than ∼ 1 % for the coastal mascons. Therefore, it is reasonable to assume α 1 and α 0 should also reflect the relationship between reality and the approximation as derived from GRACE observations. When the vector of observations y becomes the GRACE observations, we can derive an improved GRACE regional solution by applying the linear relationship to the corresponding approximationx. We will show in Sect. 4.3 that this correction yields better agreement between GRACE and the IOM. Contrary to the coastal regions, the linear relation between x andx is weak in the interior regions, where the mean correlation coefficient is ∼ 0.2. This may be due to the fact that interior regions show smaller mass change rates than the coastal regions. For simulations created within a 1σ range, the highest correlation coefficient is only 0.47 for DS7b. The strong constraint used for these regions, i.e., a priori variance of 0.1 m 2 , may cause the approximation to be more strongly determined by this constraint than in the simulation. However, even if we apply a weaker constraint, i.e., λ = 10 6 , the correlation coefficients between x andx in these regions remain below 0.5. This means that correcting the approximation errors using the same method as for the coastal regions may create larger uncertainties. Following Bonin and Chambers (2015), we choose to include the approximation errors in the uncertainties, but only for the interior regions. The uncertainties are shown in Table A2.

Results and discussions
We compare the regional mass change rate from GRACE with the IOM (Fig. 5) before and after applying the approximation error correction to GRACE, and with different discharge estimations implemented in the IOM, separately for coastal and interior regions. Note that in this figure, the time interval is January 2003 to December 2013; we extrapolate the 2013 ice discharge from Enderlin-14 D estimates. For the coastal regions, we find that the correction of the approximation errors in the GRACE solutions adjusts the mass distributions between adjacent mascons. For instance, the corrected mass loss rate in DS3a increases by 10 Gt yr −1 while it reduces the mass loss rate in the adjacent region DS4a by 15 Gt yr −1 . In DS6a, correcting for the approximation error causes a mass loss increase of 13 Gt yr −1 . One may notice that the corrected GRACE mass loss rates exceed the uncertainty range of the mass loss rates before correction; e.g., in DS1a and DS3a. This is because the uncertainty before our correction is estimated without considering the approximation error.
We only consider TMB from the IOM in order to reduce the influence of the individual uncertainties in SMB and D. We obtain two IOM solutions, using the reference D 0 by Rignot-08 (method 2) and the interpolated discharges based on the RACMO2.3 runoff (method 3). In DS1a and DS3a, we obtain lower discharge changes rate from method 3 than from method 2. In DS7a, which includes Jakobshavn glacier, method 3 results in smaller mass change than method 2. Figure 5 shows that the agreement between GRACE and the IOM improves after correcting the GRACE approximation errors and applying the runoff-based discharge estimations in DS3a, DS5a, DS6a and DS7a. The difference between GRACE and IOM estimates is also reduced in DS1a and DS2a, where the remaining difference falls within the uncertainty margins. The corrected GRACE solution in DS4a is only ∼ 3 Gt yr −1 lower than the IOM solution, while it was ∼ 10 Gt yr −1 higher before correcting for the approximation www.the-cryosphere.net/10/895/2016/ The Cryosphere, 10, 895-912, 2016 error. However, regardless of whether the approximation errors are corrected, the GRACE-inferred regional mass balance agrees with IOM mass balance in DS4a due to the large uncertainties in the GRACE solution and the RACMO2.3 model there, i.e., ±17 Gt yr −1 (see Table A2). From Fig. 5 we can also make an inference about the effect of using different methods to estimate the reference discharge. That is, a change in the reference discharge will accumulate over time and cause the same amount of change in the mass change rates. Only in DS8a, the IOM and GRACE do not agree within the uncertainties. One reason for the discrepancies could be the discharge from peripheral glaciers, which is not included in the IOM but which does affect GRACE estimates. Previous studies, e.g., Bolch et al. (2013) and Gardner et al. (2013), show that a mass rate of approximately 40 Gt yr −1 comes from the peripheral glaciers. However, this is not the reason for the difference between GRACE and the IOM. In this study we are using RACMO2.3 SMB estimates for not only the GrIS but for all of Greenland, including the majority of the mass loss from the peripheral glaciers and ice caps. Discharge from the peripheral glaciers and ice caps is expected to be small, because there are far less marine-terminating glaciers that drain the glaciers and ice caps than those that drain the ice sheet. Less than half of the glaciers and ice caps are marineterminating in Greenland (Gardner et al., 2013). Moreover, given the relationship in the discharge data found by Enderlin et al. (2014) between glacier width and area for the ice sheet's marine-terminating glaciers, we expect the discharge from these glaciers to be small. Consequently, we expect the regional mass change in these glacier areas to be dominated by changes in SMB, which are captured by RACMO2.3. The GRACE-IOM difference due to the exclusion of discharge from peripheral marine-terminating glaciers and ice caps will likely be negligible as long as the SMB for the whole of Greenland is considered, not just the ice sheet.
For the regions above 2000 m altitude, GRACE-inferred regional mass change rates agree with the estimations from IOM within their uncertainties (see Fig. 5). A noticeable mass increase appears in both the GRACE and IOM solutions in DS2b (northeastern interior). A second observation is that in the IOM the runoff dominates the regional mass balance on the edge of the southern GrIS interior, which results in a mass loss of −8 Gt yr −1 . The overall IOM uncertainties in the coastal regions are mainly influenced by the uncertainties in the SMB and D estimates, but the assumption on the flux across the 2000 m contour (F 2000 ) contributes to additional uncertainties in the GrIS interior regions. In the GRACE solution, the uncertainties are due to the errors in the GRACE coefficients which are not dependent on the altitude. Therefore, the uncertainty level is similar for the interior regions and the coastal regions.
We also compare our GRACE and IOM solutions to (1) GRACE, (2) IOM and (3) ICESat altimetry estimates from different studies, as shown in Table 1. All listed GRACE solutions agree within the uncertainty levels in DS1, DS2, DS3, DS5 and DS8. In the southeastern region DS4, there is a deceleration of the mass change after 2007, when the regional acceleration of mass loss becomes negligible (∼ −0.1 Gt yr −2 ). However, different GRACE solutions for the same time period do not result in the same mass change rates. This suggests that a large approximation error, which is associated with different approximation approaches, is likely present in the GRACE solutions for this region. As shown in Fig. 5, if we consider the time period of 2003-2013, the regional mass change is reduced by 29 % in this region after applying the correction introduced in Sect. 4.2.
The IOM is also relatively uncertain in DS4 (Sasgen et al., 2012). Even though the mass-change rates between GRACE and the IOM in this region show a relatively large difference, agreement is obtained within the large uncertainties. For ICESat-based mass loss estimates, the retrieved longterm mass loss can be very different, e.g., −75 ± 2 Gt yr −1 by Zwally et al. (2011) compared to −43 ± 11 Gt yr −1 by Sørensen et al. (2011). This may be explained by the complicated regional ice surface geometry in the coastal areas (Zwally et al., 2011) or uncertainty resulting from the conversion of height changes to mass change, e.g., different firn corrections and density conversions.

2003-2007
−14 ± 6 −1 ± 8 −49 ± 6 −60 ± 13 −18 ± 5 −6 ± 9 −25 ± 6 −32 ± 6 GRACE (IOM) IOM solutions on all the neighbor regions of DS8, no significant differences between the IOM and GRACE solutions are found. This suggests that the remaining approximation error is not the major source of the difference in DS8. The uncertainties of the GIA effect are included as part of the uncertainties of the GRACE solution for DS8 as well (see Table A3), but adding these still cannot bridge the gap between GRACE and IOM. The ICESat-and Operation IceBridge-based mass change estimates by Kjeldsen et al. (2013) yield a mass loss rates of 55 ± 8.4 Gt yr −1 from 2003 to 2010, which is consistent with the GRACE solution obtained in this study. A combination of all evidence indicates that the IOM method underestimates the mass loss rate in this drainage area by ∼ 15 Gt yr −1 . In Sasgen et al. (2012), the discharge estimations from Rignot-08 are used, in which a portion of DS8 was un-surveyed, to which they attributed the difference between GRACE and IOM (24 ± 13 Gt yr −1 ). In this study, the discharge estimation from Enderlin et al. (2014) covers the entire glacier area in this region, but only for the years after 2000. Therefore, despite observations of relatively stable terminus positions for the majority of the marine-terminating glaciers in northwestern Greenland between 1985 and 2000 (Howat and Eddy, 2011), we hypothesize that the estimated reference discharge overestimates the regional D 0 . In deriving D 0 from D-14, we have used the assumption that D from 1990 to 2000 follows Rignot-08, which contains high regional uncertainties. However, if we use the runoff-based estimate of D 0 , uncertainties are influenced by the uncertainty of the RACMO2.3 model. The SMB intercomparison study of Vernon et al. (2013) shows that the 1961-1990 reference SMB 0 of RACMO2.3 model is larger than some other SMB models, e.g., MAR or PMM5. It is interesting to see that when the cumulative TMB is calculated independently from the reference SMB 0 and D 0 (using Eq. (5) and method 1), the mass change rates agrees with the GRACE mass balance in this region within uncertainties. This indicates that the modeled SMB (as well as SMB 0 ) could have uncertainties that are larger than 18 %.

Conclusions
In this study, we implement a simulation of the GrIS mass change and show that the approximation errors caused by the least squares inversion approach can be quantified and reduced in the GRACE solution. When using the IOM, we also improve the reference discharge estimate by utilizing the modeled runoff. We show that the regional differences between our GRACE and IOM solutions are reduced and agree within their calculated confidence intervals. This is confirmed by an intercomparison with ICESat-based regional mass change rates. In the southeast, the corrections for the approximation errors in GRACE data products are especially important. We find that the IOM solutions underestimate the mass loss in the northwest compared to the GRACE and ICESat solutions, which we attribute to incorrect estimates in the reference D and/or SMB that are used to construct the IOM estimates. For the whole GrIS, we find a The GrIS ice discharge D was distributed into 34 glaciers by Rignot et al. (2008), denoted here as D-08. We define the discharge in 1996 and 2000 as D 1996 -08 and D 2000 -08, respectively. Note that the D 1996 -08 was applied as the reference D 0 in D-08 (noted D 0 -08). The deviations between D 1996 -08 and D 2000 -08 are due to discharge changes in the late 1990s (Enderlin et al., 2014). Similarly, we define Enderlin-14 as D-14, with the time series starting from the year of 2000 (D 2000 -14). In order to estimate the reference discharge D 0 for D-14, we find scaling factors between D 1996 -08 and D 2000 -08 and scale the D 2000 -14 to yield the estimation of D 1996 -14. We estimate the uncertainties of estimated D 1996 -14 via 500 pairs of randomly generated D 1996 -08, D 2000 -08 and D 2000 -14, following from a normal distribution N (D, σ D ), in which σ D is the error in the discharge estimations. For the entire GrIS, we find that the interpolated reference D 1996 -14 equals 413.8 ± 31.6 Gt yr −1 , similar to previous studies (Sasgen et al., 2012;); therefore in this study we apply D 1996 -14 as the reference D 0 when using IOM method 2).

A1 Approximation error correction
To determine the linear relationship between the simulated regional mass balances with the associated approximations after applying the least squares inversion, the linear-fit parameters α 0 and α 1 are calculated for different simulation error levels, the values of which are shown in Table A1. The values of α 0 and α 1 and their uncertainties vary slightly in all coastal regions. To determine one value for α 0 and α 1 , we assume the α 0 and α 1 follow a normal distribution in each region and draw 1000 random samples for each error level. Then we combine all the samples and fit into another normal distribution, from which the α 0 and α 1 are determined for each region (see the Table A1).

A2 The GrIS simulation
The GrIS monthly mass balance simulations that are used in Sect. 4.2 are based on the RACMO2.3 model and the discharges estimates from Enderlin et al. (2014). Note that the discharge estimates are given in the form of lumped mass change for 178 different geographical locations. To get SMB and D estimates for each drainage area we sum the discharges for all glaciers and the gridded SMB values within the same drainage area, respectively. We interpolate SMB and D onto a gridded map of EWH with a resolution of 1 • of arc for the GrIS and surrounding areas. To account for leakage from outside the GrIS, as occurs for data products obtained from GRACE, we apply the annual mass change estimates from Schrama et al. (2014) for all the major glacier areas (GrIS excluded). We convolve the gridded mass dis-tribution over the Earth's surface and obtain the potential coefficients in response to this distribution up to d/o 60. Noise in the monthly GRACE coefficients manifests mainly as north-south stripes in the spatial domain (Swenson and Wahr, 2006). In order to mimic this error in the simulation, we add randomly generated noise to the potential coefficients, as described in Bonin and Chambers (2013). The simulation model was discussed in detail by Xu et al. (2015). Note that for this study we focus on the discussion of longterm linear trend, so the linear trend of the monthly simulation is used as the simulation model.

A3 Uncertainty estimations
A summary of the uncertainties in the regional mass balance (linear trend) is shown in Table A2. In our GRACE-inferred mass balance, the uncertainties are associated with (1a) the standard deviations of the CSR RL05 GRACE spherical harmonics coefficients (including the standard deviations of the external degree l = 1 and 2 coefficients), (2a) the variations of the regional mass change due to different GIA models and (3a) the uncertainties due to the corrections of the systematic errors in the least squares inversion solutions. The uncertainties of the IOM-inferred mass balance consist of the uncertainties of the 1961-1990 reference in SMB 0 and in D 0 and (2b) the systematic error in the SMB (RACMO2.3) and (2c) the errors in the yearly D estimations (Enderlin et al., 2014;Rignot et al., 2008).

A4 Selection of the GIA model for GrIS regions
We apply the GIA correction to the GRACE data using three GIA models with a total of 11 different parametrizations before estimating the associated regional mass change in 20 GrIS and surrounding Arctic regions (see the mascon definition in Sect. 3). By comparing with one without applying GIA correction, we assume the differences are the regional GIA effects. In addition to the Paulson-07 GIA model, we use a GIA model with lateral changes in viscosity and the ICE-5G loading history (van der Wal et al., 2013). Moreover, we use another GIA model based on the ice history model from Simpson et al. (2009), provided by Glenn Milne within the scope of the IMBIE project. The upper mantle viscosity ranges from 0.3 × 10 21 to 1 × 10 21 Pa s and the lower mantle viscosity ranges from 1 × 10 21 to 10 × 10 21 Pa s. The thickness of the lithosphere is assumed to be 96 or 120 km.
In Table A3, the GIA related mass change can vary from −7 to 10 Gt yr −1 for the entire GrIS. A positive GIA effect appears in the northern GrIS while in the southern and southwestern GrIS (DS5a to DS7a), negative GIA signals prevail. In order to quantify the uncertainties of the regional GIA in the Paulson-07, since it is the GIA model we used to derive our GRACE solution, we estimate the standard deviation of www.the-cryosphere.net/10/895/2016/ The Cryosphere, 10, 895-912, 2016 all models with respect to . The uncertainties are summarized in Table A2.
The Cryosphere, 10, 895-912, 2016 www.the-cryosphere.net/10/895/2016/ Table A1. The linear-fit parameters α 0 and α 1 that describe the relationship between the regional simulated mass balance and the approximations obtained after the inversion procedure, as applied to GRACE data of the coastal regions. For the interior GrIS regions, we show the approximation errors as additional uncertainties.  Table A2. The uncertainties associated with the regional mass change rates. For the GRACE-inferred regional solutions, "coef. SD" refers to the errors due to the standard deviations in the CSR RL05 spherical coefficients, "GIA" refers to the errors obtained from comparing 11 GIA models. Note that the GIA uncertainties in the interior GrIS are all close to 0 and are therefore negligible. In the column with the header "cor" we show the uncertainties which are caused by the approximation error correction. For SMB and D trend estimations, the uncertainties consist of the reference SMB 0 and D 0 error ("SMB 0 " and "D 0 ") and the systematic errors in RACMO2.3 model and in the discharge estimations ("sys"). The column titled "cum. uncer." refers to uncertainties using the assumptions (1) and (2)