Mass balance of the Antarctic ice sheet 1992–2016: reconciling results from GRACE gravimetry with ICESat, ERS1/2 and Envisat altimetry

Abstract

GRACE and ICESat Antarctic mass-balance differences are resolved utilizing their dependencies on corrections for changes in mass and volume of the same underlying mantle material forced by ice-loading changes. Modeled gravimetry corrections are 5.22 times altimetry corrections over East Antarctica (EA) and 4.51 times over West Antarctica (WA), with inferred mantle densities 4.75 and 4.11 g cm⁻³. Derived sensitivities (S_g, S_a) to bedrock motion enable calculation of motion (δB₀) needed to equalize GRACE and ICESat mass changes during 2003–08. For EA, δB₀ is −2.2 mm a⁻¹ subsidence with mass matching at 150 Gt a⁻¹, inland WA is −3.5 mm a⁻¹ at 66 Gt a⁻¹, and coastal WA is only −0.35 mm a⁻¹ at −95 Gt a⁻¹. WA subsidence is attributed to low mantle viscosity with faster responses to post-LGM deglaciation and to ice growth during Holocene grounding-line readvance. EA subsidence is attributed to Holocene dynamic thickening. With Antarctic Peninsula loss of −26 Gt a⁻¹, the Antarctic total gain is 95 ± 25 Gt a⁻¹ during 2003–08, compared to 144 ± 61 Gt a⁻¹ from ERS1/2 during 1992–2001. Beginning in 2009, large increases in coastal WA dynamic losses overcame long-term EA and inland WA gains bringing Antarctica close to balance at −12 ± 64 Gt a⁻¹ by 2012–16.

List of symbols and units

1. Introduction

The major portion of the East Antarctic (EA) ice sheet (Fig. 1) has been dynamically stable for many millennia, as currently shown by the 800 000-year-old-basal ice at Dome C (Jouzel and others, 2007) and the million-year ice at marginal blue ice areas (Sinisalo and Moore, 2010). Surviving through major cycles of climate change with cold-glacial and warm inter-glacial periods, changes in the marginal extent and the inland thickness of the EA ice sheet have been small compared to changes in the West Antarctic (WA) and Greenland ice sheets (e.g. Denton and Hughes, 1981; Denton, 2011; Mackintosh and others, 2011; Bentley and others, 2014; Pollard and others, 2017). In contrast to EA, much of WA is grounded 1000 m below sea level, has a maximum surface elevation of 2000 m (only half of EA), may be susceptible to dynamic instabilities, and has a more uncertain and complicated long-term history, including its major retreat after the Last Glacial Maximum (LGM) and partial re-advance of the grounding lines during the Holocene (Kingslake and others, 2018).

Fig. 1. Antarctic ice sheet regions and drainage systems. East Antarctica (EA) is divided into EA1 (DS2 to DS11) and EA2 (DS12 to DS17). The Antarctic Peninsula (AP) includes DS24–27. West Antarctica (WA) is divided into WA1 (Pine Island Glacier DS22, Thwaites and Smith Glaciers DS21, and the coastal DS20) and WA2 (inland DS1, DS18, and DS19 and coastal DS23). Includes grounded ice within ice shelves and contiguous islands.

In general, variations in the total mass (M(t)) of the Antarctic ice sheet (AIS) are the sum of short-term (≲decades) accumulation-driven variations (M _a(t)) in the surface mass balance and sub-decadal to millennial dynamic variations (M _d(t)). Dynamic changes in ice velocity occur for various reasons such as changes in ice-shelf back-pressure, basal sliding or long-term changes in accumulation rate that cause changes in ice thickness and surface slope that drive long-term changes in velocity.

The mass balance of the EA ice sheet has been significantly affected by long-term changes in snowfall, as shown by the 50–200% increases in accumulation beginning after the LGM ~15 ka BP and continuing through the Holocene as derived from ice cores (Siegert, 2003). That continuing long-term accumulation increase was a key factor supporting the interpretation of the 1.59 cm a⁻¹ thickening of the EA ice sheet, derived from both ERS1/2 (1992–2001) and ICESat (2003–08) altimetry measurements, as persistent long-term dynamic thickening with a dynamic mass gain of 147 Gt a⁻¹ (Zwally and others, 2015). This Holocene ice growth in EA is also consistent with the evidence of Holocene glacier advances from the EA ice sheet through the Transantarctic mountains into the Dry Valleys (Stuiver and others, 1981; Denton and Wilson, 1982). In contrast, the most marked area of contemporary dynamic changes and coastal ice thinning in EA is on Totten glacier at 116°E (Zwally and others, 2005; Pritchard and others, 2009; Li and others, 2016).

As analysis methodologies for both satellite altimetry and gravimetry have advanced in recent years, the largest remaining difference in mass-balance estimates (Shepherd and others, 2012, 2018; Hanna and others, 2013; Zwally and others, 2015; Hanna and others, 2020) has been for the EA ice sheet (Fig. 1). The agreement has been generally better in WA. However, the behavior in the coastal portion (WA1) is dominated by dynamic losses and is markedly different from the mostly inland portion (WA2) that has significant dynamic thickening, of which some is similar to the thickening in EA (Zwally and others, 2015).

The mass balances of both EA and WA are also significantly affected by decadal variations in accumulation such as those between the 1992–2001 ERS1/2 period and the 2003–08 ICESat period: (a) the regional shift in EA of +21 Gt a⁻¹ in EA1 and −21 Gt a⁻¹ in EA2, and (b) an increase in WA snowfall that offset 50% of the increased losses of 66 Gt a⁻¹ from enhanced dynamic thinning on accelerating outlet glaciers in WA1 and the Antarctic Peninsula (AP) (Zwally and others, 2015). Therefore, determination of both the short-term accumulation-driven and the long-term dynamic-driven components of ice-sheet mass balance is critically important for understanding the causes of changes on various time scales and the ice sheet's ongoing- and future-contributions to global sea-level change.

In their Figure 3, Hanna and others (2020) show the variation in the estimates of Antarctic dM/dt from 1990 to 2018 obtained by the three principal methods (altimetry, gravimetry and input–output method), which is updated from a similar figure in Hanna and others (2013). For EA, those reviews, as well as the multi-investigator results (Shepherd and others, 2018) from the second Ice Sheet Mass Balance Inter-comparison Exercise (IMBIE), clearly show outliers on both sides of the means, indicating that disparate estimates (~100 Gt a⁻¹) of the EA mass balance have not been properly resolved.

In this paper, we focus first on the mass balance of the EA and WA ice sheets and on resolving the differences between gravimetry-based and altimetry-based estimates of the balance during the 2003–08 period of overlapping measurements. Our method is motivated by indications that a principal residual uncertainty in prior estimates was due to errors in their respective corrections (GIA_cor and dB _cor) for changes in the volume and mass of the Earth underneath the ice caused by changes in the glacial loading on the crust and mantle (Fig. 2). The process of adjusting to changes in the glacial loading is commonly called Glacial Isostatic Adjustment (GIA). For the case of full isostatic (hydrostatic) equilibrium, the vertical motion of the bedrock (dB/dt) would be zero. However, under large ice masses, the long-term isostatic state is never actually fully reached as the glacial loading continually changes and the underlying fluid mantle hydrodynamically adjusts to the changes in the gravitational forcing.

Fig. 2. Ice sheet of thickness, T, lying on Earth's crust and underlying fluid mantle. For long-term isostatic equilibrium (~10 ka) with constant ice thickness the depth of the depression would be D ≈ ρ _ice /ρ _mantle × T, which is 600 m for T = 3000 m and ρ _mantle = 4.5, and the dB/dt would be zero. As the glacial loading, T(t), on the Earth's crust continually changes, the underlying viscous mantle hydrodynamically adjusts over centuries to millennia. Illustration is for an increasing ice thickness that induces a downward motion of the crust (i.e. dB/dt < 0), outward mantle flow and mantle thinning. For this case, the GRACE senses the gravitational changes of the increasing ice mass minus the decreasing mantle mass (ΔM) under the satellite. ICESat senses the increase in ice thickness minus the downward motion of the crust and mantle caused by the change in mantle volume (ΔV).

We use the results of three GIA models (Whitehouse and others, 2012; Ivins and others, 2013; Argus and others, 2014; Peltier, 2014), which are dynamic models of the Earth's crust (lithosphere) and the underlying fluid mantle forced by changes in the glacial loading. The models provide the gravitational signal for GIA_cor and the vertical motion of the bedrock for calculation of dB _cor. We derive the respective sensitivities (S _g and S _a) of the GIA_cor and dB _cor to the bedrock motion from the models. We then use S _g and S _a to derive the bedrock motion (i.e. δB ₀) needed to match the mass changes during 2003–2008 from GRACE and ICESat without any GIA_cor and dB _cor applied, and alternatively the adjustments (i.e. δB _adj) to the modeled bedrock motions with model-based GIA_cor and dB _cor applied. We then use the adjusted GIA_cor to extend the mass-change time series using GRACE gravimetry data through to 2016.

The introduction in Whitehouse and others (2012) presents a thorough review of prior calculations of GIA corrections applied to GRACE data and the effect of residual model errors on the estimates of ice mass balance. Constraints on the models are provided by the measurements of relative sea level and GPS measurements of crustal motion, which are also used for the estimation of residual errors (Whitehouse and others, 2012). In EA where fewer constraining measurements have been made, especially inland on the vast area of the ice sheet, the errors are likely to be largest. The review by Hanna and others (2013) noted: ‘… several key challenges remain …, changes in ice (sheet) extent and thickness during the past millennium are poorly known, and typically not included in GIA models, despite the fact that they can dominate the present-day rebound signal, especially in regions of low mantle viscosity’.

In order to further establish the validity of the ICESat 2003–08 elevation and mass changes as the baseline for reconciling the GRACE GIA_cor and the ICESat dB _cor, we review the methods and corrections employed in our data analysis and derivation of elevation and mass changes in section 5 and the Appendix. We first review the compatibility and validity of our elevation and mass changes derived from ERS1/2 for 1992–2001 and ICESat 2003–08 as presented in Zwally and others (2015), showing how those results agreed with other studies. We include a new comparison of the corrected dH/dt derived from ERS1/2 and ICESat with the corrected dH/dt derived from Envisat radar altimetry from Flament and Rémy (2012). That comparison shows the essential agreement of the dH/dt measured over EA by the four satellites with differing instrumentation over 19 years from 1992 through 2010 at the level of a few mm a⁻¹. Flament and Rémy (2012) developed unique methods for correction of the highly-variable (seasonally and interannually) sub-surface radar penetration not used in other Envisat nor CryoSat radar altimeter studies, which as detailed in the Appendix is a principal reason why other studies have differed from ours.

Overall, as GIA modeling has advanced in recent years, the results remain fundamentally dependent on knowledge of the history of the glacial loading, especially in the vast inland parts of the AIS where physical constraints from measurements are not feasible and knowledge of loading history is limited. Furthermore, there has been a lag in model incorporation of new information on the glacial loading as it becomes available from paleo-rates of ice accumulation derived from ice cores (e.g. Siegert, 2003; Siegert and Payne, 2004) and radar layering (Vieli and others, 2004), from our altimetry results and conclusions on inland ice growth (Zwally and others, 2005), and from information on Antarctic glacial geology and ice modeling (e.g. Stuiver and others, 1981; Denton and Wilson, 1982; Bradley and others, 2015; Pollard and others, 2017; Kingslake and others, 2018). In our conclusions, we discuss how our regional values of δB ₀ are consistent with current knowledge and interpretation of the history of glacial loading and with new findings of a lower mantle viscosity in WA by Barletta and others (2018).

We apply the derived dB _cor corrections to the ERS1/2 results for 1992–2001 as well as the ICESat results for 2003–2009, and the corresponding GIA_cor to the GRACE results for 2003 through to the beginning of 2016, thereby showing the mass-balance variations for all of the AIS and by regions over 24 years.

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, 2012; Ivins and others, 2013; Peltier, 2014). 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 (2015) 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 (2000) proposed combining GLAS (ICESat) and GRACE measurements to slightly reduce the post-glacial rebound error in the GLAS mass-balance estimates. Shepherd and others (2012) ‘reconciled’ estimates of mass balance by taking the mean of selected estimates from three techniques (altimetry, gravimetry and input–output method). Riva and others (2009) 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 (2017) 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, 2013). 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, 2004) and subsidence rates south of the former ice sheet, for example, −1.7 mm a⁻¹ in the Chesapeake Bay region (DeJong and others, 2015) 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 (2018) 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 (2018), viscosities ~10¹⁸–10¹⁹ Pa s correspond to decadal to centennial response times as shown in the AP (Nield and others, 2014), in Southern Patagonia (Richter and others, 2016) and in Alaska (Larsen and others, 2005). Viscosities ~10²⁰–10²¹ Pa s correspond to millennial response time scales and are typically used in global and Antarctic GIA models (Whitehouse and others, 2012; Ivins and others, 2013; Argus and others, 2014; Peltier, 2014).

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, 2011). 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, 2011), 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 (2011) 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 (2011). 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, 2011) 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, 2011). Data outside of grounded-ice boundaries are excluded. Previously in Zwally and others (2015), 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, 2015) 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, 2015) 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 (2015), the calculated ρ _a = ΔM _a/Δ(H ^a−C _A) according to Eqn (7) in Zwally and others (2015) 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, 2014) 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 (2015). 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 (2015). 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 (2015), 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 (2015)

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 (2015).

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 (2013), Luthcke and others (2015) and Loomis and others (2019). 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 (2015), 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 (1998), Davis and others (2005), Zwally and others (2005) and Wingham and others (2006). In particular for 1992–2003, Davis and others (2005) 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, 2015). 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 (1997) 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 (1997) 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, 2015). 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, 2013).

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 (2015). 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 (2015), 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 (2012) 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 (2018) 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, 2013) 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 (2012) 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 (2012) 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 (2015), (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 (2012). 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 (2015). 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, 2014); 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, 2019). However, Shepherd and others (2019) 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 (2012) and alternate range retrackers (e.g. Helm and others, 2014; Nilsson and others, 2016), 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, 2018) 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, 2015), at least to the East side of the ASE. Furthermore, the preferred GIA model of Barletta and others (2018) 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 (2017) for their best-scoring model simulation. A community-based reconstruction of the AIS since the LGM (Bentley and others, 2014) 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 (2018). 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, 2015). 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 (2016) and Pollard and others (2017).

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, 2018) 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 (2007). The aforementioned reconstruction (Bentley and others, 2014) 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, 2018), 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, 2015) 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 (2015) 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 (2015) and for GRACE in Luthcke and others (2013).

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 (2000) and by drainage systems and regions in Zwally and others (2015).

^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 (2012) and Medley and others (2017) 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, 2002), and has a gain of 29 Gt a⁻¹ for 2003–08 adjusted for the new bedrock motion.

8. Discussion and conclusions

During 1992–2016, the AIS changed from a positive mass balance of over 100 Gt a⁻¹, which was reducing sea-level rise by 0.3 mm a⁻¹, to a state of balance close to zero by 2014. The mass balance successively changed from a gain of 144 ± 61 Gt a⁻¹ during 1992–2001, to 95 ± 26 Gt a⁻¹ during 2003–08, to 34 ± 85 Gt a⁻¹ during 2009–11 and to −12 ± 64 Gt a⁻¹ during 2012–2016 (Table 7a). These rates of change suggest an acceleration of −50 Gt a⁻¹ decade⁻¹ during 1992 through 2006, −138 Gt a⁻¹ decade⁻¹ during 2006 through 2010.5, and −105 Gt a⁻¹ decade⁻¹ during 2014.5 through 2016. The changes during 2003–2016, shown in Figures 12 and 13, are driven by the acceleration of outlet glaciers in the coastal WA1 with the marked increase in the dynamic loss of 119 Gt a⁻¹ beginning near the end of 2008 that reduced by 49 Gt a⁻¹ after 2012. The increased dynamic loss near the end of 2008 was enhanced by a 39 Gt a⁻¹ decrease in the gain in WA2 that was followed by an increase of 22 Gt a⁻¹ near the beginning of 2012 driven by accumulation variations. During the same periods, the mass gain in EA increased by 107 Gt a⁻¹ near the end of 2008 followed by a decrease of 124 Gt a⁻¹ after 2012 for a small net gain decrease of 16 Gt a⁻¹ during 2003–2016, also driven by contemporary accumulation variations.

Although the acceleration rates reported by Rignot and others (2019) of −48 Gt a⁻¹ decade⁻¹ during 1979–2001 and −134 Gt a⁻¹ decade⁻¹ during 2001–2017 are consistent with our acceleration estimates above, the mass-balance results from their input–output method are generally more negative throughout their analysis periods. Although their methods of interpolation or extrapolation for areas with unobserved output velocities have an insufficient description for the evaluation of associated errors, such errors in previous results (Rignot and others, 2008) caused large overestimates of the mass losses as detailed in Zwally and Giovinetto (2011).

For all of EA, our mass gains (Table 7a) are essentially unchanged at 161 ± 50 Gt a⁻¹ during 1992–2001, 150 ± 28 Gt a⁻¹ during 2003–08, 257 ± 76 Gt a⁻¹ during 2009–11 and 134 ± 58 Gt a⁻¹ during 2012–2016. The 107 Gt a⁻¹ increase during 2009–11 was driven mostly by a temporary accumulation increase in EA1. The results of Smith and others (2020) for EA, as derived from ICESat to ICESat-2 crossover differences centered, respectively, on 2006.3 and 2019.0, are a gain of 115 ± 21 Gt a⁻¹ (after adding 25 Gt a⁻¹ for the δB _adj−Iv of −2.70 mm a⁻¹ as applied to our EA results). The estimated average of our results for the same period (extrapolating the gain of 134 Gt a⁻¹ during 2012–2016 through to 2019.14) is 164 ± 55 Gt a⁻¹ in agreement with Smith within the overlapping uncertainties. In marked contrast, the IMBIE Team (2018) values for EA are much smaller during both 2007–12 at 48 ± 38 Gt a⁻¹ and 2012–17 at −3 ± 30 Gt a⁻¹ (both adjusted by +25 Gt a⁻¹) as they were during 1992–2007.

For all of WA, our mass losses (Table 7a) increased from −8 ± 20 Gt a⁻¹ during 1992–2001, to −28 ± 15 Gt a⁻¹ during 2003–2008, to −187 ± 23 Gt a⁻¹ during 2009–11, to −116 ± 24 Gt a⁻¹ during 2012–16. In coastal WA1, successively-increasing losses of −59 ± 12, −95 ± 9, −214 ± 18 and −165 ± 15 Gt a⁻¹ were partially offset by persistent gains in inland WA2 of 51 ± 14, 68 ± 9, 27 ± 15 and 49 ± 19 Gt a⁻¹. Compared to the IMBIE Team (2018), our losses for WA prior to 2009.0 are ~30 Gt a⁻¹ less negative (after adding 10.53 Gt a⁻¹ to IMBIE for the δB _adj−Iv of −6.16 mm a⁻¹ as is applied to our WA results). After 2009.0, our results for WA of −187 ± 23 Gt a⁻¹ during 2009–11 and −116 ± 24 Gt a⁻¹ during 2012–2016 are comparable to both the IMBIE (2018) results of −137 ± 27 Gt a⁻¹ during 2007–2012 and −148 ± 26 Gt a⁻¹ during 2012–2017 (adjusted by 10.53 Gt a⁻¹) and the Smith and others (2020) results of −158 ± 10 Gt a⁻¹ during 2006.3–2019.0 (adjusted by 10.53 Gt a⁻¹).

Also of interest are the 16-year mass changes in the Antarctic ice shelves from Table 1 in Smith and others (2020) showing a loss of 76 ± 49 Gt a⁻¹ from WA and 14 Gt ± 28 Gt a⁻¹ from AP while shelves in EA gained 106 ± 29 Gt a⁻¹. Similarly, during 1992–2002, the ice shelves in WA lost 57 and 38 Gt a⁻¹ in the AP, while shelves in EA gained 142 Gt a⁻¹ as obtained from ERS radar altimetry corrected for radar penetration and temperature-dependent firn compaction (Zwally and others, 2005). Although these inter-decadal changes are small (−20 Gt a⁻¹ in WA, +24 Gt a⁻¹ in AP, −36 Gt a⁻¹ in EA and −32 Gt a⁻¹ overall), they are consistent with significant changes in some drainage systems. The small changes also add support to the validity of our ERS ice-sheet results, because the same altimetry methods were used over both grounded and floating ice during each of the 1992–2001 and 2003–19 periods. Furthermore, the small magnitude of the changes suggests the lack of major inter-decadal ice-shelf thinning or thickening in Antarctica.

Significant regional mass-change rates over Antarctica ranging from tens of Gt a⁻¹ to over 100 Gt a⁻¹ occurred during 1992–2016 as shown by the deltas in Table 7b, including both regional increases in the rates of mass loss and increases in the rates of mass gain. Over all of Antarctica, the total inter-period changes are all increases in mass loss ranging from 40 to 60 Gt a⁻¹, because some regional increases in mass gains only partially or temporarily offset regional increases in losses. Over the 24 years (1992–2016), the total increase in loss is 109 Gt a⁻¹ bringing the total AIS essentially into balance at −12 ± 64 Gt a⁻¹ by 2014. As listed in Table 7b, the ratios of the changes (deltas) to the SMB provide information on the relative significance of the inter-period variations.

In both WA1 and AP, the dynamic-driven variations are more persistent and sometimes larger relative to the SMB than the sub-decadal accumulation-driven variability. In the first interval (between 1992–2001 and 2003–08), when the inter-period change in WA1 was smallest at −36 Gt a⁻¹ (−16% of SMB), the mass loss rate from DS22 with Pine Island Glacier doubled from 12 to 29 Gt a⁻¹ while the loss rate from DS21 with the Thwaites and Smith Glaciers increased from 40 to 51 Gt a⁻¹ and the loss rate from DS20 with glacier flow into the Getz Ice Shelf increased from 7 to 16 Gt a⁻¹ (Zwally and others, 2015). Studies of increases in glacier thinning and acceleration of discharge velocities on Pine Island and Thwaites glaciers in WA1 during ~1992 to the early 2000s include Rignot and others (2002), Thomas and others (2004) and Wingham and others (2009). In the second interval (2003–08 to 2009–11), the loss rate from WA1 further increased by 119 Gt a⁻¹ (−54% of SMB), likely due to continued acceleration of glacier discharge. In contrast, in the third interval (2009–11 to 2012–16.5), the loss rate from WA1 decreased by 49 Gt a⁻¹ (+22% of SMB), likely due to an unidentified slowing of glacier discharge. In the AP, the 20 Gt a⁻¹ (−10% of SMB) loss-rate increase in the first interval (1992–2001 to 2003–08) was related to the acceleration of glaciers, mainly following the collapse of the Larsen B ice shelf (Pritchard and Vaughan, 2007; Rott and others, 2011; Shuman and others, 2011). That was followed by a smaller loss-rate increase of 11 Gt a⁻¹ (−6% of SMB) between 2003–2008 and 2009–2011, which was followed by a loss-rate decrease of 6 Gt a⁻¹ (+3% of SMB) between 2009–2011 and 2012–2016.

In EA, and the EA1 and EA2 sub-regions, the inter-period variations of delta/SMB(%) (Table 7b) are mostly only a few percent, which are typical of short-term atmospheric-driven variations in accumulation rates. A marked exception is the aforementioned 117 Gt a⁻¹ increase (+25% of SMB) in EA1 between 2003–08 and 2009–11, followed by a 108 Gt a⁻¹ decrease (−23% of SMB) between 2009–11 and 2012–16. However, the net change in EA1 is only a small increase of 9 Gt a⁻¹ (+2% of SMB) during the ICESat to ICESat-2 period of 16 years. In EA2 during the same times, the rate of mass gain decreased by 10 Gt a⁻¹ (−1% of SMB) between 2003–08 and 2009–11, followed by a 16 Gt a⁻¹ decrease (−2% of SMB) between 2009–11 and 2012–16 giving a net decrease of 26 Gt a⁻¹ (−4% of SMB) during the ICESat to ICESat-2 interval. Therefore, the total accumulation-driven effect for all of EA was a decrease in the rate of 17 Gt a⁻¹ (−1.5% of SMB) from 2003 to 2016, which is most likely not the cause of the mass gain of 90 Gt a⁻¹ in EA during the ICESat to ICESat-2 interval that was reported by Smith and others (2020).

Furthermore in EA during both of the earlier periods (1992–2001 and 2003–08), there were small negative accumulation anomalies of −11.6 ± 6 Gt a⁻¹ (i.e. −1% of SMB) compared to the 27-year mean from 1982, which justified our conclusion that the mass gain in EA in those periods was dynamic ice thickening and not due to increases in contemporaneous snowfall (Zwally and others, 2015). In the last column of Table 7b, we show the net changes in the rates over the 25 years (1992–2016) by including the rates of mass change during the first interval (i.e. between 1992–2001 ERS and 2003–08 ICESat). These inter-decadal changes are very small at −17 Gt a⁻¹ (−1% of SMB) for EA, −2 Gt a⁻¹ (0% of SMB) for EA1 and −15 Gt a⁻¹ (−2% of SMB) for EA2, further supporting our conclusion that the observed mass gains in EA are from on-going dynamic thickening and not from current trends in accumulation.

The observed dynamic thickening (2003–08) extends over most of EA (Figs 16b and S1) and much of the inland WA2, as was also shown for 1992–2001 in Figure 11a of Zwally and others (2015). The average dynamic thickening from the 2003–08 ICESat analysis is 16.4 mm a⁻¹ over EA and 47.7 mm a⁻¹ over WA2 (Table 3). Comparable thickening rates were previously obtained from the average-linear-change analysis (Zwally and others, 2015), which were 17.5 mm a⁻¹ for 1992–2001 and 18.6 mm a⁻¹ for 2003–08 over EA, and 55.0 mm a⁻¹ for 1992–2001 and 62.1 mm a⁻¹ for 2003–08 over WA, after respective adjustments of δB _adj−Iv = −2.7 mm a⁻¹ for EA and −6.2 mm a⁻¹ for WA (Table 4).

Ice dynamic changes are driven by long-term changes in accumulation, but those dynamic changes remain small for long periods of time (e.g. 100–10 000 a) as changes in accumulation slowly change the ice thickness, which in turn slowly changes the gravitational forcing of the ice velocity. For decadal and sub-decadal changes that are driven by atmospheric variations in accumulation, the corresponding dynamic response is very small, including for the relatively large +25% and −23% of SMB variations in EA discussed above. Therefore, our conclusions of dynamic ice thickening in both EA and WA2 during 1992–2016 are based on the absence of persistent accumulation-driven mass changes during that time period.

However, over much longer times (e.g. >1000 a), a sustained change in accumulation significantly alters the ice velocity so that any conclusions on long-term dynamic thickening (or thinning) necessarily depend on other evidence of long-term changes. An example is the marked increase in accumulation that began in the early Holocene (~10 ka BP), with a 67–266% increase after the LGM as derived from the Vostok ice core and radar-layer linking to four other locations (Siegert, 2003). During the preceding 100 ka, accumulation had diminished from ~15–12 mm w.e. a⁻¹ as the climate cooled during the Glacial Period giving the EA ice sheet a long time to reach a dynamic balance with the low accumulation rate. During the Holocene as shown in Figure 2 of Siegert (2003), the accumulation at Vostok Station was sustained at a higher level of 22 mm w.e. a⁻¹ after rising from 12 mm w.e. a⁻¹ following the LGM. That sustained larger accumulation rate and consequent slow acceleration of the ice flow was the basis for the conclusion of long-term dynamic thickening in EA made in Zwally and others, 2015. As noted in the introduction, ice growth in EA also indicated Holocene glacier advances from the EA ice sheet through the Transantarctic mountains into the Dry Valleys (Stuiver and others, 1981; Denton and Wilson, 1982).

The effect of a long-term sustained increase in accumulation on ice thickening, and consequently the ice velocity in a large ice sheet is shown by the following basic consideration in Zwally and others (2015): a 20 mm a⁻¹ elevation increase in central EA causes only a 200 m elevation increase over 10 000 years, producing a correspondingly small ~6% increase in the driving stress under 3400 m of ice and therefore a very slow acceleration of the ice flow increasing slowly with time. In addition, a 3-D numerical model (Wang and others, 2012) of the dynamic response of the ice flow in central EA (e.g. near the ice divide in DS13 ~105°E) to a doubling of accumulation after the LGM showed the surface elevation of the grounded ice increasing at a nearly constant rate of 20 mm a⁻¹ for 10 ka, reaching a 200 m elevation increase at present, followed by a future decreasing rate of rise continuing asymptotically to a total 320 m elevation increase in another 30 ka (Wang and others, 2013).

Similar to EA, the present accumulation rate in WA at present is around twice that of the ice age rate 6400–16 000 years ago (Siegert and Payne, 2004). However, as noted in the introduction, the long-term dynamic ice history of WA with a major retreat after the LGM followed by a Holocene readvance (Bradley and others, 2015; Kingslake and others, 2018) is very different from the very long-term dynamic stability of EA. For the inland WA2, our finding of −3.48 mm a⁻¹ subsidence (δB _0−avg) is in contrast to the average uplifts ranging from 3.00 to 5.42 mm a⁻¹ from the three GIA models. As discussed in Section 6, the reasons for the difference are: (1) the differing ice-loading history in WA2 associated with the post-2014 findings of a Holocene readvance of the grounding lines of the Ross and Filchner-Ronne Ice Shelves from their maximum post-LGM inland retreats, which is not in the history of the earlier GIA models, and (2) the findings of Barletta and others (2018) of a lower mantle viscosity in WA and the consequent significantly-faster GIA response times compared to those of the higher-viscosity GIA models.

A primary glacial forcing for GIA models in WA is from the loss of ice associated with the post-LGM thinning of the ice sheet, shown in one reconstruction (Bentley and others, 2014) as a decreasing elevation at a location near the ice divide in WA2 from 200 m above the present level at 20 ka BP to the present elevation at a rate of 50 m every 5 ka. That post-LGM ice loss was a principal driver of the Antarctic contribution to global mean sea-level rise that started rising ~15 ka BP and was mostly complete by ~9 ka BP (Pollard and others, 2017).

Specific evidence for a different ice-loading history in WA2 after ~10 ka BP includes a 400 km Holocene readvance of the grounding line of the Ross Shelf from its inland retreat at 9.7 ka BP and a smaller readvance of the Filchner-Ronne Ice Shelf from its retreat at 10.2 ka BP (Kingslake and others, 2018) with associated low post-glacial rebound rates in the Weddell Sea (Bradley and others, 2015). Similar retreats and readvances are also shown (personal communication from David Pollard, 2020) in climate-driven ice-sheet modeling such as Pollard and others (2017). In Section 6, we discussed evidence for ice thickening of several hundred meters over a large area of the lower portions of DS18 and 19 during the mid-to-late Holocene causing an increase in the ice loading of several hundred meters over a rather large area of DS18 and 19 in WA2. We also noted that thickening of the lower portions can also restrain the ice flow and lead to inland thickening as is occurring in DS18 as shown in Figure 9b (and S1b). The 2003–08 mass gain in DS18 is 29 Gt a^–1 (adjusted from the 27.3 Gt a^–1 in Zwally and others, 2015).

As noted in Section 3, the response times can range from decadal to centennial for the lower viscosities found in WA1, the Antarctic Peninsula, Patagonia and elsewhere to the millennial responses for the higher viscosities used in Antarctic GIA models. In Section 6, we estimated that for the millennial response times of the high viscosity models with ~5 response times since the beginning of the Holocene 10 ka BP, the exponentially decaying uplift from post-LGM ice unloading would be reduced by a factor of ~7 × 10⁻³. In contrast, for the lower-viscosity 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.

We note again that our procedures for adjustment of the GIA_cor and dB _cor are based on the simple principle that the respective corrections are caused by the mass and volume changes of the same material in the Earth's mantle underlying the ice sheet. The matching is based on a simple linear relationship between the uncorrected GRACE and ICESat mass changes using a constant determined by the ratio of the mass change to the volume change. Although we find that the values of RatioG/dB from the GIA models give values of ρ _earth that are consistent with the knowledge of mantle densities, that physical correspondence is not essential for making the δB adjustments. However, we believe that the physical relationship implied by the consistency of the ρ _earth values strengthens the validity of our adjustments to the ICESat and GRACE mass estimates.

Finally despite the quality of GIA models, their results are very dependent on model parameters such as mantle viscosity that are estimated using model constraints from limited measurements of sea-level change and crustal motions, which are not measurable in vast ice-covered areas of Antarctica. Furthermore, the GIA models have been highly dependent on ice-sheet models and glacial-geologic evidence for their ice-loading histories that force the mantle flow. We believe our results on Antarctic dynamic thickening and our derived adjustments provide useful information that can be used for further development of the GIA models along with the recent new information on the ice loading history. Also, our attempts to calculate the spatial distribution of RatioG/dB and therefore calculate the spatial distributions of the bedrock motion adjustments for ICESat and GRACE dM/dt matching (rather than regional averages) were limited by the singularities at small values and perhaps the numerical precision of the GIA model results. Therefore, the examination of the RatioG/dB within the models, its spatial distribution and its implications regarding the density of the fluid mantle involved may provide new insights and perhaps methods for avoiding the numerical problems we encountered using the current GIA and dB/dt outputs to calculate their ratio.

Although the inter-decadal changes in total Antarctic accumulation since 1992 have been very small, future increases in accumulation with climate warming are likely to have an increasing impact on the overall Antarctic mass balance. A 200-year reconstruction of Antarctic snow accumulation (Medley and Thomas, 2019) showed that increased accumulation mitigated 20th-century sea-level rise by ~10 mm since 1901 at an average rate of 0.11 mm a⁻¹ (40 Gt a⁻¹) from 1901 to 2000 and a higher rate of 0.25 mm a⁻¹ (88 Gt a⁻¹) from 1979 to 2000, which is consistent with our mass gain of 144 ± 61 Gt a⁻¹ from ERS1/2 during 1992–2001. In that regard, the EA ice sheet is especially important because of its large area contributing 73% of the aforementioned 20th-century mass gain and sea-level rise mitigation. Estimated sensitivities of the total Antarctic mass balance to temperature change range from −0.36 to −0.80 mm a⁻¹ of global sea-level change per °C (equivalent to +130 to +290 Gt a⁻¹ of ice per °C) (Huybrechts, 2004 in Bamber and Payne, 2004). The largest estimate of −0.80 mm a⁻¹ sea-level change per °C includes the interactive effect on accumulation from changes in sea ice extent by 125 km per °C (i.e. distance to open-ocean source of moisture) from Giovinetto and Zwally (1996). A smaller estimate of −0.27 mm a⁻¹ per °C change in SMB is from a regional atmospheric climate model forced by the future climate of a global climate model (Ligtenberg and others, 2013). We also calculate −0.28 mm a⁻¹ per °C change from the temperature and precipitation data for 60° to 90°S as used in Golledge and others (2019) for several RCP climate prediction scenarios. Such accumulation-driven increases, along with the current long-term dynamic thickening in EA and WA2, might continue to offset some increases in dynamic losses such as those that have occurred in the AP and the coastal WA1.

However, the decadal-scale dynamic changes are not all causing increases in mass loss. The M(t) for the AP in Figure 14 shows reduced mass loss for the last several years. Also, as previously noted, the M(t) for WA1 in Figure 13 shows that the marked increase in dynamic loss that began around 2009 reduced some during the later years, possibly related to the solid Earth and sea-level feedbacks modeled by Larour and others (2019). Interestingly, the Kingslake model simulation does not show a post-LGM retreat to inside the present grounding line in the Amundsen Sea sector of WA1, which may have implications regarding the ongoing changes and the possible limited extent of future ice losses in WA1. Also, Barletta and others (2018) note that their finding of a lower mantle viscosity and shortening of the response time to mass changes to ‘decades up to a century … increases the potential stability of the WAIS against catastrophic collapse’, with implications for the stability of the inland WA2 as well .

Table 8. ICESat laser campaign biases determined over leads and polynyas in sea ice

DSL are the ICESat measured D corrected for changes in SSH measured concurrently by Envisat.

Table 9. Accumulation density (ρ _a) and pseudo density (ρ _pseudoI) by region

ρ _a is density associated with δA(t) anomalies.

ρ _pseudoI = dM/dt/dI/dt.

Figures S1 and S2 in alternate multi-color scales to Figures 9 and 16, as previously used in Zwally and others (2015) for example, are included in the Supplementary Material along with a link to the digital data used in those figures.