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

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.


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 MB = SMB + 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).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 (ρ ice = 917 kg/m 3 ) as an approximation of ETM model.Riva et al. [18] constructed a combined surface density distribution ρ sur f as a more appropriate choice than ρ 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 ρ ice or ρ sur f , 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.

The JIE Method to Estimate GIA
The currently available geodetic observations for AIS are mass change rates from GRACE, ∆m GRACE , surface elevation change rates from satellite altimetry, ∆h SA , and bedrock uplift rates from GPS, ∆h GPS .According to the principle of measurement and respective sensitivity of these observations, we can decompose ∆m GRACE , ∆h SA , and ∆h GPS into the following components [22,64]: where ∆m and ∆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 (∆m SMB and ∆h SMB ), ice discharges (∆m ID and ∆h ID ), GIA (∆m GI A and ∆h GI A ), firn compaction (∆h f c ), and elastic response of the lithosphere to loading by ice sheet mass change due to SMB and ice discharges (∆h ela ).
Considering Equations ( 1)-(3), we can express the rates of mass balance ∆m MB using GRACE and altimetry observations: 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: where ∆h GI A and ∆h ela in the initial estimate or zeroth iteration are equal to 0 [16], and in the following iteration are derived as: where the superscripts n and m represent degree and order in spherical harmonic domain, h ve n /k ve n 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]), ρ ave denotes the average density of the Earth, and h ela n 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, ∆h GI A and ∆h ela can be interpolated to the locations of all GPS sites as ∆h interp GI A and ∆h interp ela , and the GIA uplift rates obtained from GPS measurements are calculated by ∆h GPS GI A = ∆h GPS − ∆h interp ela , then the weight right mean square (WRMS) [66] of GIA uplift rates can be derived from the following: where q = 1 (δ∆h 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]. error propagation (∆ℎ = ∆ℎ + ∆ℎ ).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 ∆ in Equation ( 9) here should be estimated by the GRACE data with the time span as same as GPS data (2003-2013), and ∆ is calculated by the JIE method.

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 • × 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.

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 (δ∆h GPS = δ∆h 2 ela + δ∆h 2 GI A ).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 ∆m nm GRACE in Equation ( 9) here should be estimated by the GRACE data with the time span as same as GPS data (2003-2013), and ∆m nm GI A is calculated by the JIE method.

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) and compare their respective GIA results (Figure

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) 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 where ∆ can be estimated by SMB model which derived from the RACMO2 regional atmospheric climate model [83], and ∆ℎ = ∆ ,  is the density of snow in Anarctica (  = 350 kg m ), only ∆ℎ 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 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 where ∆m SMB can be estimated by SMB model which derived from the RACMO2 regional atmospheric climate model [83], and ∆h SMB = ∆m SMB ρ snow , ρ snow is the density of snow in Anarctica (ρ snow = 350 kg m −3 ), only ∆h f c 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 ∆h SMB and ∆h f c approximately, i.e., ∆h f c ≈ µ∆h SMB , then Equation ( 11) can be simplified as where 12) is equal to the suggestion of Wahr et al. [16], so ∆C can be regarded as a correction term for the initial approximation.It is easy to obtain µ = −0.619by the nonlinear least square principle ( min||WRMS||) using iterative operation, and then we can use this µ 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 µ 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 ∆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 ∆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 ∆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.
where ∆C = ∆ (1 − (1 + ) ).When  = − 1 or ∆ = 0, Equation ( 12) is equal to the suggestion of Wahr et al. [16], so ∆C can be regarded as a correction term for the initial approximation.It is easy to obtain  = −0.619by the nonlinear least square principle (‖‖) using iterative operation, and then we can use this  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  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 ∆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 ∆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 ∆ (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 Figures 3e,4a, it is clear that the large discrepancies between the uplift rates of GIA and GPS are mainly distributed in the regions where ∆ above 100 kg m −2 yr −1 , 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 (∆ > 500 kg m −2 yr −1 )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.According to the Figures 3e and 4a, it is clear that the large discrepancies between the uplift rates of GIA and GPS are mainly distributed in the regions where ∆m SMB above 100 kg m −2 yr −1 , 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 (∆m SMB > 500 kg m −2 yr −1 )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 ∆m SMB > 100 kg m −2 yr −1 in DS19), WA5 (the regions where ∆m SMB > 100 kg m −2 yr −1 in DS18), WA6 (the regions where ∆m SMB ≤ 200 kg m −2 yr −1 in DS1), and WA7 (the regions where ∆m SMB ≥ 200 kg m −2 yr −1 in DS1) (see Figure 4b) based on the difference of accumulation rates and drainage system.EA is divided into EA1 (the regions where ∆m SMB > 100 kg m −2 yr −1 in the western of DS2 and DS17), EA2 (the regions where ∆m SMB > 100 kg m −2 yr −1 in the eastern of DS17), EA3 (the regions where where ∆m SMB > 100 kg m −2 yr −1 in DS16), EA4 (the regions where ∆m SMB > 100 kg m −2 yr −1 in DS14 and DS15), EA5 (the regions where ∆m SMB > 100 kg m −2 yr −1 in DS13), EA6 (the regions where ∆m SMB > 100 kg m −2 yr −1 in DS11 and DS12), EA7 (the regions where ∆m SMB > 100 kg m −2 yr −1 in DS8 and DS9), EA8 (the regions where ∆m SMB > 100 kg m −2 yr −1 in DS5, DS6, and DS7), and EA9 (the regions where ∆m SMB > 100 kg m −2 yr −1 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 ∆m SMB ≤ 100 kg m −2 yr −1 is divided into the two subregions: WA8 (the regions where ∆m SMB ≤ 100 kg m −2 yr −1 in DS18, DS19 and the western of DS17) and EA10 (the regions where ∆m SMB ≤ 100 kg m −2 yr −1 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 µ for each subregion (see Table 3) 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 µ 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 µ 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 µ values.Thus, our GIA result is highly dependent on the accuracy of GPS observations and the reliability of the JIE method itself.

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.
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. 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]).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 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)  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].

GIA WA (7) EA (10) AP (3) All (20) WM WRMS WM WRMS WM WRMS WM WRMS
From Figures 3e and 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 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.

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 µ 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 µ values using the low rock densities ρ rock of 2000 kg/m 3 in the AP and 2200 kg/m 3 in WA by following spatial least square principle: 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 µ 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).
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.
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).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 ∆h GI A and ∆h ela can be negligible due to the effect is far smaller than δ∆h ICESat , the uncertainties of ∆C is difficult to assess due to insufficient error information for SMB model and µ values, so an standard deviation of 20% of the value for each grid point is assumed as a conservative error estimates of ∆C.Then δ∆h ICESat (Figure 1d), δ∆C and δ∆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, δ∆h ICESat has made the largest contribution to both δ∆m MB and δ∆m GI A , 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].

Regions
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  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.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.[32], and the SMB grid point values are shown in Figure 4a, all grid values are expanded into spherical harmonic coefficients of degree and order 360.The accumulated degree amplitudes [94] of n are estimated over a spectral band from 0 to n, and the total amplitudes are estimated over the spectral band from 0 to 360.The green line is the location of n = 60.

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  4a, all grid values are expanded into spherical harmonic coefficients of degree and order 360.The accumulated degree amplitudes [94] of n are estimated over a spectral band from 0 to n, and the total amplitudes are estimated over the spectral band from 0 to 360.The green line is the location of n = 60.

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.

Figure 1 .Figure 1 .
Figure 1.The Gravity Recovery and Climate Experiment (GRACE) observations for (a) mass change trends and (b) corresponding uncertainties in units of equivalent water height, and the Ice, Cloud andFigure 1.The Gravity Recovery and Climate Experiment (GRACE) observations for (a) mass change trends and (b) corresponding uncertainties in units of equivalent water height, and the Ice, Cloud and land Elevation Satellite (ICESat) observations for (c) elevation change trends and (d) corresponding uncertainties over the period February 2003 to October 2009.
Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 22 land Elevation Satellite (ICESat) observations for (c) elevation change trends and (d) corresponding uncertainties over the period February 2003 to October 2009.

Figure 2 .
Figure 2. The Global Positioning System (GPS) observations for crustal uplift rates and corresponding uncertainties in Antarctica over the time span of 2003-2013.The radius of the circle denotes the uncertainty (as the radius is larger, the uncertainty becomes smaller).Gray area denotes lowprecipitation zone (LPZ).The GPS data are from Sasgen et al. [79].

Figure 2 .
Figure 2. The Global Positioning System (GPS) observations for crustal uplift rates and corresponding uncertainties in Antarctica over the time span of 2003-2013.The radius of the circle denotes the uncertainty (as the radius is larger, the uncertainty becomes smaller).Gray area denotes low-precipitation zone (LPZ).The GPS data are from Sasgen et al. [79].

Figure 4 .
Figure 4. (a) The mean SMB (∆ ) over the time span of 2003-2009 and the Antarctic drainage systems; (b) The divided 20 subregions in Antarctica.

Figure 4 .
Figure 4. (a) The mean SMB (∆m SMB ) over the time span of 2003-2009 and the Antarctic drainage systems; (b) The divided 20 subregions in Antarctica.

Figure 6 .
Figure 6.The uncertainties for (a) ice sheet and (b) GIA-related mass change derived from joint inversion estimation, in units of equivalent water height.

Figure 6 .
Figure 6.The uncertainties for (a) ice sheet and (b) GIA-related mass change derived from joint inversion estimation, in units of equivalent water height.

Figure 7 .
Figure 7.The percentage of accumulated degree amplitudes in total for the ice discharge (ID; blue) and the SMB (red).The ID grid point values are simulated from the Antarctic ice velocity [32], and the SMB grid point values are shown in Figure 4a, all grid values are

Figure 7 .
Figure 7.The percentage of accumulated degree amplitudes in total for the ice discharge (ID; blue) and the SMB (red).The ID grid point values are simulated from the Antarctic ice velocity [32], and the SMB grid point values are shown in Figure4a, all grid values are expanded into spherical harmonic coefficients of degree and order 360.The accumulated degree amplitudes[94] of n are estimated over a spectral band from 0 to n, and the total amplitudes are estimated over the spectral band from 0 to 360.The green line is the location of n = 60.

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

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.∆h f irn and ∆m f irn represent the surface height and mass change of the firn, and ρ α is an assigned density.The unit of weight right mean square (WRMS) is mm/yr, only δ∆h GPS GI A used in the calculation of WRMS.The number in parenthesis denotes the total quantity of GPS sites within each region.

Table 3 .
[79]optimal values of µ for each subregions.The recomputed µ values using low rock densities of 2000 kg/m 3 in the Antarctic Peninsula (AP) and of 2200 kg/m 3 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.