# Reconciling GRACE and GPS estimates of long-term load deformation in southern Greenland

Reconciling GRACE and GPS estimates of long-term load deformation in southern Greenland Summary We examine vertical load deformation at four continuous Global Positioning System (GPS) sites in southern Greenland relative to Gravity Recovery and Climate Experiment (GRACE) predictions of vertical deformation over the period 2002–2016. With limited spatial resolution, GRACE predictions require adjustment before they can be compared with GPS height time series. Without adjustment, both GRACE spherical harmonic (SH) and mascon solutions predict significant vertical displacement rate differences relative to GPS. We use a scaling factor method to adjust GRACE results, based on a long-term mass rate model derived from GRACE measurements, glacial geography, and ice flow data. Adjusted GRACE estimates show significantly improved agreement with GPS, both in terms of long-term rates and interannual variations. A deceleration of mass loss is observed in southern Greenland since early 2013. The success at reconciling GPS and GRACE observations with a more detailed mass rate model demonstrates the high sensitivity to load distribution in regions surrounding GPS stations. Conversely, the value of GPS observations in constraining mass changes in surrounding regions is also demonstrated. In addition, our results are consistent with recent estimates of GIA uplift (∼4.4 mm yr−1) at the KULU site. Global change from geodesy, Loading of the Earth, Satellite geodesy, Satellite gravity, Space geodetic surveys, Time variable gravity 1 INTRODUCTION The Greenland ice sheet (GrIS) is the second largest mass of ice on Earth and would raise global mean sea level (GMSL) by ∼7 m if completely melted (Gregory et al.2004). With Arctic warming, the GrIS has undergone significant melting, and is regarded as a main contributor to GMSL rise (Shepherd et al.2012). Improved quantification of GrIS mass balance at various spatial and temporal scales and understanding of its climatological and geodynamical impacts continue to be important in understanding climate change. Since the Gravity Recovery and Climate Experiment (GRACE) satellite gravity mission was launched in March 2002 (Tapley et al.2004), investigation of ice mass changes over polar ice sheets and mountain glaciers has entered a new era (e.g. Chen et al.2006; Velicogna & Wahr 2006; Khan et al.2010). By accurately measuring Earth's gravity, GRACE provides a unique monthly measure of global large-scale mass changes (Cazenave & Chen 2010). Early GRACE studies (e.g. Chen et al.2006; Velicogna & Wahr 2006) found accelerated GrIS ice loss since 2005, mostly in East and Southeast Greenland. Khan et al. (2010) and Chen et al. (2011) reported additional losses in northwest Greenland starting in late 2005, with evidence of acceleration since 2008. A recent study (Forsberg et al.2017) found large year-to-year variations during the GRACE era, especially a large 2012 melt event (Nghiem et al.2012). The estimated average GrIS ice loss rate is 265 ± 25 Gt yr−1 (Forsberg et al.2017). Continuous Global Positioning System (GPS) measurements on bedrock adjacent to the GrIS are able to detect surface displacements caused by ice mass unloading in Greenland. A number of studies have demonstrated the value of these load deformation observations (e.g. Khan et al.2007, 2010; Jiang et al.2010; Bevis et al.2012; Nielsen et al.2012; Larson et al.2015; Liu et al.2017). In combination, GRACE and GPS data should provide additional understanding of ice mass changes. Khan et al. (2010) compared GRACE-predicted surface deformation with GPS measurements in Greenland. They compared GPS loading signals and GRACE estimates, and confirmed increasing ice mass loss towards the north. Joud et al. (2017) also compared GRACE-predicted uplift rates with GPS measurements to validate their methods for processing GRACE data in regions with a large glacial isostatic adjustment (GIA) signal. The coarse spatial resolution (several hundred kilometres) of GRACE measurements and spatial leakage effects have been a major limitation in related studies. King et al. (2012) reduced the leakage effects to refine the estimation of Antarctic sea-level contribution, and Mémin et al. (2014) also applied leakage correction in studying snow and ice height change. In load deformation studies, limited resolution, due to the range of available spherical harmonic (SH) coefficients and spatial filtering used to suppress noise, contributes to differences between load deformation estimated from GRACE and GPS data. In some earlier studies GRACE-GPS differences were attributed to errors in GPS data processing (e.g. King et al.2006, van Dam et al.2007). Khan et al. (2010) found that GRACE estimates were smaller than GPS with ratios of GPS to GRACE uplift rates at THU3, KULU and KELY stations, of ∼2.9, 2.7 and 2.3, respectively. Khan et al. (2010) assumed similar ratios for ice mass loss acceleration and trends and used these to adjust GRACE trend and acceleration values for shorter time spans. The scaling factor method has been shown to be useful for correcting spatial leakage in GRACE mass change estimates (e.g. Chen et al.2007, 2017; Landerer & Swenson 2012). Scaling factors, at either regional scales or at grid points, are commonly derived from model-predicted mass change and used to adjust GRACE estimates. However, the scaling factor method relies on accurate models, and while seasonal predictions of available climate model are often good, performance tends to degrade at longer time scales. In this study, we use an improved scaling method to compare GPS-observed and GRACE-predicted uplift rates at 4 GPS sites (KELY, QAQ1, SENU and KULU) in southern Greenland (see Fig. 1) chosen for their good data quality and availability. We examine GPS and GRACE data over the period 2002–2016, long enough to examine both long term rates and interannual variation. The improved scaling factor method is based on a mass rate model that combines independent information from GRACE, geographic (Google Earth) data, and ice dynamic data. Figure 1. View largeDownload slide Map of the Greenland ice sheet, with locations of four selected GPS sites in southern Greenland marked by red dots (revised from Vasskog et al.2015). Figure 1. View largeDownload slide Map of the Greenland ice sheet, with locations of four selected GPS sites in southern Greenland marked by red dots (revised from Vasskog et al.2015). 2 DATA PROCESSING 2.1 GRACE spherical harmonic solutions GRACE time-variable gravity data products include 153 Release 5 (RL05) GSM monthly solutions, covering the period August 2002 to August 2016, provided by the Center for Space Research (CSR) at the University of Texas at Austin (Bettadpur 2012). Each monthly solution consists of fully normalized spherical harmonic (SH) coefficients up to degree and order 60. Atmospheric and oceanic effects have been removed during GRACE data processing through the de-aliasing process (Bettadpur 2012), so that over Greenland, GRACE GSM month-to-month changes mainly reflect ice mass and terrestrial water storage (TWS) effects. Following standard practice, we replace GRACE C20 coefficients by independent estimates from satellite laser ranging (SLR) provided by CSR (Cheng et al.2013), and include degree-1 coefficients estimated using the method of Swenson et al. (2008). GIA (Post Glacial Rebound, PGR) effects are removed using the model of A et al. (2013) (based on ICE-5G ice history and VM2 viscosity profile). To compare GRACE with GPS results, atmospheric and oceanic loading effects (removed from GRACE data during de-aliasing) have to be added back. We add to GSM the GAC products, representing the combined effect of atmosphere and ocean mass variations in GRACE AOD1B de-aliasing products (Flechtner et al.2014). Elastic vertical load deformation dr(θ, ϕ), at colatitude θ and longitude ϕ, can be calculated from GRACE Stokes coefficients as follows (van Dam et al.2007):   \begin{eqnarray} dr\left( {\theta ,\phi } \right) &=& {\rm{R}}\sum\limits_{l \!=\! 0}^\infty {\sum\limits_{m = 0}^l {{W_l}{{\tilde{P}}_{l,m}}\left( {\cos \theta } \right)\left( {{C_{lm}}\cos \left( {m\phi } \right) \!+\! {S_{lm}}\sin \left( {m\phi } \right)} \right)} }\nonumber\\ &&\times \frac{{h^{^\prime}}\!_l}{{1 + {{k^{^\prime}}\!_l}}} \end{eqnarray} (1)where R is the radius of Earth, $${\tilde{P}_{l,m}}$$ are normalized Legendre functions of degree l and order m, Clm and Slm represent the SH (Stokes) coefficients of the gravity field, and $$h_l^\prime$$ and $$k_l^\prime$$ are elastic load Love numbers of degree l, whose values for a Gutenberg–Bullen (G-B) elastic Earth model are provided by Farrell (1972). Wl is the weighting value at each spherical harmonic degree to implement Gaussian smoothing (Jekeli 1981; Wahr et al.1998), with a radius of 300 km in the present study. 2.2 GRACE mascon solutions CSR RL05 GRACE 0.5° × 0.5° mascon solutions for August 2002 to August 2016 are also used (Save et al.2016). Mascon solutions are obtained using several constraints, including separation of land and ocean signals to reduce leakage, and forward modelling of annual and trend signals from ice losses in polar regions and mountain glaciers. Unlike SH solutions, mascon solutions do not require additional filtering to suppress noise (Save et al.2016). Mascon solutions are obtained with C20 coefficients replaced by SLR solutions (Cheng et al.2013), and the same degree-1 terms and GIA corrections (Swenson et al. (2008) and A et al. (2013)). Note that the mascon solutions have restored the GAD products, representing the atmospheric and oceanic effects over the ocean (in contrast to GAC, GAD data set is set to zero at the continents), so we compute the total mascon solution as (CSR mascon solution + GAC − GAD) to reflect the same load signal as SH solution (GSM + GAC). Mascon 0.5° × 0.5° solutions dh(θ, ϕ) at colatitude θ and longitude ϕ, are converted to Stokes coefficients up to degree and order 360, corresponding to about 0.5° resolution at the equator, using eqs (2)–(4),   $${\hat{C}_{lm}}\,{\rm{ = }}\,\frac{{\rm{1}}}{{\rm{R}}} \cdot \frac{{\rm{1}}}{{{\rm{4\pi }}}}\iint_{\sigma }{{{\rm{d}}h\left( {\theta ,\phi } \right){{\tilde{P}}_{l,m}}\left( {\cos \theta } \right)\cos \left( {m\phi } \right){\rm{d}}\sigma }}$$ (2)  $${\hat{S}_{lm}}\,{\rm{ = }}\,\frac{{\rm{1}}}{{\rm{R}}} \cdot \frac{{\rm{1}}}{{{\rm{4\pi }}}}\iint_\sigma {{\rm{d}}h\left( {\theta ,\phi } \right){{\tilde{P}}_{l,m}}\left( {\cos \theta } \right)\sin \left( {m\phi } \right){\rm{d}}\sigma }$$ (3)  $$\left\{ {\begin{array}{@{}*{1}{c}@{}} {{C_{lm}}}\\ {{S_{lm}}} \end{array}} \right\} = \frac{{3{\rho _w}}}{{{\rho _{{\rm{ave}}}}}}\frac{{1 + {{k^{^\prime}}\!_l}}}{{2l + 1}}\left\{ {\begin{array}{@{}*{1}{c}@{}} {{{\hat{C}}_{lm}}}\\ {{{\hat{S}}_{lm}}} \end{array}} \right\}$$ (4)where $${\hat{C}_{{\rm{lm}}}}$$ and $${\hat{S}_{{\rm{lm}}}}$$ are dimensionless coefficients defined in eq. (9) of Wahr et al. (1998), dσ is the area element, which equals to sin (θ) dθdϕ, ρw is the density of water, ρave is the average density of the Earth. Then we substitute Stokes coefficients Clm and Slm into eq. (1) to obtain vertical load deformation. Note that Wl = 1 when using mascon-based SH coefficients, because spatial smoothing is not needed. 2.3 GPS data The Jet Propulsion Laboratory's (JPL’s) GNSS-Inferred Positioning System (GIPSY) solutions of daily GPS height at the four sites (KELY, QAQ1, SENU and KULU, shown in Fig. 1) are available at ftp://sopac-ftp.ucsd.edu/pub/timeseries/measures/ats/. Data at KELY, QAQ1, and KULU span the GRACE era (August 2002 to August 2016), while SENU data are available for May 2008 to August 2016. JPL’s reanalysis orbit and clock products are determined using a set of models in a free-network reference frame, including gravity from Earth, Sun, Moon, and other planets. Models include: the DE421 planetary ephemeris (Folkner et al.2008), the GSPM10 satellite solar pressure model, the IAU06 model for precession and nutation, the IERS2010 tide model (Petit & Luzum 2010), the FES2004 ocean loading model (Lettelier et al.2004), and International GNSS Service satellite and receiver antenna phase centre models (ftp://sideshow.jpl.nasa.gov/pub/JPL_GPS_Timeseries/repro2011b). We remove GIA uplift rates using predictions from A et al. (2013), used also in correcting GRACE data. Uplift rates are −2.25 mm yr−1 at KELY, 2.15 mm yr−1 at QAQ1, 1.83 mm yr−1 at SENU, and −0.76 mm yr−1 at KULU. GIA model prediction error in Greenland is thought to be relatively small (Velicogna & Wahr 2006). However, a recent study (van Dam et al.2017) suggests that the GIA rate at KULU may differ significantly from model predictions. Using GPS observations and superconducting gravimeter measurements, van Dam et al. (2017) show that the GIA uplift rate at KULU is about 4.49 mm yr−1, compared with the (A et al.2013) prediction of −0.76 mm yr−1. This affects the GRACE-GPS uplift comparison, as discussed below. 3 RESULTS 3.1 GRACE and GPS comparison We use GRACE GSM + GAC SH coefficients (up to degree and order 60) to compute vertical displacements using eq. (1) with 300 km Gaussian smoothing at KELY, QAQ1, SENU and KULU. Fig. 2 compares monthly GRACE vertical displacements and GPS height series (blue and red curves). Upward trends at all stations reflect continued mass loss at these sites or nearby. With the exception of QAQ1 (Fig. 2b), GRACE vertical displacements are smaller than GPS, consistent with previous studies (Khan et al.2010). GRACE-GPS differences are believed to be partly attributed to the coarse spatial resolution of GRACE (Khan et al.2010). Due to the truncation of RL05 SH coefficients at degree and order 60 and the applied 300 km Gaussian smoothing, GRACE SH solutions underestimate the mass change and load deformation (Chen et al.2007; Wang et al.2016). Figure 2. View largeDownload slide Comparisons among GRACE (GRC) vertical displacements using SH solutions (green), mascon solutions (blue) and GPS heights (red). A 300 km Gaussian smoothing has been applied to the GRACE SH solutions. Model predicted GIA effects (A et al.2013) have been removed from both GRACE and GPS. Figure 2. View largeDownload slide Comparisons among GRACE (GRC) vertical displacements using SH solutions (green), mascon solutions (blue) and GPS heights (red). A 300 km Gaussian smoothing has been applied to the GRACE SH solutions. Model predicted GIA effects (A et al.2013) have been removed from both GRACE and GPS. GRACE mascon series should offer higher spatial resolution and provide more accurate load deformation estimates than the SH solutions, because they are not subject to the dominant stripe noise in SH solutions, and do not require spatial smoothing. We convert gridded mascon mass change into vertical displacements using eqs (2)–(4) and (1), and superimpose the results in Fig. 2 (green curve). SH representations of mascon solutions to degree and order 360 should retain their spatial resolution. Compared with predictions from the RL05 SH series (limited to degree and order 60), the GRACE mascon predictions are more similar to GPS height variations at KELY, SENU and KULU as shown in Figs 2(a), (c) and (d), but differences are still relatively large. However, the QAQ1 series (Fig. 2b) differs markedly from the other sites. The GPS uplift rate at QAQ1 is smaller than the GRACE prediction, and GRACE mascon agreement with GPS is poor relative to the SH series. This is discussed further below. Least square fit linear trends for the four sites are given in Table 1 (annual and semi-annual signals are fit simultaneously, but not shown). GPS trend uncertainty is estimated from GPS formal error (provided with the GPS height series). GRACE trend uncertainty is based upon least squares fitting error. Large differences among vertical rates from GPS, GRACE SH, and GRACE mascon estimates are evident in Table 1. We compare in Fig. 3 GRACE SH and mascon vertical rate predictions over Greenland. GRACE mascon rates are significantly larger and show better spatial resolution. The largest vertical displacements are along southeastern and western coastal regions. However, comparisons with GPS observations discussed above indicate that neither GRACE estimate is likely to agree very well with GPS observations, without further processing of GRACE data. Figure 3. View largeDownload slide Vertical displacement rates (in mm yr−1) computed from: (a) GRACE GSM + GAC Stokes coefficients with 300 km Gaussian smoothing; (b) GRACE mascon solutions (0.5° × 0.5°) after adding back the contributions of atmospheric at the continents. GIA correction is applied to both GRACE SH and mascon solutions. Figure 3. View largeDownload slide Vertical displacement rates (in mm yr−1) computed from: (a) GRACE GSM + GAC Stokes coefficients with 300 km Gaussian smoothing; (b) GRACE mascon solutions (0.5° × 0.5°) after adding back the contributions of atmospheric at the continents. GIA correction is applied to both GRACE SH and mascon solutions. Table 1. Predicted uplift rates using GRACE and GPS data (units: mm yr−1). UW means unweighted model, W means weighted model.   KELY  QAQ1  SENUa  KULU  GPS  7.70 ± 0.13  2.51 ± 0.06 (2.69 ± 0.15)b  8.06 ± 0.25  9.88 ± 0.06 [4.63 ± 0.06]c  GRACE SH  4.51 ± 0.08  3.23 ± 0.05 (3.39 ± 0.12) b  3.54 ± 0.13  4.15 ± 0.05  GRACE mascon  5.50 ± 0.10  4.71 ± 0.06  6.10 ± 0.18  6.54 ± 0.06  GRACE restored (UW)  6.85 ± 0.11  3.91 ± 0.06  8.04 ± 0.25  4.31 ± 0.05  GRACE restored (W)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05    KELY  QAQ1  SENUa  KULU  GPS  7.70 ± 0.13  2.51 ± 0.06 (2.69 ± 0.15)b  8.06 ± 0.25  9.88 ± 0.06 [4.63 ± 0.06]c  GRACE SH  4.51 ± 0.08  3.23 ± 0.05 (3.39 ± 0.12) b  3.54 ± 0.13  4.15 ± 0.05  GRACE mascon  5.50 ± 0.10  4.71 ± 0.06  6.10 ± 0.18  6.54 ± 0.06  GRACE restored (UW)  6.85 ± 0.11  3.91 ± 0.06  8.04 ± 0.25  4.31 ± 0.05  GRACE restored (W)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05  aThe uplift rates at SENU are computed using data from May 2008 to August 2016, and others are from August 2002 to August 2016. bThe uplift rates at QAQ1 in the parentheses are computed using data from May 2008 to August 2016. cThe GPS uplift rate at KULU in the square brackets is computed using GIA uplift rate predicted by van Dam et al. (2017). View Large 3.2 Correcting GRACE spatial leakage Due to the truncation of GRACE SH coefficients at degree and order 60 and additional spatial filtering (300 km Gaussian), GRACE predictions of vertical displacements are limited to long-wavelength features. While GRACE mascon solutions offer qualitative improvements in spatial resolution, they remain subject to leakage effects due to fundamental GRACE spatial resolution limits, and it is easier to quantitatively address resolution limits of SH solutions. If the spatial spectrum (spatial distribution) of mass change is approximately known, either from climate models or from observations, scaling factors can be estimated to adjust GRACE values to reduce bias associated spatial leakage as discussed above (Chen et al.2007, 2017; Landerer & Swenson 2012; Wang et al.2016). Recognizing that long-term ice mass changes are dominantly along the edges of the GrIS, we construct an ice mass change model based on GRACE mascon solutions, published ice flow velocity data (Khan et al.2015), and known glacier locations (Google Earth). The scaling factor at a grid point is the ratio between model-predicted load deformation and the prediction after SH truncation and spatial filtering. Scaling factors are then used to adjust GRACE SH predictions of vertical motion at the four GPS sites. The procedure includes the following steps: From GRACE mascon gridded equivalent water thickness (EWT) data, we subtract the GAD products. The residual EWT contains ice and TWS mass change effects, similar to GRACE GSM data (after GIA correction). The mass rate of residual EWT at each grid from August 2002 to August 2016 is found by least squares with the result shown in Fig. 4(a). We divide the GrIS and surrounding areas into 12 regions (Fig. 4b), and find the total mass rate within each using GRACE mascon rates (Fig. 4a). Mass rates of the 12 regions are in Table 2 (units of gigatonne per year, Gt yr−1). Note that the sum of area 1 to 10 represents the total mass loss of Greenland, that is, −252.5 Gt yr−1, which is consistent with previous publications (e.g. Velicogna & Wahr 2013; Wouters et al.2013). We construct a detailed mass rate model along the periphery of southern Greenland glaciers (to latitudes 72.5°N) by assigning mass rates in the 12 regions into narrower bands (Fig. 5b), to represent dominant coastal and central Greenland locations associated with long-term changes. Within each narrower region, mass change is uniformly distributed. The Fig. 5(b) mass rate model is on a 0.1° × 0.1° grid, with a width of 0.5° (∼25 km) along glacier edges, based on Google Earth glacier geography. In regions far from the GPS stations, a lower resolution model (0.5° × 0.5°) is used (Fig. 5a). An additional mass rate model is constructed, in which areas near the GPS stations are weighted according to observed ice flow velocity data (fig. 8 in Khan et al.2015). Weighted solutions retain the same total mass rate in each subregion. We define the weight as approximately proportional to the velocity (e.g. we define the light yellow colour region in their figure as weight 1.0), and then roughly allocate 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6 or 1.8 as gridded weight according to the colour of corresponding grid in the ice flow velocity map. The weighted (W) model is shown in Fig. 5(c), compared with the unweighted (UW) model in Fig. 5(b). The mass rate models are represented in an SH expansion to degree and order 1800, and from these vertical displacements calculated using both the full range Stokes coefficients, and Stokes coefficients limited to degree and order 60 Stokes with additional 300 km Gaussian smoothing. At each grid point, the ratio between the two displacements is the scaling factor used to adjust GRACE estimates. We apply scaling factors to vertical displacement rates predicted by GRACE SH products, and add back GAC and GIA models. These can now be compared with GPS observations for both UW and W models in Figs 6 and 7 and Table 1. Figure 4. View largeDownload slide (a) GRACE mascon mass rates (in units of cm yr−1 of equivalent water thickness change) after the GAD and GIA models are removed; (b) the 12 shaded areas used to integrate the total mass rate within each of the regions in Fig. 4(a). Figure 4. View largeDownload slide (a) GRACE mascon mass rates (in units of cm yr−1 of equivalent water thickness change) after the GAD and GIA models are removed; (b) the 12 shaded areas used to integrate the total mass rate within each of the regions in Fig. 4(a). Figure 5. View largeDownload slide (a) The 0.5° × 0.5° mass rate model approximates ice loss locations and magnitudes in areas distant from the GPS stations; (b) the unweighted 0.1° × 0.1° mass rate model that confines ice loss to locations along the periphery of the ice sheet, with uniform distribution of ice loss within each sub-region near the GPS stations; (c) the weighted 0.1° × 0.1° mass rate model. Weights are based on published ice flow velocity data. The mass rates in the far and near field mass rate models are based on Fig. 4 and Table 2. Figure 5. View largeDownload slide (a) The 0.5° × 0.5° mass rate model approximates ice loss locations and magnitudes in areas distant from the GPS stations; (b) the unweighted 0.1° × 0.1° mass rate model that confines ice loss to locations along the periphery of the ice sheet, with uniform distribution of ice loss within each sub-region near the GPS stations; (c) the weighted 0.1° × 0.1° mass rate model. Weights are based on published ice flow velocity data. The mass rates in the far and near field mass rate models are based on Fig. 4 and Table 2. Figure 6. View largeDownload slide Comparisons of monthly GPS height, GRACE spherical harmonic (SH) predicted vertical displacements and scaling factor adjusted GRACE series using the weighted mass rate model (Figs 3a and c). Uplift rates are shown using the same curve colour. Note: we use two GIA uplift rates at KULU predicted by A et al. (2013) (GIA) and van Dam et al. (2017) (GIA*). GRACE-predicted uplift rates at SENU in the parentheses are computed using the same time period as GPS (2008.05–2016.08). Figure 6. View largeDownload slide Comparisons of monthly GPS height, GRACE spherical harmonic (SH) predicted vertical displacements and scaling factor adjusted GRACE series using the weighted mass rate model (Figs 3a and c). Uplift rates are shown using the same curve colour. Note: we use two GIA uplift rates at KULU predicted by A et al. (2013) (GIA) and van Dam et al. (2017) (GIA*). GRACE-predicted uplift rates at SENU in the parentheses are computed using the same time period as GPS (2008.05–2016.08). Figure 7. View largeDownload slide Comparison of GPS heights and scaling factor adjusted GRACE vertical displacements using the weighted mass rate model after 13-month moving average filtering. The two dash lines mark dates December 2009 and March 2013. Figure 7. View largeDownload slide Comparison of GPS heights and scaling factor adjusted GRACE vertical displacements using the weighted mass rate model after 13-month moving average filtering. The two dash lines mark dates December 2009 and March 2013. Table 2. Total mass rates of the 12 regions in Fig. 4(b) from integrating observed mass rates from the mascon solutions (Fig. 4a). Area #  1  2  3  4  5  6  7  8  9  10  11  12  Gr Mass  Total  Mass Rate (Gt yr−1)  −18.2  −65.3  −7.8  −7.6  −29.2  −44.4  −57.6  −10.8  −15.8  4.2  −71.4  −7.6  −252.5  −331.5  Area #  1  2  3  4  5  6  7  8  9  10  11  12  Gr Mass  Total  Mass Rate (Gt yr−1)  −18.2  −65.3  −7.8  −7.6  −29.2  −44.4  −57.6  −10.8  −15.8  4.2  −71.4  −7.6  −252.5  −331.5  View Large 3.3 Comparing leakage-corrected GRACE with GPS vertical rates We focus on interannual and longer-term variations, because scaling factors are based on GRACE mass rates, associated mainly with mass loss along the periphery of the GrIS. Load variations at shorter time scales (e.g. seasonal) would be expected to have very different spatial distributions due to dominant effects of snowfall and TWS change. Figs 6(a) and (c) show very good agreement between adjusted GRACE and GPS uplift rates at KELY and SENU. At KELY the GPS rate of 7.70 ± 0.13 mm yr−1 is in excellent agreement with the GRACE (GRC) prediction of 7.10 ± 0.12 mm yr−1. Leakage correction steps have raised the GRACE rate from the previous value of 4.51 ± 0.08 mm yr−1. With a shorter time span at SENU (May 2008 to August 2016), GPS and GRACE uplift rates are in excellent agreement 8.06 ± 0.25 mm yr−1 and 8.05 ± 0.25 mm yr−1, respectively. The scaling factor leakage correction has raised the rate considerably, from the former value of 3.54 ± 0.13. KELY is about 40 km from the periphery of the southwest GrIS. The model includes north-south distributed mass loss along the periphery. The weighted model includes mass loss at the nearby land-terminating glacial outlets (Fig. 5c). In this region, ice loss patterns and locations are relatively well known, so the scaling factor method is expected to be especially effective. SENU is just on the periphery of the ice sheet (Fig. 1), making deformation sensitive to nearby mass losses. Consequences of SH truncation and spatial filtering are expected to be particularly large at SENU (Chen et al.2015), but with well-constrained load locations, leakage correction should be effective. At QAQ1 and KULU GRACE-GPS differences persist, especially at KULU. QAQ1 is the only site of the four whose GPS uplift rate is smaller than GRACE prediction (before or after leakage correction). SENU and QAQ1 are only about 70 km apart (Fig. 1), so it is interesting to compare the two sites. We compute uplift rates at QAQ1 for the period of SENU GPS data availability (May 2008 to August 2016) with results in parentheses in Table 1. GRACE SH solutions at these two sites have similar rates (3.39 ± 0.12 versus 3.54 ± 0.13 mm yr−1), while GPS uplift rates are quite different (2.69 ± 0.15 versus 8.06 ± 0.25 mm yr−1). This indicates that mass loss is nearest SENU causing larger elastic load deformation evident with GPS, and with a much smaller effect at QAQ1 (∼70 km away). GRACE SH solutions, with their limited resolution, apparently cause the GRACE load estimate to be distributed over a broader region that includes QAQ1. This indicates that the leakage correction was not fully effective, as discussed below. At KULU, GPS and GRACE uplift rates show the largest difference among the four stations. The GRACE rate (even with leakage correction) is less than half the GPS value (4.30 ± 0.05 versus 9.88 ± 0.06 mm yr−1). While the geographical setting around KULU is relatively complicated, the mass rate model may not be adequate, but we suspect that the large difference is not mainly related to this. As noted earlier, the GIA uplift rate at KULU has been proposed to be substantially different from model predictions (4.49 mm yr−1 by van Dam et al.2017). Khan et al. (2016) also reported a similar GIA rate at KULU (∼4.4 mm yr−1). If we adopt the van Dam et al. (2017) rate at KULU, then the GPS ice load rate becomes 4.63 ± 0.06 mm yr−1, very close to the GRACE prediction (4.30 ± 0.05 mm yr−1) (Fig. 6d). We have used a different GIA value for the GPS data correction only. For GRACE data, we still rely on GIA model prediction (A et al.2013). Furthermore, GIA correction in GRACE estimate (GIA mass rate converted to deformation) is much smaller (∼0.1 mm yr−1), which unlikely has any notable impact on our qualitative analysis. KULU is far from the outlet glaciers relative to KELY and SENU, and is unlikely to experience a much larger contemporary elastic unloading rate (9.88 mm yr−1) than the other two sites. We apply a 13-month moving average filter (with half weight at first and last values) to both GPS and GRACE monthly time series, and compare smoothed series in Fig. 7. Errors are estimated as in Table 1 (described in the 3rd paragraph of Section 3.1). Leakage-corrected GRACE and GPS series show similar trends and interannual changes. We estimate linear rates for three different time spans (August 2002 to December 2009, December 2009 to March 2013, March 2013 to August 2016). All eight series (4 GPS + 4 GRACE) show increasing uplift rates after December 2009, and decreasing rates after March 2013. This is further evidence of accelerated ice loss since 2010 (Forsberg et al.2017). Both GPS and GRACE series indicate a mass loss slow-down since 2013. Results from unweighted and weighted mass rate models are shown in the last two rows of Table 1. 4 DISCUSSION AND CONCLUSIONS Analysis of vertical deformations at four continuous GPS sites, KELY, SENU, QAQ1 and KULU in southern Greenland shows that it is possible to reconcile GPS estimates of elastic loading with GRACE time-variable gravity solutions after GRACE leakage effect is appropriately addressed. The scaling factor method used in this study is based on ice mass change model derived from GRACE observations, which is an improvement from the posteriori factors (not validated through other independent observations) applied in Khan et al. (2010), also from the scaling factors derived from deficient land surface model data (TWS model without groundwater information) by Wang et al. (2016). The good agreements between GPS and restored GRACE at KELY, SENU and QAQ1 indicate that the long-term variations in the three GPS sites are mainly caused by continuous ice mass loss in the surrounding regions. Some uncertainties may affect our comparisons, including uncertainties in GRACE estimates, GIA models, our mass rate models and GPS observations. GRACE estimates are limited by coarse spatial resolution, and subject to bias due to spatial leakage, resulting in apparently large discrepancies with GPS estimates. GRACE mascon solutions offer qualitative improvements in spatial resolution relative to SH solutions (Fig. 3), but correcting leakage bias in mascon solutions is not as straightforward as for SH solutions (Chen et al.2017). The scaling factor method is used to correct GRACE SH data, with careful attention to constructing a long-term mass rate model that combines GRACE measurements and geographical information on glacier location and ice flow velocity. Leakage-corrected GRACE vertical rate predictions then agree remarkably well with GPS observations at KELY and SENU. A relatively large difference at KULU is probably due to an erroneous GIA uplift rate, and when an estimate suggested in an earlier study (van Dam et al.2017) is used, GPS and GRACE rates agree very well. At QAQ1 differences are thought to be due to errors in both the GIA model and the mass rate model related to complex geography. GRACE-GPS results also agree well at interannual time scales, with higher uplift rates during December 2009 to March 2013 at all four sites, suggesting accelerated ice loss during this period. A decrease in uplift rates at the four sites since early 2013, suggests deceleration of mass loss in southern Greenland (Fig. 7). Our finding of interannual variations of GrIS ice mass changes are consistent with results from Forsberg et al. (2017). Excellent agreement between GPS and GRACE estimates gives us confidence in the ability to geodetically monitor ice mass change. The care required to correct for spatial leakage in the GRACE estimates illustrates the sensitivity of load deformation to details of load changes in the surrounding regions. A corollary is that GPS observations provide unique and important constraints on regional mass changes. The present study provides support, though indirect, for recent GIA uplift rate estimates at KULU (Khan et al.2016; van Dam et al.2017, both near 4.4 mm yr−1) which differ from previous GIA model predictions (e.g. Peltier 2004, 2015; A et al.2013; Stuhne & Peltier 2015). GPS uplift rates at KULU are large, near 9.88 mm yr−1, and show less seasonal and interannual variability than the other three sites. KULU is relatively far from major outlet glaciers, so it is sensible that seasonal effects are smaller and that a larger portion of the observed rate is due to GIA, rather than elastic loading. Good GIA uplift rate estimates are essential in interpreting and comparing GPS and GRACE uplift results, but uncertainty in GIA estimates is difficult to determine, mainly due to the lack of long-term in situ measurements. Besides probable uncertainties in GRACE estimates and GIA models, parametrization of the mass rate models also contributes to uncertainty in our conclusions. The current width of mass rate model is 0.5°, which is approximately 25 km at this latitude zone. We try other widths to test the sensitivity and show the results in Table 3. In the experiment, we broaden the width of the mass rate model for glaciers in southwest Greenland (as shown in Fig. 5c) eastward to 0.6° and 1°, and explore that the KELY is the most affected sites in terms of uplift rates, which demonstrates the high deformation sensitivity to near-field loading source. Indeed, Wahr et al. (2013) showed their experiment results between vertical load deformation at a point and the distance from loading disc centre in their fig. 1(a), and Argus et al. (2014) also reported the similar results in their fig. 1. Table 3. Predicted uplift rates of restored GRACE data using weighted model with different widths of the mass rate model for glaciers in southwest Greenland as shown in Fig. 5(c) (units: mm yr−1).   KELY  QAQ1  SENU  KULU  GRACE restored (0.5°)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05  GRACE restored (0.6°)  6.98 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.30 ± 0.05  GRACE restored (1°)  6.55 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.28 ± 0.05    KELY  QAQ1  SENU  KULU  GRACE restored (0.5°)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05  GRACE restored (0.6°)  6.98 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.30 ± 0.05  GRACE restored (1°)  6.55 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.28 ± 0.05  View Large Furthermore, the differences between our unweighted and weighted model results reinforce our conclusion of the effect of load location and shape of glaciers. The weighted model provides more detailed definition of the load location, especially in regions where ice mass loss concentrates in large glaciers other than distributes evenly along the edge. The uplift rate at KELY, located near the Russell Glacier, is notably affected by the weighted mass rate model, and the same is for QAQ1. The sensitivity to load location means that details of glacier location and mass change magnitudes can be important. Additionally, the ice velocity data from Khan et al. (2015) is based on observations between 2000/01 and 2005/06, which is only a first order approximation of ice velocity during GRACE period. We have also compared JPL GPS series used here with GPS solutions from the Scripps Orbit and Permanent Array Center (SOPAC) and MEaSUREs combination solutions (Bock & Webb 2012, ftp://sopac-ftp.ucsd.edu/pub/timeseries/measures/ats/). Fig. 8 plots the three GPS height series at the four selected sites (SOPAC provides KELY and QAQ1 data only), and Table 4 shows their annual amplitudes, annual phases and linear trend signals. Slight differences among the three solutions are found for both annual and linear trend signals. There are a few factors causing such discrepancies, partly due to the differences in data processing strategies between GIPSY (used in JPL) and GAMIT (used in SOPAC), and partly due to the errors during the st_filter combination (used in MEaSUREs combination solution). There will always misfit between JPL and SOPAC solutions, especially one is from Precise Point Positioning (completely depend on FIXED orbit and clock with no tie to any other station) and the other is from network solution with orbit also participating in the adjustment (Peng Fang, personal communication, 2017). The further discussion of GPS time series may be beyond the scope of this article, and need to be presented in a more specific article. Figure 8. View largeDownload slide Comparisons of monthly GPS heights from JPL, SOPAC and MEaSUREs solutions. GIA from A et al. (2013) has been applied to each series. Figure 8. View largeDownload slide Comparisons of monthly GPS heights from JPL, SOPAC and MEaSUREs solutions. GIA from A et al. (2013) has been applied to each series. Table 4. Annual amplitude, annual phase and linear trend corresponding with Fig. 8. The annual phase is defined as φ in cos (2π(t − t0) + φ), where t0 refers to the start of the year (e.g. 2004.0). The uncertainties are estimated from GPS formal error and least squares fitting error.     KELY  QAQ1  SENU  KULU  Annual amplitude (mm)    JPL  5.59 ± 0.72  4.94 ± 0.35  9.09 ± 0.86  4.60 ± 0.35    SOPAC  5.16 ± 0.73  4.81 ± 0.43  /  /    MEaSUREs  6.23 ± 0.76  5.19 ± 0.38  9.53 ± 0.95  5.55 ± 0.47  Annual phase (degree)    JPL  332 ± 7  344 ± 4  324 ± 5  330 ± 4    SOPAC  309 ± 8  316 ± 5  /  /    MEaSUREs  276 ± 7  278 ± 4  294 ± 6  280 ± 5  Trend (mm yr−1)    JPL  7.70 ± 0.13  2.51 ± 0.06  8.06 ± 0.25  9.88 ± 0.06    SOPAC  7.50 ± 0.13  2.31 ± 0.08  /  /    MEaSUREs  7.96 ± 0.14  2.32 ± 0.07  7.25 ± 0.26  8.22 ± 0.08      KELY  QAQ1  SENU  KULU  Annual amplitude (mm)    JPL  5.59 ± 0.72  4.94 ± 0.35  9.09 ± 0.86  4.60 ± 0.35    SOPAC  5.16 ± 0.73  4.81 ± 0.43  /  /    MEaSUREs  6.23 ± 0.76  5.19 ± 0.38  9.53 ± 0.95  5.55 ± 0.47  Annual phase (degree)    JPL  332 ± 7  344 ± 4  324 ± 5  330 ± 4    SOPAC  309 ± 8  316 ± 5  /  /    MEaSUREs  276 ± 7  278 ± 4  294 ± 6  280 ± 5  Trend (mm yr−1)    JPL  7.70 ± 0.13  2.51 ± 0.06  8.06 ± 0.25  9.88 ± 0.06    SOPAC  7.50 ± 0.13  2.31 ± 0.08  /  /    MEaSUREs  7.96 ± 0.14  2.32 ± 0.07  7.25 ± 0.26  8.22 ± 0.08  View Large Disagreement between GPS and GRACE results at QAQ1 may be due to a combination of errors in GIA estimates, GPS data, and mass rate model or ice velocity weighting. Similarities in GPS results at QAQ1 and SENU (two adjacent sites) suggests that errors in GIA or GRACE results may be the source of disagreement. QAQ1 is located at the juncture between southeast and southwest Greenland, so our formulation of the mass rate model from mascon rates (Figs 4b and 5b and c) may produce relatively poor results near QAQ1. Additional observations, for example from satellite altimetry, airborne altimetry, absolute gravimetry, satellite gravimetry, and surface mass balance models may be the basis for an improved mass rate model in the future. We also take into account possible insufficiencies due to the use of an outdated GIA model (A et al.2013, based on the ICE-5G ice history) as well as load Love number (Farrell 1972, based on the G-B Earth model). When using the new ICE-6G GIA model (Peltier 2015), the difference of GPS uplift rates at QAQ1 is ∼0.36 mm yr−1, as compared to the value based on the ICE-5G model (i.e. 2.15 versus 2.51 mm yr−1). When using load Love numbers based on the Preliminary Reference Earth Model (Han & Wahr 1995), the estimated uplift rate at QAQ1 is 3.67 mm yr−1, which is virtually the same as the result (3.66 mm yr−1) based on the G-B Earth model (Farrell 1972) (see Table 1). While long-term ice loss dominates GrIS mass changes, there are also significant short-term fluctuations. Thus long-term continuous observations are needed to understand long-term changes in the GrIS. Deployment of the Greenland GPS Network (GNET; Bevis et al.2009) and continuing satellite gravity observations with the GRACE Follow-On mission, are important elements of this longer-term mission to understand GrIS mass balance. Acknowledgements We are grateful to two anonymous reviewers for their insightful comments. We thank the CSR and NASA’s PODAAC for providing GRACE data products, the JPL and SOPAC for providing GPS data products. We appreciate the discussions with Dr Peng Fang at University of California at San Diego on the GPS solutions. SYW is supported by the UCAS Joint PhD Training Program (UCAS [2015]37). CRW has received support from the Geology Foundation of the University of Texas. JL is supported by the Natural Science Foundation of Shanghai (No.17ZR1435600), the Open Fund of Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University (16-01-05), and the National Key Research and Development Program of China (2016YFB0501405). REFERENCES A, G., Wahr J., Zhong S., 2013. Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to glacial isostatic adjustment in Antarctica and Canada, Geophys. J. Int. , 192( 2), 557– 572. https://doi.org/10.1093/gji/ggs030 Google Scholar CrossRef Search ADS   Argus D.F., Fu Y., Landerer F.W., 2014. Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion, Geophys. Res. Lett. , 41( 6), 1971– 1980. https://doi.org/10.1002/2014GL059570 Google Scholar CrossRef Search ADS   Bettadpur S., 2012. Gravity recovery and climate experiment, Level-2 gravity field product user handbook (Rev. 3.0, May 29, 2012), GRACE 327-734 (CSR-GR-03-01) , UTCSR, Texas, USA. Bevis M.G., Kendrick E.C., Brown A.K., Khan S.A., Knudsen P., Madsen F., Wahr J., Willis M.J., 2009. Greenland GPS Network: crustal oscillations and seasonal ice mass fluctuations, EOS, Trans. Am. geophys. Un. , 90( 52), Fall Meet. Suppl., Abstract G43B-0728. Bevis M. et al.  , 2012. Bedrock displacements in Greenland manifest ice mass variations, climate cycles and climate change, Proc. Natl. Acad. Sci. USA , 109( 30), 11 944–11 948. https://doi.org/10.1073/pnas.1204664109 Google Scholar CrossRef Search ADS   Bock Y., Webb F.H., 2012. MEaSUREs Solid Earth Science ESDR System , Digital Media. Cazenave A., Chen J., 2010. Time-variable gravity from space and present-day mass redistribution in the Earth system, Earth planet. Sci. Lett. , 298( 3–4), 263– 274. https://doi.org/10.1016/j.epsl.2010.07.035 Google Scholar CrossRef Search ADS   Chen J.L., 2006. Satellite gravity measurements confirm accelerated melting of Greenland ice sheet, Science , 313( 5795), 1958– 1960. https://doi.org/10.1126/science.1129007 Google Scholar CrossRef Search ADS PubMed  Chen J.L., Wilson C.R., Famiglietti J.S., Rodell M., 2007. Attenuation effect on seasonal basin-scale water storage changes from GRACE time-variable gravity, J. Geod. , 81( 4), 237– 245. https://doi.org/10.1007/s00190-006-0104-2 Google Scholar CrossRef Search ADS   Chen J.L., Wilson C.R., Tapley B.D., 2011. Interannual variability of Greenland ice losses from satellite gravimetry, J. geophys. Res. , 116, B07406, doi:10.1029/2010JB007789. https://doi.org/10.1029/2010JB007789 Google Scholar CrossRef Search ADS   Chen J.L., Wilson C.R., Li J., Zhang Z., 2015. Reducing leakage error in GRACE-observed long-term ice mass change: a case study in West Antarctica, J. Geod. , 89( 9), 925– 940. https://doi.org/10.1007/s00190-015-0824-2 Google Scholar CrossRef Search ADS   Chen J.L., Wilson C.R., Tapley B.D., Save H., Cretaux J.F., 2017. Long-term and seasonal Caspian Sea level change from satellite gravity and altimeter measurements, J. geophys. Res. , 192, 557– 517. Cheng M., Tapley B.D., Ries J.C., 2013. Deceleration in the Earth's oblateness, J. geophys. Res. , 118( 2), 740– 747. https://doi.org/10.1002/jgrb.50058 Google Scholar CrossRef Search ADS   Farrell W.E., 1972. Deformation of the Earth by surface loads, Rev. Geophys. , 10( 3), 761– 797. https://doi.org/10.1029/RG010i003p00761 Google Scholar CrossRef Search ADS   Flechtner F., Dobslaw H., Fagiolini E., 2014. Gravity recovery and climate experiment, AOD1B product description document for product Release 05 (Rev. 4.2, May 20, 2014), GRACE 327-750 (GR-GFZ-AOD- 0001), GFZ Potsdam, Potsdam, Germany, p. 33. Folkner W.M., Williams J.G., Boggs D.H., 2008. The Planetary and lunar ephemeris DE 421, JPL memorandum IOM 343R-08-003. Forsberg R., Sorensen L., Simonsen S., 2017. Greenland and antarctica ice sheet mass changes and effects on global sea level, Surv. Geophys. , 38( 1), 89– 104. https://doi.org/10.1007/s10712-016-9398-7 Google Scholar CrossRef Search ADS   Gregory J.M., Huybrechts P., Raper S.C.B., 2004. Climatology: threatened loss of the Greenland ice-sheet, Nature , 428( 6983), 616, doi:10.1038/428616a. https://doi.org/10.1038/428616a Google Scholar CrossRef Search ADS PubMed  Han D., Wahr J., 1995. The viscoelastic relaxation of a realistically stratified earth, and a further analysis of postglacial rebound, Geophys. J. Int. , 120( 2), 287– 311. Google Scholar CrossRef Search ADS   Jekeli C., 1981. Alternative methods to smooth the Earth's gravity field, Report No. 327, Reports of the Department of Geodetic Science and Surveying , Columbus: The Ohio State University. Jiang Y., Dixon T.H., Wdowinski S., 2010. Accelerating uplift in the North Atlantic region as an indicator of ice loss. Nat. Geosci. , 3( 6), 404– 407. https://doi.org/10.1038/ngeo845 Google Scholar CrossRef Search ADS   Joud M.S.S., Sjöberg L.E., Bagherbandi M., 2017. Use of GRACE data to detect the present land uplift rate in Fennoscandia, Geophys. J. Int. , 209( 2), 909– 922. https://doi.org/10.1093/gji/ggx063 Google Scholar CrossRef Search ADS   Khan S.A., Wahr J., Stearns L.A., Hamilton G.S., Van Dam T., Larson K.M., Francis O., 2007. Elastic uplift in southeast Greenland due to rapid ice mass loss, Geophys. Res. Lett. , 34, L21701, doi:10.1029/2007GL031468. https://doi.org/10.1029/2007GL031468 Google Scholar CrossRef Search ADS   Khan S.A., Wahr J., Bevis M., Velicogna I., Kendrick E., 2010. Spread of ice mass loss into northwest Greenland observed by GRACE and GPS, Geophys. Res. Lett. , 37, L06501, doi:10.1029/2010GL042460. https://doi.org/10.1029/2010GL042460 Google Scholar CrossRef Search ADS   Khan S.A., Aschwanden A., Bjørk A.A., Wahr J., Kjeldsen K.K., Kjær K.H., 2015. Greenland ice sheet mass balance: a review, Rep. Prog. Phys. , 78( 4), doi:10.1088/0034-4885/78/4/046801. https://doi.org/10.1088/0034-4885/78/4/046801 Khan S.A. et al., 2016. Geodetic measurements reveal similarities between post-Last Glacial Maximum and present-day mass loss from the Greenland ice sheet, Sci. Adv. , 2( 9), doi:10.1126/sciadv.1600931. https://doi.org/10.1126/sciadv.1600931 King M., Moore P., Clarke P., Lavallée D., 2006. Choice of optimal averaging radii for temporal GRACE gravity solutions, a comparison with GPS and satellite altimetry, Geophys. J. Int. , 166( 1), 1– 11. https://doi.org/10.1111/j.1365-246X.2006.03017.x Google Scholar CrossRef Search ADS   King M.A., Bingham R.J., Moore P., Whitehouse P.L., Bentley M.J., Milne G.A., 2012. Lower satellite-gravimetry estimates of Antarctic sea-level contribution, Nature , 491( 7425), 586– 589. https://doi.org/10.1038/nature11621 Google Scholar CrossRef Search ADS PubMed  Landerer F.W., Swenson S.C., 2012. Accuracy of scaled GRACE terrestrial water storage estimates, Water Resour. Res. , 48, W04531, doi:10.1029/2011WR011453. https://doi.org/10.1029/2011WR011453 Google Scholar CrossRef Search ADS   Larson K.M., Wahr J., Munneke P.K., 2015. Constraints on snow accumulation and firn density in Greenland using GPS receivers, J. Glaciol. , 61, 101– 114. https://doi.org/10.3189/2015JoG14J130 Google Scholar CrossRef Search ADS   Lettelier T., Lyard F., Lefebre F., 2004. The new global tidal solution: FES2004, in Proceeding of the Ocean Surface Topography Science Team Meeting , St. Petersburg, Florida. Liu L., Khan S.A., Van Dam T., Ma J.H.Y., Bevis M., 2017. Annual variations in GPS-measured vertical displacements near Upernavik Isstrøm (Greenland) and contributions from surface mass loading, J. geophys. Res. , 122, 677– 691. https://doi.org/10.1002/2016JB013494 Google Scholar CrossRef Search ADS   Mémin A., Flament T., Rémy F., Llubes M., 2014. Snow- and ice-height change in Antarctica from satellite gravimetry and altimetry data, Earth planet. Sci. Lett. , 404, 344– 353. https://doi.org/10.1016/j.epsl.2014.08.008 Google Scholar CrossRef Search ADS   Nghiem S.V. et al.  , 2012. The extreme melt across the Greenland ice sheet in 2012, Geophys. Res. Lett. , 39, L20502, doi:10.1029/2012GL053611. https://doi.org/10.1029/2012GL053611 Google Scholar CrossRef Search ADS   Nielsen K., Khan S.A., Korsgaard N.J., Kjær K.H., Wahr J., Bevis M., Stearns L.A., Timm L.H., 2012. Crustal uplift due to ice mass variability on Upernavik Isstrøm, west Greenland, Earth planet. Sci. Lett. , 353–354, 182– 189. https://doi.org/10.1016/j.epsl.2012.08.024 Google Scholar CrossRef Search ADS   Peltier W.R., 2004. Global glacial isostasy and the surface of the Ice-Age Earth: The ICE-5G (VM2) Model and GRACE, Annu. Rev. Earth Planet. Sci. , 32, 111– 149. https://doi.org/10.1146/annurev.earth.32.082503.144359 Google Scholar CrossRef Search ADS   Peltier W.R., Argus D.F., Drummond R., 2015. Space geodesy constrains ice age terminal deglaciation: the global ICE-6G_C (VM5a) model, J. geophys. Res. , 120, 450– 487. https://doi.org/10.1002/2014JB011176 Google Scholar CrossRef Search ADS   Petit G., Luzum B., 2010. IERS Conventions (2010), IERS Technical Note 36, Verlagdes Bundesamts für Kartographie und Geodäsie , Frankfurt am Main, Germany. Save H., Bettadpur S., Tapley B.D., 2016. High-resolution CSR GRACE RL05 mascons, J. geophys. Res. , 121, 7547– 7569. https://doi.org/10.1002/2016JB013007 Google Scholar CrossRef Search ADS   Shepherd A. et al.  , 2012. A reconciled estimate of ice-sheet mass balance, Science , 338, 1183– 1189. https://doi.org/10.1126/science.1228102 Google Scholar CrossRef Search ADS PubMed  Stuhne G.R., Peltier W.R., 2015. Reconciling the ICE-6G_C reconstruction of glacial chronology with ice sheet dynamics: the cases of Greenland and Antarctica, J. geophys. Res. , 120, 1841– 1865. https://doi.org/10.1002/2015JF003580 Google Scholar CrossRef Search ADS   Swenson S., Chambers D., Wahr J., 2008. Estimating geocenter variations from a combination of GRACE and ocean model output, J. geophys. Res. , 113, B08410, doi:10.1029/2007JB005338. https://doi.org/10.1029/2007JB005338 Google Scholar CrossRef Search ADS   Tapley B.D., Bettadpur S., Watkins M., Reigber C., 2004. The gravity recovery and climate experiment: mission overview and early results, Geophys. Res. Lett. , 31, L09607, doi:10.1029/2004GL019920. https://doi.org/10.1029/2004GL019920 Google Scholar CrossRef Search ADS   Van Dam T., Wahr J., Lavallée D., 2007. A comparison of annual vertical crustal displacements from GPS and Gravity Recovery and Climate Experiment (GRACE) over Europe, J. geophys. Res. , 112, B03404, doi:10.1029/2006JB004335. https://doi.org/10.1029/2006JB004335 Google Scholar CrossRef Search ADS   Van Dam T., Francis O., Wahr J., Khan S.A., Bevis M., Van Den Broeke M.R., 2017. Using GPS and absolute gravity observations to separate the effects of present-day and Pleistocene ice-mass changes in South East Greenland, Earth planet. Sci. Lett. , 459, 127– 135. https://doi.org/10.1016/j.epsl.2016.11.014 Google Scholar CrossRef Search ADS   Vasskog K., Langebroek P.M., Andrews J.T., Nilsen J.E.Ø., Nesje A., 2015. The Greenland Ice Sheet during the last glacial cycle: current ice loss and contribution to sea-level rise from a palaeoclimatic perspective, Earth-Sci. Rev. , 150, 45– 67. https://doi.org/10.1016/j.earscirev.2015.07.006 Google Scholar CrossRef Search ADS   Velicogna I., Wahr J., 2006. Acceleration of Greenland ice mass loss in spring 2004, Nature , 443, 329– 331. https://doi.org/10.1038/nature05168 Google Scholar CrossRef Search ADS PubMed  Velicogna I., Wahr J., 2013. Time-variable gravity observations of ice sheet mass balance: precision and limitations of the GRACE satellite data, Geophys. Res. Lett. , 40, 3055– 3063. https://doi.org/10.1002/grl.50527 Google Scholar CrossRef Search ADS   Wahr J., Molenaar M., Bryan F., 1998. Time variability of the Earth's gravity field: hydrological and oceanic effects and their possible detection using GRACE, J. geophys. Res. , 103, 30 205–30 229. https://doi.org/10.1029/98JB02844 Google Scholar CrossRef Search ADS   Wahr J., Khan S.A., Van Dam T., Liu L., Van Angelen J.H., Van Den Broeke M.R., Meertens C.M., 2013. The use of GPS horizontals for loading studies, with applications to northern California and southeast Greenland, J. geophys. Res. , 118, 1795– 1806. https://doi.org/10.1002/jgrb.50104 Google Scholar CrossRef Search ADS   Wang S.Y., Chen J.L., Li J., Hu X.G., Ni S., 2016. Geophysical interpretation of GPS loading deformation over western Europe using GRACE measurements, Ann. Geophys. , 59( 5), S0538, doi:10.4401/ag-7058. Wouters B., Bamber J.L., Van Den Broeke M.R., Lenaerts J.T.M., Sasgen I., 2013. Limits in detecting acceleration of ice sheet mass loss due to climate variability, Nat. Geosci. , 6, 613– 616. https://doi.org/10.1038/ngeo1874 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Reconciling GRACE and GPS estimates of long-term load deformation in southern Greenland

, Volume 212 (2) – Feb 1, 2018
12 pages

Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx473
Publisher site
See Article on Publisher Site

### Abstract

Summary We examine vertical load deformation at four continuous Global Positioning System (GPS) sites in southern Greenland relative to Gravity Recovery and Climate Experiment (GRACE) predictions of vertical deformation over the period 2002–2016. With limited spatial resolution, GRACE predictions require adjustment before they can be compared with GPS height time series. Without adjustment, both GRACE spherical harmonic (SH) and mascon solutions predict significant vertical displacement rate differences relative to GPS. We use a scaling factor method to adjust GRACE results, based on a long-term mass rate model derived from GRACE measurements, glacial geography, and ice flow data. Adjusted GRACE estimates show significantly improved agreement with GPS, both in terms of long-term rates and interannual variations. A deceleration of mass loss is observed in southern Greenland since early 2013. The success at reconciling GPS and GRACE observations with a more detailed mass rate model demonstrates the high sensitivity to load distribution in regions surrounding GPS stations. Conversely, the value of GPS observations in constraining mass changes in surrounding regions is also demonstrated. In addition, our results are consistent with recent estimates of GIA uplift (∼4.4 mm yr−1) at the KULU site. Global change from geodesy, Loading of the Earth, Satellite geodesy, Satellite gravity, Space geodetic surveys, Time variable gravity 1 INTRODUCTION The Greenland ice sheet (GrIS) is the second largest mass of ice on Earth and would raise global mean sea level (GMSL) by ∼7 m if completely melted (Gregory et al.2004). With Arctic warming, the GrIS has undergone significant melting, and is regarded as a main contributor to GMSL rise (Shepherd et al.2012). Improved quantification of GrIS mass balance at various spatial and temporal scales and understanding of its climatological and geodynamical impacts continue to be important in understanding climate change. Since the Gravity Recovery and Climate Experiment (GRACE) satellite gravity mission was launched in March 2002 (Tapley et al.2004), investigation of ice mass changes over polar ice sheets and mountain glaciers has entered a new era (e.g. Chen et al.2006; Velicogna & Wahr 2006; Khan et al.2010). By accurately measuring Earth's gravity, GRACE provides a unique monthly measure of global large-scale mass changes (Cazenave & Chen 2010). Early GRACE studies (e.g. Chen et al.2006; Velicogna & Wahr 2006) found accelerated GrIS ice loss since 2005, mostly in East and Southeast Greenland. Khan et al. (2010) and Chen et al. (2011) reported additional losses in northwest Greenland starting in late 2005, with evidence of acceleration since 2008. A recent study (Forsberg et al.2017) found large year-to-year variations during the GRACE era, especially a large 2012 melt event (Nghiem et al.2012). The estimated average GrIS ice loss rate is 265 ± 25 Gt yr−1 (Forsberg et al.2017). Continuous Global Positioning System (GPS) measurements on bedrock adjacent to the GrIS are able to detect surface displacements caused by ice mass unloading in Greenland. A number of studies have demonstrated the value of these load deformation observations (e.g. Khan et al.2007, 2010; Jiang et al.2010; Bevis et al.2012; Nielsen et al.2012; Larson et al.2015; Liu et al.2017). In combination, GRACE and GPS data should provide additional understanding of ice mass changes. Khan et al. (2010) compared GRACE-predicted surface deformation with GPS measurements in Greenland. They compared GPS loading signals and GRACE estimates, and confirmed increasing ice mass loss towards the north. Joud et al. (2017) also compared GRACE-predicted uplift rates with GPS measurements to validate their methods for processing GRACE data in regions with a large glacial isostatic adjustment (GIA) signal. The coarse spatial resolution (several hundred kilometres) of GRACE measurements and spatial leakage effects have been a major limitation in related studies. King et al. (2012) reduced the leakage effects to refine the estimation of Antarctic sea-level contribution, and Mémin et al. (2014) also applied leakage correction in studying snow and ice height change. In load deformation studies, limited resolution, due to the range of available spherical harmonic (SH) coefficients and spatial filtering used to suppress noise, contributes to differences between load deformation estimated from GRACE and GPS data. In some earlier studies GRACE-GPS differences were attributed to errors in GPS data processing (e.g. King et al.2006, van Dam et al.2007). Khan et al. (2010) found that GRACE estimates were smaller than GPS with ratios of GPS to GRACE uplift rates at THU3, KULU and KELY stations, of ∼2.9, 2.7 and 2.3, respectively. Khan et al. (2010) assumed similar ratios for ice mass loss acceleration and trends and used these to adjust GRACE trend and acceleration values for shorter time spans. The scaling factor method has been shown to be useful for correcting spatial leakage in GRACE mass change estimates (e.g. Chen et al.2007, 2017; Landerer & Swenson 2012). Scaling factors, at either regional scales or at grid points, are commonly derived from model-predicted mass change and used to adjust GRACE estimates. However, the scaling factor method relies on accurate models, and while seasonal predictions of available climate model are often good, performance tends to degrade at longer time scales. In this study, we use an improved scaling method to compare GPS-observed and GRACE-predicted uplift rates at 4 GPS sites (KELY, QAQ1, SENU and KULU) in southern Greenland (see Fig. 1) chosen for their good data quality and availability. We examine GPS and GRACE data over the period 2002–2016, long enough to examine both long term rates and interannual variation. The improved scaling factor method is based on a mass rate model that combines independent information from GRACE, geographic (Google Earth) data, and ice dynamic data. Figure 1. View largeDownload slide Map of the Greenland ice sheet, with locations of four selected GPS sites in southern Greenland marked by red dots (revised from Vasskog et al.2015). Figure 1. View largeDownload slide Map of the Greenland ice sheet, with locations of four selected GPS sites in southern Greenland marked by red dots (revised from Vasskog et al.2015). 2 DATA PROCESSING 2.1 GRACE spherical harmonic solutions GRACE time-variable gravity data products include 153 Release 5 (RL05) GSM monthly solutions, covering the period August 2002 to August 2016, provided by the Center for Space Research (CSR) at the University of Texas at Austin (Bettadpur 2012). Each monthly solution consists of fully normalized spherical harmonic (SH) coefficients up to degree and order 60. Atmospheric and oceanic effects have been removed during GRACE data processing through the de-aliasing process (Bettadpur 2012), so that over Greenland, GRACE GSM month-to-month changes mainly reflect ice mass and terrestrial water storage (TWS) effects. Following standard practice, we replace GRACE C20 coefficients by independent estimates from satellite laser ranging (SLR) provided by CSR (Cheng et al.2013), and include degree-1 coefficients estimated using the method of Swenson et al. (2008). GIA (Post Glacial Rebound, PGR) effects are removed using the model of A et al. (2013) (based on ICE-5G ice history and VM2 viscosity profile). To compare GRACE with GPS results, atmospheric and oceanic loading effects (removed from GRACE data during de-aliasing) have to be added back. We add to GSM the GAC products, representing the combined effect of atmosphere and ocean mass variations in GRACE AOD1B de-aliasing products (Flechtner et al.2014). Elastic vertical load deformation dr(θ, ϕ), at colatitude θ and longitude ϕ, can be calculated from GRACE Stokes coefficients as follows (van Dam et al.2007):   \begin{eqnarray} dr\left( {\theta ,\phi } \right) &=& {\rm{R}}\sum\limits_{l \!=\! 0}^\infty {\sum\limits_{m = 0}^l {{W_l}{{\tilde{P}}_{l,m}}\left( {\cos \theta } \right)\left( {{C_{lm}}\cos \left( {m\phi } \right) \!+\! {S_{lm}}\sin \left( {m\phi } \right)} \right)} }\nonumber\\ &&\times \frac{{h^{^\prime}}\!_l}{{1 + {{k^{^\prime}}\!_l}}} \end{eqnarray} (1)where R is the radius of Earth, $${\tilde{P}_{l,m}}$$ are normalized Legendre functions of degree l and order m, Clm and Slm represent the SH (Stokes) coefficients of the gravity field, and $$h_l^\prime$$ and $$k_l^\prime$$ are elastic load Love numbers of degree l, whose values for a Gutenberg–Bullen (G-B) elastic Earth model are provided by Farrell (1972). Wl is the weighting value at each spherical harmonic degree to implement Gaussian smoothing (Jekeli 1981; Wahr et al.1998), with a radius of 300 km in the present study. 2.2 GRACE mascon solutions CSR RL05 GRACE 0.5° × 0.5° mascon solutions for August 2002 to August 2016 are also used (Save et al.2016). Mascon solutions are obtained using several constraints, including separation of land and ocean signals to reduce leakage, and forward modelling of annual and trend signals from ice losses in polar regions and mountain glaciers. Unlike SH solutions, mascon solutions do not require additional filtering to suppress noise (Save et al.2016). Mascon solutions are obtained with C20 coefficients replaced by SLR solutions (Cheng et al.2013), and the same degree-1 terms and GIA corrections (Swenson et al. (2008) and A et al. (2013)). Note that the mascon solutions have restored the GAD products, representing the atmospheric and oceanic effects over the ocean (in contrast to GAC, GAD data set is set to zero at the continents), so we compute the total mascon solution as (CSR mascon solution + GAC − GAD) to reflect the same load signal as SH solution (GSM + GAC). Mascon 0.5° × 0.5° solutions dh(θ, ϕ) at colatitude θ and longitude ϕ, are converted to Stokes coefficients up to degree and order 360, corresponding to about 0.5° resolution at the equator, using eqs (2)–(4),   $${\hat{C}_{lm}}\,{\rm{ = }}\,\frac{{\rm{1}}}{{\rm{R}}} \cdot \frac{{\rm{1}}}{{{\rm{4\pi }}}}\iint_{\sigma }{{{\rm{d}}h\left( {\theta ,\phi } \right){{\tilde{P}}_{l,m}}\left( {\cos \theta } \right)\cos \left( {m\phi } \right){\rm{d}}\sigma }}$$ (2)  $${\hat{S}_{lm}}\,{\rm{ = }}\,\frac{{\rm{1}}}{{\rm{R}}} \cdot \frac{{\rm{1}}}{{{\rm{4\pi }}}}\iint_\sigma {{\rm{d}}h\left( {\theta ,\phi } \right){{\tilde{P}}_{l,m}}\left( {\cos \theta } \right)\sin \left( {m\phi } \right){\rm{d}}\sigma }$$ (3)  $$\left\{ {\begin{array}{@{}*{1}{c}@{}} {{C_{lm}}}\\ {{S_{lm}}} \end{array}} \right\} = \frac{{3{\rho _w}}}{{{\rho _{{\rm{ave}}}}}}\frac{{1 + {{k^{^\prime}}\!_l}}}{{2l + 1}}\left\{ {\begin{array}{@{}*{1}{c}@{}} {{{\hat{C}}_{lm}}}\\ {{{\hat{S}}_{lm}}} \end{array}} \right\}$$ (4)where $${\hat{C}_{{\rm{lm}}}}$$ and $${\hat{S}_{{\rm{lm}}}}$$ are dimensionless coefficients defined in eq. (9) of Wahr et al. (1998), dσ is the area element, which equals to sin (θ) dθdϕ, ρw is the density of water, ρave is the average density of the Earth. Then we substitute Stokes coefficients Clm and Slm into eq. (1) to obtain vertical load deformation. Note that Wl = 1 when using mascon-based SH coefficients, because spatial smoothing is not needed. 2.3 GPS data The Jet Propulsion Laboratory's (JPL’s) GNSS-Inferred Positioning System (GIPSY) solutions of daily GPS height at the four sites (KELY, QAQ1, SENU and KULU, shown in Fig. 1) are available at ftp://sopac-ftp.ucsd.edu/pub/timeseries/measures/ats/. Data at KELY, QAQ1, and KULU span the GRACE era (August 2002 to August 2016), while SENU data are available for May 2008 to August 2016. JPL’s reanalysis orbit and clock products are determined using a set of models in a free-network reference frame, including gravity from Earth, Sun, Moon, and other planets. Models include: the DE421 planetary ephemeris (Folkner et al.2008), the GSPM10 satellite solar pressure model, the IAU06 model for precession and nutation, the IERS2010 tide model (Petit & Luzum 2010), the FES2004 ocean loading model (Lettelier et al.2004), and International GNSS Service satellite and receiver antenna phase centre models (ftp://sideshow.jpl.nasa.gov/pub/JPL_GPS_Timeseries/repro2011b). We remove GIA uplift rates using predictions from A et al. (2013), used also in correcting GRACE data. Uplift rates are −2.25 mm yr−1 at KELY, 2.15 mm yr−1 at QAQ1, 1.83 mm yr−1 at SENU, and −0.76 mm yr−1 at KULU. GIA model prediction error in Greenland is thought to be relatively small (Velicogna & Wahr 2006). However, a recent study (van Dam et al.2017) suggests that the GIA rate at KULU may differ significantly from model predictions. Using GPS observations and superconducting gravimeter measurements, van Dam et al. (2017) show that the GIA uplift rate at KULU is about 4.49 mm yr−1, compared with the (A et al.2013) prediction of −0.76 mm yr−1. This affects the GRACE-GPS uplift comparison, as discussed below. 3 RESULTS 3.1 GRACE and GPS comparison We use GRACE GSM + GAC SH coefficients (up to degree and order 60) to compute vertical displacements using eq. (1) with 300 km Gaussian smoothing at KELY, QAQ1, SENU and KULU. Fig. 2 compares monthly GRACE vertical displacements and GPS height series (blue and red curves). Upward trends at all stations reflect continued mass loss at these sites or nearby. With the exception of QAQ1 (Fig. 2b), GRACE vertical displacements are smaller than GPS, consistent with previous studies (Khan et al.2010). GRACE-GPS differences are believed to be partly attributed to the coarse spatial resolution of GRACE (Khan et al.2010). Due to the truncation of RL05 SH coefficients at degree and order 60 and the applied 300 km Gaussian smoothing, GRACE SH solutions underestimate the mass change and load deformation (Chen et al.2007; Wang et al.2016). Figure 2. View largeDownload slide Comparisons among GRACE (GRC) vertical displacements using SH solutions (green), mascon solutions (blue) and GPS heights (red). A 300 km Gaussian smoothing has been applied to the GRACE SH solutions. Model predicted GIA effects (A et al.2013) have been removed from both GRACE and GPS. Figure 2. View largeDownload slide Comparisons among GRACE (GRC) vertical displacements using SH solutions (green), mascon solutions (blue) and GPS heights (red). A 300 km Gaussian smoothing has been applied to the GRACE SH solutions. Model predicted GIA effects (A et al.2013) have been removed from both GRACE and GPS. GRACE mascon series should offer higher spatial resolution and provide more accurate load deformation estimates than the SH solutions, because they are not subject to the dominant stripe noise in SH solutions, and do not require spatial smoothing. We convert gridded mascon mass change into vertical displacements using eqs (2)–(4) and (1), and superimpose the results in Fig. 2 (green curve). SH representations of mascon solutions to degree and order 360 should retain their spatial resolution. Compared with predictions from the RL05 SH series (limited to degree and order 60), the GRACE mascon predictions are more similar to GPS height variations at KELY, SENU and KULU as shown in Figs 2(a), (c) and (d), but differences are still relatively large. However, the QAQ1 series (Fig. 2b) differs markedly from the other sites. The GPS uplift rate at QAQ1 is smaller than the GRACE prediction, and GRACE mascon agreement with GPS is poor relative to the SH series. This is discussed further below. Least square fit linear trends for the four sites are given in Table 1 (annual and semi-annual signals are fit simultaneously, but not shown). GPS trend uncertainty is estimated from GPS formal error (provided with the GPS height series). GRACE trend uncertainty is based upon least squares fitting error. Large differences among vertical rates from GPS, GRACE SH, and GRACE mascon estimates are evident in Table 1. We compare in Fig. 3 GRACE SH and mascon vertical rate predictions over Greenland. GRACE mascon rates are significantly larger and show better spatial resolution. The largest vertical displacements are along southeastern and western coastal regions. However, comparisons with GPS observations discussed above indicate that neither GRACE estimate is likely to agree very well with GPS observations, without further processing of GRACE data. Figure 3. View largeDownload slide Vertical displacement rates (in mm yr−1) computed from: (a) GRACE GSM + GAC Stokes coefficients with 300 km Gaussian smoothing; (b) GRACE mascon solutions (0.5° × 0.5°) after adding back the contributions of atmospheric at the continents. GIA correction is applied to both GRACE SH and mascon solutions. Figure 3. View largeDownload slide Vertical displacement rates (in mm yr−1) computed from: (a) GRACE GSM + GAC Stokes coefficients with 300 km Gaussian smoothing; (b) GRACE mascon solutions (0.5° × 0.5°) after adding back the contributions of atmospheric at the continents. GIA correction is applied to both GRACE SH and mascon solutions. Table 1. Predicted uplift rates using GRACE and GPS data (units: mm yr−1). UW means unweighted model, W means weighted model.   KELY  QAQ1  SENUa  KULU  GPS  7.70 ± 0.13  2.51 ± 0.06 (2.69 ± 0.15)b  8.06 ± 0.25  9.88 ± 0.06 [4.63 ± 0.06]c  GRACE SH  4.51 ± 0.08  3.23 ± 0.05 (3.39 ± 0.12) b  3.54 ± 0.13  4.15 ± 0.05  GRACE mascon  5.50 ± 0.10  4.71 ± 0.06  6.10 ± 0.18  6.54 ± 0.06  GRACE restored (UW)  6.85 ± 0.11  3.91 ± 0.06  8.04 ± 0.25  4.31 ± 0.05  GRACE restored (W)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05    KELY  QAQ1  SENUa  KULU  GPS  7.70 ± 0.13  2.51 ± 0.06 (2.69 ± 0.15)b  8.06 ± 0.25  9.88 ± 0.06 [4.63 ± 0.06]c  GRACE SH  4.51 ± 0.08  3.23 ± 0.05 (3.39 ± 0.12) b  3.54 ± 0.13  4.15 ± 0.05  GRACE mascon  5.50 ± 0.10  4.71 ± 0.06  6.10 ± 0.18  6.54 ± 0.06  GRACE restored (UW)  6.85 ± 0.11  3.91 ± 0.06  8.04 ± 0.25  4.31 ± 0.05  GRACE restored (W)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05  aThe uplift rates at SENU are computed using data from May 2008 to August 2016, and others are from August 2002 to August 2016. bThe uplift rates at QAQ1 in the parentheses are computed using data from May 2008 to August 2016. cThe GPS uplift rate at KULU in the square brackets is computed using GIA uplift rate predicted by van Dam et al. (2017). View Large 3.2 Correcting GRACE spatial leakage Due to the truncation of GRACE SH coefficients at degree and order 60 and additional spatial filtering (300 km Gaussian), GRACE predictions of vertical displacements are limited to long-wavelength features. While GRACE mascon solutions offer qualitative improvements in spatial resolution, they remain subject to leakage effects due to fundamental GRACE spatial resolution limits, and it is easier to quantitatively address resolution limits of SH solutions. If the spatial spectrum (spatial distribution) of mass change is approximately known, either from climate models or from observations, scaling factors can be estimated to adjust GRACE values to reduce bias associated spatial leakage as discussed above (Chen et al.2007, 2017; Landerer & Swenson 2012; Wang et al.2016). Recognizing that long-term ice mass changes are dominantly along the edges of the GrIS, we construct an ice mass change model based on GRACE mascon solutions, published ice flow velocity data (Khan et al.2015), and known glacier locations (Google Earth). The scaling factor at a grid point is the ratio between model-predicted load deformation and the prediction after SH truncation and spatial filtering. Scaling factors are then used to adjust GRACE SH predictions of vertical motion at the four GPS sites. The procedure includes the following steps: From GRACE mascon gridded equivalent water thickness (EWT) data, we subtract the GAD products. The residual EWT contains ice and TWS mass change effects, similar to GRACE GSM data (after GIA correction). The mass rate of residual EWT at each grid from August 2002 to August 2016 is found by least squares with the result shown in Fig. 4(a). We divide the GrIS and surrounding areas into 12 regions (Fig. 4b), and find the total mass rate within each using GRACE mascon rates (Fig. 4a). Mass rates of the 12 regions are in Table 2 (units of gigatonne per year, Gt yr−1). Note that the sum of area 1 to 10 represents the total mass loss of Greenland, that is, −252.5 Gt yr−1, which is consistent with previous publications (e.g. Velicogna & Wahr 2013; Wouters et al.2013). We construct a detailed mass rate model along the periphery of southern Greenland glaciers (to latitudes 72.5°N) by assigning mass rates in the 12 regions into narrower bands (Fig. 5b), to represent dominant coastal and central Greenland locations associated with long-term changes. Within each narrower region, mass change is uniformly distributed. The Fig. 5(b) mass rate model is on a 0.1° × 0.1° grid, with a width of 0.5° (∼25 km) along glacier edges, based on Google Earth glacier geography. In regions far from the GPS stations, a lower resolution model (0.5° × 0.5°) is used (Fig. 5a). An additional mass rate model is constructed, in which areas near the GPS stations are weighted according to observed ice flow velocity data (fig. 8 in Khan et al.2015). Weighted solutions retain the same total mass rate in each subregion. We define the weight as approximately proportional to the velocity (e.g. we define the light yellow colour region in their figure as weight 1.0), and then roughly allocate 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6 or 1.8 as gridded weight according to the colour of corresponding grid in the ice flow velocity map. The weighted (W) model is shown in Fig. 5(c), compared with the unweighted (UW) model in Fig. 5(b). The mass rate models are represented in an SH expansion to degree and order 1800, and from these vertical displacements calculated using both the full range Stokes coefficients, and Stokes coefficients limited to degree and order 60 Stokes with additional 300 km Gaussian smoothing. At each grid point, the ratio between the two displacements is the scaling factor used to adjust GRACE estimates. We apply scaling factors to vertical displacement rates predicted by GRACE SH products, and add back GAC and GIA models. These can now be compared with GPS observations for both UW and W models in Figs 6 and 7 and Table 1. Figure 4. View largeDownload slide (a) GRACE mascon mass rates (in units of cm yr−1 of equivalent water thickness change) after the GAD and GIA models are removed; (b) the 12 shaded areas used to integrate the total mass rate within each of the regions in Fig. 4(a). Figure 4. View largeDownload slide (a) GRACE mascon mass rates (in units of cm yr−1 of equivalent water thickness change) after the GAD and GIA models are removed; (b) the 12 shaded areas used to integrate the total mass rate within each of the regions in Fig. 4(a). Figure 5. View largeDownload slide (a) The 0.5° × 0.5° mass rate model approximates ice loss locations and magnitudes in areas distant from the GPS stations; (b) the unweighted 0.1° × 0.1° mass rate model that confines ice loss to locations along the periphery of the ice sheet, with uniform distribution of ice loss within each sub-region near the GPS stations; (c) the weighted 0.1° × 0.1° mass rate model. Weights are based on published ice flow velocity data. The mass rates in the far and near field mass rate models are based on Fig. 4 and Table 2. Figure 5. View largeDownload slide (a) The 0.5° × 0.5° mass rate model approximates ice loss locations and magnitudes in areas distant from the GPS stations; (b) the unweighted 0.1° × 0.1° mass rate model that confines ice loss to locations along the periphery of the ice sheet, with uniform distribution of ice loss within each sub-region near the GPS stations; (c) the weighted 0.1° × 0.1° mass rate model. Weights are based on published ice flow velocity data. The mass rates in the far and near field mass rate models are based on Fig. 4 and Table 2. Figure 6. View largeDownload slide Comparisons of monthly GPS height, GRACE spherical harmonic (SH) predicted vertical displacements and scaling factor adjusted GRACE series using the weighted mass rate model (Figs 3a and c). Uplift rates are shown using the same curve colour. Note: we use two GIA uplift rates at KULU predicted by A et al. (2013) (GIA) and van Dam et al. (2017) (GIA*). GRACE-predicted uplift rates at SENU in the parentheses are computed using the same time period as GPS (2008.05–2016.08). Figure 6. View largeDownload slide Comparisons of monthly GPS height, GRACE spherical harmonic (SH) predicted vertical displacements and scaling factor adjusted GRACE series using the weighted mass rate model (Figs 3a and c). Uplift rates are shown using the same curve colour. Note: we use two GIA uplift rates at KULU predicted by A et al. (2013) (GIA) and van Dam et al. (2017) (GIA*). GRACE-predicted uplift rates at SENU in the parentheses are computed using the same time period as GPS (2008.05–2016.08). Figure 7. View largeDownload slide Comparison of GPS heights and scaling factor adjusted GRACE vertical displacements using the weighted mass rate model after 13-month moving average filtering. The two dash lines mark dates December 2009 and March 2013. Figure 7. View largeDownload slide Comparison of GPS heights and scaling factor adjusted GRACE vertical displacements using the weighted mass rate model after 13-month moving average filtering. The two dash lines mark dates December 2009 and March 2013. Table 2. Total mass rates of the 12 regions in Fig. 4(b) from integrating observed mass rates from the mascon solutions (Fig. 4a). Area #  1  2  3  4  5  6  7  8  9  10  11  12  Gr Mass  Total  Mass Rate (Gt yr−1)  −18.2  −65.3  −7.8  −7.6  −29.2  −44.4  −57.6  −10.8  −15.8  4.2  −71.4  −7.6  −252.5  −331.5  Area #  1  2  3  4  5  6  7  8  9  10  11  12  Gr Mass  Total  Mass Rate (Gt yr−1)  −18.2  −65.3  −7.8  −7.6  −29.2  −44.4  −57.6  −10.8  −15.8  4.2  −71.4  −7.6  −252.5  −331.5  View Large 3.3 Comparing leakage-corrected GRACE with GPS vertical rates We focus on interannual and longer-term variations, because scaling factors are based on GRACE mass rates, associated mainly with mass loss along the periphery of the GrIS. Load variations at shorter time scales (e.g. seasonal) would be expected to have very different spatial distributions due to dominant effects of snowfall and TWS change. Figs 6(a) and (c) show very good agreement between adjusted GRACE and GPS uplift rates at KELY and SENU. At KELY the GPS rate of 7.70 ± 0.13 mm yr−1 is in excellent agreement with the GRACE (GRC) prediction of 7.10 ± 0.12 mm yr−1. Leakage correction steps have raised the GRACE rate from the previous value of 4.51 ± 0.08 mm yr−1. With a shorter time span at SENU (May 2008 to August 2016), GPS and GRACE uplift rates are in excellent agreement 8.06 ± 0.25 mm yr−1 and 8.05 ± 0.25 mm yr−1, respectively. The scaling factor leakage correction has raised the rate considerably, from the former value of 3.54 ± 0.13. KELY is about 40 km from the periphery of the southwest GrIS. The model includes north-south distributed mass loss along the periphery. The weighted model includes mass loss at the nearby land-terminating glacial outlets (Fig. 5c). In this region, ice loss patterns and locations are relatively well known, so the scaling factor method is expected to be especially effective. SENU is just on the periphery of the ice sheet (Fig. 1), making deformation sensitive to nearby mass losses. Consequences of SH truncation and spatial filtering are expected to be particularly large at SENU (Chen et al.2015), but with well-constrained load locations, leakage correction should be effective. At QAQ1 and KULU GRACE-GPS differences persist, especially at KULU. QAQ1 is the only site of the four whose GPS uplift rate is smaller than GRACE prediction (before or after leakage correction). SENU and QAQ1 are only about 70 km apart (Fig. 1), so it is interesting to compare the two sites. We compute uplift rates at QAQ1 for the period of SENU GPS data availability (May 2008 to August 2016) with results in parentheses in Table 1. GRACE SH solutions at these two sites have similar rates (3.39 ± 0.12 versus 3.54 ± 0.13 mm yr−1), while GPS uplift rates are quite different (2.69 ± 0.15 versus 8.06 ± 0.25 mm yr−1). This indicates that mass loss is nearest SENU causing larger elastic load deformation evident with GPS, and with a much smaller effect at QAQ1 (∼70 km away). GRACE SH solutions, with their limited resolution, apparently cause the GRACE load estimate to be distributed over a broader region that includes QAQ1. This indicates that the leakage correction was not fully effective, as discussed below. At KULU, GPS and GRACE uplift rates show the largest difference among the four stations. The GRACE rate (even with leakage correction) is less than half the GPS value (4.30 ± 0.05 versus 9.88 ± 0.06 mm yr−1). While the geographical setting around KULU is relatively complicated, the mass rate model may not be adequate, but we suspect that the large difference is not mainly related to this. As noted earlier, the GIA uplift rate at KULU has been proposed to be substantially different from model predictions (4.49 mm yr−1 by van Dam et al.2017). Khan et al. (2016) also reported a similar GIA rate at KULU (∼4.4 mm yr−1). If we adopt the van Dam et al. (2017) rate at KULU, then the GPS ice load rate becomes 4.63 ± 0.06 mm yr−1, very close to the GRACE prediction (4.30 ± 0.05 mm yr−1) (Fig. 6d). We have used a different GIA value for the GPS data correction only. For GRACE data, we still rely on GIA model prediction (A et al.2013). Furthermore, GIA correction in GRACE estimate (GIA mass rate converted to deformation) is much smaller (∼0.1 mm yr−1), which unlikely has any notable impact on our qualitative analysis. KULU is far from the outlet glaciers relative to KELY and SENU, and is unlikely to experience a much larger contemporary elastic unloading rate (9.88 mm yr−1) than the other two sites. We apply a 13-month moving average filter (with half weight at first and last values) to both GPS and GRACE monthly time series, and compare smoothed series in Fig. 7. Errors are estimated as in Table 1 (described in the 3rd paragraph of Section 3.1). Leakage-corrected GRACE and GPS series show similar trends and interannual changes. We estimate linear rates for three different time spans (August 2002 to December 2009, December 2009 to March 2013, March 2013 to August 2016). All eight series (4 GPS + 4 GRACE) show increasing uplift rates after December 2009, and decreasing rates after March 2013. This is further evidence of accelerated ice loss since 2010 (Forsberg et al.2017). Both GPS and GRACE series indicate a mass loss slow-down since 2013. Results from unweighted and weighted mass rate models are shown in the last two rows of Table 1. 4 DISCUSSION AND CONCLUSIONS Analysis of vertical deformations at four continuous GPS sites, KELY, SENU, QAQ1 and KULU in southern Greenland shows that it is possible to reconcile GPS estimates of elastic loading with GRACE time-variable gravity solutions after GRACE leakage effect is appropriately addressed. The scaling factor method used in this study is based on ice mass change model derived from GRACE observations, which is an improvement from the posteriori factors (not validated through other independent observations) applied in Khan et al. (2010), also from the scaling factors derived from deficient land surface model data (TWS model without groundwater information) by Wang et al. (2016). The good agreements between GPS and restored GRACE at KELY, SENU and QAQ1 indicate that the long-term variations in the three GPS sites are mainly caused by continuous ice mass loss in the surrounding regions. Some uncertainties may affect our comparisons, including uncertainties in GRACE estimates, GIA models, our mass rate models and GPS observations. GRACE estimates are limited by coarse spatial resolution, and subject to bias due to spatial leakage, resulting in apparently large discrepancies with GPS estimates. GRACE mascon solutions offer qualitative improvements in spatial resolution relative to SH solutions (Fig. 3), but correcting leakage bias in mascon solutions is not as straightforward as for SH solutions (Chen et al.2017). The scaling factor method is used to correct GRACE SH data, with careful attention to constructing a long-term mass rate model that combines GRACE measurements and geographical information on glacier location and ice flow velocity. Leakage-corrected GRACE vertical rate predictions then agree remarkably well with GPS observations at KELY and SENU. A relatively large difference at KULU is probably due to an erroneous GIA uplift rate, and when an estimate suggested in an earlier study (van Dam et al.2017) is used, GPS and GRACE rates agree very well. At QAQ1 differences are thought to be due to errors in both the GIA model and the mass rate model related to complex geography. GRACE-GPS results also agree well at interannual time scales, with higher uplift rates during December 2009 to March 2013 at all four sites, suggesting accelerated ice loss during this period. A decrease in uplift rates at the four sites since early 2013, suggests deceleration of mass loss in southern Greenland (Fig. 7). Our finding of interannual variations of GrIS ice mass changes are consistent with results from Forsberg et al. (2017). Excellent agreement between GPS and GRACE estimates gives us confidence in the ability to geodetically monitor ice mass change. The care required to correct for spatial leakage in the GRACE estimates illustrates the sensitivity of load deformation to details of load changes in the surrounding regions. A corollary is that GPS observations provide unique and important constraints on regional mass changes. The present study provides support, though indirect, for recent GIA uplift rate estimates at KULU (Khan et al.2016; van Dam et al.2017, both near 4.4 mm yr−1) which differ from previous GIA model predictions (e.g. Peltier 2004, 2015; A et al.2013; Stuhne & Peltier 2015). GPS uplift rates at KULU are large, near 9.88 mm yr−1, and show less seasonal and interannual variability than the other three sites. KULU is relatively far from major outlet glaciers, so it is sensible that seasonal effects are smaller and that a larger portion of the observed rate is due to GIA, rather than elastic loading. Good GIA uplift rate estimates are essential in interpreting and comparing GPS and GRACE uplift results, but uncertainty in GIA estimates is difficult to determine, mainly due to the lack of long-term in situ measurements. Besides probable uncertainties in GRACE estimates and GIA models, parametrization of the mass rate models also contributes to uncertainty in our conclusions. The current width of mass rate model is 0.5°, which is approximately 25 km at this latitude zone. We try other widths to test the sensitivity and show the results in Table 3. In the experiment, we broaden the width of the mass rate model for glaciers in southwest Greenland (as shown in Fig. 5c) eastward to 0.6° and 1°, and explore that the KELY is the most affected sites in terms of uplift rates, which demonstrates the high deformation sensitivity to near-field loading source. Indeed, Wahr et al. (2013) showed their experiment results between vertical load deformation at a point and the distance from loading disc centre in their fig. 1(a), and Argus et al. (2014) also reported the similar results in their fig. 1. Table 3. Predicted uplift rates of restored GRACE data using weighted model with different widths of the mass rate model for glaciers in southwest Greenland as shown in Fig. 5(c) (units: mm yr−1).   KELY  QAQ1  SENU  KULU  GRACE restored (0.5°)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05  GRACE restored (0.6°)  6.98 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.30 ± 0.05  GRACE restored (1°)  6.55 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.28 ± 0.05    KELY  QAQ1  SENU  KULU  GRACE restored (0.5°)  7.10 ± 0.12  3.66 ± 0.05  8.05 ± 0.25  4.30 ± 0.05  GRACE restored (0.6°)  6.98 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.30 ± 0.05  GRACE restored (1°)  6.55 ± 0.11  3.62 ± 0.05  7.75 ± 0.24  4.28 ± 0.05  View Large Furthermore, the differences between our unweighted and weighted model results reinforce our conclusion of the effect of load location and shape of glaciers. The weighted model provides more detailed definition of the load location, especially in regions where ice mass loss concentrates in large glaciers other than distributes evenly along the edge. The uplift rate at KELY, located near the Russell Glacier, is notably affected by the weighted mass rate model, and the same is for QAQ1. The sensitivity to load location means that details of glacier location and mass change magnitudes can be important. Additionally, the ice velocity data from Khan et al. (2015) is based on observations between 2000/01 and 2005/06, which is only a first order approximation of ice velocity during GRACE period. We have also compared JPL GPS series used here with GPS solutions from the Scripps Orbit and Permanent Array Center (SOPAC) and MEaSUREs combination solutions (Bock & Webb 2012, ftp://sopac-ftp.ucsd.edu/pub/timeseries/measures/ats/). Fig. 8 plots the three GPS height series at the four selected sites (SOPAC provides KELY and QAQ1 data only), and Table 4 shows their annual amplitudes, annual phases and linear trend signals. Slight differences among the three solutions are found for both annual and linear trend signals. There are a few factors causing such discrepancies, partly due to the differences in data processing strategies between GIPSY (used in JPL) and GAMIT (used in SOPAC), and partly due to the errors during the st_filter combination (used in MEaSUREs combination solution). There will always misfit between JPL and SOPAC solutions, especially one is from Precise Point Positioning (completely depend on FIXED orbit and clock with no tie to any other station) and the other is from network solution with orbit also participating in the adjustment (Peng Fang, personal communication, 2017). The further discussion of GPS time series may be beyond the scope of this article, and need to be presented in a more specific article. Figure 8. View largeDownload slide Comparisons of monthly GPS heights from JPL, SOPAC and MEaSUREs solutions. GIA from A et al. (2013) has been applied to each series. Figure 8. View largeDownload slide Comparisons of monthly GPS heights from JPL, SOPAC and MEaSUREs solutions. GIA from A et al. (2013) has been applied to each series. Table 4. Annual amplitude, annual phase and linear trend corresponding with Fig. 8. The annual phase is defined as φ in cos (2π(t − t0) + φ), where t0 refers to the start of the year (e.g. 2004.0). The uncertainties are estimated from GPS formal error and least squares fitting error.     KELY  QAQ1  SENU  KULU  Annual amplitude (mm)    JPL  5.59 ± 0.72  4.94 ± 0.35  9.09 ± 0.86  4.60 ± 0.35    SOPAC  5.16 ± 0.73  4.81 ± 0.43  /  /    MEaSUREs  6.23 ± 0.76  5.19 ± 0.38  9.53 ± 0.95  5.55 ± 0.47  Annual phase (degree)    JPL  332 ± 7  344 ± 4  324 ± 5  330 ± 4    SOPAC  309 ± 8  316 ± 5  /  /    MEaSUREs  276 ± 7  278 ± 4  294 ± 6  280 ± 5  Trend (mm yr−1)    JPL  7.70 ± 0.13  2.51 ± 0.06  8.06 ± 0.25  9.88 ± 0.06    SOPAC  7.50 ± 0.13  2.31 ± 0.08  /  /    MEaSUREs  7.96 ± 0.14  2.32 ± 0.07  7.25 ± 0.26  8.22 ± 0.08      KELY  QAQ1  SENU  KULU  Annual amplitude (mm)    JPL  5.59 ± 0.72  4.94 ± 0.35  9.09 ± 0.86  4.60 ± 0.35    SOPAC  5.16 ± 0.73  4.81 ± 0.43  /  /    MEaSUREs  6.23 ± 0.76  5.19 ± 0.38  9.53 ± 0.95  5.55 ± 0.47  Annual phase (degree)    JPL  332 ± 7  344 ± 4  324 ± 5  330 ± 4    SOPAC  309 ± 8  316 ± 5  /  /    MEaSUREs  276 ± 7  278 ± 4  294 ± 6  280 ± 5  Trend (mm yr−1)    JPL  7.70 ± 0.13  2.51 ± 0.06  8.06 ± 0.25  9.88 ± 0.06    SOPAC  7.50 ± 0.13  2.31 ± 0.08  /  /    MEaSUREs  7.96 ± 0.14  2.32 ± 0.07  7.25 ± 0.26  8.22 ± 0.08  View Large Disagreement between GPS and GRACE results at QAQ1 may be due to a combination of errors in GIA estimates, GPS data, and mass rate model or ice velocity weighting. Similarities in GPS results at QAQ1 and SENU (two adjacent sites) suggests that errors in GIA or GRACE results may be the source of disagreement. QAQ1 is located at the juncture between southeast and southwest Greenland, so our formulation of the mass rate model from mascon rates (Figs 4b and 5b and c) may produce relatively poor results near QAQ1. Additional observations, for example from satellite altimetry, airborne altimetry, absolute gravimetry, satellite gravimetry, and surface mass balance models may be the basis for an improved mass rate model in the future. We also take into account possible insufficiencies due to the use of an outdated GIA model (A et al.2013, based on the ICE-5G ice history) as well as load Love number (Farrell 1972, based on the G-B Earth model). When using the new ICE-6G GIA model (Peltier 2015), the difference of GPS uplift rates at QAQ1 is ∼0.36 mm yr−1, as compared to the value based on the ICE-5G model (i.e. 2.15 versus 2.51 mm yr−1). When using load Love numbers based on the Preliminary Reference Earth Model (Han & Wahr 1995), the estimated uplift rate at QAQ1 is 3.67 mm yr−1, which is virtually the same as the result (3.66 mm yr−1) based on the G-B Earth model (Farrell 1972) (see Table 1). While long-term ice loss dominates GrIS mass changes, there are also significant short-term fluctuations. Thus long-term continuous observations are needed to understand long-term changes in the GrIS. Deployment of the Greenland GPS Network (GNET; Bevis et al.2009) and continuing satellite gravity observations with the GRACE Follow-On mission, are important elements of this longer-term mission to understand GrIS mass balance. Acknowledgements We are grateful to two anonymous reviewers for their insightful comments. We thank the CSR and NASA’s PODAAC for providing GRACE data products, the JPL and SOPAC for providing GPS data products. We appreciate the discussions with Dr Peng Fang at University of California at San Diego on the GPS solutions. SYW is supported by the UCAS Joint PhD Training Program (UCAS [2015]37). CRW has received support from the Geology Foundation of the University of Texas. JL is supported by the Natural Science Foundation of Shanghai (No.17ZR1435600), the Open Fund of Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University (16-01-05), and the National Key Research and Development Program of China (2016YFB0501405). REFERENCES A, G., Wahr J., Zhong S., 2013. Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to glacial isostatic adjustment in Antarctica and Canada, Geophys. J. Int. , 192( 2), 557– 572. https://doi.org/10.1093/gji/ggs030 Google Scholar CrossRef Search ADS   Argus D.F., Fu Y., Landerer F.W., 2014. Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion, Geophys. Res. Lett. , 41( 6), 1971– 1980. https://doi.org/10.1002/2014GL059570 Google Scholar CrossRef Search ADS   Bettadpur S., 2012. Gravity recovery and climate experiment, Level-2 gravity field product user handbook (Rev. 3.0, May 29, 2012), GRACE 327-734 (CSR-GR-03-01) , UTCSR, Texas, USA. Bevis M.G., Kendrick E.C., Brown A.K., Khan S.A., Knudsen P., Madsen F., Wahr J., Willis M.J., 2009. Greenland GPS Network: crustal oscillations and seasonal ice mass fluctuations, EOS, Trans. Am. geophys. Un. , 90( 52), Fall Meet. Suppl., Abstract G43B-0728. Bevis M. et al.  , 2012. Bedrock displacements in Greenland manifest ice mass variations, climate cycles and climate change, Proc. Natl. Acad. Sci. USA , 109( 30), 11 944–11 948. https://doi.org/10.1073/pnas.1204664109 Google Scholar CrossRef Search ADS   Bock Y., Webb F.H., 2012. MEaSUREs Solid Earth Science ESDR System , Digital Media. Cazenave A., Chen J., 2010. Time-variable gravity from space and present-day mass redistribution in the Earth system, Earth planet. Sci. Lett. , 298( 3–4), 263– 274. https://doi.org/10.1016/j.epsl.2010.07.035 Google Scholar CrossRef Search ADS   Chen J.L., 2006. Satellite gravity measurements confirm accelerated melting of Greenland ice sheet, Science , 313( 5795), 1958– 1960. https://doi.org/10.1126/science.1129007 Google Scholar CrossRef Search ADS PubMed  Chen J.L., Wilson C.R., Famiglietti J.S., Rodell M., 2007. Attenuation effect on seasonal basin-scale water storage changes from GRACE time-variable gravity, J. Geod. , 81( 4), 237– 245. https://doi.org/10.1007/s00190-006-0104-2 Google Scholar CrossRef Search ADS   Chen J.L., Wilson C.R., Tapley B.D., 2011. Interannual variability of Greenland ice losses from satellite gravimetry, J. geophys. Res. , 116, B07406, doi:10.1029/2010JB007789. https://doi.org/10.1029/2010JB007789 Google Scholar CrossRef Search ADS   Chen J.L., Wilson C.R., Li J., Zhang Z., 2015. Reducing leakage error in GRACE-observed long-term ice mass change: a case study in West Antarctica, J. Geod. , 89( 9), 925– 940. https://doi.org/10.1007/s00190-015-0824-2 Google Scholar CrossRef Search ADS   Chen J.L., Wilson C.R., Tapley B.D., Save H., Cretaux J.F., 2017. Long-term and seasonal Caspian Sea level change from satellite gravity and altimeter measurements, J. geophys. Res. , 192, 557– 517. Cheng M., Tapley B.D., Ries J.C., 2013. Deceleration in the Earth's oblateness, J. geophys. Res. , 118( 2), 740– 747. https://doi.org/10.1002/jgrb.50058 Google Scholar CrossRef Search ADS   Farrell W.E., 1972. Deformation of the Earth by surface loads, Rev. Geophys. , 10( 3), 761– 797. https://doi.org/10.1029/RG010i003p00761 Google Scholar CrossRef Search ADS   Flechtner F., Dobslaw H., Fagiolini E., 2014. Gravity recovery and climate experiment, AOD1B product description document for product Release 05 (Rev. 4.2, May 20, 2014), GRACE 327-750 (GR-GFZ-AOD- 0001), GFZ Potsdam, Potsdam, Germany, p. 33. Folkner W.M., Williams J.G., Boggs D.H., 2008. The Planetary and lunar ephemeris DE 421, JPL memorandum IOM 343R-08-003. Forsberg R., Sorensen L., Simonsen S., 2017. Greenland and antarctica ice sheet mass changes and effects on global sea level, Surv. Geophys. , 38( 1), 89– 104. https://doi.org/10.1007/s10712-016-9398-7 Google Scholar CrossRef Search ADS   Gregory J.M., Huybrechts P., Raper S.C.B., 2004. Climatology: threatened loss of the Greenland ice-sheet, Nature , 428( 6983), 616, doi:10.1038/428616a. https://doi.org/10.1038/428616a Google Scholar CrossRef Search ADS PubMed  Han D., Wahr J., 1995. The viscoelastic relaxation of a realistically stratified earth, and a further analysis of postglacial rebound, Geophys. J. Int. , 120( 2), 287– 311. Google Scholar CrossRef Search ADS   Jekeli C., 1981. Alternative methods to smooth the Earth's gravity field, Report No. 327, Reports of the Department of Geodetic Science and Surveying , Columbus: The Ohio State University. Jiang Y., Dixon T.H., Wdowinski S., 2010. Accelerating uplift in the North Atlantic region as an indicator of ice loss. Nat. Geosci. , 3( 6), 404– 407. https://doi.org/10.1038/ngeo845 Google Scholar CrossRef Search ADS   Joud M.S.S., Sjöberg L.E., Bagherbandi M., 2017. Use of GRACE data to detect the present land uplift rate in Fennoscandia, Geophys. J. Int. , 209( 2), 909– 922. https://doi.org/10.1093/gji/ggx063 Google Scholar CrossRef Search ADS   Khan S.A., Wahr J., Stearns L.A., Hamilton G.S., Van Dam T., Larson K.M., Francis O., 2007. Elastic uplift in southeast Greenland due to rapid ice mass loss, Geophys. Res. Lett. , 34, L21701, doi:10.1029/2007GL031468. https://doi.org/10.1029/2007GL031468 Google Scholar CrossRef Search ADS   Khan S.A., Wahr J., Bevis M., Velicogna I., Kendrick E., 2010. Spread of ice mass loss into northwest Greenland observed by GRACE and GPS, Geophys. Res. Lett. , 37, L06501, doi:10.1029/2010GL042460. https://doi.org/10.1029/2010GL042460 Google Scholar CrossRef Search ADS   Khan S.A., Aschwanden A., Bjørk A.A., Wahr J., Kjeldsen K.K., Kjær K.H., 2015. Greenland ice sheet mass balance: a review, Rep. Prog. Phys. , 78( 4), doi:10.1088/0034-4885/78/4/046801. https://doi.org/10.1088/0034-4885/78/4/046801 Khan S.A. et al., 2016. Geodetic measurements reveal similarities between post-Last Glacial Maximum and present-day mass loss from the Greenland ice sheet, Sci. Adv. , 2( 9), doi:10.1126/sciadv.1600931. https://doi.org/10.1126/sciadv.1600931 King M., Moore P., Clarke P., Lavallée D., 2006. Choice of optimal averaging radii for temporal GRACE gravity solutions, a comparison with GPS and satellite altimetry, Geophys. J. Int. , 166( 1), 1– 11. https://doi.org/10.1111/j.1365-246X.2006.03017.x Google Scholar CrossRef Search ADS   King M.A., Bingham R.J., Moore P., Whitehouse P.L., Bentley M.J., Milne G.A., 2012. Lower satellite-gravimetry estimates of Antarctic sea-level contribution, Nature , 491( 7425), 586– 589. https://doi.org/10.1038/nature11621 Google Scholar CrossRef Search ADS PubMed  Landerer F.W., Swenson S.C., 2012. Accuracy of scaled GRACE terrestrial water storage estimates, Water Resour. Res. , 48, W04531, doi:10.1029/2011WR011453. https://doi.org/10.1029/2011WR011453 Google Scholar CrossRef Search ADS   Larson K.M., Wahr J., Munneke P.K., 2015. Constraints on snow accumulation and firn density in Greenland using GPS receivers, J. Glaciol. , 61, 101– 114. https://doi.org/10.3189/2015JoG14J130 Google Scholar CrossRef Search ADS   Lettelier T., Lyard F., Lefebre F., 2004. The new global tidal solution: FES2004, in Proceeding of the Ocean Surface Topography Science Team Meeting , St. Petersburg, Florida. Liu L., Khan S.A., Van Dam T., Ma J.H.Y., Bevis M., 2017. Annual variations in GPS-measured vertical displacements near Upernavik Isstrøm (Greenland) and contributions from surface mass loading, J. geophys. Res. , 122, 677– 691. https://doi.org/10.1002/2016JB013494 Google Scholar CrossRef Search ADS   Mémin A., Flament T., Rémy F., Llubes M., 2014. Snow- and ice-height change in Antarctica from satellite gravimetry and altimetry data, Earth planet. Sci. Lett. , 404, 344– 353. https://doi.org/10.1016/j.epsl.2014.08.008 Google Scholar CrossRef Search ADS   Nghiem S.V. et al.  , 2012. The extreme melt across the Greenland ice sheet in 2012, Geophys. Res. Lett. , 39, L20502, doi:10.1029/2012GL053611. https://doi.org/10.1029/2012GL053611 Google Scholar CrossRef Search ADS   Nielsen K., Khan S.A., Korsgaard N.J., Kjær K.H., Wahr J., Bevis M., Stearns L.A., Timm L.H., 2012. Crustal uplift due to ice mass variability on Upernavik Isstrøm, west Greenland, Earth planet. Sci. Lett. , 353–354, 182– 189. https://doi.org/10.1016/j.epsl.2012.08.024 Google Scholar CrossRef Search ADS   Peltier W.R., 2004. Global glacial isostasy and the surface of the Ice-Age Earth: The ICE-5G (VM2) Model and GRACE, Annu. Rev. Earth Planet. Sci. , 32, 111– 149. https://doi.org/10.1146/annurev.earth.32.082503.144359 Google Scholar CrossRef Search ADS   Peltier W.R., Argus D.F., Drummond R., 2015. Space geodesy constrains ice age terminal deglaciation: the global ICE-6G_C (VM5a) model, J. geophys. Res. , 120, 450– 487. https://doi.org/10.1002/2014JB011176 Google Scholar CrossRef Search ADS   Petit G., Luzum B., 2010. IERS Conventions (2010), IERS Technical Note 36, Verlagdes Bundesamts für Kartographie und Geodäsie , Frankfurt am Main, Germany. Save H., Bettadpur S., Tapley B.D., 2016. High-resolution CSR GRACE RL05 mascons, J. geophys. Res. , 121, 7547– 7569. https://doi.org/10.1002/2016JB013007 Google Scholar CrossRef Search ADS   Shepherd A. et al.  , 2012. A reconciled estimate of ice-sheet mass balance, Science , 338, 1183– 1189. https://doi.org/10.1126/science.1228102 Google Scholar CrossRef Search ADS PubMed  Stuhne G.R., Peltier W.R., 2015. Reconciling the ICE-6G_C reconstruction of glacial chronology with ice sheet dynamics: the cases of Greenland and Antarctica, J. geophys. Res. , 120, 1841– 1865. https://doi.org/10.1002/2015JF003580 Google Scholar CrossRef Search ADS   Swenson S., Chambers D., Wahr J., 2008. Estimating geocenter variations from a combination of GRACE and ocean model output, J. geophys. Res. , 113, B08410, doi:10.1029/2007JB005338. https://doi.org/10.1029/2007JB005338 Google Scholar CrossRef Search ADS   Tapley B.D., Bettadpur S., Watkins M., Reigber C., 2004. The gravity recovery and climate experiment: mission overview and early results, Geophys. Res. Lett. , 31, L09607, doi:10.1029/2004GL019920. https://doi.org/10.1029/2004GL019920 Google Scholar CrossRef Search ADS   Van Dam T., Wahr J., Lavallée D., 2007. A comparison of annual vertical crustal displacements from GPS and Gravity Recovery and Climate Experiment (GRACE) over Europe, J. geophys. Res. , 112, B03404, doi:10.1029/2006JB004335. https://doi.org/10.1029/2006JB004335 Google Scholar CrossRef Search ADS   Van Dam T., Francis O., Wahr J., Khan S.A., Bevis M., Van Den Broeke M.R., 2017. Using GPS and absolute gravity observations to separate the effects of present-day and Pleistocene ice-mass changes in South East Greenland, Earth planet. Sci. Lett. , 459, 127– 135. https://doi.org/10.1016/j.epsl.2016.11.014 Google Scholar CrossRef Search ADS   Vasskog K., Langebroek P.M., Andrews J.T., Nilsen J.E.Ø., Nesje A., 2015. The Greenland Ice Sheet during the last glacial cycle: current ice loss and contribution to sea-level rise from a palaeoclimatic perspective, Earth-Sci. Rev. , 150, 45– 67. https://doi.org/10.1016/j.earscirev.2015.07.006 Google Scholar CrossRef Search ADS   Velicogna I., Wahr J., 2006. Acceleration of Greenland ice mass loss in spring 2004, Nature , 443, 329– 331. https://doi.org/10.1038/nature05168 Google Scholar CrossRef Search ADS PubMed  Velicogna I., Wahr J., 2013. Time-variable gravity observations of ice sheet mass balance: precision and limitations of the GRACE satellite data, Geophys. Res. Lett. , 40, 3055– 3063. https://doi.org/10.1002/grl.50527 Google Scholar CrossRef Search ADS   Wahr J., Molenaar M., Bryan F., 1998. Time variability of the Earth's gravity field: hydrological and oceanic effects and their possible detection using GRACE, J. geophys. Res. , 103, 30 205–30 229. https://doi.org/10.1029/98JB02844 Google Scholar CrossRef Search ADS   Wahr J., Khan S.A., Van Dam T., Liu L., Van Angelen J.H., Van Den Broeke M.R., Meertens C.M., 2013. The use of GPS horizontals for loading studies, with applications to northern California and southeast Greenland, J. geophys. Res. , 118, 1795– 1806. https://doi.org/10.1002/jgrb.50104 Google Scholar CrossRef Search ADS   Wang S.Y., Chen J.L., Li J., Hu X.G., Ni S., 2016. Geophysical interpretation of GPS loading deformation over western Europe using GRACE measurements, Ann. Geophys. , 59( 5), S0538, doi:10.4401/ag-7058. Wouters B., Bamber J.L., Van Den Broeke M.R., Lenaerts J.T.M., Sasgen I., 2013. Limits in detecting acceleration of ice sheet mass loss due to climate variability, Nat. Geosci. , 6, 613– 616. https://doi.org/10.1038/ngeo1874 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations