List of symbols and units

2. Summary of approach to reconciling altimetry and gravimetry mass-changes

In the same way that satellite gravimetry measures changes in the ice mass on the Earth's crust and altimetry measures changes in the ice volume, the respective measurements include the effects of ongoing changes in the mass and volume (ΔM, ΔV) of the Earth under the ice. The fundamental concept of our approach for resolving the difference between GRACE- and ICESat-based estimates of ice mass changes is based on the realization that the respective mass and volume corrections are for the ΔM and ΔV of the same underlying material. The changing Earth material is illustrated schematically in Figure 2 as a distinct element (ΔM, ΔV) of the mantle, even though the actual material involved is spatially distributed in three dimensions within the mantle and may include elastic deformation of the crust. Furthermore, the required mass and volume corrections are both provided by the same dynamical models of the motion within the Earth caused by changes in the glacial loading (e.g. either Whitehouse, Ivins or Peltier GIA models). These models calculate the change in gravity caused by the ΔM and the vertical motion of the bedrock, dB/dt, caused by the ΔV. Although our concept emphasizes the mass and volume changes of the fluid mantle that force motion of the crust, it also includes the more rapid elastic deformation of the crust noted in Section 3, as do our bedrock motions calculated for the matching of GRACE and ICESat mass changes in Section 6.

For gravimetry, the correction (GIA_cor) is for the rate of change in gravity caused by the ΔM/Δt mass-change underneath the ice in units of the rate of mass change, which is essentially

$${\rm GI}{\rm A}_{{\rm cor}} = \Delta M/\Delta t = \rho _{{\rm earth}}\;\Delta V/\Delta t\comma \;$$

where ρ _earth is the relative density of mantle material involved in the ΔM/Δt change. For altimetry, the correction (dB _cor) to the mass changes calculated from changes in ice-sheet surface elevation (dH/dt) is made for the vertical motion of the bedrock (dB/dt) caused by the ΔV/Δt. The dB _cor is equal to ρ _ice ΔV/Δt, where ρ _ice is the relative density of ice, 0.91, that is typical of the density in deep ice cores rather than 0.917. (Throughout the paper, we use densities relative to that of water, 1 g cm⁻³, so densities are used as dimensionless quantities, being 1 for water). The use of ρ _ice is appropriate, because basal motion displaces solid ice and does not affect the density nor volume of the firn column. GIA_cor and dB _cor are used here as rates of mass change per unit area, so using ΔV/Δt = dB/dt × area, the GIA_cor and dB _cor per unit area are:

(1a)$${\rm GI}{\rm A}_{{\rm cor}} = \rho _{{\rm earth}}\;{\rm d}B/{\rm d}t$$

(1b)$${\rm d}B_{{\rm cor}} = \rho _{{\rm ice}}\;{\rm d}B/{\rm d}t$$

where dB/dt is positive upward. The GIA_cor and dB _cor corrections are always subtracted from the uncorrected observations. For example, positive values of GIA_cor and dB _cor corrections reduce mass gains or increase mass losses.

Although GIA_cor and dB _cor are defined as rates of mass change per unit area, as are the dM/dt rates of mass change, both are often written in units of rates of vertical mass change such as mm w.e. a⁻¹ without an explicitly associated area that requires multiplication by an area to get mass change per unit area. Examples are the modeled GIA given in units of mm w.e. a⁻¹ in Figure 3 and the cm w.e. a⁻¹ scale for the dM/dt in Figure 14, for which the implicit area for the latter is 1.0 cm² and the rate of mass change is 1.0 g a⁻¹ cm⁻² that is equivalent to 0.1 Gt a⁻¹ (100 km)⁻² as shown in the color scale in Figure 16 (and S2). We also use units of mm w.e. a⁻¹ for the average values of GIA_cor over specific regional areas as in column 2 of Table 1 with the regional rates of mass change in units of Gt a⁻¹ in column 3.

Fig. 3. Glacial Isostatic Adjustment (GIA) in mm w.e. a⁻¹, basal uplift (dB/dt) in mm a⁻¹, and RatioG/db equal to GIA/(0.91 × dB/dt) derived by three Earth models labeled Ivins, Whitehouse and Peltier (Whitehouse and others, Reference Whitehouse, Bentley, Milne, King and Thomas2012; Ivins and others, Reference Ivins2013; Peltier, Reference Peltier2014). Subsidence rate from glacial loading in the central part of EA ice sheet is largest in Whitehouse model and smallest in Ivins.

Table 1. Glacial Isostatic Adjustment (GIA) and uplift (dB/dt) from Ivins, Whitehouse and Peltier Earth models

GIA (mm w.e. a⁻¹) and dB/dt (mm a⁻¹) are regional-average rates and GIAcor (Gt a⁻¹) and dB _cor(Gt a⁻¹) = 0.91 dB/dt are area-integrated regional rates of the corrections. S _g and S _a sensitivities and RatioG/dB are also regional averages defined in the text.

^a See text about the calculation of S _g and RatioG/dB.

We define RatioG/dB as

(2)$${\rm Ratio}G/{\rm d}B\equiv {\rm GI}{\rm A}_{{\rm cor}}/{\rm d}B_{{\rm cor}} = \rho _{{\rm earth}}/\rho _{{\rm ice}} = \rho _{{\rm earth}}/0.91\comma \;$$

the gravimetry sensitivity (S _g)_md to bedrock motion as

(3)$$\lpar {S_{\rm g}} \rpar _{{\rm md}}\equiv{-}{\rm GI}{\rm A}_{{\rm cor}}/\lpar {{\rm d}B/{\rm d}t} \rpar \comma \;$$

where the subscript (md) indicates the GIA model used (Iv, Pe or Wh), and the altimetry sensitivity (S _a) to bedrock motion as

(4)$$S_{\rm a}\;\equiv{-}{\rm d}B_{{\rm cor}}/\lpar {{\rm d}B/{\rm d}t} \rpar .$$

The units of dB/dt are mm a⁻¹, the units of GIA_cor and dB _cor are Gt a⁻¹, the units of S _g and S _a are both Gt a⁻¹/mm a⁻¹, and

(5)$${\rm Ratio}G/{\rm d}B = S_{\rm g}/S_{\rm a}.$$

We include minus signs in the sensitivity definitions so a positive change in dB/dt (i.e. more uplift) causes the derived mass change to decrease and a negative change (i.e. more subsidence) causes it to increase. Whereas S _a is a geometric factor depending only on ρ _ice, area and the dB/dt from the GIA models, the S _g includes additional dependencies on the characteristics of the models.

S _g and S _a provide a straightforward linear relation for reconciling the differences in the GRACE and ICESat mass estimates by calculating the rate of uplift or subsidence (δB _0-md) needed to provide the GIA_cor and dB _cor corrections that bring the respective mass estimates into full agreement (i.e. [(dM/dt)_GRACE]_eq = [(dM/dt)_ICESat]_eq). The required uplift or subsidence, δB _0−md, relative to zero is given by:

(6)$$\eqalign{& \left[ {{\lpar {d{\rm M}/d{\rm t}} \rpar }_{{\rm GRACE}}} \right] _0 \, + \, \lpar {S_{\rm g}} \rpar _{{\rm md}}\;\delta B_{0-{\rm md}} = \left[ {{\lpar {{\rm d}M/{\rm d}t} \rpar }_{{\rm ICESat}}} \right] _0 \cr &+ S_{\rm a}\;\delta B_{0-{\rm md}} \cr & \delta B_{{\rm 0}-{\rm md}} = {{ { \big [ {{\lpar {{\rm d}M/{\rm d}t} \rpar }_{{\rm GRACE}}} \big ] }_{{\rm 0}}-{\big [ {{\lpar {{\rm d}M/{\rm d}t} \rpar }_{{\rm ICESat}}} \big ] }_{{\rm 0}}} \over {S_{\rm a}- {\lpar {S_{\rm g}} \rpar }_{\rm md}}}} $$

where [(dM/dt)_GRACE]₀ and [(dM/dt)_ICESat]₀ are the respective GRACE and ICESat measurements with zero GIA_cor and zero dB _cor applied as indicated by the subscript (0). The second subscript (md) indicates the GIA model used to calculate the gravity sensitivity, i.e. (S _g)_Iv, (S _g)_Wh or (S _g)_Pe for Ivins, Whitehouse or Peltier model. For example, δB _0-Iv indicates that (S_g)_Iv derived from the Ivins model of GIA and dB/dt was used with no dB _cor nor GIA_cor applied to the measured dM/dt.

The required uplift or subsidence (δB _adj-md) is also calculated relative to the GIA modeled

(7)$$\delta B_{{\rm adj}-{\rm md}} = {{ { \big [ {{\lpar {{\rm d}M/{\rm d}t} \rpar }_{{\rm GRACE}}} \big ] }_{{\rm md}}-{\big [ {{\lpar {{\rm d}M/{\rm d}t} \rpar }_{{\rm ICESat}}} \big ] }_{{\rm md}}} \over {S_{\rm a}- {\lpar {S_{\rm g}} \rpar }_{{\rm md}}}}\comma \;$$

where md is either Iv, Pe or Wh. The resulting GRACE and ICESat equalized mass changes using either δB _0-md or δB _adj-md are denoted

(8)$$\lpar {{\rm d}M/{\rm d}t} \rpar _{{\rm eq}-{\rm md}} = \big [ {{\lpar {{\rm d}M/{\rm d}t} \rpar }_{{\rm GRACE}}} \big ] _{{\rm eq}-{\rm md}} = \big [ {{\lpar {{\rm d}M/{\rm d}t} \rpar }_{{\rm ICESat}}} \big ] _{{\rm eq}-{\rm md}}.$$

As shown in Section 6, the differences among the three (dM/dt)_eq-md are small, and therefore the mass-change adjustment is largely independent of the particular GIA model used, even though relative differences among the modeled dB/dt are large.

Previously, Zwally and others (Reference Zwally2015) used a preliminary estimate of RatioG/dB = 6 with S _g = −55.7 Gt mm⁻¹ and S _a = −9.3 Gt mm⁻¹ for EA. For EA, the uncorrected GRACE and ICESat dM/dt of 61 Gt a⁻¹ and 136 t a⁻¹, respectively, came into agreement at 150 Gt a⁻¹ after adjusting the uplift by δB _adj-Ivx = −1.6 mm a⁻¹. (The subscript Ivx indicates that some parameters of the Ivins model run previously used for the calculation of dB _cor were not exactly the same as those for GIA_cor and S _g.)

The RatioG/dB also provides a basis for estimating the incremental long-term effect (δB′) on the rate of bedrock motion of a long-term dynamic ice thickening, (dH _d/dt)_obs, using

(9)$$\delta {B}^{\prime} = {-}\lpar {{\rm d}H_{\rm d}/{\rm d}t} \rpar _{{\rm obs}}/{\rm Ratio}\,G/{\rm d}B.$$

Eqn (9) is based on the hypothesis that the long-term dynamic response of the Earth's mantle to a continued long-term ice loading produces a corresponding downward flow of mantle material with mass and ice-volume changes in the ratio of RatioG/dB with respect to the ice loading. As noted in the introduction, the 15.9 mm a⁻¹ ice thickening observed in EA was interpreted as commencing at the beginning of the Holocene. Therefore, the corresponding estimated change in the long-term compensation rate was δB′ = −15.9/6 = −2.65 mm a⁻¹. This δB′ is 1.7 times larger than the δB _adj-Ivx = −1.6 mm a⁻¹ required for the mass-matching adjustment, which suggests that some but not all of the observed thickening may have been included in the model's ice loading history. In the following, we derive more accurate values of RatioG/dB and related parameters from the results of the three GIA models.

Finally, we note that our approach to resolving differences in the GRACE- and ICESat-based estimates of ice mass changes is fundamentally different from those proposed or applied by others. Wahr and others (Reference Wahr, Wingham and Bentley2000) proposed combining GLAS (ICESat) and GRACE measurements to slightly reduce the post-glacial rebound error in the GLAS mass-balance estimates. Shepherd and others (Reference Shepherd2012) ‘reconciled’ estimates of mass balance by taking the mean of selected estimates from three techniques (altimetry, gravimetry and input–output method). Riva and others (Reference Riva2009) combined ICESat and GRACE measurements using: (1) for ICESat data a surface snow density, ρ _surf, ranging from 0.32 to 0.45 for some ice areas, an intermediate (between firn and ice) density of 0.60 in other ice areas, and the density of pure ice (0.92) in areas where rapid changes in ice velocity have been documented; and (2) for GRACE data a rock density, ρ _rock, under grounded ice ranging from 3.4 to 4.0 in order to obtain the GIA impact on GRACE-derived estimates of mass balance of 100 ± 67 Gt a⁻¹. Martin-Español and others (Reference Martin-Español, Bamber and Zammit-Mangion2017) performed a statistical analysis combining satellite altimetry, gravimetry and GPS with prior assumptions characterizing the underlying geophysical processes and concluded that gains in EA are smaller than losses in WA, although we show in the Appendix that their use of a single density for estimating mass changes from elevation changes is not correct.

3. GIA and bedrock vertical motion (dB/dt)

The fundamental physical process involved in GIA is glacial loading/unloading that bends the Earth's crust and forces 3-D viscous flow in the underlying fluid mantle, as illustrated in Figure 2. Part of the elastic bending is relatively rapid, for example, as shown by GPS-measured seasonal vertical motions of the crust in response to the seasonal cycle of summer surface melting and water runoff from the ablation zone of Greenland (Nielsen and others, Reference Nielsen2013). In contrast, another part of the crustal bending occurs along with the viscous flow of the mantle with uplift and subsidence rates that decay exponentially, with response times depending on the viscosity, following changes in the glacial loading or unloading. For example, adjustments following the relatively abrupt demise of the Laurentide ice sheet ~10 K years ago are continuing with current uplift rates on the order of +15 mm a⁻¹ in central Canada (Peltier, Reference Peltier2004) and subsidence rates south of the former ice sheet, for example, −1.7 mm a⁻¹ in the Chesapeake Bay region (DeJong and others, Reference DeJong2015) from the hinge effect in the crustal bending. However, the response time strongly depends on the viscosity of the mantle, which is a principal parameter typically varied in the models to improve agreement with the constraining information available on uplift rates. For example, the analysis of Barletta and others (Reference Barletta2018) indicated a lower viscosity and faster uplift rate in the Amundsen Sea Embayment (ASE) in WA than previous studies. As noted in Barletta and others (Reference Barletta2018), viscosities ~10¹⁸–10¹⁹ Pa s correspond to decadal to centennial response times as shown in the AP (Nield and others, Reference Nield2014), in Southern Patagonia (Richter and others, Reference Richter2016) and in Alaska (Larsen and others, Reference Larsen, Motyka, Freymueller, Echelmeyer and Ivins2005). Viscosities ~10²⁰–10²¹ Pa s correspond to millennial response time scales and are typically used in global and Antarctic GIA models (Whitehouse and others, Reference Whitehouse, Bentley, Milne, King and Thomas2012; Ivins and others, Reference Ivins2013; Argus and others, Reference Argus, Peltier, Drummond and Moore2014; Peltier, Reference Peltier2014).

The density of the mantle ranges from ~3.4 to 4.4 in the upper mantle and from 4.4 to 5.6 in the lower mantle (e.g. Robinson, Reference Robinson2011). In contrast, the density of the crust is generally lighter, ranging from 2.2 to 2.9 similar to surface rocks such as granite, basalt and quartz. A somewhat common misconception is that the material involved in the GIA correction has the density of the surface or crustal rocks (e.g. ρ ≥ 2.7 in Zwally and Giovinetto, Reference Zwally and Giovinetto2011), rather than primarily the greater densities of the underlying fluid mantle.

In our analysis, we use the GIA and dB/dt uplift results provided by three GIA models labeled Ivins, Whitehouse and Peltier with the maps of the modeled data given in Figure 3. GIA models may have variations in model characteristics, parameters (e.g. mantle viscosities, mantle densities and crustal thickness), ice-loading histories, and their use of GPS and other data to constrain the model results, details of which are given in the references. In EA, the models generally show crustal subsidence in the central portions of the ice sheet with an uplift in the coastal regions and along the boundary with WA. This pattern of subsidence and uplift implies a radial outflow of mantle material from the central region, and inflow at the outer regions from both the central region and southward from the Southern Ocean. Over many millennia, the spatial and temporal variability of the glacial loading history produces a complex 3-D flow of the mantle, which on a continental scale at any given time can have flow in multiple directions at different depths with regions of convergence and divergence. During the short decadal times of satellite measurements, temporal variations in the mantle flow and the resulting uplift and subsidence rates are small.

The regional average values of GIA_cor (mm w.e. a⁻¹) and dB/dt (mm a⁻¹) for the Ivins, Whitehouse and Peltier models are given in Table 1 along with the regional GIA_cor (Gt a⁻¹) and dB _cor (Gt a⁻¹) mass corrections. [Note: total regional values are calculated as GIA_cor (Gt a⁻¹) = GIA_corr (mm w.e. a⁻¹) area (km²) 10⁻⁶ and dB _cor (Gt a⁻¹) = 0.91 dB/dt (mm a⁻¹) area (km²) 10⁻⁶.] The GIA_corr and dB _cor are both positive for positive dB/dt (i.e. uplift) and are subtracted from the measured gravity and altimetry mass changes (i.e. conventional usage). For the three models, the GIA_cor and dB _cor mass corrections for EA and WA are mostly comparable in magnitude, with the smaller area of WA (18% as large as EA) offset by its seven times greater average uplift.

For EA, the area of subsidence inland is largest in the Whitehouse model (Fig. 3) with subsidence more than −2 mm a⁻¹ in three locations and an area-averaged value of −0.19 mm a⁻¹ (subsidence), in contrast to average uplift rates of 0.42 mm a⁻¹ for Ivins and 0.60 mm a⁻¹ for Peltier (Table 1). For EA, the Ivins average GIA_cor is 1.9 mm w.e. a⁻¹ uplift and the regional dM/dt adjustment is −19.9 Gt a⁻¹. The GIA_cor is largest for the Peltier model at 3.1 mm w.e. a⁻¹ with a regional dM/dt adjustment of −31.5 Gt a⁻¹. For Whitehouse, the average GIA_corr is −0.9 mm w.e. a⁻¹ with a regional dM/dt adjustment of +8.8 Gt a⁻¹. Differences among the modeled GIA_corr are as large as 40 Gt a⁻¹ between the Peltier and Whitehouse models for EA (mostly EA1), 15 Gt a⁻¹ between Peltier and Ivins for WA (mostly WA2), and 45 Gt a⁻¹ between Peltier and Whitehouse for AIS.

The spatial variations of the model results and variations among the models are also illustrated by the profiles of GIA_cor, dB/dt and RatioG/db along longitudes 90°W and 90°E across WA and EA shown in Figure 4. Local-scale dB/dt differences along the transect are up to 4 mm a⁻¹ in the coastal WA1 for Peltier minus Ivins, up to 6 mm a⁻¹ in WA2 between Whitehouse minus Peltier, and up to 2 mm a⁻¹ in EA between both Ivins minus Whitehouse and between Peltier minus Whitehouse. For WA, the regional-average dB/dt difference among the models is largest for Peltier minus Ivins at 2 mm a⁻¹, and for EA, the difference is largest for Peltier minus Whitehouse at 0.79 mm a⁻¹ as shown in Table 5 in Section 6, along with our δB adjustments for comparison.

Fig. 4. Profiles of GIA, dB/dt and RatioG/dB from three dynamic Earth models Ivins (red), Peltier (green) and Whitehouse (blue) along 90°W across West Antarctica and along 90°E across East Antarctica extending into oceans. Singularities in RatioG/dB are avoided by calculating regional averages. Extent of continental ice is indicated by red lines.

A singularity in the RatioG/dB occurs where dB/dt approaches zero, changing from uplift to subsidence (or the reverse) ~70°S, 90°E in Whitehouse and Peltier models, ~76°S, 90°E in the Ivins model, and ~90°S in all three models. The location of the singularity forms an oval-shaped ring mostly in EA in the three models (Fig. 3). For the RatioG/db in Table 1, we avoid the effect of the singularity by using regional averages to calculate

(10)$$ {\lt \!\rm Ratio}G/{\rm d}B \! \gt _{{\rm regavg}} = \lt \!\! {\rm GIA} \!\! \,\gt _{{\rm regavg}}/ 0.91 \lt \! \! {\rm d}B/{\rm d}t \! \gt _{{\rm regavg}}\!\! \,.$$

The S _g and S _a sensitivities to bedrock motion in Table 1 are also regional averages. A special case occurs for the calculation of the RatioG/dB and S _g for the Whitehouse model data for EA, because the small values of the regional averages cause anomalous ratios in two drainage systems labeled DS16 and DS17. For these, we use area-weighted averages of <GIA>_DSavg and <dB/dt>_DSavg by drainage system (DS) to calculate the <RatioG/dB>_regavg, excluding DS16 and DS17 from the EA2 calculation. The regional average S _g are calculated using <RatioG/dB>_regavg S _a and <GIA>_regavg = <dB/dt>_regavg S _g.

The regional-average sensitivities to bedrock motion S _g and S _a in Table 1 are the only two parameters used in Section 6 to calculate the bedrock motions (i.e. the δB _0-md in Eqn (6) and the δB _adj-md in Eqn (7)) and the GIA_cor and dB _cor corrections needed for equalization of the ICESat and GRACE dM/dt. As previously noted in Section 2, whereas S _a is a geometric factor independent of the GIA model, S _g is model dependent. However, it is important to note that the differences (≈10%) among the three models in the values of S _g (and similarly for RatioG/dB) are relatively small compared to the relative differences in both the modeled dB/dt and the resulting GIA_cor in Table 1, which means that the results of the dM/dt equalization shown in Table 4 are not very dependent on which modeled S _g is used. For EA, EA1, EA2 and WA1, the equalized dM/dt are the same to two or three significant figures after rounding, and for WA2 differ by only 4% for some numerical reason.

Also shown in Table 1 are the inferred ρ _earth derived from the three GIA models by region according to Eqn (2), with values mostly in the range of 4–5. These ρ _earth are consistent with those in the fluid upper- to mid-mantle from Robinson (Reference Robinson2011) in our concept (Fig. 2), and larger than the typical densities of crustal rocks. Beyond showing this consistency in support of our concept, it is important to emphasize that these ρ _earth densities are not used for calculation of the δB for dM/dt equalization in Section 6.

4. Time-series of elevation and mass changes from ICESat and grace data

For ICESat, elevation time series, H_j,k(t_i), in 50 km grid cells (j, k) are created by a second stage of analysis following the along-track solution method described in Zwally and others (Reference Zwally and Giovinetto2011). In the first stage, the ICESat elevation measurements h (x_i, y, t_i), which are made at 172 m along-track spacings in the y-direction on repeat tracks lying within ±100 m (1σ) in the cross-track x-direction (c.f. Fig. 1 in Zwally and others, Reference Zwally2011) during 16 laser campaigns (Table 8) from Fall 2003 to Fall 2009, are first interpolated to equally-spaced reference points along track. The measured elevations depend on the cross-track position, x_i, and cross-track slope, α_r, as well as on real elevation variations with time according to

(11)$$h\lpar {x_i\comma \;y_r\comma \;t_i} \rpar = x_i\tan \alpha _r + t_i\lpar {{\rm d}h/{\rm d}t} \rpar _r + h_0\lpar {y_r\comma \;t_0} \rpar \comma \;$$

where h ₀ is the elevation at the position y_r on the reference track at t ₀. The use of constant (dh/dt)_r assumes that height changes at each reference point are a linear function of time over the period of measurement (e.g. 2003–2009). Equation (11) is solved by least-squares methods for the three parameters α_r, (dh/dt)_r and h ₀ at each reference point and other procedures (e.g. a seven-reference point solution using a calculated quadratic along-track slope) (Zwally and others, Reference Zwally and Giovinetto2011). Data outside of grounded-ice boundaries are excluded. Previously in Zwally and others (Reference Zwally2015), the (dh/dt)_r were averaged in 50 km cells creating multiyear-average [dH/dt]_j.k by cell, but those dH/dt are not used here.

In the second stage, a time series h_r(t_i) = h ₀(t ₀), h ₁(t ₁), h ₂(t ₂), … h ₁₆(t ₁₆) is created for each reference point using the cross-track slope α_r and x_i to correct each height for the cross-track displacement. Very importantly, any non-linear height variations with time (such as a seasonal cycle) relative to the constant (dh/dt)_r are retained in derived time series. The h_i(t_i) terms of the series in the 50 km cells are then averaged and 16 grid maps of the terms are created. Cells with any missing terms (i.e. 1, 2, … 16) are filled by interpolation creating a complete [H(t)]_j,k time series for each cell. The [H(t)]_j,k are then averaged (weighted by cell area) over drainage systems (DS) and ice-sheet regions creating H(t) for DS and for regions accounting for the splitting of partial cells at DS boundaries.

Calculation of mass changes (M(t)) from measured surface elevation changes (H(t)) requires correction for the elevation changes that do not involve changes in ice mass caused by variations in the rate of firn compaction (FC) as well as by the bedrock motion (e.g. Zwally and others, Reference Zwally2015) according to:

(12)$$\eqalign{H\lpar t \rpar = \, & H_{\rm d}\lpar t \rpar + H^{\rm a}\lpar t \rpar + C_{\rm A}\lpar t \rpar + C_{\rm T}\lpar t \rpar + \lpar {{\rm d}B/{\rm d}t} \rpar \;t \cr H_{\rm d}\lpar t \rpar = \, & H\lpar t \rpar -H^{\rm a}\lpar t \rpar -C_{\rm A}\lpar t \rpar -C_{\rm T}\lpar t \rpar -\lpar {{\rm d}B/{\rm d}t} \rpar \;t\comma \;} $$

where H _d(t) and H ^a(t) are the elevation components driven, respectively, by ice dynamics and by contemporary accumulation variations. The dB/dt used are the adjusted δB _0-Iv derived in Section 6 (Table 4). The C _A(t) and C _T(t) are the changes in surface elevation resulting from changes in the rate of FC driven by variations in accumulation rate (δA(t) = A(t)−<A(t)>) and by variations in firn temperature (T(t)). The C _A(t), C _T(t) and C _AT(t) ≡ C _A(t) + C _T(t) are calculated with an FC model (Li and Zwally, Reference Li and Zwally2015) using satellite measured surface temperatures and ERA-Interim re-analysis data for δA(t). The term H ^a_CA(t) = H ^a(t) + C _A(t) combines the direct height change from accumulation variations and the resulting accumulation-driven change in FC.

The δA(t) are also used to calculate the cumulative accumulation-driven mass change

(13)$$M_{\rm a}\lpar t \rpar = \int ^{\rm t}\delta A\lpar t \rpar \;{\rm d}t$$

and the cumulative accumulation-driven height change

(14)$$H^{\rm a}\lpar t \rpar = 1/\rho _{\rm s} \int ^{\rm t}\;\delta A\lpar t \rpar \;{\rm d}t\comma \;$$

where ρ _s = 0.3 is the density of new surface firn. Separation of the H _d(t) and H ^a(t) components of elevation change is essential for proper calculation of the total mass change, as well as the respective components of mass change caused by ice dynamics and by the δA(t) variations in the surface mass balance (SMB). The dynamic mass change is

(15)$$M_{\rm d}\lpar t \rpar = \rho _{{\rm ice}}H_{\rm d}\lpar t \rpar $$

using the well-defined ρ _ice = 0.91 and the total mass change is

(16)$$M\lpar t \rpar = M_{\rm d}\lpar t \rpar + M_{\rm a}\lpar t \rpar .$$

Calculation of M _a(t) (Eqn 13) very importantly avoids the need to use a firn density (ρ _a) that can only be known by first calculating M _a(t). As shown in Figure 8 in Zwally and others (Reference Zwally2015), the calculated ρ _a = ΔM _a/Δ(H ^a−C _A) according to Eqn (7) in Zwally and others (Reference Zwally2015) has a wide distribution over Antarctica from 0.2 to 0.9 with an average of 0.39 (see also maps of ρ _a variability in Fig. 17 and discussion in the Appendix). Therefore, a priori selection of appropriate single or multiple firn/ice densities (e.g. McMillan and others, Reference McMillan2014) is not possible due to the extensive spatial and temporal variabilities of the actual ρ _a, and because H _a and H _d have differing spatial variations in both magnitude and sign.

The ICESat measured H(t) for EA1, EA2 and EA and the other components of elevation change according to Eqn (12) are in Figure 5 including B(t) using the adjusted model values of dB/dt derived in Section 6 (Table 4). The corresponding M(t), M _d(t) and M_a(t) are in Figure 6. The series are fitted to a linear-quadratic-sinusoidal function [y(t) = A + B t + C t ² + D sin(ω t) + E cos(ω t)] with an annual period representing a seasonal cycle with the phase and amplitude selected by the fit. The derived values of most interest here are the linear terms, which we evaluate at the midpoint of the time period (at year 2006.0 for the period 2003 through 2008 and at year 2006.5 for 2003 through 2009). A clear seasonal-cycle is evident in the C _AT(t) firn-compaction term that is mainly driven by the seasonal cycle in temperature as shown in Figures 6b, c, h and i in Li and Zwally (Reference Li and Zwally2015). The seasonal cycles in both A(t) and H ^a(t) are very small even at the specific locations of South Pole and Law Dome, but their multi-year variability is large locally as shown in Figures 6a, d, g and j in Li and Zwally (Reference Li and Zwally2015). Significant multi-year to decadal scale variations in the regional averages for EA, EA1 and EA2 are evident in the H ^a(t) in Figure 5. Similarly, height and mass time series for WA1, WA2 and WA are given in Figures 7 and 8 also using our adjusted values of dB/dt.

Fig. 5. Components of elevation change from ICESat for EA, EA1 and EA2 from H _d(t) = H(t)−H ^a(t)−C _AT(t)−(dB/dt) × t with LQS fit through 2008 data only. Linear trends and the adjusted dB/dt used for B(t) are in Table 3. The dynamic H _d(t) is more linear than other elevation terms.

For the purpose of comparing the rates of change derived from the current time series method with those previously derived from the average-linear-change method in Zwally and others (Reference Zwally2015), the respective linear rates for 2003–2008 are given in Table 2 using the previous Ivins (dB/dt)₂₀₁₅ for both. The rates of change from the two methods are all in good agreement with the exception of those for the AP, for which the average method gave a loss of 28.8 Gt a⁻¹ versus only 10.3 Gt a⁻¹ for the time series method. Less significant are some of the differences in the accumulation-driven rates that may be due to the difference between using the LQS fit to the time series versus the linear-only fit for the previous method due to the more non-linear variation of the H ^a(t) and M ^a(t) as shown in the figures. Table 2 also shows the relation between the previously-used combined parameter, dH ^a_CAT/dt, which combines the FC and direct accumulation-driven height changes, and the separate FC parameter, dC _AT(t)/dt, and direct accumulation-driven height change, dH ^a/dt.

Table 2. ICESat elevation and mass change components from time series analysis for 2003–2008 using the Ivins (dB/dt) 2015 in Zwally and others (Reference Zwally2015)

Slopes are linear term at midpoint of time period (2006.0) from LQS fitting. Terms in italics with 2015 subscript are from the average-linear-change analysis for 2003–2008 from Zwally and others (Reference Zwally2015).

The linear trends of the time series in Figures 5–8 using the adjusted model values of dB/dt derived in the next section are in Table 3. These time series along with the values of their trends clearly illustrate: (1) the importance of the C _AT(t) correction for FC that does not involve changes in mass, and (2) the need to separate the elevation changes driven by the accumulation variations in surface mass balance to obtain the dynamic ice changes. In particular, the dynamic elevation, H _d(t), and dynamic mass series, M _d(t), are more linear than the total H(t) and M(t), especially in EA2, consistent with the expectation that decadal-scale dynamic changes are small in EA, while the M _a(t) varies on shorter time scales. In EA1 and EA, some non-linearity in the last year might be caused by errors in the non-linear accumulation term from the method used for interpolating monthly-accumulation-rates for the laser campaigns from annual averages. For this reason, the more linear 2003–08 period will be used for the adjustments of the GIA_cor and dB _cor in the next section, rather than the full 2003–2009.

Fig. 6. Components of mass change from ICESat for EA, EA1 and EA2 from M _d(t) = ρ _ice × H _d(t) from Figure 5 and $M_{\rm a}\lpar t \rpar = \int ^t\delta A\lpar t \rpar \times {\rm d}t$ with LQS fit through 2008 data only. Linear trends and the adjusted dB_cor applied are in Table 3. The dynamic M _d(t) is more linear than the total M(t).

Fig. 7. Components of elevation change from ICESat for WA, WA1 and WA2 from H _d(t) = H(t)−H ^a(t)−C _AT(t)−(dB/dt) × t with LQS fit through 2008 data only. Linear trends and the adjusted dB/dt used for B(t) are in Table 3. The dynamic H _d(t) is more linear than H(t) and other elevation terms.

Fig. 8. Components of mass change from ICESat for WA, WA1 and WA2 from M _d(t) = ρ _ice × H _d(t) from Figure 7 and $M_{\rm a}\lpar t \rpar = \int ^t\delta A\lpar t \rpar \times {\rm d}t$ with LQS fit through 2008 data only. Linear trends and the adjusted dB_cor and GIA_cor applied are in Table 3. The dynamic M _d(t) is more linear than the total M(t).

Table 3. ICESat elevation and mass change components for 2003–2008 and 2003–2009 from time series analysis using dB/dt equal δB _0-Iv (Table 4) and corresponding dB _cor from the matching of ICESat and GRACE dM/dt during 2003–2008 as described in Section 6

Slopes are linear term at midpoint of time period from LQS fitting (2006.0 for 2003–2008 and 2006.5 for 2003–2009).

For GRACE, time series are created using the mascon-solution methods described in Luthcke and others (Reference Luthcke2013), Luthcke and others (Reference Luthcke and Tedesco2015) and Loomis and others (Reference Loomis, Luthcke and Sabaka2019). Information on the GRACE Mascons and our data used is at https://earth.gsfc.nasa.gov/geo/data/grace-mascons.

5. Consistency of elevation changes from ERS1/2 1992–2001, ICESAT 2003–2008 and ENVISAT 2002–2010

In this section, we first review the compatibility and validity of our elevation and mass changes derived from ERS1/2 for 1992–2001 and ICESat 2003–08 presented in Zwally and others (Reference Zwally2015), including comparisons of the corrections for firn compaction and the accumulation-driven and dynamic-driven changes. Our results are in essential agreement with other studies that show an increasing mass loss in the Antarctic Peninsula and the coastal WA1, where large changes are observed over relatively small areas. In the interior WA2 and in EA, where the changes are small over large areas, our results are in agreement with some studies, but differ from others.

Previous unrefuted results showing ice-sheet growth in EA based on ERS/1/2 include Wingham and others (Reference Wingham, Ridout, Scharroo, Arthern and Shum1998), Davis and others (Reference Davis, Li, McConnell, Frey and Hanna2005), Zwally and others (Reference Zwally2005) and Wingham and others (Reference Wingham, Shepherd, Muir and Marshall2006). In particular for 1992–2003, Davis and others (Reference Davis, Li, McConnell, Frey and Hanna2005) found: ‘Using a near-surface snow density of 350 kg m⁻³, an average elevation change of 18 ± 3 mm a⁻¹ over an area of 7.1 million km² for the EA interior … corresponds to a mass gain of 45 ± 8 Gt a⁻¹’. However, the density of ice is the more appropriate density, because the increase in elevation has been shown to not be from contemporaneous increasing snowfall (Zwally and others, Reference Zwally2015). Therefore, the corrected result for their observed area would be a mass gain of 117 ± 18 Gt a⁻¹. For all of EA, their gain would be ~168 Gt a⁻¹, since the average elevation change south of the ERS coverage is similar to the northern area as shown in Figure 9 (and S1).

Fig. 9. Maps of dH/dt: (a) for 1992–2001 from ERS1/2, (b) for 2003–2008 from ICESat, and (c) for 2002.7–2010.7 from Envisat showing regional dH/dt for areas of common coverage. (Areas south of 81.6° coverage of ERS and Envisat and south of 86° of ICESat are interpolated in images.)

In comparison to Davis's (Reference Davis1997) 18 mm a⁻¹, our EA elevation changes for our calculated ERS coverage of 8.13 × 10⁶ km² are 10.7 mm a⁻¹ for 1992–2001 and 13.1 mm a⁻¹ for 2003–08 from ICESat, which are both smaller than Davis's (Reference Davis1997) 18 mm a⁻¹. For all EA (10.2 × 10⁶ km²), our changes are 11.1 mm a⁻¹ for 1992–2001 and 13.0 m a⁻¹ for 2003–08 (Table 2 in Zwally and others, Reference Zwally2015). Over Lake Vostok, the respective ERS and ICESat dH/dt of 20.3 and 20.2 mm a⁻¹ are in close agreement as shown in Table 1 and Figure 7 in Zwally and others (2015), which demonstrates the compatibility of the radar and laser altimetry over a flat area where the radar is not affected by its slope-induced errors. Furthermore, the accuracy of ERS altimetry for constructing time series is demonstrated by its measurement of global sea-level rise in good agreement with TOPEX and other ocean radar altimeters at the rate of 2.7 mm a⁻¹ (Scharroo and others, Reference Scharroo2013).

Our mass changes for EA from ERS1/2 and ICESat are also in very close agreement with each other at (dM/dt)₂₀₁₅ of 147 Gt a⁻¹, (dM _a/dt)₂₀₁₅ of −11 Gt a⁻¹ and (dM _d/dt)₂₀₁₅ of 136 Gt a⁻¹ as in Table 2 and in Table 5 of Zwally and others (Reference Zwally2015). Even though the respective measured dH/dt over EA differed by 1.9 mm a⁻¹, the long-term dynamic changes (dH _d/dt) were essentially the same at 15.8 and 15.9 mm a⁻¹ after correction for the FC and direct accumulation-driven changes (dC _T/dt and dH ^a_CA/dt) as shown in Table 2 of Zwally and others (Reference Zwally2015), which is consistent with the long-term dynamic stability of EA. In the EA1 and EA2 sub-regions, the elevation-change differences between periods are larger, likely due to variability in the accumulation-driven dH _a/dt. Overall, there is no apparent bias of the ICESat measurements compared to the ERS1/2 measurements.

In order to further establish the validity of the ICESat 2003–08 elevation and mass changes as the baseline for reconciling the GRACE and ICESat mass changes and the GIA_cor and dB _cor, we further review the methods and corrections employed in our data analysis and derivation of mass changes from elevation changes in the Appendix. We provide reasons why our results agree with some studies and differ from others. Among other things for ICESat laser altimetry, we review our ICESat inter-campaign biases and the Gaussian-Centroid (G-C) error correction including (1) the critical importance of our use of an independent determination of the motion of reference surface for bias determinations, and (2) the critical importance of using bias corrections determined using altimeter data with the G-C error applied (or vice versa) and the consequent substantial dH/dt error of 1.29 cm a⁻¹ if that compatibility is not maintained as noted on NSIDC ICESat-data website in 2013 (see Appendix). For example, Shepherd and others (Reference Shepherd2012) IMBIE-1 included (see Table S8 and data contributors in SOM) mass gain estimates from ICESat for EA of 118 ± 56 Gt a⁻¹ by Sorensen and Forsberg, 126 ± 60 Gt a⁻¹ by Smith, and a smaller gain of 86 ± 55 Gt a⁻¹ by Yi and Zwally, all of which were done before the G-C laser error correction was discovered, and therefore with campaign bias corrections consistently determined and when the effect of the biases was small as noted in the Appendix. In contrast, Shepherd and others (Reference Shepherd2018) IMBIE-2 did not include ICESat results from Forsberg nor Smith and at least some of the ICESat results from other data contributors (other than Zwally, Reference Zwally2013) had laser biases determined with the G-C inconsistency causing a significant dH/dt bias.

For radar altimetry, we review the major problem of the highly-variable (seasonally and interannually) penetration and backscatter depth and the correction methods used (or not used) by various investigators that are a likely source of residual errors. Whereas successful penetration-backscatter corrections were developed and applied for ERS1/2 radar altimetry by several investigators (as detailed in the Appendix), the problem became substantially more complex for Envisat and CryoSat data, because the linearly-polarized radar signals (oriented across-track on Envisat at 120° and CryoSat at 90°) interact with firn properties related to the direction of the surface slope and the relative directions differ significantly at track crossings. However, a successful radar penetration-correction method was developed for Envisat data by Flament and Rémy (Reference Flament and Rémy2012) using repeat-track analysis and waveform-dependent correction parameters, but has not been adopted in other studies. Specifically, Figure 1 in Flament and Rémy (Reference Flament and Rémy2012) for Envisat (2002.7–10.7) shows significant elevation increases over EA that are consistent with our ERS and ICESat increases.

In Figure 9 (and S1), we compare the dH/dt from: (a) ERS1/2 (1992–2001) from Figure 6a of Zwally and others (Reference Zwally2015), (b) ICESat (2003–2008) from Figure 6b (Zwally and others, 2015) and (c) from Envisat 2002.7-10.7 as mapped from data presented in Figure 1 of Flament and Rémy (Reference Flament and Rémy2012). We added a correction of +2.06 mm a⁻¹ to the Envisat dH/dt for the Point Target Response calibration that changed the derived MSL (mean sea level) trend from 0.463 to 2.52 mm a⁻¹ for mid-2002 to 2012 (Fig. 1 in https://earth.esa.int/web/guest/content/-/article/improvement-of-envisat-ra-2-reprocessed-data-v2-1).

In Figure 10, we compare the dH/dt averaged by DS from the three satellites in EA for their common coverage north of 81.6°S and in four DS in WA completely covered by all three. In WA1, the increasing ice loss from the coastal DS20, 21 and 22 is shown by the average dH/dt of −110, −151 and −177 mm a⁻¹ from the successive satellites. In WA, some of the features evident in Figure 9 (and S1) are: (1) the more extensive thinning extending inland in DS20, 21 and 22 during the later ICESat and Envisat periods compared to ERS1/2 and (2) thickening in the western part of DS21 and over much of DS19 draining into the Ross Ice Shelf in the later periods compared to thinning during ERS, which is likely due to the increased accumulation extending over the base of the AP and into WA as shown by the dM _a/dt from ERS and ICESat in Figures 10a and b of Zwally and others (Reference Zwally2015). That strong inter-period increase in accumulation also extended over WA1 offsetting part of the increase in dynamic thinning in the coastal DS20, 21 and 22.

Fig. 10. Average dH/dt from ERS1/2 1992–2001 (dashed red), ICESat 2003–2008 (solid blue) and Envisat 2002.7–20010.7 (dotted green) by DS and sub-regions for areas of common coverage. DS20, 21, 22, 19 and 4 to 16 are completely covered.

In EA, the large average dH/dt of −63 mm a⁻¹ from ICESat in the small part of DS2, for most of their common coverage over Berkner Island, is due to the large negative values on the southern point of the Island that apparently are not resolved in the radar altimetry. Similarly in the small coastal DS15 of EA, which has numerous alpine-like glaciers, the average dH/dt from ICESat is also notably more negative than from ERS and Envisat.

Over much of EA, the variability among the periods is also driven by accumulation variations as shown by the aforementioned dM _a/dt from ERS and ICESat. In EA1, supporting examples of this accumulation variability with corresponding variations in dH/dt between ERS and ICESat are: (1) the increase in dM _a/dt in the coastal DS4 following the ERS period, (2) the marked decrease in dM _a/dt in the adjacent coastal DS5 extending into the western part of DS6, (3) the increase in dM _a/dt in the eastern part of the coastal DS6 and in the adjacent coastal DS7, and (4) the increase in dM _a/dt in the mostly inland DS3 and the inland DS10.

For DS4, the increase in the average dH/dt from ERS to ICESat continued into Envisat as shown by the successive dH/dt of 27, 59 and 58 mm a⁻¹, and similarly for the decrease in DS5 with successive 66, 15 and 29 mm a⁻¹. In contrast, in DS8 the dH/dt of 68 mm a⁻¹ during ERS lowered to 24 mm a⁻¹ during ICESat and raised to 53 mm a⁻¹ during Envisat. Also, in DS10 the dH/dt of −3 mm a⁻¹ during ERS increased to 28 mm a⁻¹ during ICESat and decreased to 3 mm a⁻¹ during Envisat.

Overall of EA1 (DS2 to DS11), the successive average dH/dt are 13, 24 and 20 mm a⁻¹. Over EA2 (DS12 to DS17), the successive average dH/dt of 8, 1 and −4 mm showed a progressive decrease, which is mostly over the inland portions as shown in Figure 9 (and S1) and is likely due to a progressive shift in accumulation continuing the aforementioned increase of 21 Gt a⁻¹ in EA1 and decrease of 21 Gt a⁻¹ in EA2 between ERS1/2 1992–2001 and ICESat 2003–08. That is also consistent with the increasing mass gain in EA1 for several years after 2008 and the decreasing mass gain in EA2 after 2008 as shown by the M(t) from ICESat and GRACE in Figure 12 beginning around 2007 and continuing through 2010. For all of EA, the ERS to ICESat to Envisat variation is from 11 to 13 to 8 mm a⁻¹.

Fig. 11. ICESat and GRACE dM/dt for EA with no dB _cor or GIA_cor corrections (•) and with corrections from models of Ivins (), Peltier () and Whitehouse (). ICESat and GRACE equalized dM/dt mass changes range from 148 Gt a⁻¹ () using S _g = −52.6 Gt a⁻¹/mm a⁻¹ and δB ₀ = −1.99 mm a⁻¹ from Peltier model, to 151 Gt a⁻¹ () using S _g = −47.1 Gt a⁻¹/mm a⁻¹ and δB ₀ = −2.28 mm a⁻¹ from Ivins model, to 151 Gt a⁻¹ () using S _g = −45.7 Gt a⁻¹/mm a⁻¹ and δB ₀ = −2.36 mm a⁻¹ from Whitehouse model.

Fig. 12. M(t) time series for East Antarctica from ICESat (blue) and GRACE (red) using the equalizing dB _cor and GIA_cor listed in Table 3. The linear trends from LQS fits at the midpoints of 2003–2009, 2009–2012 and 2012–2016.3 also in Table 7a.

Considering the accumulation variability and the differing time periods, these dH/dt for EA from ERS1/2, ICESat and Envisat are consistent at the level of a few mm a⁻¹, and are all significantly more positive than the results of other studies. For example, the result from CryoSat data for 2010–13 for EA was only 1 ± 2 mm a⁻¹, from which they calculated a mass loss of 3 ± 36 Gt a⁻¹ for an area of 9 499 900 km² (McMillan and others, Reference McMillan2014); and from ERS, Envisat and CryoSat data for 1992–2017 was 6 ± 1 mm a⁻¹, from which they calculated a mass gain of only 16.3 ± 5.5 Gt a⁻¹ for an area of 9 909 800 km² (Shepherd and others, Reference Shepherd2019). However, Shepherd and others (Reference Shepherd2019) did not make nor assess the impact of corrections based on parameters not included in the satellite Level-2 data products, including the time-variable penetration corrections as made by Flament and Rémy (Reference Flament and Rémy2012) and alternate range retrackers (e.g. Helm and others, Reference Helm, Humbert and Miller2014; Nilsson and others, Reference Nilsson, Gardner, Sørensen and Forsberg2016), which along with their binary choice of firn or ice density affect their conclusions on ice-sheet elevation and mass changes, especially in EA.

6. Equalization of GRACE and ICESat mass change (dM/dt) determinations 2003–08

The required uplift or subsidence to bring the GRACE and ICESat dM/dt into agreement is calculated both relative to zero, giving δB _0-md according to Eqn (6), and relative to the modeled dB/dt, giving δB _adj-md according to Eqn (7). The resulting δB _0-md, δB _adj-md and the corresponding (dM/dt)_eq-md are given in Table 4 for the WA and EA regions and sub-regions. The linear solution for EA is also illustrated graphically in Figure 11. Corrections for increasing uplift linearly decrease the ice mass change according to the respective sensitivities: S _a = −9.29 Gt a⁻¹ per mm a⁻¹ for altimetry and S _g-Iv = −47.1, S _g-Pe = −52.6, or S _g-Wh = −45.7 Gt a⁻¹ per mm a⁻¹ for gravimetry (Table 1). For EA, the derived uplift adjustments are δB _adj-Iv = −2.70 mm a⁻¹, δB _adj-Pe = −2.59 mm a⁻¹ and δB _adj-Wh = −2.17 mm a⁻¹ with an average of −2.49 mm; and δB _0-Iv = −2.28 mm a⁻¹, δB _0-Pe = −1.99 mm a⁻¹ and δB _0-Wh = −2.36 mm a⁻¹ with an average of −2.21 mm a⁻¹. The corresponding (dM/dt)_eq-Iv is 150.5 Gt a⁻¹, the (dM/dt)_eq-Pe is 147.8 Gt a⁻¹ and the (dM/dt)_eq-Wh is 151.3 Gt a⁻¹ with an average of 149.9 Gt a⁻¹. The required δB _0-avg and δB _md-avg bedrock motions for mass matching from Table 4 for the EA and WA regions and their sub-regions are summarized in Table 5 along with their corresponding dB _cor and the dB/dt from the three GIA models for comparison.

Table 4. Values of adjustments to rate of uplift/subsidence needed to bring the ICESat and GRACE rates of mass change into agreement at [(dM/dt)_eq]_md

The δB _0-md is relative to zero uplift using dM/dt with no dB _cor nor GIA_cor applied and (δB _adj-md) is relative to the modeled dB/dt using dM/dt with the corresponding dB _cor and GIA_cor applied using S _a and (S _g)_md given in Table 1 in both cases.

dM/dt* is the linear term at year 2006.0 from LQS fit to regional M(t) series obtained using dB _cor for Ivins dB/dt (Table 1) + δB _adj-Iv.

Table 5. Bedrock motions δB _0-avg and δB _md-avg with their corresponding dB _cor that bring ICESat and GRACE dM/dt into agreement, dB/dt from Ivins, Peltier and Whitehouse models, maximum difference, δ(dB/dt)_max, among models

For EA, the range of (dM/dt)_eq-md among the three models is only 2.3% of their mean compared to the larger ranges of 17% in the δB _0-md and δB _adj-md. The smaller fractional difference among the (dM/dt)_eq-md occurs because of its primary sensitivity to the slope S _a, compared to the primary sensitivity of δB _adj-md and δB _0-md to the five times larger S _g (c.f. the solution in Fig. 11). Similarly for the sub-regions of EA, and for WA and its sub-regions, the differences among the (dM/dt)_eq-md are also small (≤2.5% range). Therefore, the (dM/dt)_eq-md vary <2.5% among the models used to equalize the GRACE and ICESat dM/dt's. Furthermore, for regional averages, it makes no difference whether δB ₀ and its corresponding dB _cor and GIA_cor are applied to [(dM/dt)_ICESat]₀ and [(dM/dt)_GRACE]₀ with no dB _cor nor GIA_cor applied, or whether δB _adj and its corresponding dB _cor and GIA_cor are applied to the [(dM/dt)_ICESat]_md and [(dM/dt)_GRACE]_md with their GIA modeled dB _cor and GIA_cor already applied.

Importantly, in the coastal WA1, the ICESat and GRACE measurements give nearly the same dM/dt of 95.2 and 96.0 Gt a⁻¹ with neither dB _cor nor GIA_cor corrections for volume and mass changes beneath the ice sheet. The required δB _0-avg for exact mass matching at 95.0 Gt a⁻¹ is only −0.35 mm a⁻¹, which is consistent with uplift in the part of WA1 nearest the coast and subsidence in the inner portions toward the ice divide with WA2. That spatial response is consistent with the differing histories of ice unloading in the coastal part of WA1 compared to the inner portion of WA1 that should be more similar to the inland WA2 with a different history of ice loading during the Holocene as discussed below.

Recent GPS measurements (Barletta and others, Reference Barletta2018) of land motion in the ASE of WA1 gave strong uplift rates (15–41 mm a⁻¹) at four locations that are much larger than the Peltier modeled dB/dt as shown in their Figure 1c. Those results imply errors in all other altimetry and gravimetry estimates of mass changes that necessarily use dB _cor and GIA_cor corrections from GIA models. However, the new GPS measurements (6 to −2 mm a⁻¹) at two locations a few hundred km to the Northeast outside of the ASE are small and closer to modeled values, suggesting that the strong uplift is confined to the ASE where recent grounding-retreat and ice thinning near the coast has occurred on Pine Island, Thwaites and Smith glaciers (e.g. Figure 4 in Zwally and others, Reference Zwally2015), at least to the East side of the ASE. Furthermore, the preferred GIA model of Barletta and others (Reference Barletta2018) in their Figure S12 shows that the associated gravitational uplift from the recent ASE ice changes is confined to an area of 300 km North-South and 800 km East-West.

The large uplifts measured in the ASE essentially have little or no effect on our results, because we do not use the GIA models' dB _cor or GIA_cor in our calculation of the δB _0-md adjustments that are made relative to the measured ICESat and GRACE dM/dt without any dB _cor or GIA_cor applied. Although the S _g sensitivities used in the adjustment Eqn (6) are calculated from the GIA-models, the differences among the modeled S _g are small (Table 1) and cause little differences among the resulting δB _0-md (Table 4). In WA1, with the aforementioned near equality of the uncorrected ICESat and GRACE dM/dt, the δB _0-md have a small range from −0.33 to −0.38 mm a⁻¹ (Table 4) due to small range in the S _g from −2.6 to −2.9 Gt a⁻¹/mm a⁻¹ (Table 1). The same comments apply when the δB _adj-md are calculated with Eqn (7), because the corrected mass changes are essentially the same for both calculations (Table 4). Similarly, additional uplift measurements in the other regions (WA2, EA, EA1 and EA2) would have little or no effect for the same reasons.

For WA2, the required δB _0-avg is −3.48 mm a⁻¹ (subsidence) in contrast to the average uplifts ranging from 3.00 to 5.42 mm a⁻¹ from the three GIA models due to their histories of post-LGM ice unloading over WA and the models' high mantle viscosities and millennial response times. That post-LGM ice loss was the principal driver of the Antarctic contribution to global mean sea-level rise that started ~15 ka BP and was mostly complete by ~9 ka BP, as shown in Figure 2 of Pollard and others (Reference Pollard, Gomez and DeConto2017) for their best-scoring model simulation. A community-based reconstruction of the AIS since the LGM (Bentley and others, Reference Bentley2014) shows the elevation at an inland site near the ice divide in WA2 as 200 m above present at 20 ka BP, decreasing to 150 m above present at 10 ka BP and to 50 m above present at 5 ka BP indicating a sustained history of ice unloading at a rate bringing it to the current elevation.

Evidence for a different ice-loading history in at least the lower elevations of WA2 after ~10 ka BP includes a 400 km Holocene readvance of the grounding line of the Ross Shelf from its simulated maximum inland retreat at 9.7 ka BP in Figure 3 of Kingslake and others (Reference Kingslake2018). That readvance is supported by their analyses of sediment cores. The ice-loading history also includes a smaller Holocene readvance from the simulated maximum-inland retreat of the Filchner-Ronne Ice Shelf grounding line at 10.2 ka BP, which is consistent with the low post-glacial rebound rates in the Weddell Sea that were attributed to a late Holocene ice-sheet readvance (Bradley and others, Reference Bradley, Hindmarsh, Whitehouse, Bentley and King2015). These grounding-line retreats and readvances on both sides of WA2 are also visible to a smaller extent (personal communication from David Pollard, 2020) in climate-driven ice-sheet modeling such as in Pollard and others (Reference Pollard, Chang, Haran, Applegate and DeConto2016) and Pollard and others (Reference Pollard, Gomez and DeConto2017).

When the grounding line in the Ross Ice Shelf retreated to its maximum inland position at 9.7 ka BP, it was ~400 km inland of its current position (Fig. 3 of Kingslake and others, Reference Kingslake2018) and at a location where the surface elevation is now ~500 m above sea level, indicating a thinning there of ~400 m at the maximum post-LGM retreat and subsequent thickening to its current elevation. Also, the ice at Siple Dome located ~70 km inside the present Ross-Ice-Shelf grounding line at 81.5°S, thinned by 350 m between 15 and 14 ka BP and its ice-divide advance began 2.5 ka BP, as derived from ice core data by Price and others (Reference Price, Conway and Waddington2007). The aforementioned reconstruction (Bentley and others, Reference Bentley2014) also shows the elevation at a location near Siple Dome to have been 350 m higher than present at 15 ka BP.

Therefore, during the Holocene readvance of the ice sheet to the current grounding-line position, the ice sheet thickened by ~300 m or more over a rather large area of the lower portions of DS18 and 19 in WA2. This suggests a mid-to-late Holocene increase in the ice loading of several hundred meters over a rather large area of DS18 and 19 in WA2. The thickening of the lower portions can also restrain the ice flow of the inland ice and lead to inland thickening as is now apparently occurring in DS18 as shown in Figure 9b (and S1b).

Therefore, one fundamental reason explaining why our results show subsidence in WA2 rather than the uplift of the three GIA models is the differing ice-loading history in WA2 associated with the Holocene readvance based on the above post-2013 findings not used in the ice load history of the earlier GIA models. A second reason is the findings of a significantly lower viscosity of the mantle in the ASE with implications for all of WA (Barletta and others, Reference Barletta2018), due in particular to the significantly faster GIA response times (as noted in Section 3) compared to those for the higher-viscosity GIA models. Since it has been ≈15 ka since the maximum post-LGM retreat and subsequent initiation of the readvance, the number of response times for the millennial response of the uplift from the post-LGM ice unloading to decay in the high-viscosity GIA models is ~5, which would reduce the exponentially decaying uplift from post-LGM ice unloading by a factor of ~7 × 10⁻³. However, for the lower viscosities of ~10¹⁸–10¹⁹ Pa s and their decadal to centennial response times, the corresponding reductions would be by a much-larger factor of <~2 × 10⁻²². Therefore, the primary on-going response should be subsidence from the later Holocene readvance that has been driven by the associated thickening of the grounded ice sheet. Subsidence is also consistent with our currently observed dynamic thickening in WA2.

As noted in Section 2, the RatioG/dB also provides a basis for estimating the incremental long-term effect on the rate of bedrock motion (δB′) of a long-term dynamic ice thickening using Eqn (9). Values of δB′ (calculated using the RatioG/dB from Table 1) for the EA, EA1, EA2 and WA2 regions with currently-observed dynamic thickening are −3.23, −2.50, −3.97 and −10.76 mm a⁻¹ as listed in column 4 of Table 6. These δB′ are compared to the δB _0-avg adjustments (i.e. relative to zero dB/dt) for the averages of the three δB _0-md (listed in column 5, taken from column 3 of Table 4) showing that the δB′ estimated from the regional dynamic thickenings are 1.2–3.1 times larger (column 6) than the required δB _0-avg and 1.1–1.7 times larger than the δB _adj-Iv (last column). For both comparisons, the estimated long-term response (subsidence) to ice thickening of the magnitudes observed is larger than the required bedrock motion adjustments (subsidence) for mass matching.

Table 6. Estimated bedrock motion, δB′, caused by the observed dynamic thickening

The δB′ equal to −(dH _d/dt)_obs/RatioG/dB is larger than the bedrock motion (both δB _0-avg and δB _adj-Iv) needed to bring ICESat and GRACE dM/dt into agreement.

For the ICESat analysis, we use the Ivins dB/dt grid plus the regional δB _adj-Iv in the calculation of the dynamic H _d(t) in 50 km cells using Eqn (12); this retains the spatial variation of the modeled dB/dt to which the regional average δB _adj-md are added. The dB _cor = −(S _a) [(dB/dt)_Iv + δB _adj-Iv] are listed in Table 3. The spatial variation is included in the ICESat grid maps of dM/dt and dM _d/dt (Figs 16 and S2), but is not distinguishable. The adjusted ICESat M _d(t) and M(t) are calculated using Eqns (15) and (16) for the height and mass series shown in Figures 5–8. To obtain the adjusted GRACE M(t), we calculated the regional GIA_cor = −(S_g)_Iv (δB ₀)_Iv for the EA1, EA2, WA1 and WA2 sub-regions using values from column 4 of Table 1 and column 3 of Table 4. The GIA_cor for EA and WA are the sums of their respective sub-regions. The GIA_cor listed in Table 3 are applied (subtracted) to the GRACE M(t) that had no correction already applied. The corrected M(t) for ICESat 2003–2009 and GRACE 2003–2016.5 for the EA and WA regions and sub-regions are shown in Figures 12 and 13, and for AP and AIS in Figures 14 and 15.

Fig. 13. M(t) time series for West Antarctica from ICESat (blue) and GRACE (red) using the equalizing dB _cor and GIA_cor listed in Table 3. The linear trends from LQS fits at the midpoints of 2003–2009, 2009–2012 and 2012–2016.3 also in Table 7a.

Fig. 14. M(t) time-series for Antarctic Peninsula from ICESat (blue) and GRACE (red) using dB _cor = −0.5 a⁻¹ and GIA_cor = −2.3 Gt a⁻¹ from Ivins2. The linear trends from LQS fits at the midpoints of 2003–2009, 2009–2012 and 2012–2016.3 are also in Table 7a. *The −10 Gt a⁻¹ from LQS is replaced by −29 Gt a⁻¹ from average-linear change analysis in AIS sum in Figure 15 and Table 7a.

Fig. 15. M(t) time series for Antarctica from ICESat (blue) and GRACE (red). The linear trends from LQS fits at the midpoints of 2003–2009, 2009–2012 and 2012–2016.3 are also in Table 7a.

Fig. 16. ICESat maps for 2003–2008, (a) dM/dt, (b) dM _d/dt and (c) dM_a/dt using dB/dt equal to IvinsdB/dt + δB _adj. Rates are linear terms of LQS fits at year 2006.0. *Rates for AP from average-linear-change analysis.

Fig. 17. Maps of the calculated firn density ρ _a = ΔM _a/Δ(H ^a−C _A) (see text following Eqn (16)) associated with the accumulation driven dM _a/dt mass changes for (a) 1992–2001 and (b) 2003–08, showing the large spatial and temporal variations.

7. Antarctic regional changes 1992–2016

The regional changes during 1992 through 2016 are examined for four periods as labeled in Table 7a: (1) the first is the 1992–2001 period of ERS1/2 measurements, (2) the second is the 2003–08 period of ICESat and GRACE measurements and mass-change matching, (3) the third is the 2009–11 period of GRACE measurements, and (4) the fourth is the 2012–16 period of GRACE measurements. The second, third and fourth periods are chosen for the analysis of the linear trends in the ICESat and GRACE M(t) series, because (a) the 2003–08 period has near linear trends and is used for ICESat GRACE mass change matching and (b) there are discernable changes in the slopes of the M(t) series around 2009.0 and around 2012.0 in both the EA and WA regions as well as their sub-regions.

Table 7a. Summary of linear rates of mass change (dM/dt) from ERS1/2, ICESat and GRACE for select periods during 1992–2016

^a AP ICESat and all ERS from average-linear-change analysis (Zwally and others, Reference Zwally2015) with ERS using adjusted dB _cor from Table 3.

The linear mass trends from LQS fits at the midpoints of the 2003–08, 2009–11 and 2012–2016 periods are in Table 7a and discussed in the next section. The ERS1/2 dM/dt for 1992–2001 are from Zwally and others (Reference Zwally2015) with the (dB _cor)₂₀₁₅ in Table 2 replaced with the dB _cor in Table 3. The differences between successive periods are given as the deltas in Table 7b along with a comparison of the deltas as fractions of the average-annual SMB. The ICESat and ERS1/2 estimates of uncertainties are made using the methods detailed in the Appendix of Zwally and others (Reference Zwally2015) and for GRACE in Luthcke and others (Reference Luthcke2013).

Table 7b. Summary of changes (delta) in the linear rates of mass change between periods compared to the annual SMB

^a SMB from Giovinetto and Zwally (Reference Giovinetto and Zwally2000) and by drainage systems and regions in Zwally and others (Reference Zwally2015).

^b These delta are adjusted by 11 Gt a⁻¹ to account for the difference between the 11 Gt a⁻¹ larger ICESat dM/dt from the prior average-linear-change analysis (see Table 2) as was also used for ERS1/2.

In the EA1 sub-region, the rate of mass gain more than doubled from 79 Gt a⁻¹ during 2003–08 to 196 Gt a⁻¹ beginning around the 2009.0. That increased gain of 117 Gt a⁻¹ occurred mostly in the Queen Maud Land portion of EA1, where Shepherd and others (Reference Shepherd2012) and Medley and others (Reference Medley2017) reported mass gains and accumulation increases, but it did not persist after 2012 when the EA1 gain reduced to 88 Gt a⁻¹, close to the prior rate of 79 Gt a⁻¹. In the EA2 sub-region, successive decreases of 10 and 16 Gt a⁻¹ helped to reduce the overall gain in EA from a high of 257 Gt a⁻¹ during 2009–11 to 134 G a⁻¹ during 2012–16, which is similar to the prior rates of 150 Gt a⁻¹ during 2003–08 and 161 Gt a⁻¹ during 1992–2001.

As the mass gain doubled in EA, the mass loss in the coastal WA1 doubled from 95 Gt a⁻¹ during 2003–08 to 214 Gt a⁻¹ during 2009–11. WA1 includes DS22 with the Pine Island Glacier, DS21 with the Thwaites and Smith Glaciers, and DS20 with grounded ice discharging into Getz ice shelf along the coast of Marie Byrd Land. The increased loss of 119 Gt a⁻¹ in WA1 was enhanced by a 39 Gt a⁻¹ reduction in the mass gain in the mostly inland WA2 bringing the total WA loss rate to 187 Gt a⁻¹ during 2009–11. In the last period, the loss from WA1 reduced by 49 Gt a⁻¹ as the gain in WA2 increased by 22 Gt a⁻¹, which together reduced the overall loss from WA to 116 Gt a⁻¹ during 2012–16. This reduced loss is still significantly greater than the 8 Gt a⁻¹ loss rates during 1992–2001 and 29–26 Gt a⁻¹ during 2003–08 from ICESat and GRACE. In the Antarctic Peninsula, the rate of loss increased from 9 Gt a⁻¹ during 1992–2001, to 29–24 Gt a⁻¹ from ICESat and GRACE during 2003–08, followed by losses of 36 Gt a⁻¹ during 2009–11 and 30 Gt a⁻¹ during 2012–16.

The spatial distributions of the rates of dynamic-driven mass changes (dM _d/dt), the accumulation-driven changes (dM _a/dt) and the total mass changes (dM/dt) during 2003–08 are shown in Figures 16a–c (and S2a, S2b and S2c). The magnitude and spatial distribution of the dM/dt and dM _d/dt are very similar and differ from the dM _a/dt that are generally smaller and more spatially variable. Areas of significant dynamic thinning are mostly in the coastal areas of WA1, parts of the AP and on the Totten Glacier at 115°E in DS13 of EA2. In DS22 of WA1 with the Pine Island outlet glacier, the dynamic thinning and negative dM/dt both extend inland close to the ice divide except for an area of positive rates in the Southeast corner. Similarly in DS21, dynamic thinning and negative dM/dt extend inland to the ice divide, except for an area of small positive rates in the Southwest corner (see Figs S2a and S2b). Inland dynamic thinning is also inland of the Mercer and Whillans Ice Streams in the Eastern part of DS17 of EA2 and the Western part of DS18 in WA2 inland of the Ross Ice Shelf.

As shown in Figure 16b (and S2), dynamic thickening (discussed further in the next section) extends over most of EA, WA2 and DS27 in the AP. A marked area of dynamic thickening is in DS18 of WA2, inland from the Kamb Ice Stream that stagnated 150 years ago (Joughin and others, Reference Joughin, Tulaczyk, Bindschadler and Price2002), and has a gain of 29 Gt a⁻¹ for 2003–08 adjusted for the new bedrock motion.