A Joint Inversion Estimate of Antarctic Ice Sheet Mass Balance Using Multi-Geodetic Data Sets

Abstract

Many recent mass balance estimates using the Gravity Recovery and Climate Experiment (GRACE) and satellite altimetry (including two kinds of sensors of radar and laser) show that the ice mass of the Antarctic ice sheet (AIS) is in overall decline. However, there are still large differences among previously published estimates of the total mass change, even in the same observed periods. The considerable error sources mainly arise from the forward models (e.g., glacial isostatic adjustment [GIA] and firn compaction) that may be uncertain but indispensable to simulate some processes not directly measured or obtained by these observations. To minimize the use of these forward models, we estimate the mass change of ice sheet and present-day GIA using multi-geodetic observations, including GRACE and Ice, Cloud and land Elevation Satellite (ICESat), as well as Global Positioning System (GPS), by an improved method of joint inversion estimate (JIE), which enables us to solve simultaneously for the Antarctic GIA and ice mass trends. The GIA uplift rates generated from our JIE method show a good agreement with the elastic-corrected GPS uplift rates, and the total GIA-induced mass change estimate for the AIS is 54 ± 27 Gt/yr, which is in line with many recent GPS calibrated GIA estimates. Our GIA result displays the presence of significant uplift rates in the Amundsen Sea Embayment of West Antarctica, where strong uplift has been observed by GPS. Over the period February 2003 to October 2009, the entire AIS changed in mass by −84 ± 31 Gt/yr (West Antarctica: −69 ± 24, East Antarctica: 12 ± 16 and the Antarctic Peninsula: −27 ± 8), greater than the GRACE-only estimates obtained from three Mascon solutions (CSR: −50 ± 30, JPL: −71 ± 30, and GSFC: −51 ± 33 Gt/yr) for the same period. This may imply that single GRACE data tend to underestimate ice mass loss due to the signal leakage and attenuation errors of ice discharge are often worse than that of surface mass balance over the AIS.

1. Introduction

As the largest single mass of ice on Earth, the Antarctic Ice Sheet (AIS) contains about 30 million cubic kilometers of ice [1] that makes that even a small relative change in its mass balance can have a noticeable effect on global sea level. To quantify the mass balance (MB) of an ice sheet, traditionally, we firstly have to determine mass input at its surface and mass output at its edges, surface mass balance (SMB) and ice discharge (ID), respectively [2], then

$$\mathrm{MB}=\mathrm{SMB}+\mathrm{ID}$$

(1)

where the value of SMB is usually positive, and ID is the opposite. This way to calculate mass balance is called the input-output method (IOM) or the mass budget method [3]. However, both surface mass balance and ice discharge in the AIS are large quantities with large uncertainties, making it difficult to accurately capture the much smaller difference between the two [2]. In the past two decades, some space geodetic techniques are available to measure the changes in mass and volume of ice sheet, such as the Gravity Recovery and Climate Experiment (GRACE) dual-satellite mission [4,5,6] and satellite altimetry includes two kinds of sensors of radar [7,8] and laser [9,10].

Most of estimates based on the three methods show the presence of significant loss of ice mass in the AIS over the past 40 years (see

Table 1. The trends of mass change over the Antarctic ice sheet (AIS) published in the last two decades.

Method	Time Period	Trend (Gt/yr)	Reference
IOM	1995–2000	−26 ± 27	Rignot and Thomas [30]
IOM	1996	−112 ± 91	Rignot et al. [3]
	2000	−138 ± 92
	2006	−196 ± 92
IOM	1992–2009	−83 ± 91	Rignot et al. [31,32]
IOM	2003–2009	−105 ± 59	Gardner et al. [33]
IOM	1979–1990	−40 ± 9	Rignot et al. [34]
	1989–2000	−50 ± 14
	1999–2009	−166 ± 18
	2009–2017	−252 ± 26
Radar	1992–1996	−60 ± 76	Wingham et al. [7]
Radar	1992–2001	−31 ± 12	Zwally et al. [24]
Radar	1992–2003	27 ± 29	Wingham et al. [35]
Radar	1995–2003	25 ± 14	Helsen et al. [28]
Radar	1992–2001	−31 ± 12	Zwally and Giovinetto [36]
Radar	2010–2013	−159 ± 48	McMillan et al. [8]
Radar	2011–2014	−128 ± 83	Helm et al. [37]
Radar	1992–2001	112 ± 61	Zwally et al. [25]
Laser	2003–2007	−98	Gunter et al. [9]
Laser	2003–2008	−78 ± 5	Shi et al. [38]
Laser	2003–2009	−6 ± 10	Gardner et al. [33]
Laser	2003–2009	−60 ± 44	Helm et al. [37]
Laser	2003–2009	−126 ± 20	Groh et al. [39]
Laser	2003–2008	82 ± 25	Zwally et al. [25]
Laser	2003–2008	−44 ± 21	Li et al. [10]
Radar and Laser	1992–2017	−85 ± 16	Schröder et al. [26]
	2010–2017	−137 ± 25
GRACE	2002–2005	−139 ± 73	Velicogna and Wahr [4]
GRACE	2002–2005	−40 ± 36	Ramillien et al. [40]
GRACE	2002–2006	−150 ± 73	Moore and King [41]
GRACE	2003–2007	−74	Gunter et al. [9]
GRACE	2002–2009	−143 ± 73	Velicogna [42]
	2002–2006	−104
	2006–2009	−246
GRACE	2002–2009	−190 ± 77	Chen et al. [5]
GRACE	2002–2008	−109 ± 48	Horwath and Dietrich [43]
GRACE	2002–2007	−78 ± 37	E et al. [44]
GRACE	2003–2008	−198 ± 22	Cazenave et al. [45]
GRACE	2002–2008	−91	Zhu et al. [46]
GRACE	2002–2010	−82 ± 29	Jia et al. [47]
GRACE	2002–2010	−69 ± 18	King et al. [6]
GRACE	2003–2010	−165 ± 72	Jacob et al. [48]
GRACE	2002–2010	−80	Luo et al. [49]
GRACE	2006–2011	−211 ± 75	Tang et al. [50]
GRACE	2003–2011	−83 ± 36	Barletta et al. [51]
GRACE	2003–2012	−57 ± 34	Ivins et al. [52]
GRACE	2003–2010	−81 ± 26	Luthcke et al. [29]
GRACE	2003–2012	−83 ± 49	Velicogna and Wahr [14]
		−147 ± 80
GRACE	2004–2012	−196 ± 21	Ju et al. [53]
GRACE	2003–2012	−107 ± 30	Groh et al. [39]
GRACE	2003–2012	−58 ± 16	Williams et al. [54]
GRACE	2003–2013	−171 ± 22	Schrama et al. [55]
	2003–2010	−156 ± 24
GRACE	2003–2014	−92 ± 10	Harig and Simons [56]
GRACE	2003–2013	−81 ± 27	Gao et al. [15]
GRACE	2003–2012	−107 ± 34	Mu et al. [57]
JIE	2002–2008	−87 ± 43	Wu et al. [19]
JIE	2002–2013	−114 ± 23	Sasgen et al. [58]
JIE	2003–2010	−47 ± 35	Mémin et al. [59]
JIE	2003–2009	−100 ± 44	Gunter et al. [20]
JIE	2003–2009	−67 ± 55	Gao et al. [21]
JIE	2003–2013	−84 ± 22	Martín-Español et al. [12]
JIE	2002–2015	−95 ± 50	Forsberg et al. [60]
JIE	1993–2000	−56 ± 28	Talpe et al. [61]
	2000–2005	20 ± 41
	2005–2013	−103 ± 20
JIE	2003–2009	−46 ± 43	Zhang et al. [22]
JIE	2003–2016	−141 ± 27	Sasgen et al. [23]
RE	1992–2011	−71 ± 53	Shepherd et al. [13]
	1992–2000	−48 ± 65
	1993–2003	−71 ± 61
	2000–2011	−87 ± 43
	2005–2010	−81 ± 37
RE	1992–1997	−49 ± 67	Shepherd et al. [11]
	1997–2002	−38 ± 64
	2007–2012	−160 ± 50
	2012–2017	−219 ± 43
	1992–2011	−76 ± 59
	1992–2017	−109 ± 56

). However, there are large differences among published estimates with a range of the total mass change from −252 to 112 Gt/yr (see Table 1). Even within the same observation periods, the differences are still considerable [11]. Some studies have suggested that the wide range is mainly caused by the uncertain forward models (e.g., snow accumulation model for IOM, glacial isostatic adjustment [GIA] model for GRACE, and firn compaction and surface density models for altimetry) [12], which are indispensable to simulate some processes that are unmeasured or unknown for these observations.

Of which, GIA is the ongoing deformation of the solid Earth due to past changes in ice-ocean surface loading, which causes an apparent and secular change in the Earth’s gravitational field. We have to remove the GIA effect from GRACE measurements before estimating ice sheet mass changes. Many studies confirmed that GRACE estimates of the AIS mass change are greatly affected by the use of different GIA forward models [13,14,15]. Mass change signals in Antarctica associated with GIA are poorly known in the forward models due to the lack of climatological and geophysical data to constrain the glacial history and viscoelastic Earth structure. To decrease the uncertainty of GIA, an alternative joint inversion estimate (JIE) approach was developed by combining multiple geodetic data [12,16,17,18,19,20,21,22,23]. Table 1 shows the disparity of mass change trends estimated by JIE or reconciled estimates (RE) using multi-geodetic data is evidently reduced comparing with that using only single method.

As to the altimetry, the unknown time-varying compaction effects of the firn can complicates the conversion of elevation-to-mass (ETM) changes (we call that “ETM model” in the following), e.g., the significant difference (more than 140 Gt/yr) can be found between Zwally et al. [24] and Zwally et al. [25] with same altimetry data but different ETM model. Wahr et al. [16] suggested multiplying by the density of ice (${\rho }_{ice}$ = 917 kg/m³) as an approximation of ETM model. Riva et al. [18] constructed a combined surface density distribution ${\rho }_{surf}$ as a more appropriate choice than ${\rho }_{ice}$. Gunter et al. [20], Gao et al. [21], Martín-Español et al. [12], Zhang et al. [22], Sasgen et al. [23], and Schröder et al. [26] used a regional climate model and associated firn densification model (FDM) [27] as a key constraint for ETM model. However, the density of conversion, whether ${\rho }_{ice}$ or ${\rho }_{surf}$, is only exactly correct in an ideal condition that obviously deviates from the actual in many areas over the AIS, especially in the marginal areas where the mass fluctuations are noticeable. This will introduce some uncertainty into the estimation of ice sheet mass balance. As to FDM, its reference period (1979 to present) is too short to accurately represent long-term (about 100-year cycle or longer) accumulation variability [28], which is an important assumption when using FDM data to estimate surface elevation changes [27].

Besides, some studies suggested using Global Positioning System (GPS) measurement of uplift rates as additional constraints to improve the JIE results [12,17,21,22,23]. The simulation results presented by Velicogna and Wahr [17] show the addition of GPS data is capable ofreducing the GIA errors in combined estimates of GRACE and Ice, Cloud and land Elevation Satellite (ICESat). Gao et al. [21] demonstrated the use of 35 GPS sites could help to reduce the GIA uplift rate of mm-level biases in Antarctica. Martín-Español et al. [12] developed a JIE method via the Bayesian hierarchical framework based on multiple data sets including GRACE, altimetry, GPS, FDM and so on, which would be able to obtain the components of mass change, GIA and firn compaction simultaneously. Sasgen et al. [23] proposed a new JIE method on the basis of the viscoelastic response function, of which the GPS data as key given parameters of function to solve the unknown component of GIA-related change. However, the GPS data used in these studies is mainly for optimizing the signal of GIA, it is not sensitive enough to improve the mass balance estimates directly. It is not doubtful that the improvement of GIA would be of great help for estimating the mass balance from GRACE, but, aside from GIA, the significant leakage (including “leakage-in” and “leakage-out”) and attenuation errors also make a mess of GRACE estimates due to the coarse spatial resolution [29].

In this paper, we present an improved JIE method, which evolved from the approach of Wahr et al. [16], to determine the mass change and GIA using multi-geodetic data sets including GRACE, ICESat and GPS. A new contribution of our method is to improve the ETM model based on the GPS uplift rates instead of the previously proposed FDM. The GPS data is used as a main controller to drive the estimations of ice sheet mass change and GIA based on regional clustering and blocked processing. In addition, we assume a linear relationship between the elevation change of SMB and firn compaction to simulate the ETM model of ICESat. Our JIE method allows to avoid the use of uncertain forward models as much as possible and to improve the precision and reliability of estimating ice mass change and GIA over the AIS.

In addition, as was done in Martín-Español et al. [62], we present an assessment of our GIA estimate and eight recent forward and inverse GIA solutions based on a comparison with vertical height displacement measurements from GPS data. We analyze similarities and differences of these GIA solutions on regional scale. Finally, we discuss our mass change estimates in 27 drainage systems [63] and the major regions of West Antarctica (WA), East Antarctica (EA) and the Antarctic Peninsula (AP) and compare them against the estimates from single ICESat and GRACE data in certain drainage systems and regions.

2. Method and Data

2.1. The JIE Method to Estimate GIA

The currently available geodetic observations for AIS are mass change rates from GRACE, $\Delta m_{GRACE}$, surface elevation change rates from satellite altimetry, $\Delta h_{SA}$, and bedrock uplift rates from GPS, $\Delta h_{GPS}$. According to the principle of measurement and respective sensitivity of these observations, we can decompose $\Delta m_{GRACE}$, $\Delta h_{SA}$, and $\Delta h_{GPS}$ into the following components [22,64]:

$$\Delta m_{GRACE}=\Delta m_{SMB}+\Delta m_{ID}+\Delta m_{GIA}$$

(2)

$$\Delta h_{SA}=\Delta h_{SMB}+\Delta h_{fc}+\Delta h_{ID}+\Delta h_{GIA}+\Delta h_{ela}$$

(3)

$$\Delta h_{GPS}=\Delta h_{GIA}+\Delta h_{ela}$$

(4)

where $\Delta m$ and $\Delta h$ are rates of mass and elevation change, respectively, with the individual components on the right of Equations (2)–(4) representing the rates related to SMB ($\Delta m_{SMB}$ and $\Delta h_{SMB}$), ice discharges ($\Delta m_{ID}$ and $\Delta h_{ID}$), GIA ($\Delta m_{GIA}$ and $\Delta h_{GIA}$), firn compaction ($\Delta h_{fc}$), and elastic response of the lithosphere to loading by ice sheet mass change due to SMB and ice discharges ($\Delta h_{ela}$).

Considering Equations (1)–(3), we can express the rates of mass balance $\Delta m_{MB}$ using GRACE and altimetry observations:

$$\Delta m_{MB}=\Delta m_{GRACE}-\Delta m_{GIA}$$

(5)

$$\Delta m_{MB}=\mathrm{ETM}(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela})$$

(6)

where ETM is the conversion model of elevation-to-mass changes for altimetry data (see Section 2.3 below for more discussions).

From Equations (5) and (6), GIA-related mass change rates can be estimated by combining GRACE and altimetry observations:

$$\Delta m_{GIA}=\Delta m_{GRACE}-\mathrm{ETM}(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela})$$

(7)

where $\Delta h_{GIA}$ and $\Delta h_{ela}$ in the initial estimate or zeroth iteration are equal to 0 [16], and in the following iteration are derived as:

$$\Delta h_{GIA}^{nm}=\frac{h_n^{ve}}{k_n^{ve}}\frac{1}{2n+1}\frac{3}{{\rho }_{ave}}\Delta m_{GIA}^{nm}\left(n\ge 2\right)$$

(8)

$$\Delta h_{ela}^{nm}=\frac{h_n^{ela}}{2n+1}\frac{3}{{\rho }_{ave}}\left(\Delta m_{GRACE}^{nm}-\Delta m_{GIA}^{nm}\right)\left(n\ge 2\right)$$

(9)

where the superscripts $n$ and $m$ represent degree and order in spherical harmonic domain, $h_n^{ve}/k_n^{ve}$ is the ratio of viscoelastic Love loading numbers (here we use the improved ratio values as provided by Table S3 of Purcell et al. [65]), ${\rho }_{ave}$ denotes the average density of the Earth, and $h_n^{ela}$ is the elastic load Love number of degree.

Here we can obtain the optimal result of GIA by iterating Equations (7)–(9) until the improvement of GIA is negligible. To evaluate the accuracy of GIA, $\Delta h_{GIA}$ and $\Delta h_{ela}$ can be interpolated to the locations of all GPS sites as $\Delta h_{GIA}^{interp}$ and $\Delta h_{ela}^{interp}$, and the GIA uplift rates obtained from GPS measurements are calculated by $\Delta h_{GIA}^{GPS}=\Delta h_{GPS}-\Delta h_{ela}^{interp}$, then the weight right mean square (WRMS) [66] of GIA uplift rates can be derived from the following:

$$\mathrm{WRMS}=\sqrt{\frac{\sum q{\left(\Delta h_{GIA}^{GPS}-\Delta h_{GIA}^{interp}\right)}^2}{\sum q}}$$

(10)

where $q=\frac{1}{{(\delta \Delta h_{GIA}^{GPS})}^2+{(\delta \Delta h_{GIA}^{interp})}^2}$, of which $\delta \Delta h_{GIA}^{GPS}$ and $\delta \Delta h_{GIA}^{interp}$ are the uncertainties of $\Delta h_{GIA}^{GPS}$ and $\Delta h_{GIA}^{interp}$ respectively.

2.2. Input Datasets

2.2.1. GRACE

Here, we use the Release 05 (RL05) monthly GRACE gravity field solutions of the Center for Space Research (CSR) at the University of Texas, Austin. Each solution is represented by fully normalized Stokes potential coefficients complete to spherical-harmonic degree and order 60. The coefficients C20 are replaced with values from satellite laser ranging [67] and the values of missing degree-one coefficients are inserted from Swenson et al. [68]. The rates of mass change in terms of mm/yr of equivalent water height are determined from GRACE gravity field models using the classical transformation method as described by Wahr et al. [69] and the least squares adjustment of eight-parameter function consisting of a constant, a linear trend, annual and semi-annual periodic signals, and 161-day period of S2 tidal alias [70]. The error of typical north–south oriented stripes is significant in the GRACE monthly solutions, but it is not apparent in the fitting trend especially for the polar regions. The rates of trend are what we need in the JIE method, so the de-striping filter [71] is not used here. The higher order noise is suppressed using the Fan filter of 300 km [72], which is chosen as the optimal half-width by considering the spatial resolution of GRACE [73] and the sensitive wavelength of GIA in Antarctica [17]. The leakage-in errors are reduced using the method proposed by Velicogna and Wahr [14]. In Figure 1a,b we show the linear trend rates and corresponding uncertainties of mass change measured by GRACE over only the period between February of 2003 and October of 2009, which is the same time for the entire ICESat mission period. The uncertainties are estimated using the approach described by Gao [74], which are propagated from the error of measurement, leakage, atmosphere, and tidal alias [14].

2.2.2. ICESat

We use satellite altimetry measurements from the latest version (V34) of ICESat level−2 products [75], providing surface elevations for ice sheets (GLAH12) over the period February 2003 to October 2009. The surface elevation change rates are estimated using a near repeat-track analysis based on triangular irregular networks [76]. The inter-campaign biases are corrected using the value of 1.04 ± 0.48 cm/yr from Hofton et al. [77]. The uncertainties of elevation change rates are derived from the propagation of the residuals after fitting, and more details for ICESat data processing can be found in Shi [78]. In Figure 1c,d, we show the linear trend rates and corresponding uncertainties of elevation changes estimated by ICESat, which are gridded into a ${0.125}^{°}\times {0.125}^{°}$ grid. To match with the spatial resolution of GRACE, the mass changes estimated by ICESat should be expanded into spherical harmonic coefficients of degree and order 60 with the same Fan filter of 300 km radii as are used on GRACE (this process is referred to simply as T60F300 in the following) when applied in the JIE method.

2.2.3. GPS

The GPS uplift rates from a total of 118 Antarctic sites covered the time span from 1995 to 2013 are listed in Sasgen et al. [79]. We use 86 of them with the temporal coverage limited to 2003–2013 and remove of which errors more than 5 mm/yr that exceed the uncertainties in altimetry/GRACE recovery without GPS data [17]. Regional clusters are sorted using the K-means algorithm [23] with a pre-defined threshold value of 400 km, which is analogous to the length scale recovered with GRACE after smoothing, to reduce stochastic and geophysical noise of neighboring stations. In addition, we use a low-precipitation zone (see grey area in Figure 2) analogous to that of Gunter et al. [20] to constrain the inland of EA where GPS sites are scarce. The rates of GIA vertical deformation and the surface height change in low-precipitation zone are expected to be very small, so the GIA uplift rates can be assumed to be zero [20] and the GPS uplift rates in that region can be determined by Equation (9) from GRACE data directly. Considering that the fluctuation of most GIA results in low-precipitation zone is less than ±2 mm/yr, the uncertainties of GIA uplift rates in the region can be assumed to be ±2 mm/yr conservatively, then the errors of GPS uplift rates can be computed by error propagation ($\delta \Delta h_{GPS}=\sqrt{\delta \Delta h_{ela}^2+\delta \Delta h_{GIA}^2}$). In Figure 2, we show the uplift rates and their uncertainties of 22 GPS sites after clustering in Antarctica, and four assumed sites that are arranged uniformly in low-precipitation region. The elastic response, which can be determined by Equation (9), should be deducted from GPS uplift rates when evaluating the GIA using Equation (10). It is noted that $\Delta m_{GRACE}^{nm}$ in Equation (9) here should be estimated by the GRACE data with the time span as same as GPS data (2003–2013), and $\Delta m_{GIA}^{nm}$ is calculated by the JIE method.

2.3. The Method to Improve ETM Model

As described in Section 1, there are two kinds of ETM model, one is the use of the ice or assigned surface density directly, and the other is to add the FDM as a constraint. Here we test various ETM models suggested by previous authors (see

Table 2. Comparison of the estimated glacial isostatic adjustment (GIA) uplift rates derived from different elevation-to-mass (ETM) models and other GIA solutions with 26 GPS sites over Antarctica. $\Delta h_{firn}$ and $\Delta m_{firn}$ represent the surface height and mass change of the firn, and ${\rho }_{\alpha }$ is an assigned density. The unit of weight right mean square (WRMS) is mm/yr, only $\delta \Delta h_{GIA}^{GPS}$ used in the calculation of WRMS. The number in parenthesis denotes the total quantity of GPS sites within each region.

Reference	ETM or GIA Model	WRMS
		WA (5)	EA (18)	AP (3)	All (26)
Wahr et al. [16]	$(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela}){\rho }_{ice}$	6.0	2.8	8.9	3.5
Riva et al. [18]	$(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela}){\rho }_{surf}$	6.3	3.0	8.9	3.7
Gunter et al. [20]	$(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela}-\Delta h_{firn}){\rho }_{\alpha }+\Delta m_{firn}$	6.8	3.2	7.1	3.8
Ligtenberg (personal communication, 2017)	$(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela}-\Delta h_{firn}){\rho }_{ice}+\Delta m_{firn}$	6.1	3.0	11.4	3.8
This study	Equation (12) (using $\mu =-0.619$)	6.0	2.8	9.0	3.5
This study	Equation (12) (using $\mu$ values as shown in Table 3)	3.0	0.4	0.8	1.0
Wang et al. [80]	RF3L20	3.9	2.9	4.9	3.1
Whitehouse et al. [66]	W12a	3.6	1.3	1.7	1.7
A et al. [81]	A13	3.7	1.5	3.2	1.9
Peltier et al. [82]	ICE6G	3.4	0.7	1.6	1.3
Gunter et al. [20]	G14	3.9	1.5	2.5	1.9
Gao et al. [21]	G16	3.5	1.2	1.9	1.6
Martín-Español et al. [12]	RATES	3.5	0.6	1.9	1.3
Sasgen et al. [23]	REGINA	3.6	1.1	2.2	1.6

) and compare their respective GIA results (Figure 3a–d) estimated by the JIE method as mentioned in Section 2.1 with the 26 elastic-corrected GPS uplift rates. In Table 2, we show all the GIA results have large WRMS values in WA and the AP where the key areas for estimating mass balance of the AIS are. Test result displays the application of FDM is not of much help to improve the inversion of GIA. This proves that the spatial feature of surface height change in marginal area of the AIS revealed by FDM (see Figure 6a of Gunter et al. [20]) is not supported by multi-geodetic observations.

Based on the above test, all of the existing ETM models have significant deviations that can contaminate our GIA estimates. To resolve this issue, we first use the suggestion of Wahr et al. [16] as the initial approximation for ETM model and the completed ETM function can be expressed as

$$\Delta m_{MB}=(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela}){\rho }_{ice}+\Delta m_{SMB}-(\Delta h_{SMB}+\Delta h_{fc}){\rho }_{ice}$$

(11)

where $\Delta m_{SMB}$ can be estimated by SMB model which derived from the RACMO2 regional atmospheric climate model [83], and $\Delta h_{SMB}=\frac{\Delta m_{SMB}}{{\rho }_{snow}}$, ${\rho }_{snow}$ is the density of snow in Anarctica (${\rho }_{snow}=350{\mathrm{kg}\mathrm{m}}^{-3}$), only $\Delta h_{fc}$ as an unknown quantity if the FDM is not used. We noted that both of SMB and firn compaction effects are dominated by the accumulation rates in Antarctica [28], so we can assume a simple linear correlation is existed between $\Delta h_{SMB}$ and $\Delta h_{fc}$ approximately, i.e., $\Delta h_{fc}\approx \mu \Delta h_{SMB}$, then Equation (11) can be simplified as

$$\Delta m_{MB}=(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela}){\rho }_{ice}+\Delta \mathrm{C}$$

(12)

where $\Delta \mathrm{C}=\Delta m_{SMB}(1-\left(1+\mu \right)\frac{{\rho }_{ice}}{{\rho }_{snow}})$. When $\mu =\frac{{\rho }_{snow}}{{\rho }_{ice}}-1$ or $\Delta m_{SMB}=0$, Equation (12) is equal to the suggestion of Wahr et al. [16], so $\Delta \mathrm{C}$ can be regarded as a correction term for the initial approximation. It is easy to obtain $\mu =-0.619$ by the nonlinear least square principle ($min\left|\right|WRMS\left|\right|$) using iterative operation, and then we can use this $\mu$ value in Equation (12) as the further approximation for ETM model to determine associated GIA estimates (see Figure 3e) using the JIE method. Table 2 and Figure 3 show the uplift rates of GIA estimated by the further approximation are very close to that by the initial approximation of ETM model. That means a uniform $\mu$ value in entire Antarctica is unsuitable due to the significant difference between WA and EA in the characteristic of snow accumulation and firn compaction. Coupled with the various characteristics for ices discharges and the errors of GRACE and ICESat observations in different regions, it is difficult to find the reliable $\Delta \mathrm{C}$ with the same linear relation in entire Antarctica.

The quasi-uniform spatial distribution of GPS uplift rates over the entire Antarctica (see Figure 2) has made it possible to construct the correction component of $\Delta \mathrm{C}$ with the various values in different subregions. Therefore, we divide the entire Antarctica into 20 subregions according to the comprehensive consideration of the mean (2003–2009) surface mass balance $\Delta m_{SMB}$ (see Figure 4a) that is closely associated with firn compaction in the AIS [27], surface temperature, the Antarctic drainage systems (DS) [63] that are influenced by ice dynamic, and the distribution of GPS sites.

According to the Figure 3e and Figure 4a, it is clear that the large discrepancies between the uplift rates of GIA and GPS are mainly distributed in the regions where $\Delta m_{SMB}$ above 100 kg m⁻² yr⁻¹, maybe due to the large errors of firn compaction, so we take these areas as main reconstructed areas in the first step. Of which, the AP, where high accumulation values lead to a relatively high ice sheet integrated surface mass balance ($\Delta m_{SMB}$ > 500 kg m⁻² yr⁻¹)and large errors of firn compaction, is divided into AP1 (DS25 and DS26) and AP2 (DS24 and DS27) (see Figure 4b) in consideration of the difference between the two areas in surface temperature [84], drainage system, and the position of GPS sites. WA is divided into WA1 (coastal DS23 where large accumulation amounts occur), WA2 (DS21 and DS22 where large accumulation amounts occur), WA3 (coastal DS20 where large accumulation amounts occur), WA4 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS19), WA5 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS18), WA6 (the regions where $\Delta m_{SMB}$ ≤ 200 kg m⁻² yr⁻¹ in DS1), and WA7 (the regions where $\Delta m_{SMB}$ ≥ 200 kg m⁻² yr⁻¹ in DS1) (see Figure 4b) based on the difference of accumulation rates and drainage system. EA is divided into EA1 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in the western of DS2 and DS17), EA2 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in the eastern of DS17), EA3 (the regions where where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS16), EA4 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS14 and DS15), EA5 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS13), EA6 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS11 and DS12), EA7 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS8 and DS9), EA8 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS5, DS6, and DS7), and EA9 (the regions where $\Delta m_{SMB}$ > 100 kg m⁻² yr⁻¹ in DS3, DS4, and the northeast of DS2) (see Figure 4b) in view of the different accumulation rates, drainage systems and the position of GPS sites. For the next step, the region where $\Delta m_{SMB}$ ≤ 100 kg m⁻² yr⁻¹ is divided into the two subregions: WA8 (the regions where $\Delta m_{SMB}$ ≤ 100 kg m⁻² yr⁻¹ in DS18, DS19 and the western of DS17) and EA10 (the regions where $\Delta m_{SMB}$ ≤ 100 kg m⁻² yr⁻¹ in DS2–16 and the eastern of DS17). The distributions of all subregions are shown in Figure 4b.

Based on the above division, we can calculate the optimal values of $\mu$ for each subregion (see

Table 3. The optimal values of $\mu$ for each subregions. The recomputed $\mu$ values using low rock densities of 2000 kg/m³ in the Antarctic Peninsula (AP) and of 2200 kg/m³ in West Antarctica (WA), which are computed by the viscoelastic response function [79] in consideration of a thin lithosphere and low viscosity in these regions.

Subregions	$\mathit{\mu }$	The Recomputed $\mathit{\mu }$	Subregions	$\mathit{\mu }$
AP1	−0.530	−0.542	EA1	−0.656
AP2	−0.603	−0.611	EA2	−0.563
WA1	−0.550	−0.579	EA3	−1.000
WA2	−0.522	−0.575	EA4	−0.595
WA3	−0.579	−0.575	EA5	−0.613
WA4	−0.793	−0.814	EA6	−0.660
WA5	−0.384	−0.453	EA7	−0.672
WA6	−0.652	−0.708	EA8	−0.664
WA7	−0.684	−0.672	EA9	−0.595
WA8	−0.474	−0.502	EA10	−0.583

) using the nonlinear least square principle as used in the further approximation. Substituting these values into Equation (12) as the optimal approximation of ETM model, then we can obtain a new GIA estimate (see Figure 3f) by the JIE method. Comparing with GPS data, we show the new GIA result with just 1.0 mm/yr of WRMS error in entire Antarctica (Table 2), there is a significant improvement relative to that determined by other ETM models (3.5 to 3.8 mm/yr).

It is notable that the magnitude of $\mu$ values for each subregions are not only sensitive to the correlation between accumulation and compaction, but also to the errors or biases of GRACE, GPS, and in particular, the ICESat dataset. For instance, there are some abnormal $\mu$ values occurred in WA5, WA8, WA4, and EA3 (see Table 3). In WA5 and WA8, the subnormal values may be related with the considerable ICESat errors of surface height changes (shown in Figure 1d). In WA4, however, it is likely associated with the significant leakage-in errors from adjacent WA2 of strong mass change. In EA3, it is possibly from the obvious stripe errors that can be found in all inversed GIA estimates (see Figure 3). This means not only the compaction effects but also the observed deviations of altimetry or GRACE can be corrected by adjusting $\mu$ values. Thus, our GIA result is highly dependent on the accuracy of GPS observations and the reliability of the JIE method itself.

3. Assessment of the GIA Estimates

In order to gain more insight into the performance of our GIA uplift rates, we also calculate the WRMS of the residuals between the 26 GPS uplift rates and the uplift rates derived from the recent forward and inverse solutions available for Antarctica: RF3L20 [80], W12a [66], A13 [81], ICE6G [82], G14 [20], G16 [21], RATES [12], and REGINA [23]. The WRMS comparisons are shown in Table 2. To examine the differences of spatial patterns, all above existing GIA uplift rates and their discrepancies with GPS observations computed at each station are plotted in Figure 5. Furthermore, we use another 67 elastic-corrected GPS uplift rates to re-evaluate all above GIA solutions through the weighted mean (WM) and WRMS as proposed by Martín-Español et al. [62] to avoid the use of same GPS data in the processes of both estimation and assessment that may make the evaluation biased. Also we use the same regional clustering as described in Section 2.2.3 to reduce the bias of spatial correlation in neighboring GPS stations, only 20 sites left after clustering. The WM and WRMS errors for all GIA solutions over WA, EA, AP, and entire Antarctica estimated by new GPS data as shown in

Table 4. The weighted mean (WM) and WRMS residuals, in mm/yr, between the estimated and observed GIA uplift rates at each GPS location. Note that the elastic-corrected GPS uplift rates proposed by Martín-Español et al. [62] and the way of WRMS calculation here is slightly different from Equation (10) (see Equation (2) in Martín-Español et al. [62]).

GIA	WA (7)		EA (10)		AP (3)		All (20)
	WM	WRMS	WM	WRMS	WM	WRMS	WM	WRMS
RF3L20	1.4	4.4	2.5	3.7	6.4	7.3	2.5	4.1
W12a	−0.2	4.4	1.4	2.0	−0.2	0.7	1.0	2.5
A13	0.4	4.2	1.0	2.3	4.4	6.0	1.0	3.0
ICE6G	2.0	4.3	0.5	0.9	1.9	3.0	0.8	2.1
G14	−0.3	3.8	0.7	1.5	−1.4	2.7	0.4	2.2
G16	−0.8	4.2	0.7	1.3	2.5	3.3	0.5	2.2
RATES	0.1	3.8	0.1	0.6	1.8	2.2	0.2	1.8
REGINA	−2.1	5.3	−0.2	1.1	1.7	2.9	−0.4	2.5
This study	−0.8	4.1	0.2	0.8	1.2	2.1	0.1	1.9

.

In Table 2, we show our GIA uplift rates derived from JIE method using ETM model of the optimal approximation have the smallest WRMS value, whether in WA, EA, the AP or entire AIS, among all GIA estimates. Over entire Antarctica, the WRMS value from our estimate of 1.0 mm/yr (WA: 3.0, EA: 0.4, the AP: 0.8) is smaller than other solutions with the range of 1.3–3.1 mm/yr (WA: 3.4–3.9, EA: 0.6–2.9, the AP: 1.6–4.9). In Table 4 we show our GIA estimate still has good overall agreement with an independent set of GPS-derived uplift rates with 1.9 mm/yr (WA: 4.1, EA: 0.8, the AP: 2.1) of WRMS error over entire Antarctica. It is closely followed by the best result of 1.8 mm/yr (WA: 3.8, EA: 0.6, the AP: 2.2) from RATES that has been corrected by this GPS dataset [12].

From Figure 3e and Figure 5, we can find that the largest misfits between GIA solutions and GPS observations are located in the Amundsen Sea Embayment of WA (DS21 and DS22) where the glaciers of Pine Island and Thwaites have thinned at an accelerating rates [34]. The consistent positive biases occurred in this region for most of GIA solutions, meaning that these GIA solutions underestimate the uplift rates of bedrock. Several studies suggest a very low viscosity in the upper mantle and thin elastic lithosphere [85] over there, which would lead to large uplift rates due to a rapid viscoelastic response to more recent ice loss in this area [23,62,86].

Most of the GIA solutions, especially for the inverse solutions, underestimate the uplift rates over Ellsworth Mountains (WA5). This may be because the leakage-in errors from neighboring Kamb Ice Stream that is thickening revealed by ICESat observation (see Figure 1c), or because the rapid viscoelastic response to more recent ice loss in this area with the thin elastic lithosphere [23].

In the AP and the area surrounding the ice shelves of Filchner-Ronne and Ross, most of GIA forward models overestimate the uplift rates, but all inversed GIA results show good agreement with the GPS observations. In the coastal DS23 of WA, a systematic overestimation of the uplift rates across most GIA solutions likely due to the leakage errors from the neighboring Amundsen Sea Embayment where the strong ice loss and GIA uplift are taking place. In Wilkes Land of EA, a significant uplift rate is found in our GIA result, which also can be found in REGINA. This could be because the negative load change in the last 16 kyr in that region [87], or because an insufficient correction for the JIE method under the lack of nearby GPS sites to constrain the ETM model.

4. Results and Discussion

Sasgen et al. [23] found that a strong uplift rate of GIA can be associated with a small gravity perturbation in the AP and WA due to the lithosphere is thin in these regions. This means that the $\mu$ values of these regions will be incorrectly estimated due to averaged ratio values of viscoelastic Love loading numbers used in our JIE method, resulting in significant errors into subsequent estimations of mass changes of ice sheet and GIA. Thus, we recalculate $\mu$ values using the low rock densities ${\rho }_{rock}$ of 2000 kg/m³ in the AP and 2200 kg/m³ in WA by following spatial least square principle:

$$min||\Delta m_{GRACE}-\Delta h_{GIA}·{\rho }_{rock}-\mathrm{T}60\mathrm{F}300[(\Delta h_{SA}-\Delta h_{GIA}-\Delta h_{ela}){\rho }_{ice}+\Delta \mathrm{C}]||$$

(13)

The low rock densities are obtained from the viscoelastic response function [79] under an assumed equilibrium state between load forcing and deformational response based on a thin elastic lithosphere [23,88] without ductile layer and low asthenosphere viscosity of 1 × 1018 Pa s [79]. The recomputed $\mu$ values are shown in Table 3. Then we can estimate the ice sheet mass balance using Equation (12) for time period February 2003 to October 2009 and GIA-related mass change using Equation (7) over the 27 Antarctic drainage systems (basin locations as shown in Figure 4a, without ice shelves), WA, EA, AP, and the entire Antarctica (see

Table 5. The mean trends of 27 Antarctic drainage systems, the AP, WA, East Antarctica (EA), and the entire AIS for GIA-induced mass change and ice sheet mass balance derived from different methods over the period February 2003 to October 2009, in Gt/yr. The GIA effects for all ICESat-only and GRACE-only estimates are corrected using our GIA result. To reduce leakage-out errors across land/ocean boundaries, the regions of ice shelf are included in the Antarctic drainage systems when estimating from the Mascon solutions of Center for Space Research (CSR) and Goddard Space Flight Center (GSFC) (no need for the solution of Jet Propulsion Laboratory [JPL] due to that has been used a CRI filter).

Regions	The JIE Estimate in This Study			ICESat-Only Estimates Using ETM Models as Suggested by		GRACE-Only Estimates Using Mascon Products from
	$\mathbf{\Delta }{\mathbf{m}}_{\mathbf{MB}}$	$\mathbf{\Delta }\mathbf{C}$	GIA	Riva et al. [18]	Gunter et al. [20]	CSR	JPL	GSFC
DS1	33.7 ± 6.8	21.8 ± 2.8	3.3 ± 6.0	4.5 ± 2.4	−0.1 ± 0.9	19.8 ± 8.3	15.6 ± 8.4	20.1 ± 9.4
DS2	0.2 ± 1.2	−1.2 ± 0.4	4.1 ± 2.2	0.2 ± 0.3	−3.3 ± 0.6	2.3 ± 2.9	0.1 ± 2.6	2.6 ± 3.2
DS3	11.7 ± 7.1	−7.6 ± 0.7	0.9 ± 5.5	7.6 ± 4.8	−3.0 ± 1.4	12.3 ± 5.9	13.4 ± 5.7	13.2 ± 5.8
DS4	6.1 ± 1.5	−3.2 ± 0.4	1.1 ± 2.2	4.6 ± 0.5	0.0 ± 0.8	7.4 ± 2.7	9.1 ± 2.7	8.7 ± 2.7
DS5	3.7 ± 1.3	3.0 ± 0.5	0.4 ± 1.4	0.3 ± 0.5	−0.6 ± 0.3	3.9 ± 1.8	3.6 ± 1.9	3.3 ± 1.8
DS6	0.6 ± 3.8	4.2 ± 1.3	0.9 ± 3.1	−1.6 ± 1.4	−2.0 ± 1.0	3.7 ± 3.4	1.3 ± 3.6	1.6 ± 3.4
DS7	8.2 ± 1.6	6.9 ± 1.2	1.0 ± 2.8	0.9 ± 0.5	0.3 ± 0.5	10.6 ± 2.9	11.8 ± 3.1	9.1 ± 2.9
DS8	6.9 ± 1.2	4.7 ± 0.8	−0.3 ± 1.1	1.0 ± 0.3	0.7 ± 0.4	5.0 ± 1.1	6.4 ± 1.2	4.8 ± 1.1
DS9	7.5 ± 2.0	2.3 ± 0.6	−0.4 ± 1.2	2.2 ± 0.7	0.8 ± 0.4	2.8 ± 1.5	3.0 ± 1.5	3.1 ± 1.5
DS10	2.3 ± 3.4	−4.6 ± 0.6	2.6 ± 3.2	2.4 ± 1.2	−0.1 ± 0.5	3.3 ± 3.6	3.4 ± 3.5	1.8 ± 3.4
DS11	1.7 ± 1.5	−1.3 ± 0.2	0.8 ± 1.4	4.2 ± 0.6	1.1 ± 0.1	−0.4 ± 1.6	−0.5 ± 1.6	−1.2 ± 1.6
DS12	12.1 ± 7.0	13.5 ± 2.5	1.0 ± 3.6	−1.2 ± 2.9	0.9 ± 1.4	7.1 ± 3.9	8.2 ± 4.0	6.1 ± 4.0
DS13	−17.3 ± 2.2	−5.1 ± 0.5	9.4 ± 4.4	−10.3 ± 1.1	−4.0 ± 1.3	−12.0 ± 4.9	−13.6 ± 5.1	−11.2 ± 5.3
DS14	−12.8 ± 4.1	−8.5 ± 1.2	2.7 ± 3.8	−2.9 ± 2.8	−2.7 ± 1.7	−7.8 ± 3.9	−11.4 ± 4.0	−6.0 ± 4.0
DS15	−6.3 ± 2.0	−1.5 ± 0.5	0.4 ± 1.5	−3.0 ± 0.7	−5.3 ± 1.1	−3.4 ± 1.6	−5.2 ± 1.6	−3.8 ± 1.6
DS16	−3.7 ± 1.2	1.6 ± 0.1	1.5 ± 1.6	−4.1 ± 0.4	−5.3 ± 0.2	−1.9 ± 1.7	−2.9 ± 1.7	−1.9 ± 1.7
DS17	−9.1 ± 4.1	−3.5 ± 0.7	8.2 ± 3.6	−1.6 ± 1.7	−2.3 ± 0.3	−6.3 ± 4.5	−7.6 ± 4.1	−5.6 ± 4.3
DS18	6.1 ± 4.7	−12.5 ± 2.2	4.7 ± 4.3	15.6 ± 1.5	6.7 ± 1.1	12.0 ± 4.6	11.8 ± 4.6	12.0 ± 4.9
DS19	18.4 ± 6.6	18.5 ± 3.5	0.4 ± 6.7	−0.3 ± 2.6	0.2 ± 1.5	7.3 ± 7.2	4.6 ± 7.2	8.6 ± 7.6
DS20	−25.3 ± 12.3	−14.9 ± 3.2	1.0 ± 3.4	−9.7 ± 10.3	−9.6 ± 3.4	−21.4 ± 3.5	−21.4 ± 3.5	−20.9 ± 3.7
DS21	−59.9 ± 4.3	−16.1 ± 2.3	5.5 ± 6.6	−39.6 ± 2.4	−43.6 ± 3.2	−47.8 ± 7.7	−43.0 ± 7.8	−49.5 ± 8.9
DS22	−38.2 ± 6.4	−13.1 ± 1.9	3.8 ± 5.8	−23.8 ± 4.5	−31.0 ± 5.2	−31.5 ± 5.8	−36.1 ± 5.9	−31.1 ± 6.4
DS23	−3.4 ± 1.3	−7.8 ± 0.8	0.4 ± 0.7	2.1 ± 0.3	−4.2 ± 1.2	−1.7 ± 2.2	−2.6 ± 2.2	−5.0 ± 2.7
DS24	3.0 ± 5.1	−2.3 ± 0.4	0.0 ± 2.3	3.3 ± 2.5	−1.7 ± 1.7	1.2 ± 3.3	−0.2 ± 3.4	1.1 ± 3.4
DS25	−14.8 ± 2.1	−11.4 ± 2.0	0.2 ± 0.5	−3.1 ± 0.3	−6.6 ± 0.5	−6.5 ± 0.5	−9.9 ± 0.5	−5.4 ± 0.6
DS26	−14.9 ± 1.8	−6.6 ± 1.1	0.3 ± 0.6	−7.9 ± 0.7	−6.9 ± 0.5	−10.6 ± 1.3	−10.1 ± 1.3	−7.9 ± 1.4
DS27	−0.4 ± 0.8	−0.3 ± 0.1	0.4 ± 0.7	0.1 ± 0.3	−0.8 ± 0.5	2.3 ± 0.9	1.8 ± 0.9	2.0 ± 1.0
WA	−68.5 ± 24.2	−24.1 ± 10.0	19.1 ± 20.9	−51.3 ± 14.3	−81.7 ± 10.1	−63.3 ± 23.5	−71.1 ± 23.7	−65.8 ± 26.2
EA	11.7 ± 16.0	−0.5 ± 4.4	34.4 ± 15.1	−1.2 ± 7.7	−24.8 ± 4.2	26.5 ± 16.8	19.0 ± 16.8	24.8 ± 16.9
AP	−27.1 ± 7.8	−20.7 ± 2.5	0.9 ± 2.7	−7.6 ± 2.7	−16.0 ± 3.6	−13.6 ± 4.6	−18.4 ± 4.7	−10.3 ± 4.9
AIS	−83.9 ± 31.0	−45.4 ± 11.3	54.4 ± 26.5	−60.1 ± 16.7	−122.5 ± 11.8	−50.4 ± 30.2	−70.5 ± 30.4	−51.3 ± 32.6

).

As a comparison, we also use ICESat data with two ETM models as suggested by Riva et al. [18], Gunter et al. [20] and three GRACE Mascon solutions to estimate the ice mass trends for same period February 2003 to October 2009. Three Mascon products respectively are: (i) the RL05 solution from CSR [89]; (ii) the RL05 solution employs a Coastal Resolution Improvement (CRI) filter from Jet Propulsion Laboratory (JPL) [90]; (iii) the Version 02.4 solution from Goddard Space Flight Center (GSFC) [29]. The GIA effects for all ICESat-only and GRACE-only estimates are corrected using our GIA result. The estimates are listed in Table 5.

The uncertainties of our JIE result are computed by formal error propagation from the individual input sources. For the right side of Equation (12), the uncertainties of $\Delta h_{GIA}$ and $\Delta h_{ela}$ can be negligible due to the effect is far smaller than $\delta \Delta h_{ICESat}$, the uncertainties of $\Delta \mathrm{C}$ is difficult to assess due to insufficient error information for SMB model and $\mu$ values, so an standard deviation of 20% of the value for each grid point is assumed as a conservative error estimates of $\Delta \mathrm{C}$. Then $\delta \Delta h_{ICESat}$ (Figure 1d), $\delta \Delta \mathrm{C}$ and $\delta \Delta m_{GRACE}$ (Figure 1b) can be formally propagated using Equations (7) and (12) to generate total uncertainties for ice sheet and GIA-related mass change as shown in Figure 6. Among the three error sources, $\delta \Delta h_{ICESat}$ has made the largest contribution to both $\delta \Delta m_{MB}$ and $\delta \Delta m_{GIA}$, so the uncertainties for both ice sheet and GIA-related mass change are largest in the Amundsen Sea Embayment of WA (shown in Figure 6) where the errors of ICESat are most notable (Figure 1d). The overall errors in per drainage basin are computed from the error covariance of all grid points within a given region (see Table 5). And, for more accuracy, we consider that the errors in the gridded data are spatially correlated, a Gaussian window of decorrelation-length scale of 300 km is used to approximate the covariance that is a function of the distance between grid points [91].

Our JIE result shows WA, EA and the AP changed in mass by −69 ± 24, 12 ± 16, and −27 ± 8 Gt/yr, respectively, over the period February 2003 to October 2009. The total mass loss over the entire AIS of −84 ± 31 Gt/yr, which falls within the middle range of the previously published altimetry, GRACE and combined estimates for the similar period (see Table 1). The total GIA mass change estimate in Antarctica of 54 ± 27 Gt/yr, which is very close to recent inversed GIA solutions: RATES (55 ± 15 Gt/yr) and REGINA (55 ± 23 Gt/yr).

The largest signal of ice loss occurred in the Amundsen Sea Embayment (DS21 and DS22) and adjacent Marie Byrd Land (DS20) of the WA with a rate of −123 ± 15 Gt/yr over the period February 2003 to October 2009. Meanwhile, the largest signal of ice gain occurred in the area surrounding the Filchner-Ronne ice shelf (DS1) of 34 ± 7 Gt/yr, which is considerable higher than previous solutions of GRACE, altimetry and combined estimates in the similar observed periods [6]. In addition, the Siple Coast (DS18 and DS19) of WA exhibits a clear positive trend (25 ± 8 Gt/yr), mainly driven by the thickening Kamb Ice Stream that can be reflected by the observed positive elevation change from ICESat (see Figure 1c). In the AP, the Northern Antarctic Peninsula (DS25 and DS26) shows a large mass loss of −29 ± 3 Gt/yr, which is in agreement with some recent studies based on different observations that revealed an ongoing rapid mass loss among the large glacier tributaries in this areas after the disintegration of the Larsen B Ice Shelf in 2002 [25,56,92]. The Southern Antarctic Peninsula (DS24 and DS27) is close to balance with a rate of 3 ± 5 Gt/yr, which is different from the previous estimates that suggested a moderate positive trend for the same period [6,25]. Mass change in the East Antarctic Ice Sheet displays two contrasted patterns: Coats Land, Dronning Maud Land, Enderby Land, and Princess Elizabeth Land (DS3-12) exhibit a total positive trend (61 ± 12 Gt/yr), and Wilkes Land, Victoria Land, and Oates Land (DS13-17), conversely, show a similar scale of mass loss with a trend of −49 ± 7 Gt/yr.

The comparison with the estimates derived from the single ICESat data with two different ETM models shows that the total ice mass change in the AIS with trends of −60 ± 17 and −123 ± 12 Gt/yr, respectively (see Table 5). There is an obvious difference, not only in magnitude, but also in spatial pattern, when using different ETM models to constrain altimetry data. This further confirms that the large differences among altimetry estimates from the uncertain ETM model.

The GRACE Mascon estimates show a good agreement with our JIE result in spatial pattern. However, there are some differences in the order of magnitude over the key regions where the ice mass loss are significant, such as the Amundsen Sea Embayment sector (DS20–22), DS13-17 of EA, and the Northern Antarctic Peninsula (DS25 and DS26). Table 5 shows that the mass loss rates in the Amundsen Sea Embayment (CSR: −101 ± 10, JPL: −101 ± 10, and GSFC: −101 ± 12 Gt/yr), DS13-17 of EA (CSR: −31 ± 8, JPL: −41 ± 8, and GSFC: −28 ± 8 Gt/yr), and Northern Antarctic Peninsula (CSR: −17 ± 1, JPL: −20 ± 1, and GSFC: −13 ± 2 Gt/yr) estimated by three GRACE Mascon solutions are all less than that by our JIE method (−123 ± 15, −49 ± 7, and −30 ± 3 Gt/yr, for the three regions, respectively). In the entire AIS, the mass loss rates from GRACE Mascon estimates (CSR: −50 ± 30, JPL: −71 ± 30, and GSFC: −51 ± 33 Gt/yr) also less than that from our JIE estimates (−84 ± 31 Gt/yr). While the Mascon solution has been shown to have better spatial resolution and more optimal removal of correlated error [90], the local estimates of mass flux still can be contaminated by the leakage-out errors across land/ocean boundaries [93]. We note that the JPL Mascon estimate, which employs a CRI filter to reduce leakage-out errors across coastlines [93], is much closer to our JIE estimate than CSR and GSFC estimates.

Due to the ice discharges all distributed in the edges of an ice sheet, the leakage-out errors of ID are much greater than that of SMB. Beyond the effect of leakage-out errors, signal attenuation errors due to incomplete power spectrum of GRACE gravity field also complicates GRACE estimates. Identically, the signal attenuation for the high-frequency signals of gravity field change (that is, those signals that the dramatic mass changes occurred in narrow and small region, e.g., rapid dynamic thinning of glacier or ice stream) are often more serious than the low-frequency signals (e.g., the signal of mass gain due to wide-range snowing). As shown in Figure 7, the simulated ID shows the higher degree terms, i.e., the shorter wavelengths in the signal spectra, compared to the signal of SMB. When spherical harmonic coefficients are limited to degree and order 60, the percentage of accumulated degree amplitudes in total for the SMB is 79%, but the ID is only 56%. This implies that the signal attenuation errors of ID are also worse than that of SMB. Thus, when using the limited spatial resolution of GRACE data to estimate the mass balance of AIS, the rates of ice mass loss are usually underestimated. This may be the reason for the difference between JIE and GRACE-only estimates.

5. Conclusions

In this study, we developed an improved JIE method for estimating GIA and ice mass change in Antarctica by a combination of multi-geodetic datasets including satellite gravimetry, altimetry, and GPS. The new JIE method uses the GPS uplift rates instead of FDM to constrain the conversion of elevation-to-mass changes that play a key role in the process of ICESat estimate. The various linear relationships between the elevation change of SMB and firn compaction were assumed based on a blocked processing. Through the improved JIE method, we could obtain present-day ice sheet mass change and GIA simultaneously. When assessing our GIA estimate and eight recent forward and inverse GIA solutions based on a comparison with two sets of elastic-corrected GPS-derived uplift rate data, we found that our GIA result is very close agreement with the GPS observations that shows good precision of WRMS values of 1.0 (Table 2) and 1.9 mm/yr (Table 4), respectively. The total GIA-induced mass change estimates for the AIS is 54 ± 27 Gt/yr (WA: 19 ± 21, EA: 34 ± 15, the AP: 1 ± 3) (Table 5), which is in line with many recent GPS calibrated GIA estimates (RATES: 55 ± 15 and REGINA: 55 ± 23 Gt/yr). The comparison results show some large differences in magnitude and spatial pattern exist between our GIA uplift rates and other GIA solutions, such as stronger uplift in the Amundsen Sea Embayment and weaker in the area surrounding the Filchner-Ronne ice shelf.

Over the period February 2003 to October 2009, the total ice sheet mass change in Antarctica estimated by our JIE method is −84 ± 31 Gt/yr (WA: −69 ± 24, EA: −2 ± 21, the AP: −27 ± 8), which falls within the middle range of two ICESat-only estimates (−60 ± 17 and −123 ± 12 Gt/yr), and greater than three GRACE-only estimates (CSR: −50 ± 30, JPL: −71 ± 30, and GSFC: −51 ± 33 Gt/yr) (see Table 5). Some of the disagreement between our JIE result and these single-data estimates may be due to the significant leakage and attenuation errors in GRACE estimates and ETM errors in altimetry estimates. While our mass balance estimate lies in the middle range of aggregated estimates from satellite altimetry (−43 ± 21 Gt/yr), gravimetry (−76 ± 20 Gt/yr) and the input–output method (−201 ± 82 Gt/yr) provided by the Ice sheet Mass Balance Inter-comparison Exercise (IMBIE) team [11] for the similar period 2003–2010 and it is less positive, but consistent within error bounds, compared with their averaged estimate (−105 ± 51 Gt/yr). The improved JIE method using multiple space-geodetic observations in this study shows its potential to reduce the uncertainties of ice sheet mass balance and GIA estimates in Antarctica.