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

H. Jay Zwally, John W. Robbins, Scott B. Luthcke, Bryant D. Loomis, Frederique Remy

► To cite this version:

H. Jay Zwally, John W. Robbins, Scott B. Luthcke, Bryant D. Loomis, Frederique Remy. Mass balance of the Antarctic ice sheet 1992-2016: reconciling results from GRACE gravimetry with ICESat, ERS1/2 and Envisat altimetry. Journal of Glaciology, 2021, 10.1017/jog.2021.8. hal-04720517

HAL Id: hal-04720517

https://hal.science/hal-04720517v1

Submitted on 7 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

**Cite this article:** Zwally HJ, Robbins JW, Luthcke SB, Loomis BD, Rémy F (2021). Mass balance of the Antarctic ice sheet 1992–2016: reconciling results from GRACE gravimetry with ICESat, ERS1/2 and Envisat altimetry.

Journal of Glaciology 67(263), 533–559. https://doi.org/10.1017/jog.2021.8

Received: 31 December 2019

Revised: 8 December 2020

Accepted: 23 December 2020

First published online: 29 March 2021

**Key words:**

Antarctic glaciology; climate change; ice dynamics; ice-sheet mass balance; laser altimetry

**Author for correspondence:**

H. Jay Zwally, E-mail: jayzwallyice@verizon.net

© The Author(s), 2021. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

H. Jay Zwally¹,² ID, John W. Robbins³, Scott B. Luthcke⁴, Bryant D. Loomis⁴ ID and Frédérique Rémy⁵

¹Cryospheric Sciences Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA; ²Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD, USA; ³Craig Technologies, NASA Goddard Space Flight Center, Greenbelt, MD, USA; ⁴Geodesy and Geophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA and ⁵Laboratoire d'Etudes en Geophysique et Oceanographie Spatiale (Legos), Toulouse, France

## 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₉, Sₐ) 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

<table><thead><tr><th>Symbol</th><th>Meaning</th><th>Usual units</th></tr></thead><tbody><tr><td>ρ<sub>x</sub></td><td>Density of material: ρ<sub>ice</sub> (ice below the firn), ρ<sub>earth</sub> (Earth material involved in GIA correction), ρ<sub>a</sub> (new firn from accumulation variations distributed over range of depths)</td><td>dimensionless, relative to water (1 g cm<sup>-3</sup>)</td></tr><tr><td>M(t), dM/dt</td><td>Total mass time series, rate of change</td><td>Gt, Gt a<sup>-1</sup> or mm w.e., mm w.e. a<sup>-1</sup></td></tr><tr><td>M<sub>a</sub>(t), dM<sub>a</sub>/dt</td><td>Component of mass change from variations in accumulation rate: time series, rate of change</td><td>Gt, Gt a<sup>-1</sup> or mm w.e., mm w.e. a<sup>-1</sup></td></tr><tr><td>A(t)</td><td>Accumulation rate</td><td>mm w.e. a<sup>-1</sup></td></tr><tr><td>δA(t)</td><td>Variations in accumulation rate (A(t) - &lt;A(t)&gt;)</td><td>mm w.e. a<sup>-1</sup></td></tr><tr><td>M<sub>d</sub>(t), dM<sub>d</sub>/dt</td><td>Component of mass-change from ice dynamics: time series, rate of change</td><td>Gt, Gt a<sup>-1</sup> or mm w.e. mm w.e. a<sup>-1</sup></td></tr><tr><td>dB/dt</td><td>Bedrock elevation-change rate</td><td>mm a<sup>-1</sup></td></tr><tr><td>H(t), dH/dt</td><td>Time-series, rate of change of surface elevation</td><td>mm, mm a<sup>-1</sup></td></tr><tr><td>GIA<sub>cor</sub></td><td>Glacial Isostatic Adjustment correction</td><td>Gt a<sup>-1</sup> or mm w.e. a<sup>-1</sup></td></tr><tr><td>dB<sub>cor</sub></td><td>Correction to dH/dt for bedrock motion (dB/dt)</td><td>Gt a<sup>-1</sup> or mm w.e. a<sup>-1</sup></td></tr><tr><td>RatioG/dB</td><td>Ratio of GIA<sub>cor</sub> to dB<sub>cor</sub></td><td>dimensionless</td></tr><tr><td>(S<sub>g</sub>)<sub>md</sub></td><td>Sensitivity of gravimetry to bedrock motion for model md (Iv=Ivins, Pe=Peltier, Wh=Whitehouse)</td><td>Gt mm<sup>-1</sup></td></tr><tr><td>S<sub>a</sub></td><td>Sensitivity of altimetry to bedrock motion</td><td>Gt mm<sup>-1</sup></td></tr><tr><td>δB<sub>0</sub>-md</td><td>Rate of bedrock uplift or subsidence needed to provide the GIA<sub>cor</sub> and dB<sub>cor</sub> needed to bring gravimetry and altimetry measured dM/dt (with no GIA<sub>cor</sub> nor dB<sub>cor</sub> applied) into agreement using (S<sub>g</sub>)<sub>md</sub> and S<sub>a</sub></td><td>mm a<sup>-1</sup></td></tr><tr><td>δB<sub>adj</sub>-md</td><td>Same as above, but using dM/dt with modeled GIA<sub>cor</sub> and dB<sub>cor</sub> already applied</td><td>mm a<sup>-1</sup></td></tr><tr><td>δB'</td><td>Estimated long-term effect on the rate of bedrock motion caused by a long-term dynamic thickening/thinning</td><td>mm a<sup>-1</sup></td></tr><tr><td>h(x<sub>i</sub>, y<sub>r</sub>, t<sub>i</sub>)</td><td>ICESat-based surface elevation at cross-track positions x<sub>i</sub> and y<sub>r</sub> (172 m spacings) on the reference track at time t<sub>i</sub> derived from repeat-track analysis</td><td>length units</td></tr><tr><td>α<sub>d</sub></td><td>Cross-track slope</td><td>dimensionless</td></tr><tr><td>h<sub>r</sub>(t<sub>i</sub>)</td><td>Surface-elevation time series at reference point</td><td>length units</td></tr><tr><td>H<sub>j,k</sub>(t<sub>i</sub>)</td><td>ICESat surface-elevation time series for grid cell (j,k)</td><td>length units</td></tr><tr><td>H(t)</td><td>Time series (t<sub>i</sub>) of surface elevation from area integration of h<sub>r</sub>(t<sub>i</sub>) over 50 km grid cells, drainage systems, or regions</td><td>length units</td></tr></tbody></table>

<figure><img src="image_3.png" /></figure>
<table><tr><td>H<sub>d</sub>(t)</td><td>Time series of surface elevation driven by ice dynamics</td><td>length units</td></tr><tr><td>H<sup>a</sup>(t)</td><td>Surface elevation changes driven directly by contemporary accumulation variations</td><td>length units</td></tr><tr><td>C<sub>A</sub>(t)</td><td>Surface-elevation change from changes in firn-compaction rate driven by variations in accumulation rate</td><td>length units</td></tr><tr><td>C<sub>T</sub>(t)</td><td>Surface-elevation change from changes in firn-compaction rate driven by variations in near-surface firn temperature</td><td>length units</td></tr><tr><td>C<sub>AT</sub>(t)</td><td>Surface-elevation change from changes in firn-compaction rate driven by the combined effect of variations in accumulation and firn temperature</td><td>length units</td></tr><tr><td>H<sup>a</sup><sub>CA</sub>(t)</td><td>Combination of elevation changes directly by H<sup>a</sup>(t) accumulation variations and accumulation-driven changes in firn compaction rate (H<sup>a</sup><sub>CA</sub>(t) = H<sup>a</sup>(t) + C<sub>A</sub>(t))</td><td>length units</td></tr><tr><td>(dM/dt)<sub>eq-md</sub></td><td>Mass change rate after bringing gravimetry and altimetry estimates into agreement using (S<sub>g</sub>)<sub>md</sub> and S<sub>a</sub> with either δB<sub>0-md</sub> or δB<sub>adj-md</sub></td><td>Gt a<sup>-1</sup> or mm w.e. a<sup>-1</sup></td></tr></table>

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

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 glacial 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<sub>cor</sub> and dB<sub>cor</sub>) 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.

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
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

for GIAcor and the vertical motion of the bedrock for calculation of dBcor. We derive the respective sensitivities (Sg and Sa) of the GIAcor and dBcor to the bedrock motion from the models. We then use Sg and Sa to derive the bedrock motion (i.e. δB0) needed to match the mass changes during 2003–2008 from GRACE and ICESat without any GIAcor and dBcor applied, and alternatively the adjustments (i.e. δBadj) to the modeled bedrock motions with model-based GIAcor and dBcor applied. We then use the adjusted GIAcor 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 GIAcor and the ICESat dBcor, 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) subsurface 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 δB0 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 dBcor corrections to the ERS1/2 results for 1992–2001 as well as the ICESat results for 2003–2009, and the corresponding GIAcor 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 (GIAcor) 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

$$GIA_{cor} = \Delta M / \Delta t = \rho_{earth} \Delta V / \Delta t,$$

where ρearth is the relative density of mantle material involved in the ΔM/Δt change. For altimetry, the correction (dBcor) to the mass changes calculated from changes in ice-sheet surface
<figure><img src="image_5.png" /><figcaption>**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 \approx \rho_{ice} / \rho_{mantle} \times T$, which is 600 m for T = 3000 m and $\rho_{ice} = 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).</figcaption></figure>

elevation (dH/dt) is made for the vertical motion of the bedrock (dB/dt) caused by the ΔV/Δt. The dBcor is equal to $\rho_{ice} \Delta V/\Delta t$, where $\rho_{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 $\rho_{ice}$ is appropriate, because basal motion displaces solid ice and does not affect the density nor volume of the firn column. GIAcor and dBcor are used here as rates of mass change per unit area, so using $\Delta V/\Delta t = dB/dt \times$ area, the GIAcor and dBcor per unit area are:

$$GIA_{cor} = \rho_{earth} \, dB/dt \qquad (1a)$$

$$dB_{cor} = \rho_{ice} \, dB/dt \qquad (1b)$$

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

Although GIAcor and dBcor 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 GIAcor 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.

We define RatioG/dB as

$$\text{RatioG/dB} \equiv \text{GIA}_{\text{cor}}/\text{dB}_{\text{cor}} = \rho_{\text{earth}}/\rho_{\text{ice}} = \rho_{\text{earth}}/0.91, \qquad (2)$$

the gravimetry sensitivity ($S_g$)$_{md}$ to bedrock motion as

$$(S_g)_{md} \equiv -\text{GIA}_{\text{cor}}/(\text{dB/dt}), \qquad (3)$$

where the subscript (md) indicates the GIA model used (Iv, Pe

or Wh), and the altimetry sensitivity ($S_a$) to bedrock motion as

$$S_a \equiv -\frac{dB_{cor}}{(dB/dt)}. \qquad (4)$$

The units of dB/dt are mm a⁻¹, the units of GIAcor and dBcor are Gt a⁻¹, the units of Sg and Sa are both Gt a⁻¹/mm a⁻¹, and

$$\text{RatioG/dB} = S_g/S_a. \qquad (5)$$

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 $\rho_{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 ($\delta B_{0-\text{md}}$) needed to provide the GIAcor and dBcor corrections that bring the respective mass estimates into full agreement (i.e. $[(dM/dt)_{\text{GRACE}}]_{\text{eq}} = [(dM/dt)_{\text{ICESat}}]_{\text{eq}}$). The required uplift or subsidence, $\delta B_{0-\text{md}}$, relative to zero is given by:

$$[(dM/dt)_{\text{GRACE}}]_0 + (S_g)_{\text{md}} \delta B_{0-\text{md}} = [(dM/dt)_{\text{ICESat}}]_0 + S_a \delta B_{0-\text{md}} \\ \delta B_{0-\text{md}} = \frac{[(dM/dt)_{\text{GRACE}}]_0 - [(dM/dt)_{\text{ICESat}}]_0}{S_a - (S_g)_{\text{md}}} \qquad (6)$$

where $[(dM/dt)_{\text{GRACE}}]_0$ and $[(dM/dt)_{\text{ICESat}}]_0$ are the respective GRACE and ICESat measurements with zero GIAcor and zero dBcor 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)_{\text{Iv}}$, $(S_g)_{\text{Wh}}$ or $(S_g)_{\text{Pe}}$ for Ivins, Whitehouse or Peltier model. For example, $\delta B_{0-\text{Iv}}$ indicates that $(S_g)_{\text{Iv}}$ derived from the Ivins model of GIA and dB/dt was used with no dBcor nor GIAcor applied to the measured dM/dt.

The required uplift or subsidence ($\delta B_{\text{adj}-\text{md}}$) is also calculated relative to the GIA modeled

$$\delta B_{\text{adj}-\text{md}} = \frac{[(dM/dt)_{\text{GRACE}}]_{\text{md}} - [(dM/dt)_{\text{ICESat}}]_{\text{md}}}{S_a - (S_g)_{\text{md}}}, \qquad (7)$$
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

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

$$ (dM/dt)_{\text{eq-md}} = [(dM/dt)_{\text{GRACE}}]_{\text{eq-md}} \\ = [(dM/dt)_{\text{ICESat}}]_{\text{eq-md}} \qquad (8) $$

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

$$ \delta B' = -(dH_d/dt)_{\text{obs}}/\text{Ratio G/dB.} \qquad (9) $$

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
<table><caption>**Table 1.** Glacial Isostatic Adjustment (GIA) and uplift ($dB/dt$) from Ivins, Whitehouse and Peltier Earth models</caption><thead><tr><th rowspan="2">Model/region</th><th colspan="2">GIA<sub>cor</sub></th><th rowspan="2">S<sub>g</sub><br>(Gt a<sup>-1</sup>)</th><th rowspan="2">dB/dt<br>(mm a<sup>-1</sup>)</th><th rowspan="2">dB<sub>cor</sub><br>(Gt a<sup>-1</sup>)</th><th rowspan="2">S<sub>a</sub><br>(Gt a<sup>-1</sup>)</th><th rowspan="2">RatioG/dB = S<sub>g</sub>/S<sub>a</sub></th><th rowspan="2">ρ<sub>earth</sub></th></tr><tr><th>(mm w.e. a<sup>-1</sup>)</th><th>(Gt a<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="9"><strong>Ivins/</strong></td></tr><tr><td>EA total</td><td>1.9</td><td>19.9</td><td>-47.1</td><td>0.422</td><td>3.9</td><td>-9.29</td><td>5.07</td><td>4.61</td></tr><tr><td>EA1 (2-11)</td><td>2.3</td><td>12.4</td><td>-24.4</td><td>0.508</td><td>2.5</td><td>-4.93</td><td>4.95</td><td>4.50</td></tr><tr><td>EA2 (12-17)</td><td>1.6</td><td>7.5</td><td>-23.0</td><td>0.325</td><td>1.4</td><td>-4.36</td><td>5.28</td><td>4.80</td></tr><tr><td>WA total</td><td>11.1</td><td>20.9</td><td>-7.6</td><td>2.758</td><td>4.7</td><td>-1.71</td><td>4.44</td><td>4.04</td></tr><tr><td>WA 1</td><td>9.3</td><td>5.9</td><td>-2.6</td><td>2.278</td><td>1.3</td><td>-0.58</td><td>4.47</td><td>4.07</td></tr><tr><td>WA 2</td><td>12.1</td><td>15.1</td><td>-5.0</td><td>3.002</td><td>3.4</td><td>-1.13</td><td>4.43</td><td>4.03</td></tr><tr><td>AP</td><td>8.0</td><td>2.3</td><td>-1.2</td><td>1.983</td><td>0.5</td><td>-0.27</td><td>4.43</td><td>4.03</td></tr><tr><td>AIS (all)</td><td>3.5</td><td>43.1</td><td>-53.0</td><td>0.814</td><td>9.2</td><td>-11.26</td><td>4.71</td><td>4.29</td></tr><tr><td colspan="9"><strong>Peltier/</strong></td></tr><tr><td>EA total</td><td>3.1</td><td>31.5</td><td>-52.6</td><td>0.599</td><td>5.6</td><td>-9.29</td><td>5.67</td><td>5.16</td></tr><tr><td>EA1 (2-11)</td><td>3.5</td><td>19.0</td><td>-29.1</td><td>0.655</td><td>3.2</td><td>-4.93</td><td>5.90</td><td>5.37</td></tr><tr><td>EA2 (12-17)</td><td>2.6</td><td>12.5</td><td>-23.7</td><td>0.527</td><td>2.3</td><td>-4.36</td><td>5.44</td><td>4.95</td></tr><tr><td>WA total</td><td>19.0</td><td>35.6</td><td>-7.5</td><td>4.758</td><td>8.1</td><td>-1.71</td><td>4.38</td><td>3.99</td></tr><tr><td>WA 1</td><td>15.3</td><td>9.7</td><td>-2.8</td><td>3.451</td><td>2.0</td><td>-0.58</td><td>4.87</td><td>4.43</td></tr><tr><td>WA 2</td><td>20.8</td><td>25.9</td><td>-4.8</td><td>5.423</td><td>6.1</td><td>-1.13</td><td>4.22</td><td>3.84</td></tr><tr><td>AP</td><td>16.4</td><td>4.8</td><td>-1.2</td><td>4.030</td><td>1.1</td><td>-0.27</td><td>4.48</td><td>4.08</td></tr><tr><td>AIS (all)</td><td>5.8</td><td>72.0</td><td>-55.0</td><td>1.308</td><td>14.7</td><td>-11.26</td><td>4.88</td><td>4.44</td></tr><tr><td colspan="9"><strong>Whitehouse/</strong></td></tr><tr><td>EA total<sup>a</sup></td><td>-0.9</td><td>-8.8</td><td>-45.7</td><td>-0.192</td><td>-1.8</td><td>-9.29</td><td>4.92</td><td>4.48</td></tr><tr><td>EA1 (2-11)<sup>a</sup></td><td>-2.2</td><td>-11.7</td><td>-24.0</td><td>-0.487</td><td>-2.4</td><td>-4.93</td><td>4.86</td><td>4.43</td></tr><tr><td>EA2 (12-17)<sup>a</sup></td><td>0.6</td><td>3.1</td><td>-21.7</td><td>0.143</td><td>0.6</td><td>-4.36</td><td>4.99</td><td>4.54</td></tr><tr><td>WA total</td><td>17.5</td><td>32.9</td><td>-8.1</td><td>4.074</td><td>7.0</td><td>-1.71</td><td>4.72</td><td>4.30</td></tr><tr><td>WA 1</td><td>14.6</td><td>9.2</td><td>-2.9</td><td>3.231</td><td>1.9</td><td>-0.58</td><td>4.95</td><td>4.50</td></tr><tr><td>WA 2</td><td>19.0</td><td>23.7</td><td>-5.2</td><td>4.504</td><td>5.1</td><td>-1.13</td><td>4.64</td><td>4.22</td></tr><tr><td>AP</td><td>6.5</td><td>1.9</td><td>-2.1</td><td>0.909</td><td>0.2</td><td>-0.27</td><td>7.90</td><td>7.19</td></tr><tr><td>AIS (all)</td><td>2.1</td><td>26.9</td><td>-55.9</td><td>0.482</td><td>5.4</td><td>-11.26</td><td>4.79</td><td>4.36</td></tr></tbody></table>

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

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 GIAcor (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 GIAcor (Gt a⁻¹) and dBcor (Gt a⁻¹) mass corrections. [Note: total regional values are calculated as GIAcor (Gt a⁻¹) = GIAcorr (mm w.e. a⁻¹) area (km²) 10⁻⁶ and dBcor (Gt a⁻¹) = 0.91 dB/dt (mm a⁻¹) area (km²) 10⁻⁶.] The GIAcorr and dBcor 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 GIAcor and dBcor 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 GIAcor is 1.9 mm w.e. a⁻¹ uplift and the regional dM/dt adjustment is -19.9 Gt a⁻¹. The GIAcor 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 GIAcorr is -0.9 mm w.e. a⁻¹ with a regional dM/dt adjustment of +8.8 Gt a⁻¹. Differences among the modeled GIAcorr are as large as 40 Gt a⁻¹ between the Peltier and Whitehouse models for EA (mostly EA1), 15 Gt a⁻¹

<footnote>ªSee text about the calculation of Sₗ and RatioG/dB.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Components of mass change from ICESat for EA, EA1 and EA2 from $M_d(t) = \rho_{ice} \times H_d(t)$ from Figure 5 and $M_a(t) = \int \delta A(t) \times dt$ 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)$.</figcaption></figure>

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 GIAcor, 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.

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

$$<\text{RatioG}/\text{dB}>_{\text{regavg}} = <\text{GIA}>_{\text{regavg}} / 0.91 <\text{dB}/\text{dt}>_{\text{regavg}}. \quad (10)$$

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 $<\text{GIA}>_{\text{DSavg}}$ and $<\text{dB}/\text{dt}>_{\text{DSavg}}$ by drainage system (DS) to calculate the $<\text{RatioG}/\text{dB}>_{\text{regavg}}$, excluding DS16 and DS17 from the EA2 calculation. The regional average $S_g$ are calculated using $<\text{RatioG}/\text{dB}>_{\text{regavg}}$ $S_a$ and $<\text{GIA}>_{\text{regavg}} = <\text{dB}/\text{dt}>_{\text{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 $\delta B_{0-md}$ in Eqn (6) and the $\delta B_{adj-md}$ in Eqn (7)) and the GIAcor 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 GIAcor 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 $\rho_{earth}$ derived from the three GIA models by region according to Eqn (2), with values mostly in the range of 4–5. These $\rho_{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 $\rho_{earth}$ densities are not used for calculation of the $\delta 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, $\alpha_r$, as well as on real elevation variations with time according to

$$h(x_i, y_r, t_i) = x_i \tan \alpha_r + t_i(\text{dh}/\text{dt})_r + h_0(y_r, t_0), \quad (11)$$

where $h_0$ is the elevation at the position $y_r$ on the reference track at
<table><caption>**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)</caption><thead><tr><th></th><th colspan="8">2003–2008</th></tr><tr><th></th><th>EA</th><th>EA1</th><th>EA2</th><th>WA</th><th>WA1</th><th>WA2</th><th>AP</th><th>AIS</th></tr></thead><tbody><tr><td>(mm a<sup>-1</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>dH/dt</strong></td><td>12.85</td><td>22.66</td><td>1.85</td><td>-7.08</td><td>-146.70</td><td>65.30</td><td>-55.50</td><td>8.38</td></tr><tr><td>(dH/dt)<sub>2015</sub></td><td>13.00</td><td>23.10</td><td>1.50</td><td>-5.10</td><td>-150.80</td><td>69.20</td><td>-97.50</td><td>7.60</td></tr><tr><td><strong>dH<sub>a</sub>/dt</strong></td><td>13.69</td><td>10.45</td><td>17.39</td><td>-39.93</td><td>-196.10</td><td>40.96</td><td>-47.93</td><td>4.38</td></tr><tr><td>(dH<sub>a</sub>/dt)<sub>2015</sub></td><td>15.90</td><td>14.70</td><td>17.30</td><td>-24.80</td><td>-183.30</td><td>55.90</td><td>-102.00</td><td>6.90</td></tr><tr><td><strong>dH<sub>a</sub>/dt</strong></td><td>-0.11</td><td>11.03</td><td>-12.69</td><td>55.40</td><td>80.10</td><td>42.67</td><td>-2.39</td><td>8.10</td></tr><tr><td><strong>dC<sub>AT</sub>/dt</strong></td><td>-1.14</td><td>0.65</td><td>-3.13</td><td>-25.50</td><td>-32.70</td><td>-21.73</td><td>-8.44</td><td>-4.91</td></tr><tr><td><strong>dH<sub>a</sub>C<sub>AT</sub>/dt</strong></td><td>-1.25</td><td>11.68</td><td>-15.82</td><td>29.90</td><td>47.40</td><td>20.94</td><td>-10.83</td><td>3.19</td></tr><tr><td>(dH<sup>o</sup>C<sub>AT</sub>/dt)<sub>2015</sub></td><td>-3.50</td><td>8.40</td><td>-16.20</td><td>17.30</td><td>30.80</td><td>10.50</td><td>2.90</td><td>0.00</td></tr><tr><td>(dB/dt)<sub>2015</sub></td><td>0.42</td><td>0.49</td><td>0.33</td><td>2.57</td><td>2.04</td><td>2.84</td><td>1.70</td><td>0.77</td></tr><tr><td>(Gt a<sup>-1</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>dM/dt</strong></td><td>125.5</td><td>69.6</td><td>55.9</td><td>-34.3</td><td>-96.4</td><td>62.0</td><td>-10.3</td><td>80.9</td></tr><tr><td>(dM/dt)<sub>2015</sub></td><td>136.1</td><td>85.5</td><td>50.6</td><td>-25.0</td><td>-96.8</td><td>71.8</td><td>-28.8</td><td>82.4</td></tr><tr><td><strong>dM<sub>a</sub>/dt</strong></td><td>125.4</td><td>51.4</td><td>74.0</td><td>-65.2</td><td>-111.6</td><td>46.4</td><td>-10.4</td><td>49.8</td></tr><tr><td>(dM<sub>a</sub>/dt)<sub>2015</sub></td><td>147.4</td><td>72.2</td><td>75.1</td><td>-42.4</td><td>-105.7</td><td>63.4</td><td>-27.1</td><td>77.9</td></tr><tr><td><strong>dM<sub>a</sub>/dt</strong></td><td>0.1</td><td>18.2</td><td>-18.1</td><td>30.9</td><td>15.2</td><td>15.6</td><td>0.1</td><td>31.1</td></tr><tr><td>(dM<sub>a</sub>/dt)<sub>2015</sub></td><td>-11.2</td><td>13.3</td><td>-24.5</td><td>17.4</td><td>8.9</td><td>8.5</td><td>-2.8</td><td>4.5</td></tr><tr><td>(dB<sub>cor</sub>)<sub>2015</sub></td><td>3.9</td><td>2.4</td><td>1.4</td><td>4.4</td><td>1.2</td><td>3.2</td><td>0.5</td><td>8.7</td></tr></tbody></table>

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

$t_0$. 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 $\alpha_r$, $(dh/dt)_r$ and $h_0$ 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_0(t_0)$, $h_1(t_1)$, $h_2(t_2)$, ... $h_{16}(t_{16})$ is created for each reference point using the cross-track slope $\alpha_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

<table><caption>**Table 3.** ICESat elevation and mass change components for 2003–2008 and 2003–2009 from time series analysis using dB/dt equal $\delta 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</caption><thead><tr><th></th><th colspan="6">2003–2008</th></tr><tr><th></th><th>EA<br>(mm a<sup>-1</sup>)</th><th>EA1<br>(mm a<sup>-1</sup>)</th><th>EA2<br>(mm a<sup>-1</sup>)</th><th>WA<br>(mm a<sup>-1</sup>)</th><th>WA1<br>(mm a<sup>-1</sup>)</th><th>WA2<br>(mm a<sup>-1</sup>)</th></tr></thead><tbody><tr><td>dH/dt</td><td>12.85</td><td>22.66</td><td>1.85</td><td>-7.08</td><td>-146.72</td><td>65.30</td></tr><tr><td>dH<sub>a</sub>/dt</td><td>16.38</td><td>12.37</td><td>20.98</td><td>-34.67</td><td>-193.62</td><td>47.66</td></tr><tr><td>dH<sub>a</sub>/dt</td><td>-0.11</td><td>11.03</td><td>-12.69</td><td>55.40</td><td>80.07</td><td>42.67</td></tr><tr><td>dC<sub>AT</sub>/dt</td><td>-1.14</td><td>0.65</td><td>-3.13</td><td>-25.50</td><td>-32.66</td><td>-21.73</td></tr><tr><td>dB/dt</td><td>-2.28</td><td>-1.32</td><td>-3.25</td><td>-3.40</td><td>-0.38</td><td>-3.47</td></tr><tr><td></td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td></tr><tr><td>dM<sub>a</sub>/dt</td><td>150.0</td><td>60.5</td><td>89.5</td><td>-57.0</td><td>-110.2</td><td>53.2</td></tr><tr><td>dM<sub>a</sub>/dt</td><td>0.1</td><td>18.2</td><td>-18.1</td><td>30.9</td><td>15.2</td><td>15.6</td></tr><tr><td>dM/dt</td><td>150.1</td><td>78.7</td><td>71.4</td><td>-26.1</td><td>-94.9</td><td>68.8</td></tr><tr><td>dB<sub>cor</sub></td><td>-20.6</td><td>-6.5</td><td>-14.1</td><td>-4.1</td><td>-0.2</td><td>-3.9</td></tr><tr><td>GlA<sub>cor</sub></td><td>-106.8</td><td>-32.2</td><td>-74.6</td><td>-18.4</td><td>-1.0</td><td>-17.4</td></tr></tbody><thead><tr><th></th><th colspan="6">2003–2009</th></tr><tr><th></th><th>(mm a<sup>-1</sup>)</th><th>(mm a<sup>-1</sup>)</th><th>(mm a<sup>-1</sup>)</th><th>(mm a<sup>-1</sup>)</th><th>(mm a<sup>-1</sup>)</th><th>(mm a<sup>-1</sup>)</th></tr></thead><tbody><tr><td>dH/dt</td><td>18.15</td><td>32.12</td><td>2.36</td><td>-25.36</td><td>-182.42</td><td>55.98</td></tr><tr><td>dH<sub>a</sub>/dt</td><td>21.35</td><td>20.38</td><td>22.26</td><td>-33.08</td><td>-201.43</td><td>54.03</td></tr><tr><td>dH<sub>a</sub>/dt</td><td>-0.53</td><td>13.54</td><td>-16.26</td><td>32.46</td><td>47.10</td><td>24.84</td></tr><tr><td>dC<sub>AT</sub>/dt</td><td>-0.48</td><td>-0.52</td><td>-0.34</td><td>-22.45</td><td>-26.60</td><td>-19.79</td></tr><tr><td></td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td><td>(Gt a<sup>-1</sup>)</td></tr><tr><td>dM<sub>a</sub>/dt</td><td>195.2</td><td>99.7</td><td>95.5</td><td>-54.2</td><td>-114.5</td><td>60.3</td></tr><tr><td>dM<sub>a</sub>/dt</td><td>-0.9</td><td>22.2</td><td>-23.1</td><td>18.0</td><td>9.0</td><td>9.1</td></tr><tr><td>dM/dt</td><td>194.3</td><td>121.9</td><td>72.4</td><td>-36.2</td><td>-105.6</td><td>69.4</td></tr></tbody></table>

Slopes are linear term at midpoint of time period from LQS fitting (2006.0 for 2003–2008 and 2006.5 for 2003–2009).
<table><caption>**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}$</caption><thead><tr><th>Region</th><th>Model</th><th>δB<sub>0-md</sub><br>(mm a<sup>-1</sup>)</th><th>δB<sub>adj-md</sub><br>(mm a<sup>-1</sup>)</th><th>[(dM/dt)<sub>eq</sub>]md<br>(Gt a<sup>-1</sup>)</th><th>(dM/dt)*<br>(Gt a<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="4">EA</td><td>Ivins</td><td>-2.28</td><td>-2.70</td><td>150.5</td><td rowspan="4">150.1</td></tr><tr><td>Peltier</td><td>-1.99</td><td>-2.59</td><td>147.8</td></tr><tr><td>Whitehouse</td><td>-2.36</td><td>-2.17</td><td>151.3</td></tr><tr><td>Average</td><td>-2.21</td><td>-2.49</td><td>149.9</td></tr><tr><td rowspan="4">EA1</td><td>Ivins</td><td>-1.32</td><td>-1.83</td><td>78.5</td><td rowspan="4">78.7</td></tr><tr><td>Peltier</td><td>-1.06</td><td>-1.72</td><td>77.3</td></tr><tr><td>Whitehouse</td><td>-1.35</td><td>-0.86</td><td>78.7</td></tr><tr><td>Average</td><td>-1.24</td><td>-1.47</td><td>78.2</td></tr><tr><td rowspan="4">EA2</td><td>Ivins</td><td>-3.25</td><td>-3.57</td><td>71.5</td><td rowspan="4">71.4</td></tr><tr><td>Peltier</td><td>-3.13</td><td>-3.66</td><td>71.0</td></tr><tr><td>Whitehouse</td><td>-3.48</td><td>-3.62</td><td>72.5</td></tr><tr><td>Average</td><td>-3.29</td><td>-3.62</td><td>71.7</td></tr><tr><td rowspan="4">WA</td><td>Ivins</td><td>-3.40</td><td>-6.16</td><td>-24.1</td><td rowspan="4">-26.1</td></tr><tr><td>Peltier</td><td>-3.47</td><td>-8.22</td><td>-24.0</td></tr><tr><td>Whitehouse</td><td>-3.15</td><td>-7.22</td><td>-24.6</td></tr><tr><td>Average</td><td>-3.34</td><td>-7.20</td><td>-24.2</td></tr><tr><td rowspan="4">WA1</td><td>Ivins</td><td>-0.38</td><td>-2.66</td><td>-95.0</td><td rowspan="4">-94.9</td></tr><tr><td>Peltier</td><td>-0.34</td><td>-3.79</td><td>-95.0</td></tr><tr><td>Whitehouse</td><td>-0.33</td><td>-3.56</td><td>-95.0</td></tr><tr><td>Average</td><td>-0.35</td><td>-3.34</td><td>-95.0</td></tr><tr><td rowspan="4">WA2</td><td>Ivins</td><td>-3.47</td><td>-6.47</td><td>66.0</td><td rowspan="4">68.8</td></tr><tr><td>Peltier</td><td>-3.70</td><td>-9.12</td><td>66.2</td></tr><tr><td>Whitehouse</td><td>-3.27</td><td>-7.78</td><td>65.7</td></tr><tr><td>Average</td><td>-3.48</td><td>-7.79</td><td>66.0</td></tr></tbody></table>

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_s$ 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}$.

$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:

$$ 
\begin{aligned}
H(t) &= H_d(t) + H^a(t) + C_A(t) + C_T(t) + (\text{dB/dt}) t \\
H_d(t) &= H(t) - H^a(t) - C_A(t) - C_T(t) - (\text{dB/dt}) t,
\end{aligned} 
\quad (12) $$

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) \equiv 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

$$ M_a(t) = \int^t \delta A(t) \, dt \quad (13) $$

and the cumulative accumulation-driven height change

$$ H^a(t) = 1/\rho_s \int^t \delta A(t) \, dt, \quad (14) $$

where $ρ_s = 0.3$ is the density of new surface firn. Separation of the

<figure><img src="image_10.png" /><figcaption>**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) - (\text{dB/dt}) \times 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.</figcaption></figure>

$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

$$ M_d(t) = \rho_{\text{ice}} H_d(t) \quad (15) $$

using the well-defined $ρ_{\text{ice}} = 0.91$ and the total mass change is

$$ M(t) = M_d(t) + M_a(t). \quad (16) $$

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 = \Delta M_a / (\Delta H^a - C_A)$ according to Eqn (7) in Zwally
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Components of mass change from ICESat for WA, WA1 and WA2 from $M_d(t) = \rho_{ice} \times H_d(t)$ from Figure 7 and $M_a(t) = \int^1 \delta A(t) \times dt$ with LQS fit through 2008 data only. Linear trends and the adjusted $\text{dB}_{cor}$ and $\text{GIA}_{cor}$ applied are in Table 3. The dynamic $M_d(t)$ is more linear than the total $M(t)$.</figcaption></figure>

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 $\rho_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 $\rho_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 $\text{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^2 + D \ \sin(\omega \ t) + E \ \cos(\omega \ t)$] with an annual period representing a seasonal cycle with the phase and amplitude selected by the fit.

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Maps of $\text{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 $\text{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.)</figcaption></figure>

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 the period 2009 through 2015).
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**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 $\delta B_0 = -1.99$ mm a⁻¹ from Peltier model, to 151 Gt a⁻¹ (●) using $S_g = -47.1$ Gt a⁻¹/mm a⁻¹ and $\delta B_0 = -2.28$ mm a⁻¹ from Ivins model, to 151 Gt a⁻¹ (●) using $S_g = -45.7$ Gt a⁻¹/mm a⁻¹ and $\delta B_0 = -2.36$ mm a⁻¹ from Whitehouse model.</figcaption></figure>

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$.

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$)$_{2015}$ 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$.

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.

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 ERS1/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).

In comparison to Davis’s (1997) 18 mm a⁻¹, our EA elevation changes for our calculated ERS coverage of 8.13 × 10⁶ km² are
<table><caption>**Table 5.** Bedrock motions $\delta B_{0-avg}$ and $\delta 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, $\delta(dB/dt)_{max}$, among models</caption><thead><tr><th rowspan="2">Region</th><th colspan="2">δB<sub>0-avg</sub> and dB<sub>cor</sub></th><th colspan="2">δB<sub>md-avg</sub> and dB<sub>cor</sub></th><th colspan="3">dB/dt (mm a<sup>-1</sup>)</th><th rowspan="2">δ(dB/dt)<sub>max</sub></th></tr><tr><th>(mm a<sup>-1</sup>)</th><th>(Gt a<sup>-1</sup>)</th><th>(mm a<sup>-1</sup>)</th><th>(Gt a<sup>-1</sup>)</th><th>Ivins</th><th>Peltier</th><th>Wthse</th></tr></thead><tbody><tr><td>EA</td><td>-2.21</td><td>-20.5</td><td>-2.49</td><td>-23.1</td><td>0.42</td><td>0.60</td><td>-0.19</td><td>0.79</td></tr><tr><td>EA1</td><td>-1.24</td><td>-6.1</td><td>-0.86</td><td>-4.2</td><td>0.51</td><td>0.66</td><td>-0.49</td><td>1.14</td></tr><tr><td>EA2</td><td>-3.29</td><td>-14.3</td><td>-3.62</td><td>-15.8</td><td>0.33</td><td>0.53</td><td>0.14</td><td>0.38</td></tr><tr><td>WA</td><td>-3.34</td><td>-5.7</td><td>-7.20</td><td>-12.3</td><td>2.76</td><td>4.76</td><td>4.07</td><td>2.00</td></tr><tr><td>WA1</td><td>-0.35</td><td>-0.2</td><td>-3.34</td><td>-1.9</td><td>2.28</td><td>3.45</td><td>3.23</td><td>1.17</td></tr><tr><td>WA2</td><td>-3.48</td><td>-3.9</td><td>-7.79</td><td>-8.8</td><td>3.00</td><td>5.42</td><td>4.50</td><td>2.42</td></tr></tbody></table>

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 mm 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)_{2015}$ of 147 Gt a⁻¹, $(dM_d/dt)_{2015}$ of -11 Gt a⁻¹ and $(dM_d/dt)_{2015}$ 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_{CA}^a/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_d/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<sub>cor</sub> 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

<table><caption>**Table 6.** Estimated bedrock motion, $\delta B'$, caused by the observed dynamic thickening</caption><thead><tr><th>Region</th><th>(dH<sub>d</sub>/dt)<sub>obs</sub><br>(mm a<sup>-1</sup>)</th><th>RatioG/dB</th><th>δB'<br>(mm a<sup>-1</sup>)</th><th>δB<sub>0-avg</sub><br>(mm a<sup>-1</sup>)</th><th>δB''/<br>δB<sub>0-avg</sub></th><th>δB'<br>δB<sub>adj-lv</sub></th></tr></thead><tbody><tr><td>EA</td><td>16.38</td><td>5.07</td><td>-3.23</td><td>-2.21</td><td>1.46</td><td>1.20</td></tr><tr><td>EA1</td><td>12.37</td><td>4.95</td><td>-2.50</td><td>-1.24</td><td>2.02</td><td>1.37</td></tr><tr><td>EA2</td><td>20.98</td><td>5.28</td><td>-3.97</td><td>-3.29</td><td>1.21</td><td>1.11</td></tr><tr><td>WA2</td><td>47.66</td><td>4.43</td><td>-10.76</td><td>-3.48</td><td>3.09</td><td>1.66</td></tr></tbody></table>

The $\delta B'$ equal to -(dH_d/dt)<sub>obs</sub>/RatioG/dB is larger than the bedrock motion (both $\delta B_{0-avg}$ and $\delta B_{adj-lv}$) needed to bring ICESat and GRACE $dM/dt$ into agreement.

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_d/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.

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
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** *M(t)* time series for East Antarctica from ICESat (blue) and GRACE (red) using the equalizing *db*<sub>cor</sub> and GIA<sub>cor</sub> 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.</figcaption></figure>

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

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** *M(t)* time series for West Antarctica from ICESat (blue) and GRACE (red) using the equalizing *db*<sub>cor</sub> and GIA<sub>cor</sub> 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.</figcaption></figure>

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⁻¹.

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
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** *M(t)* time-series for Antarctic Peninsula from ICESat (blue) and GRACE (red) using $dB_{cor} = -0.5 \text{ a}^{-1}$ and $GIA_{cor} = -2.3 \text{ Gt a}^{-1}$ 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 \text{ Gt a}^{-1}$ from LQS is replaced by $-29 \text{ Gt a}^{-1}$ from average-linear change analysis in AIS sum in Figure 15 and Table 7a.*</figcaption></figure>

example, the result from CryoSat data for 2010–13 for EA was only $1 \pm 2 \text{ mm a}^{-1}$, from which they calculated a mass loss of $3 \pm 36 \text{ Gt a}^{-1}$ for an area of 9 499 900 km$^2$ (McMillan and others, 2014); and from ERS, Envisat and CryoSat data for 1992–2017 was $6 \pm 1 \text{ mm a}^{-1}$, from which they calculated a mass gain of only $16.3 \pm 5.5 \text{ Gt a}^{-1}$ for an area of 9 909 800 km$^2$ (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 $\delta B_{0\text{-md}}$ according to Eqn (6), and relative to the modeled $dB/dt$, giving $\delta B_{\text{adj\text{-md}}}$ according to Eqn (7). The resulting $\delta B_{0\text{-md}}$, $\delta B_{\text{adj\text{-md}}}$ and the corresponding $(dM/dt)_{\text{eq\text{-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 \text{ Gt a}^{-1}$ per mm $\text{ a}^{-1}$ for altimetry and $S_{g\text{-lv}} = -47.1$, $S_{g\text{-Pe}} = -52.6$, or $S_{g\text{-Wh}} = -45.7 \text{ Gt a}^{-1}$ per mm $\text{ a}^{-1}$ for gravimetry (Table 1). For EA, the derived uplift adjustments are $\delta B_{\text{adj\text{-lv}}}$ = $-2.70 \text{ mm a}^{-1}$, $\delta B_{\text{adj\text{-Pe}}}$ = $-2.59 \text{ mm a}^{-1}$ and $\delta B_{\text{adj\text{-Wh}}}$ = $-2.17 \text{ mm a}^{-1}$ with an average of $-2.49 \text{ mm}$; and $\delta B_{0\text{-lv}}$ = $-2.28 \text{ mm a}^{-1}$, $\delta B_{0\text{-Pe}}$ = $-1.99 \text{ mm a}^{-1}$ and $\delta B_{0\text{-Wh}}$ = $-2.36 \text{ mm a}^{-1}$ with an average of $-2.21 \text{ mm a}^{-1}$. The corresponding $(dM/dt)_{\text{eq\text{-lv}}}$ is 150.5 Gt $\text{ a}^{-1}$, the $(dM/dt)_{\text{eq\text{-Pe}}}$ is 147.8 Gt $\text{ a}^{-1}$ and the $(dM/dt)_{\text{eq\text{-Wh}}}$ is 151.3 Gt $\text{ a}^{-1}$ with an average of 149.9 Gt $\text{ a}^{-1}$. The required $\delta B_{0\text{-avg}}$ and $\delta B_{\text{md\text{-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.

For EA, the range of $(dM/dt)_{\text{eq\text{-md}}}$ among the three models is only 2.3% of their mean compared to the larger ranges of 17% in the $\delta B_{0\text{-md}}$ and $\delta B_{\text{adj\text{-md}}}$. The smaller fractional difference among the $(dM/dt)_{\text{eq\text{-md}}}$ occurs because of its primary sensitivity to the slope $S_a$, compared to the primary sensitivity of $\delta B_{\text{adj\text{-md}}}$ and $\delta B_{0\text{-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)_{\text{eq\text{-md}}}$ are also small ($\le 2.5\%$ range). Therefore, the $(dM/dt)_{\text{eq\text{-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 $\delta B_0$ and its

<figure><img src="image_18.png" /><figcaption>**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.*</figcaption></figure>

corresponding $dB_{cor}$ and $GIA_{cor}$ are applied to $[(dM/dt)_{\text{ICESat}}]_0$ and $[(dM/dt)_{\text{GRACE}}]_0$ with no $dB_{cor}$ nor $GIA_{cor}$ applied, or whether $\delta B_{\text{adj}}$ and its corresponding $dB_{cor}$ and $GIA_{cor}$ are applied to the $[(dM/dt)_{\text{ICESat}}]_{\text{md}}$ and $[(dM/dt)_{\text{GRACE}}]_{\text{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 $\text{ a}^{-1}$ with neither $dB_{cor}$ nor $GIA_{cor}$ corrections for volume and mass changes beneath the ice sheet. The required $\delta B_{0\text{-avg}}$ for exact mass matching at 95.0 Gt $\text{ a}^{-1}$ is only $-0.35 \text{ mm a}^{-1}$, 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 $\text{ a}^{-1}$) 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 \text{ mm a}^{-1}$) 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 $\delta B_{0\text{-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 $\delta B_{0\text{-md}}$ (Table 4). In WA1, with the aforementioned near equality of the uncorrected ICESat and GRACE $dM/dt$, the $\delta B_{0\text{-md}}$ have a small range from $-0.33$ to $-0.38 \text{ mm a}^{-1}$ (Table 4) due to small range in the $S_g$ from $-2.6$ to $-2.9 \text{ Gt a}^{-1}/\text{mm a}^{-1}$ (Table 1). The same comments apply when the $\delta B_{\text{adj\text{-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 $\delta B_{0-avg}$ is $-3.48 \text{ mm a}^{-1}$ (subsidence) in contrast to the average uplifts ranging from 3.00 to $5.42 \text{ mm a}^{-1}$ 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 $\sim$15 ka BP and was mostly complete by $\sim$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 $\sim$10 ka BP includes a 400 km Holocene read advance 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 read advance is supported by their analyses of sediment cores. The ice-loading history also includes a smaller Holocene read advance 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 read advance (Bradley and others, 2015). These grounding-line retreats and read advances 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 $\sim$400 km inland of its current position (Fig. 3 of Kingslake and others, 2018) and at a location where the surface elevation is now $\sim$500 m above sea level, indicating a thinning there of $\sim$400 m at the maximum post-LGM retreat and subsequent thickening to its current elevation. Also, the ice at Siple Dome located $\sim$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 read advance of the ice sheet to the current grounding-line position, the ice sheet thickened by $\sim$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 read advance 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 $\approx$15 ka since the maximum post-LGM retreat and subsequent initiation of the read advance, 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 $\sim$5, which would reduce the exponentially decaying uplift from post-LGM ice unloading by a factor of $\sim$7 $\times 10^{-3}$. However, for the lower viscosities of $\sim 10^{18}$–$10^{19}$ Pa s and their decadal to centennial response times, the corresponding reductions would be by a much-larger factor of $\ll$2 $\times 10^{-22}$. Therefore, the primary on-going response should be subsidence from the later Holocene read advance 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 ($\delta B'$) of a long-term dynamic ice thickening using Eqn (9). Values of $\delta 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$^{-1}$ as listed in column 4 of Table 6. These $\delta B'$ are compared to the $\delta B_{0-avg}$ adjustments (i.e. relative to zero dB/dt) for the averages of the three $\delta B_{0-md}$ (listed in column 5, taken from column 3 of Table 4) showing that the $\delta B'$ estimated from the regional dynamic thickenings are 1.2–3.1 times larger (column 6) than the required $\delta B_{0-avg}$ and 1.1–1.7 times larger than the $\delta 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.

For the ICESat analysis, we use the Ivins dB/dt grid plus the regional $\delta 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 $\delta B_{adj-md}$ are added. The $\delta B_{cor} = -(S_a) [(\delta B/dt)_{IV} + \delta 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} (\delta B_0)_{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.

## 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.

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 $(\delta B_{cor})_{2015}$ in Table 2 replaced with the $\delta B_{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).

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 Gt 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

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** ICESat maps for 2003–2008, (a) $dM/dt$, (b) $dM_d/dt$ and (c) $dM_a/dt$ using $dB/dt$ equal to $\text{IvinsdB/dt} + \delta B_{\text{adj}}$. Rates are linear terms of LQS fits at year 2006.0. \*Rates for AP from average-linear-change analysis.</figcaption></figure>

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
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Maps of the calculated firn density $\rho_a = \Delta M_a / \Delta(H^3 - 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.</figcaption></figure>

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 $\delta 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 $\delta 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
<table><caption>**Table 7a.** Summary of linear rates of mass change (dM/dt) from ERS1/2, ICESat and GRACE for select periods during 1992–2016</caption><thead><tr><th rowspan="2">Period<br>Mid-point</th><th>ERS1/2<sup>a</sup></th><th>ICESat</th><th colspan="3">GRACE</th></tr><tr><th>1992–2001<br>1996.5<br>Gt a<sup>-1</sup></th><th>2003–2008<br>2006.0<br>Gt a<sup>-1</sup></th><th>2003–2008<br>2006.0<br>Gt a<sup>-1</sup></th><th>2009–2011<br>2010.5<br>Gt a<sup>-1</sup></th><th>2012–2016<br>2014.5<br>Gt a<sup>-1</sup></th></tr></thead><tbody><tr><td>EA</td><td>161 ± 50</td><td>150 ± 28</td><td>150 ± 21</td><td>257 ± 76</td><td>134 ± 58</td></tr><tr><td>EA1</td><td>101 ± 33</td><td>79 ± 23</td><td>79 ± 15</td><td>196 ± 46</td><td>88 ± 35</td></tr><tr><td>EA2</td><td>60 ± 21</td><td>71 ± 26</td><td>71 ± 16</td><td>61 ± 60</td><td>45 ± 47</td></tr><tr><td>WA</td><td>-8 ± 20</td><td>-26 ± 15</td><td>-29 ± 12</td><td>-187 ± 23</td><td>-116 ± 24</td></tr><tr><td>WA1</td><td>-59 ± 12</td><td>-95 ± 6</td><td>-95 ± 9</td><td>-214 ± 18</td><td>-165 ± 15</td></tr><tr><td>WA2</td><td>51 ± 14</td><td>69 ± 9</td><td>66 ± 9</td><td>27 ± 15</td><td>49 ± 19</td></tr><tr><td>AP</td><td>-9 ± 10</td><td>-29<sup>a</sup> ± 2</td><td>-24 ± 9</td><td>-36 ± 15</td><td>-30 ± 19</td></tr><tr><td>AIS</td><td>144 ± 61</td><td>95 ± 25</td><td>97 ± 26</td><td>34 ± 85</td><td>-12 ± 64</td></tr></tbody></table>

<table><caption>**Table 7b.** Summary of changes (delta) in the linear rates of mass change between periods compared to the annual SMB</caption><thead><tr><th rowspan="2"></th><th rowspan="2">SMB<sup>a</sup><br>Gt a<sup>-1</sup></th><th colspan="2">Change from (1992–2001)<br>to (2003–2008)</th><th colspan="2">Change from (2003–2008)<br>to (2009–2011)</th><th colspan="2">Change from (2009–2011)<br>to (2012–2016)</th><th colspan="2">Change from (1992–2001)<br>to (2012–2016)</th></tr><tr><th>delta<br>Gt a<sup>-1</sup></th><th>delta/SMB<br>(%)</th><th>delta<br>Gt a<sup>-1</sup></th><th>delta/SMB<br>(%)</th><th>delta<br>Gt a<sup>-1</sup></th><th>delta/SMB<br>(%)</th><th>delta<br>Gt a<sup>-1</sup></th><th>delta/SMB<br>(%)</th></tr></thead><tbody><tr><td>EA</td><td>1145</td><td>0<sup>b</sup></td><td>0</td><td>107</td><td>9</td><td>-124</td><td>-11</td><td>-17</td><td>-1</td></tr><tr><td>EA1</td><td>463</td><td>-11<sup>b</sup></td><td>-2</td><td>117</td><td>25</td><td>-108</td><td>-23</td><td>-2</td><td>0</td></tr><tr><td>EA2</td><td>683</td><td>11</td><td>2</td><td>-10</td><td>-1</td><td>-16</td><td>-2</td><td>-15</td><td>-2</td></tr><tr><td>WA</td><td>501</td><td>-18</td><td>-4</td><td>-158</td><td>-32</td><td>71</td><td>14</td><td>-105</td><td>-21</td></tr><tr><td>WA1</td><td>221</td><td>-36</td><td>-16</td><td>-119</td><td>-54</td><td>49</td><td>22</td><td>-106</td><td>-48</td></tr><tr><td>WA2</td><td>281</td><td>18</td><td>6</td><td>-39</td><td>-14</td><td>22</td><td>8</td><td>1</td><td>0</td></tr><tr><td>AP</td><td>196</td><td>-20</td><td>-10</td><td>-11</td><td>-6</td><td>6</td><td>3</td><td>-25</td><td>-13</td></tr><tr><td>AIS</td><td>1843</td><td>-39<sup>b</sup></td><td>-2</td><td>-62</td><td>-3</td><td>-47</td><td>-3</td><td>-109</td><td>-6</td></tr></tbody></table>

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<sup>-1</sup> (-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<sup>-1</sup> (+22% of SMB), likely due to an unidentified slowing of glacier discharge. In the AP, the 20 Gt a<sup>-1</sup> (-10% of SMB) loss-rate increase in the first interval (1992–2001 to 2003–08) was

<table><caption>**Table 8.** ICESat laser campaign biases determined over leads and polynya in sea ice</caption><thead><tr><th>Campaign</th><th>Boreal season/year</th><th>DSL (m)</th></tr></thead><tbody><tr><td>L2a</td><td>F03</td><td>-0.258</td></tr><tr><td>L2b</td><td>W04</td><td>-0.233</td></tr><tr><td>L2c</td><td>S04</td><td>-0.230</td></tr><tr><td>L3a</td><td>F04</td><td>-0.248</td></tr><tr><td>L3b</td><td>W05</td><td>-0.303</td></tr><tr><td>L3c</td><td>S05</td><td>-0.281</td></tr><tr><td>L3d</td><td>F05</td><td>-0.302</td></tr><tr><td>L3e</td><td>W06</td><td>-0.301</td></tr><tr><td>L3f</td><td>S06</td><td>-0.331</td></tr><tr><td>L3g</td><td>F06</td><td>-0.302</td></tr><tr><td>L3h</td><td>W07</td><td>-0.290</td></tr><tr><td>L3i</td><td>F07</td><td>-0.305</td></tr><tr><td>L3j</td><td>W08</td><td>-0.327</td></tr><tr><td>L3k & L2d</td><td>F08</td><td>-0.287</td></tr><tr><td>L2e</td><td>S09</td><td>-0.262</td></tr><tr><td>L2f</td><td>F09</td><td>-0.310</td></tr></tbody></table>

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

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<sup>-1</sup> (-6% of SMB) between 2003–2008 and 2009–2011, which was followed by a loss-rate decrease of 6 Gt a<sup>-1</sup> (+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<sup>-1</sup> increase (+25% of SMB) in EA1 between 2003–08 and 2009–11, followed by a 108 Gt a<sup>-1</sup> 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<sup>-1</sup> (+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<sup>-1</sup> (-1% of SMB) between 2003–08 and 2009–11, followed by a 16 Gt a<sup>-1</sup> decrease (-2% of SMB) between 2009–11 and 2012–16 giving a net decrease of 26 Gt a<sup>-1</sup> (-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<sup>-1</sup> (-1.5% of SMB) from 2003 to 2016, which is most likely not the cause of the mass gain of 90 Gt a<sup>-1</sup> 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<sup>-1</sup> (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

<footnote><sup>a</sup>AP ICESat and all ERS from average-linear-change analysis (Zwally and others, 2015) with ERS using adjusted dB<sub>cor</sub> from Table 3.</footnote>

<footnote><sup>a</sup>SMB from Giovinetto and Zwally (2000) and by drainage systems and regions in Zwally and others (2015).</footnote>

<footnote><sup>b</sup>These delta are adjusted by 11 Gt a<sup>-1</sup> to account for the difference between the 11 Gt a<sup>-1</sup> larger ICESat dM/dt from the prior average-linear-change analysis (see Table 2) as was also used for ERS1/2.</footnote>
(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<sub>adj–IV</sub> = −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 read advance (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₀–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 read advance 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 read advance of the grounding line of the Ross Shelf from its inland retreat at 9.7 ka BP and a smaller read advance 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 read advances 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⁻¹ (adjusted from the 27.3 Gt a⁻¹ 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-millennial 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 read advance 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 GIAcor and dBcor are based on the simple principle that the
<table><caption>**Table 9.** Accumulation density ($\rho_a$) and pseudo density ($\rho_{\text{pseudol}}$) by region</caption><thead><tr><th rowspan="2">Region</th><th colspan="2">$\rho_a$</th><th colspan="2">$\rho_{\text{pseudol}}$</th></tr><tr><th>1992-2001</th><th>2003-2008</th><th>1992-2001</th><th>2003-2008</th></tr></thead><tbody><tr><td>WA1</td><td>0.51</td><td>0.49</td><td>0.93</td><td>1.01</td></tr><tr><td>WA2</td><td>0.44</td><td>0.45</td><td>1.24</td><td>0.79</td></tr><tr><td>WA</td><td>0.46</td><td>0.46</td><td>0.55</td><td>5.78</td></tr><tr><td>EA1</td><td>0.35</td><td>0.36</td><td>1.07</td><td>0.67</td></tr><tr><td>EA2</td><td>0.41</td><td>0.39</td><td>1.22</td><td>2.86</td></tr><tr><td>EA</td><td>0.37</td><td>0.37</td><td>1.12</td><td>0.93</td></tr><tr><td>AP</td><td>0.61</td><td>0.59</td><td>1.53</td><td>1.17</td></tr><tr><td>AIS</td><td>0.39</td><td>0.39</td><td>1.27</td><td>0.70</td></tr></tbody></table>

$\rho_a$ is density associated with $\delta A(t)$ anomalies.

$\rho_{\text{pseudol}} = dM/dt/dl/dt$.

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 $\rho_{\text{earth}}$ that are consistent with the knowledge of mantle densities, that physical correspondence is not essential for making the $\delta B$ adjustments. However, we believe that the physical relationship implied by the consistency of the $\rho_{\text{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.

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.

**Supplementary material.** The supplementary material for this article can be found at https://doi.org/10.1017/jog.2021.8.

**Acknowledgements.** We appreciate the GIA and $dB/dt$ model data provided by E. Ivins, P. Whitehouse and R. Peltier. NASA GSFC mascon solutions were developed partly under the NASA GRACE and GRACE Follow-On Science Team Grant NNH15ZDA0-01N. This research was also supported by NASA’s ICESat project science funding. We appreciate the reviewers’ efforts and their helpful comments, especially those of the fourth reviewer with expertise in GIA processes and an anonymous colleague for thorough reading of the manuscript and helpful suggestions. Following our original submission, the Editors’ requested addressing reviewers’ concerns about why our ICESat results on Antarctic mass gains differ from others in the literature. That request resulted in additional documentation of the validity of our ICESat and ERS1/2 results and our methodologies and reasons why we believe the results of others had some significant errors, now in Section 5 and the Appendix. It also led to the inclusion of the supportive Envisat results and co-author F. Rémy. We also greatly appreciate the dedicated efforts of so many scientists, engineers, analysts, and managers in NASA and supporting companies that led to the successes of the ICESat missions and similarly of the many individuals involved in the successes of GRACE and the ESA missions of ERS and Envisat.

## References

*Adodo FI, Rémy F and Picard G* (2018) Seasonal variations of the backscattering coefficient measured by radar altimeters over the Antarctic ice sheet. *The Cryosphere* **12**, 1767–1778. doi: 10.5194/tc-12-1767-2018

*Argus DF, Peltier WR, Drummond R and Moore AW* (2014) The Antarctica component of postglacial rebound model ICE-6G_C (VM5a) based upon GPS positioning, exposure age dating of ice thicknesses, and relative sea level histories. *Geophysical Journal International* **198**(1), 537–563. doi: 10.1093/gji/ggu140
Arthern R, Wingham D and Ridout A (2001) Controls on ERS altimeter measurements over ice sheets: footprint-scale topography, backscatter fluctuations, and the dependence of microwave penetration depth on satellite orientation. *Journal of Geophysical Research: Atmospheres* **106**, 33471–33484. doi: 10.1029/2001JD000498.

Bamber JL (1994) Ice sheet altimeter processing scheme. *International Journal of Remote Sensing* **15**(4), 925–938.

Bamber JL and Payne AJ (eds) (2004) *Mass Balance of the Cryosphere: Observations and Modelling of Contemporary and Future Changes*. Cambridge: Cambridge University Press.

Bamber JL, Westaway RM, Marzeion B and Wouters B (2018) The land ice contribution to sea level during the satellite era. *Environmental Research Letters* **13**, 063008. doi: 10.1088/1748-9326/aac2f0.

Barletta VR and 16 others (2018) Observed rapid bedrock uplift in Amundsen Sea Embayment promotes ice-sheet stability. *Science* (New York, N.Y.) **360**(6395), 1335–1339. doi: 10.1126/science.aa01447.

Bentley MJ and 87 others (2014) A community-based geological reconstruction of Antarctic ice sheet deglaciation since the last glacial maximum. *Quaternary Science Reviews* **100**, 1–9. doi: 10.1016/j.quascirev.2014.06.025.

Bradley SL, Hindmarsh RCA, Whitehouse PL, Bentley MJ and King MA (2015) Low post-glacial rebound rates in the Weddell Sea due to Late Holocene ice-sheet readvance. *Earth and Planetary Science Letters* **413**, 79–89. doi: 10.1016/j.epsl.2014.12.039.

Davis CH (1997) A robust threshold retracking algorithm for measuring ice sheet surface elevation change from satellite radar altimeters. *IEEE Transactions on Geoscience and Remote Sensing* **35**(4), 974–979.

Davis CH and Ferguson AC (2004) Elevation change of the Antarctic ice sheet, 1995–2000, from ERS-2 satellite radar altimetry. *IEEE Transactions on Geoscience and Remote Sensing* **42**(11), 2437–2445.

Davis CH, Li Y, McConnell JR, Frey MM and Hanna E (2005) Snowfall-driven growth in East Antarctic ice sheet mitigates recent sea-level rise. *Science* (New York, N.Y.) **308**.

Dejong BD and 5 others (2015) Pleistocene relative sea levels in the Chesapeake Bay region and their implications for the next century. *GSA Today* **25**(8). doi: 10.1130/GSATG223A.1.

Denton GH (2011) East Antarctic retreat. *Nature Geoscience* **4**, 135–136.

Denton GH and Hughes TJ (eds) (1981) *The Last Great Ice Sheets*. New York: Wiley, pp. 319–436.

Denton GH and Wilson SC (1982) Late quaternary geology of the Rennick Glacier area, northern Victoria Land. *Antarctic Journal of the United States* **17**, 49–52.

Farrell SL and 6 others (2012) Mean dynamic topography of the Arctic ocean. *Geophysical Research Letters* **39**, L01601. doi: 10.1029/2011GL05005.

Farrell SL, Laxon SW, McAdoo DC, Yi D and Zwally HJ (2009) Five years of Arctic sea ice freeboard measurements from the ice, cloud and land elevation satellite. *Journal of Geophysical Research* **114**, C04008. doi: 10.1029/2008JC005074.

Femenias P (1996) ERS QLOPR and OPR range processing. ESA/ESRIN Tech. Note ER-TN-RS-RA-0022. Frascati, European Space Agency/European Space Research Institute.

Flament T and Rémy F (2012) Dynamic thinning of Antarctic glaciers from along-track repeat radar altimetry. *Journal of Glaciology* **58**(211), 830–840. doi: 10.3189/2012JoG11J11.

Giles KA, Laxon SW, Ridout AL, Wingham DJ and Bacon S (2012) Western Arctic Ocean freshwater storage increased by wind-driven spin-up of the Beaufort Gyre. *Nature Geoscience* **5**(3), 194–197. doi: 10.1038/NGEO1379.

Giovinetto MB and Zwally HJ (1996) Annual changes in sea ice extent and accumulation on ice sheets: implications for sea level change. *Zietschrift fur Gletscherkunde und Glazialgeologie* **31**, 39–49.

Giovinetto MB and Zwally HJ (2000) Spatial distribution of net surface accumulation on the Antarctic ice sheet. *Annals of Glaciology* **31**, 171–178.

Golledge NA and 5 others (2019) The multi-millennial Antarctic commitment to future sea-level rise. *Nature* **9**, 34–39. doi: 10.1038/nature15706.

Hanna E and 11 others (2013) Ice sheet mass balance and climate change. *Nature* **498**, 51–59. doi: 10.1038/nature1223.

Hanna E and 10 others (2020) Mass balance of the ice sheets and glaciers – progress since AR5 and challenges. *Earth-Science Reviews* **201**, 102976. doi: 10.1016/j.earscirev.2019.102976.

Helm V, Humbert A and Miller H (2014) Elevation and elevation change of Greenland and Antarctica derived from CryoSat-2. *The Cryosphere* **8**, 1539–1559. doi: 10.5194/tc-8-1539-2014.

Huybrechts P (2004) Antarctica: modelling. In Bamber JL and Payne AJ eds. *Mass Balance of the Cryosphere: Observations and Modelling of Contemporary and Future Changes*. Cambridge: Cambridge University Press.

Ivins ER and 5 others (2013) Antarctic contribution to sea-level rise observed by GRACE with improved GIA correction. *Journal of Geophysical Research: Solid Earth* **118**(6), 3126–3141. doi: 10.1002/jgrb.50208.

Joughin I, Tulaczyk S, Bindschadler R and Price SF (2002) Changes in West Antarctic ice stream velocities: observation and analysis. *Journal of Geophysical Research* **107**(B11). doi: 10.1029/2001JB001029.

Jouzel J and 31 others (2007) Orbital and millennial Antarctic climate variability over the past 800,000 years. *Science* (New York, N.Y.) **317**. 793–795. doi: 10.1126/science.1141038.

Khvorostovsky KS (2012) Merging and analysis of elevation time series over Greenland ice sheet from satellite radar altimetry. *IEEE Transactions on Geoscience and Remote Sensing* **50**(1), 23–36. doi: 10.1109/TGRS.2011.2160071.

Kingslake J and 9 others (2018) Extensive retreat and re-advance of the West Antarctic ice sheet during the Holocene. *Nature*. doi: 10.1038/s41586-018-0208-x.

Larour E and 6 others (2019) Slowdown in Antarctic mass loss from solid Earth and sea-level feedbacks. *Science* (New York, N.Y.) **364**, eaav7908. doi: 10.1126/science.aav7908.

Larsen CF, Motyka RJ, Freymueller JT, Echelmeyer KA and Ivins ER (2005) Rapid viscoelastic uplift in southeast Alaska caused by post-Little Ice Age glacial retreat. *Earth and Planetary Science Letters* **237**, 548–560. doi: 10.1016/j.epsl.2005.06.032.

Legresy B, Rémy F and Schaeffer P (1999) Different ERS altimeter measurements between ascending and descending tracks caused by wind induced features over ice sheets. *Geophysical Research Letters* **26**(15), 2231–2234.

Li X, Rignot E, Mouginot J and Scheuchl B (2016) Ice flow dynamics and mass loss of Totten Glacier; East Antarctica from 1989 to 2015. *Geophysical Research Letters* **43**, 6366–6373. doi: 10.1002/2016GL069173.

Li J and Zwally HJ (2011) Modeling of firn compaction for estimating ice-sheet mass change from observed ice-sheet elevation change. *Annals of Glaciology* **52**, 1–7.

Li J and Zwally HJ (2015) Response times of ice-sheet surface heights to changes in the rate of Antarctic firn compaction caused by accumulation and temperature variations. *Journal of Glaciology* **61**(230), 1037–1047. doi: 10.3189/2015JoG14J082.

Ligtenberg SRM, van de Berg MJ, van den Broeke MR, Rae JGL and van Meijgaard E (2013) Future surface mass balance of the Antarctic ice sheet and its influence on sea level change, simulated by a regional atmospheric climate model. *Climate Dynamics* **41**, 867–884. doi: 10.1007/s00382-013-1749-1.

Loomis BD, Luthcke SB and Sabaka TJ (2019) Regularization and error characterization of GRACE mascons. *Journal of Geodesy*. doi: 10.1007/s00190-019-01252-y].

Luthcke SB and 5 others (2013) Antarctica, Greenland and Gulf of Alaska land-ice evolution from an iterated GRACE global mascon solution. *Journal of Glaciology* **59**(216), 613–631. doi: 10.3189/2013JoG12J147.

Luthcke SB and 5 others (2015) Chapter 10 in Gravimetry measurements from space. In Tedesco M (ed.), *Remote Sensing of the Cryosphere*. Oxford: Wiley Blackwell, pp. 231–247. ISBN 9781118368855.

Mackintosh A and 11 others (2011) Retreat of the East Antarctic ice sheet during the last glacial termination. *Nature Geoscience* **4**, 195–202. doi: 10.1038/NGEO1061.

Martin-Español A, Bamber JL and Zammit-Mangion A (2017) Constraining the mass balance of East Antarctica. *Geophysical Research Letters* **44**, 4168–4175. doi: 10.1002/2017GL072937.

Martin TV, Zwally HJ, Brenner AC and Bindschadler RA (1983) Analysis and retracking of continental ice sheet radar altimeter waveforms. *Journal of Geophysical Research* **88**, 1608–1616.

McMillan M and 7 others (2014) Increased ice losses from Antarctica detected by CryoSat-2. *Geophysical Research Letters* **41**. doi: 10.1002/2014GL060111.

McMillan M and 12 others (2016) A high-resolution record of Greenland mass balance. *Geophysical Research Letters* **43**, 7002–7010.

Medley B and 12 others (2013) Airborne-radar and ice-core observations of annual snow accumulation over Thwaites Glacier, West Antarctica confirm the spatiotemporal variability of global and regional atmospheric models.
*Geophysical Research Letters* **40**, 3649–3654. doi: 10.1002/grl.50706.

**Medley B and 6 others** (2017) Temperature and snowfall in Western Queen Maud Land increasing faster than climate model projections. *Geophysical Research Letters* **45**, 1472–1480. doi: 10.1002/2017GL075992.

**Medley B and Thomas ER** (2019) Increased snowfall over the Antarctic ice sheet mitigated twentieth-century sea-level rise. *Nature Climate Change* **9**, 34–39. doi: 10.1038/s41558-018-0356-x.

**Moholdt G, Nuth C, Hagen JO and Kohler J** (2010) Recent elevation changes of Svalbard glaciers derived from ICESat laser altimetry. *Remote Sensing of Environment* **114**(11), 2756–2767. doi: 10.1016/j.rse.2010.06.008.

**Nield GA and 9 others** (2014) Rapid bedrock uplift in the Antarctic Peninsula explained by viscoelastic response to recent ice unloading. *Earth and Planetary Science Letters* **397**, 32–41. doi: 10.1016/j.epsl.2014.04.019.

**Nielsen K and 6 others** (2013) Vertical and horizontal surface displacements near Jakobshavn Isbræ driven by melt-induced and dynamic ice loss. *Journal of Geophysical Research: Atmospheres* **118**(4), 1837–1844. doi: 10.1002/jgb.50145.

**Nilsson J, Gardner A, Sørensen SL and Forsberg R** (2016) Improved retrieval of land ice topography from CryoSat-2 data and its impact for volume-change estimation of the Greenland ice sheet. *The Cryosphere* **10**, 2953–2969. doi: 10.5194/tc-10-2953-2016.

**Partington KC, Ridley JK, Rapley CG and Zwally HJ** (1989) Observations of the surface properties of the ice sheets by satellite radar altimetry. *Journal of Glaciology* **35**(120), 267–275.

**Peltier WR** (2004) Global glacial isostasy and the surface of the ice age earth: the ICE-5G (VM2) model and GRACE. *Annual Review of Earth and Planetary Sciences* **32**, 111–149. doi: 10.1146/annurev.earth.32.082503.144359.

**Peltier WR** (2014) http://www.atmosp-physics.utoronto.ca/~peltier/data.php.

**Pollard D, Chang W, Haran W, Applegate M and DeConto R** (2016) Large ensemble modeling of last deglacial retreat of the West Antarctic ice sheet: comparison of simple and advanced statistical techniques. *Geoscientific Model Development* **9**, 1697–1723. doi: 10.5194/gmd-9-1697-2016.

**Pollard D, Gomez N and DeConto RM** (2017) Variations of the Antarctic ice sheet in a coupled ice sheet-Earth-sea level model: sensitivity to viscoelastic Earth properties. *Journal of Geophysical Research: Earth Surface* **122**, 2124–2138. doi: 10.1002/2017JF004371.

**Price SF, Conway H and Waddington ED** (2007) Evidence for late Pleistocene thinning of Siple Dome, West Antarctica. *Journal of Geophysical Research: Earth Surface* **112**, F03021.

**Pritchard HD, Arthern RJ, Vaughan DG and Edwards LA** (2009) Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets. *Nature* **461**, 971–975. doi: 10.1038/nature08471.

**Pritchard H and Vaughan D** (2007) Widespread acceleration of tidewater glaciers on the Antarctic Peninsula. *Journal of Geophysical Research* **112**, doi: 10.1029/2006JF000597.

**Rémy F, Flament T, Blarel F and Benveniste J** (2012) Radar altimetry measurements over Antarctic ice sheet: a focus on antenna polarization and change in backscatter problems. *Advances in Space Research* **50**.

**Richter A and 11 others** (2016) Crustal deformation across the Southern Patagonian Icefield observed by GNSS. *Earth and Planetary Science Letters* **452**, 206–215. doi: 10.1016/j.epsl.2016.07.042.

**Rignot EJ and 6 others** (2008) Recent Antarctic ice mass loss from radar interferometry and regional climate modelling. *Nature Geoscience* **1**, 106–110.

**Rignot E and 5 others** (2019) Four decades of Antarctic ice sheet mass balance from 1979–2017. *Proceedings of the National Academy of Sciences of the USA* **116**(4), 1095–1103.

**Rignot E, Vaughan DG, Schmeltz M, Dupont T and MacAyeal D** (2002) Acceleration of Pine Island and Thwaites Glaciers, West Antarctica. *Annals of Glaciology* **34**, 189–194.

**Riva REM and 8 others** (2009) Glacial isostatic adjustment over Antarctica from combined ICESat and GRACE satellite data. *Earth and Planetary Science Letters* **288**, 516–523. doi: 10.1016/j.epsl.2009.10.013.

**Robinson EC** (2011) The Interior of the Earth. https://pubs.usgs.gov/gip/interior/.

**Rott H, Müller F, Nagler T and Floricioiu D** (2011) The imbalance of glaciers after disintegration of Larsen-B ice shelf, Antarctic Peninsula. *The Cryosphere* **5**, doi: 10.5194/tc-5-125-2011.

**Scambos T and Shuman C** (2016) Comment on mass gains of the Antarctic ice sheet exceed losses by H. J. Zwally and others. *Journal of Glaciology* **62**(233), 599–603. doi: 10.1017/jog.2016.59.

**Scharroo R and 5 others** (2013) RADS: Consistent multi-mission products. *Proc. Symp. on 20 Years of Progress in Radar Altimetry, Venice, 20–28 September 2012*, European Space Agency Specification Publ., ESA SP-710, p. 4.

**Schroder L and 5 others** (2019) Four decades of Antarctic surface elevation changes from multi-mission satellite altimetry. *Cryosphere* **13**, 427–449.

**Shepherd A and 46 others** (2012) A reconciled estimate of ice-sheet mass balance. *Science* **338**, 1183–1189. doi: 10.1126/science.1228102.

**Shepherd A and 68 others** (2018) Mass balance of the Antarctic ice sheet from 1992 to 2017. *Nature* **558**, 219–222. doi: 10.1038/s41586-018-0179-y.

**Shepherd A and 5 others** (2019) Trends in Antarctic ice sheet elevation and mass. *Geophysical Research Letters* **46**, 8174–8183. doi: 10.1029/2019GL082182.

**Shuman CA, Berthier E and Scambos TA** (2011) 2001–2009 Elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula. *Journal of Glaciology* **57**(204), 737–754. doi: 10.3189/002214311797409811.

**Siegert MJ** (2003) Glacial-interglacial variations in central East Antarctic ice accumulation rates. *Quaternary Science Reviews* **22**.

**Siegert MJ and Payne AJ** (2004) Past rates of accumulation in central West Antarctica. *Geophysical Research Letters* **31**, L12403. doi: 10.1029/2004GL020290.

**Sinisalo A and Moore J** (2010) Antarctic blue ice areas – towards extracting palaeoclimate information. *Antarctic Science* **22**(2), 99–115. doi: 10.1017/S0954102009990691.

**Smith BE and 14 others** (2020) Pervasive ice sheet mass loss reflects competing ocean and atmosphere processes. *Science* **368**(6496), 1239–1242. doi: 10.1126/science.aaz5845.

**Smith BE, Fricker HA, Joughin IR and Tulaczyk S** (2009) An inventory of active subglacial lakes in Antarctica detected by ICESat (2003–2008). *Journal of Glaciology* **55**(192), 573–595. doi: 10.3189/002214309789470879.

**Stuiver M, Denton GH, Hughes TJ and Fastook JL** (1981) History of the marine ice sheet in West Antarctica during the last glaciation: a working hypothesis. *In The Last Great Ice Sheets* (eds Denton GH and Hughes TJ). New York: Wiley

**Thomas R and 17 others** (2004) Accelerated sea-level rise from West Antarctica. *Science* **306**, 255–258. doi: 10.1126/science.1099650.

**Vieli GJ-MC, Siegert MJ and Payne AJ** (2004) Reconstructing ice sheet accumulation rates at ridge B, East Antarctica. *Annals of Glaciology* **39**, 326–328. doi: 10.3189/172756404781814519.

**Wahr J, Wingham D and Bentley C** (2000) A method of combining ICESat and GRACE satellite data to constrain Antarctic mass balance. *Journal of Geophysical Research* **105**(B7), 16279–16294.

**Wang W, Li J and Zwally HJ** (2012) Dynamic inland propagation of thinning due to ice loss at the margins of the Greenland ice sheet. *Journal of Glaciology* **58**(210), 734–740.

**Wang W, Li J and Zwally HJ** (2013) Modeling dynamic thickening in East Antarctica as observed from ICESat, Abstract C53B-0571 presented at 2013 Fall Meeting, AGU, San Francisco, California, 9–13 December.

**Whitehouse PL, Bentley MJ, Milne G, King M and Thomas I** (2012) A new glacial isostatic adjustment model for Antarctica: calibrated and tested using observations of relative sea-level change and present-day uplift rates. *Geophysical Journal International* **190**(3), 1464–1482. doi: 10.1111/j.1365-246X.2012.05557.x.

**Wingham DJ, Ridout AL, Scharroo R, Arthern RJ and Shum CK** (1998) Antarctic elevation change 1992 to 1996. *Science* (New York, N.Y.) **282** (5388), 456–458.

**Wingham DJ, Shepherd A, Muir A and Marshall GJ** (2006) Mass balance of the Antarctic ice sheet. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **364**, 1627–1635. doi: 10.1098/rsta.2006.1792.

**Wingham DJ, Wallis DW and Shepherd A** (2009) Spatial and temporal evolution of Pine Island Glacier thinning, 1995–2006. *Geophysical Research Letters* **36**, 189–194. doi: 10.1029/2009GL039126.

**Yi D, Zwally HJ, Cornejo HG, Barbieri KA and DiMarzio JP** (2011) Sensitivity of elevations observed by satellite radar altimeter over ice sheets to variations in backscatter power and derived corrections. *CryoSat Validation Workshop, Frascati, ESA SP-693.

**Zwally HJ** (2013) Correction to the ICESat data product surface elevations due to an error in the range determination from transmit-pulse reference-point selection (Centroid vs Gaussian). (Tech. rep.) National Snow and Ice Data Center, Boulder, CO. http://nsidc.org/data/icesat/correction-to-product-surface-elevations.html.
**Zwally HJ and 7 others** (2005) Mass changes of the Greenland and Antarctic ice sheets and shelves and contributions to sea-level rise: 1992–2002. *Journal of Glaciology* **51**(175), 509–527. doi: **10.3189/172756505781829007**.

**Zwally HJ and 11 others** (2011) Greenland Ice sheet mass balance: distribution of increased mass loss with climate warming. *Journal of Glaciology* **57**(201), 88–102. doi: **10.3189/002214311795306682**.

**Zwally HJ and 5 others** (2015) Mass gains of the Antarctic ice sheet exceed losses. *Journal of Glaciology* **61**(230), 1019–1035. doi: **10.3189/2015JoG15J071**.

**Zwally HJ and 7 others** (2016) Response to comment by T. SCAMBOS and C. SHUMAN (2016) on ‘mass gains of the Antarctic ice sheet exceed losses’ by H. J. Zwally and others (2015). *Jounal of Glaciology* **62**(235), 990–992. doi: **10.1017/jog.2016.91**.

**Zwally HJ and Brenner AC** (2001) Ice sheet dynamics and mass balance. In Fu L-L and Cazanave A (eds), *Satellite Altimetry and Earth Sciences*. New York: Academic Press Inc., pp. 351–369. (International Geophysical Series 69.).

**Zwally HJ, Brenner AC, Major JA, Bindschadler RA and Marsh JG** (1989) Growth of Greenland ice sheet: measurement. *Science* (New York, N.Y.) **246**, 1587–1589.

**Zwally HJ and Giovinetto MB** (2011) Overview and assessment of Antarctic ice-sheet mass balance estimates: 1992–2009. *Surveys in Geophysics* **32**, 351. doi: **10.1007/s10712-011-9123-5**.

**Zwally HJ and Li J** (2002) Seasonal and interannual variations of firn densification and ice-sheet surface elevation at Greenland summit. *Journal of Glaciology* **48**(161), 199–207.

**Zwally HJ, Yi D, Kwok R and Zhao Y** (2008) ICESat measurements of sea ice freeboard and estimates of sea ice thickness in the Weddell Sea). *Journal of Geophysical Research* **113**(C2), C02S15. doi: **10.1029/2007J0C004284**.

# Appendix

## Introduction

We examine the compatibility of elevation changes derived from satellite altimeters including basic corrections made to the data, the methods to obtain valid ice-sheet elevation changes, and the methods to derive mass changes from the elevation changes. We review our methods and provide reasons why our results differ from some studies and agree with others. The first type of reason includes differences in the various corrections and calibrations applied in the data processing and those that may be developed later by investigators. For radar altimetry, a second reason is differences in the methods of correcting for the highly-variable penetration of the radar signal into the firn and the depth of the backscatter signal detected by the altimeter, from which the range to the surface is derived, thereby affecting the derived $H(t)$ and $dH(t)/dt$. The third reason is differences in the methods of deriving mass changes from the measured elevation changes, which includes (1) accounting for the densities of the firn and ice that are associated with the elevation changes, (2) corrections for firn compaction (FC), and (3) correction for the $dB/dt$ bedrock motion, thereby affecting the $M(t)$ and $dM(t)/dt$.

## Basic corrections and elevation-change analysis

An example of the first reason from Zwally and others (2005) is: 'Instrument corrections include subtraction of a 40.9 cm bias from ERS-1 elevations to account for a different instrument parameter used for ERS-2 (Femenias, 1996) and corrections for drifts in the ultra-stable oscillator and bias changes in the scanning point target response that are obtained from the European Space Agency'. Those corrections required application by the data users and are not necessarily applied nor noted in publications. A second example is the correction for the ERS-1/ERS-2 inter-satellite elevation bias that was discovered and empirically-determined during 13 months of simultaneous operation; from Zwally and others (2005): 'The bias correction lowers the ERS-2 elevations by an average of ... 17.5 cm ... over Antarctic grounded ice and by 12.0 cm... over Antarctic floating ice. ... the correction lowers the average $dH/dt$ by 2.4 cm $a^{-1}$... on grounded ice and by 1.6 cm $a^{-1}$... on floating ice. The effects ... on calculations of mass change ($dM/dt$) for the ERS gridpoints are roughly ... -205 Gt $a^{-1}$ for Antarctica ... indicating the importance of this correction. Davis and others (2005) in effect apply a bias correction by calculating separate $H(t)$ series for ERS-1 and ERS-2 and adjusting them together during the 12 month overlap period, but do not state the magnitude of their adjustments'. This elevation bias was

very spatially variable over the ice sheet and at least partially related to the surface slope.

Another important factor is our use of ERS ice-mode data only, because we found that ocean-mode only and mixed-mode data had differing biases that were also spatially variable and difficult to determine. Davis and others (2005) also used ice-mode data only that were obtained with corrections from our reprocessing of ESA provided data. At this level, it is possible to inter-compare results from some studies, but not all.

Another factor affecting the accuracy of the derived elevation changes is the methods used for crossover analysis and construction of elevation time series from which $dH/dt$ is derived. Our methodology (Zwally and Brenner, 2001; Zwally and others, 2005) includes two important features that affect the accuracy: (1) the averaging of elevation differences at ascending-descending crossover with those at descending-ascending crossover differences according to Eqn (20) in Zwally and Brenner (2001) [a method first used in Zwally and others (1989) to remove orbital biases but also removes the effects of penetration (Arthern and others, 2001)], and (2) the construction of time-series from crossover differences that uses not only crossovers between the first repeat cycle and all successive repeat cycles, which gives N terms for N repeat cycles including N pairings of crossover differences (e.g. Wingham and others, 1998), but also uses crossovers between the second repeat cycle and all successive cycles, plus between the third repeat and all successive repeats, and so forth constructing a series also with N terms but includes $N^2/2$ pairings of independent crossover differences. The quality of the time series in select 50 km squares from which the $dH/dt$ are calculated was shown in Figures 3 and 4 in Zwally and others (2005).

## ICESat inter-campaign biases and G-C error correction

As described in Zwally and others (2015): 'We use methods ... used in ... mapping of the level of open water and thin ice in leads and polynyas in sea ice by ICESat in the Antarctic (Zwally and others, 2008) and the Arctic (Farrell and others, 2009), in the joint mapping by ICESat and Envisat of the mean dynamic topography in the Arctic Ocean (Farrell and others, 2012), and in the analysis of temporal changes in the ocean dynamic topography ... by Envisat in the western Arctic Ocean (Giles and others, 2012). Advantages of our method compared to other studies of campaign biases ... include: (1) smooth surfaces in leads and polynyas that do not require a sea-state bias ... correction, (2) measured laser reflectivity of 0.42 that is closer to the 0.53 reflectivity of the adjacent sea ice and of ice sheets compared to the measured low reflectivity of 0.12 over open ocean, (3) availability of independent Envisat measurements of the vertical motion of the sea surface reference level, and (4) coverage over the reference surface by most of the laser tracks during each campaign.'

'As of December 2012, the ranges for ICESat/GLAS ... ice-sheet data products had been incorrectly calculated from the centroid (amplitude-weighted center of leading and trailing edge thresholds) of the transmit laser pulse to the center of a Gaussian fit of the return pulse (Zwally, 2013). Applying the range correction for the transmit Gaussian to centroid (G-C) offset improved the range precision by 1.7 to <2 cm, and changed (but did not remove) the laser campaign biases (Zwally, 2013). Our current analysis uses elevation data with the G-C correction applied and compatible bias corrections determined with data with the G-C correction also applied. Before the G-C correction was applied, the G-C offset had been in both the data for the ice-sheet $dH/dt$ along-track solutions and in our bias calculations, so the effect of the offsets cancelled. We confirmed that cancellation by comparing our previous and current analyses of $dH/dt$. The average $dH/dt$ for the AIS changed by only +0.01 cm $a^{-1}$, and the average $dH/dt$ error reduced from 0.024 to 0.012 cm $a^{-1}$, reflecting the improved range accuracy. The corresponding $dM/dt$ for the AIS changed by only +1 Gt $a^{-1}$. Therefore, although the net effect of using ice-sheet data without the G-C correction applied is very small if compensurate bias corrections are applied, the error is significant ($-1.29$ cm $a^{-1}$) if the G-C correction is only applied to the data and not to the bias determinations (i.e. incorrectly causing a less positive or more negative $dH/dt$). The error is similar if the G-C correction is applied, but ... [earlier, before G-C corrected] bias adjustments are applied as in Helm and others (2014) in which the volume change obtained from ICESat for 2003-09 for the AIS is [consequently too] negative at $-60 \pm 44$ km$^3$ $a^{-1}$. Helm and others (2014) value of $-23 \pm 36$ km$^3$ $a^{-1}$ (ICESat 2003-09) for EA would adjust to $+109 \pm 36$ km$^3$ $a^{-1}$ if their laser biases had been estimated using data with the G-C correction applied. Scambos and Shuman (2016) also compared an incompatible mixture of biases estimated using data with or without the G-C correction applied.
Importantly, before the G-C error was discovered, the trend in the estimated biases determined without the G-C correction was small, so that applying those bias corrections improved the relative accuracy of the laser campaigns (https://nsidc.org/data/icesat/laser_op periods.html) but made only a small change in trends derived from the data. Specifically, using biases determined over open-water and thin ice in the Arctic Ocean from Zwally and others (2011): ‘We reduce the time variation of these d values [biases] by 0.003 m a⁻¹ to account for the current rate of sea-level rise, and then subtract the reduced d values from the measured elevations. The linear trend in the reduced d is 0.006 m a⁻¹, which averaged over all of Greenland increases the overall mass loss by 9 Gt a⁻¹ compared with data without the d correction applied.’

Shepherd and others (2012) IMBIE-1 included mass gain estimates from ICESat for EA (in Table S8) 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 were done with campaign bias corrections consistently determined. As noted above, trends in the bias corrections were small before the G-C correction, but changed significantly afterward. Shepherd and others (2018) IMBIE-2 did not include ICESat results from Forsberg nor Smith and at least some of the included ICESat results from other investigators (other than Zwally) had laser biases determined with the G-C inconsistency causing a significant dH/dt bias as noted on NSIDC ICESat-data website in 2013.

The bias corrections used in this paper in Table 8 are the same as those in Zwally and others (2015), except for the addition of values for campaigns L2d and L2F in 2009 and the removal of a sinusoidal component with a peak-to-peak amplitude of 4.3 cm and with maxima at day 123 of the annual cycles. These and other bias estimates are available at https://nsidc.org/data/icesat/correction-to-product-surface-elevations.html along with the evaluation criteria such as whether a correction was made for an independently determined vertical motion of the reference surfaces. The NSIDC website includes the recommendation: ‘Applying the per-shot G-C changes, but does not remove all the inter-campaign biases. Any new “campaign level” bias adjustments should be determined with compatible (corrected) data and applied only to analysis of corrected data’.

## Variable radar penetration and backscatter depth

Ice-sheet surface elevations measured by radar altimeters are seriously affected by the strengths of the surface reflection and the sub-surface volume scattering and reflection from internal layers, which were modeled and analyzed in altimeter waveform data over Greenland and Antarctica (Partington and others, 1989). Numerous other papers also addressed the spatial variability of the penetration and its effects on various waveform retracking algorithms, and therefore on the calculated ‘surface’ elevation. In general, altimeter waveforms as depicted in Figures 4–6 in Partington and others (1989) have an initial rise (return vs time) with a slope that is dependent on surface roughness (on the scale of sastrugi) as the pulse-limited footprint expands over the surface, followed by a decreasing return from the radar penetrating into the firn and the consequent volume scattering and reflection from internal ice layers. The three principal waveform-retracking algorithms differ mainly in their points selected on the waveform for the range correction, and therefore differ in the level of their derived surface or near-surface elevation. The threshold tracker (Davis, 1997), which selects a point on the leading edge at 20% of the waveform peak, is least sensitive to sub-surface returns, as is the similar threshold first maximum retracker (TFMRA) (Helms and others, 2014). The multi-parameter waveform fitting tracker (Martin and others, 1983) selects the mid-point of the leading edge corresponding to the mean surface elevation and is also relatively insensitive to volume scattering. In contrast, the Offset-Center-of-Gravity (OCOG) (Bamber, 1994), used by Wingham and others (1998) and by ESA for one of the CryoSat data products, uses the whole waveform and is therefore more sensitive to the sub-surface backscatter and its variability.

While retracking algorithms give different surface or sub-surface elevations, and may have differing accuracies and precisions, those differences would not be a major problem for the measurement of elevation changes if the strengths of the surface reflection and the sub-surface reflections and scattering were constant in time. However, the penetration/reflection depth and the backscatter power are highly-variable seasonally and have multi-year trends, as clearly shown in Figure 3 of Yi and others (2011). Adodo and others (2018) provide a detailed analysis of the seasonal variations of the backscatter-ing over the Antarctic ice sheet including the theoretical dependence on firn

properties and analysis of multi-frequency radar-altimeter measurements made by Envisat and SARAL/AltiKa.

The first elevation correction for the temporal variability of the penetration depth as a function of radar backscatter used the ‘gradient’ of the observed elevation to strength the backscatter derived from the waveforms (note 10 in Wingham and others, 1998). The gradient was called ‘sensitivity’ in Zwally and others (2005), who used the altimeter AGC as a measure of the backscatter and applied other correlation criteria for its application as shown in Figure 6 of Yi and others (2011), thereby improving the correction. Yi and others (2011) also considered alternate methods (short-term, mixed-term and long-term) of calculating the sensitivity that give different sensitivities and correlation coefficients. Successful corrections for ERS1/ERS2 were also made by Davis and Ferguson (2004) and Khvorostovsky (2012).

Unfortunately for Envisat and CryoSat data, the correction for the time-variable penetration depth became substantially more difficult. The linearly-polarized radar signals, which were oriented across-track on Envisat at 120° and CryoSat at 90°, interact with firn properties related to the direction of the surface slope (sometimes called surface anisotropy) and the relative directions (polarization vs surface slope) differ significantly at track crossings (e.g. Legresy and others, 1999; Arthern and others, 2001; Rémy and others, 2012). In contrast, the orientation of the polarization along-track (at 0°) on ERS1/ERS2 tended to be more oriented in the direction of maximum surface slope at high-latitude crossovers rather than across-slope, especially at the steeper ice-sheet margins, which may have enabled the more successful penetration corrections for ERS crossover analysis.

For Envisat data, a successful correction was developed using repeat-track analysis and an advanced correction algorithm (Flament and Rémy, 2012). Repeat-track analysis significantly mitigates the variable penetration problem, because the polarization orientation relative to the surface slope is essentially identical on the repeating tracks. A critical point is that their solution makes a time-dependent backscatter correction for the variable depth penetration, and also uses time-variable waveform parameters. They used 84 of the 35-day repeat cycles from September 2002 to October 2010 and computed ‘the elevation trend every kilometer along-track’ using ‘All available measurements within a 500 m radius of a point on the mean ground track’. ‘... In the central part of the East Antarctica, the height and the leading edge width fluctuations vary together while elsewhere, height fluctuations may occur with no variations in the waveform shape, mostly during winter. As a consequence, these induced errors cannot be corrected with solely the help of the backscatter: waveform shape parameters are also needed. They are however not enough to fully correct these two errors. We propose an empirical correction for these effects. ... In terms of volume change, the estimation may vary up to 4 cm a⁻¹ at cross-overs depending on the correction used and is reduced in average to 2.3 cm a⁻¹ with our correction. The difference between the height trends estimated with both corrections is weak in average but may locally reach 5 cm a⁻¹ with a clear geographical pattern.’

Consistency of the dH/dt from Envisat radar altimetry after correction for the variable-radar penetration is shown in the comparisons in Figures 9, 10 and S1 with the corrected dH/dt from ERS1/ERS2 and the dH/dt from the ICESat laser altimetry. The three average dH/dt over EA from ERS1/ERS2, ICESat and Envisat are 10.7, 13.1 and 8.3 mm a⁻¹ showing agreement at the level of only a few mm a⁻¹.

The method of McMillan and others (2014) for CryoSat is: ‘To compute changes in ... elevation, we adapted a repeat-track method (Smith and others, 2009; Moholdt and others, 2010; Flament and Rémy, 2012) to suit the Cryosat-2 dataset, ...’ and ‘...Elevation measurements are accumulated in 469 451 regularly spaced (5 by 5 km) geographical regions, and within each region, we solve, simultaneously, for spatial and temporal fluctuations in elevation and for a fixed contribution due to the impact of surface anisotropy on the tracked range (see supporting information).... ...and a correction is applied to account for temporal fluctuations in backscatter that cause spurious fluctuations in range (Davis and Ferguson, 2004; Khvorostovsky, 2012; Wingham and others, 1998)’. Their solution is complicated because: (1) their 5 km covers a 100× larger area with more variable surface conditions than that used by Flament and Rémy (2012) and the long 365-day near-repeat cycle includes few near-repeat orbits, (2) their ‘contribution’ for the ‘impact of surface anisotropy’ is very large (+1 to -1 m in their Supplementary Material Figs 1 and 3); their separation into fixed and time-varying fluctuations is of dubious validity. Their range measurements are ‘corrected for the lag of the leading edge tracker’ (Wingham and others, 2006), which used ERS ‘WAP v. 3 altimeter data’ and presumably the OCOG retracker that is more sensitive to sub-surface penetration.
In contrast, Helm and others (2014) stated: ‘… our study show(s) that a correction for the static “Antarctic pattern” in $dH/dt$ estimates as applied in McMillan and others (2014) (for penetration) can be avoided when using the TFMRA re-tracker’. Table 4 in Helm and others (2014) for EA shows volume changes of +78 ± 19 km³ a⁻¹ (IMBIE 2003–2008) and +59 ± 63 km³ a⁻¹ (CryoSat 2011–14), compared to the -2.7 ± 33 km³ a⁻¹ (CryoSat 2010–2013) from McMillan and others (2014), giving a 62 ± 71 km³ a⁻¹ difference between CryoSat investigators for EA.

For Greenland, Nilsson and others (2016) showed that an improved leading-edge retracker for CryoSat-2, which changes the sensitivity to depth penetration, can cause a very-large 50 cm a⁻¹ difference in the derived surface elevation in the normally dry snow zone of Northern Greenland and significant differences in the volume change estimates compared to ESA’s public data product.

## Deriving mass changes from elevation changes

Our methods of deriving mass changes, as applied to Greenland (Zwally and others, 2011) and to Antarctica (Zwally and others, 2015) and followed in this paper for the $dM/dt$, $dM/dt$ and $Ma/dt$ in Figures 16 and S2, have distinct advantages not employed in other studies. The advantages are: (1) correction for accumulation-driven and temperature-driven changes in surface elevation that do not involve changes in mass using a state-of-art FC model (Li and Zwally, 2015); and (2) separation of accumulation-driven and dynamic-driven mass changes and the assignment of proper ice ($\rho_i$) and near-surface firn ($\rho_a$) densities to each, even though $\rho_a$ is not necessarily calculated (see text following Eqn (13)).

Initially, investigators used a single density $\rho$ to estimate $dM/dt = \rho \times dH/dt$ (with $dH/dt$ corrected for bedrock motion and perhaps FC), even though it was known that elevation changes were likely due to a combination of accumulation-driven changes with a density of $\rho_a$ and dynamic-driven changes with the density of $\rho_i$. For example, Zwally and others (2005) calculated a mass change $dF/dt$ using $\rho_a = 0.4$, which ‘is a typical mean density for the top strata corresponding to 10 years of accumulation’, and $dM/dt$ using $\rho_i = 0.91$, which provided their preferred estimate. Clearly, choosing either $\rho_a$ or $\rho_i$ makes a factor of 2.3 or more difference causing significant errors in mass estimates one way or the other.

More recently, users of the old method (e.g. McMillan and others, 2014; McMillan and others, 2016; Martín-Español and others, 2017; Schroder and others, 2019) take $dM/dt$ to be equal to $\rho_{\text{firn/ice}} \times dH/dt$, where H is corrected for bedrock motion and perhaps FC, and $\rho_{\text{firn/ice}}$ is chosen/assumed to be either $\rho_{\text{firn}}$ equal to ~0.350 or $\rho_{\text{ice}}$ equal to 0.917, sometimes based on a limited spatial mask as in McMillan and others (2014) and Schroder and others (2019). From McMillan and others (2016): “To convert the resulting altimeter rates of change to mass, we constructed a density model that accounted for both surface and dynamic processes. In regions where high rates of elevation change and ice flow suggested a state of dynamic imbalance, we used an ice density of 917 kg m⁻³ (see Text S8). Elsewhere, detected elevation changes were assumed to be driven by SMB processes, and we used an ice density within the ablation zone and the density of the IMAU-FDM firn layers gained or lost across the remaining areas’, for which use of the density of firn layers instead of their former 350 kg m⁻³ made a small improvement. However, the method maintains the critical flaw of not actually accounting ‘for both surface and dynamic processes’ where surface and dynamic processes occur in the same location, which is mostly everywhere in the accumulation zone.

As we noted following Eqn (16), ‘a priori selection of appropriate single or multiple firn/ice densities … is not possible due to the extensive spatial and temporal variabilities of the actual $\rho_a$, and because $H_a$ and $H_d$ have differing spatial variations in magnitude and sign’. This is further illustrated for Greenland in Zwally and others (2011) in their Figure 7 ‘Maps for the 2003–07 period. (a) Accumulation-driven elevation change, $dH_{\text{CA}}^a/dt$. (b) Ablation- and dynamic-driven elevation change, $dH_{\text{bd}}/dt$. (c) Relative density, $\rho_a$, of the firn for the $dH_{\text{CA}}^a/dt$ component’. Their Figure 7b clearly shows the extensive area of dynamic thickening over much of the higher elevations of the accumulation zone, and in their Figures 7a and b, the mixture of surface and dynamic processes everywhere. The large variability of the density for the surface processes is shown in their Figure 7c. Furthermore, the surface processes (i.e. $H_{\text{CA}}^a(t)$) are more variable with time on decadal and sub-decadal time scales, and therefore vary in sign from the more constant dynamic processes, both of which contribute to the measured $H(t)$ according to Eqn (12).

Similarly for Antarctica, the large spatial and temporal variations of the accumulation-driven mass change, $dM_a/dt$, are shown in Zwally and others (2015) in their Figure 10a for 1992–2001 and 10b for 2003–2008, and are

also evident in the measured $dH/dt$ in their Figures 6a and b. In contrast, the minimal temporal variations of the dynamic-driven changes are shown in their Figures 11a and b, with the exception of the increases in dynamic thinning in WA1. For the ICESat period, the large spatial variability of the $dM_a/dt$ is also shown in our Figures 16c and S2c, compared to the mostly small spatial variations in the dynamic thickening in EA and the large variations in dynamic thinning in WA1 and thickening in WA2 shown in Figures 16b and S2b.

The difficulty of choosing a correct density for the firn changes is further illustrated by the calculated spatial distributions of $\rho_a = \Delta M_a / \Delta(H^a - C_A)$ in Figure 17 for 1992–2001 and 2003–2008. The $\rho_a$ represent the firn distributed over a range of depths depending on the time history of the accumulation anomalies as they propagate into the firn, and do not represent the density of a particular firn layer at a specific depth. The regional average $\rho_a$ are listed in Table 9, adapted from Table 4 in Zwally and others (2015). Also in Table 9 are the $\rho_{\text{pseudol}} \equiv dM/dt/(dI/dt \times \text{Area})$ using the derived $dM/dt$ and $dI/dt$, which is the rate of ice thickness change corrected for temperature-driven FC and bedrock motion (i.e. $dI/dt \equiv dH/dt - dC_T/dt - dB/dt$). The range of $\rho_{\text{pseudol}}$ from 0.55 to 5.78, with 12 out of 16 values outside the range of 0.2–0.92 firn/ice densities, demonstrates the impossibility of selecting a single value of $\rho_{\text{firn/ice}}$ to calculate correct mass changes. This result is relevant to the critiques (Martín-Español and others, 2017; Bamber and others, 2018) of our results that are at least partially based on the premise that a single density can be used to derive accurate mass changes from elevation changes.

Finally, we note that although many altimeter studies use some form of FC modeling in their analysis, there are major differences in the validity of the models and their specific applications to altimeter data. Furthermore, quantitative evaluation of those differences is typically not possible because of the lack of details provided in various papers such as the time series of the modeled compaction parameters $C_A(t)$ and $C_T(t)$, for example, as we show combined as $C_{\text{AT}}(t)$ in Figures 5 and 7. Although the FC models mostly have a common heritage based on the semi-empirical formulation of Heron and Langway (1998), which as used in Zwally and Li (2002) included the important innovation of a greater sensitivity of the compaction rate to firn temperature based on laboratory measurements of ice creep. However, several differing temperature sensitivities have been used by other investigators giving differing temperature-driven trends in elevation.

A critically important advance not used in other FC models is the time-dependent formulation of the compaction equations on the accumulation rate $A(t)$, which was first introduced in Li and Zwally (2011) and in Eqn (9) in Li and Zwally (2015). For example, in the often-used model of Ligtenberg and others (2011), the accumulation rate appears as a constant in their Eqns (5), (8) and (9), as it was initially in Heron and Langway (1988). As detailed in Li and Zwally (2015), the time-dependent treatment of the $A(t)$ is essential for determining the proper time response of the firn to accumulation variations and for calculating the resulting accumulation-driven trends in surface elevation. Proper time-dependence of the FC modeling is critically important because the rate of FC and the consequent rate of change of the surface elevation at any given time for correction of the measured $dH/dt$ depend on the time history of both accumulation and temperature for decades (Li and Zwally, 2015) prior to the measurement.

The accumulation and temperature datasets chosen to drive the FC models are also very important and contribute to significant differences. In Zwally and others (2015), we justified and used the ERA-Interim re-analysis data on accumulation rates, $A(t)$, instead of other models partially based on the more realistic spatial distribution of the temporal variability, particularly in coastal regions. Further support was provided by a detailed analysis (Medley and others, 2013) of the spatial and temporal correlations from 1980 through 2009 in WA between $A(t)$ derived from layering shown by an airborne snow radar. Correlations among (1) four re-analyses (including ERA-Interim and RACMO) and (2) ice cores gave a temporal correlation for ERA-Interim of 0.93 compared to only 0.68 for RACMO, 0.91 and 0.92 for the other two re-analyses, and 0.80 for the ice cores. Also, we believe our use of the satellite AVHRR-measured temperature is preferred to modeled temperatures used by others because the trends in the modeled temperature vary widely among models and differ significantly from the measured temperatures.

After long-term FC model spinup with a constant mean $A$, it is extremely important to drive the models with the variability in accumulation variations ($\delta A(t) = A(t) - \langle A(t) \rangle_{>27}$) with respect to the long-term (e.g. 27-year model mean) rather than with $A(t)$ for two reasons. First, the $\delta A(t)$ are mostly more accurate than the model mean ($\langle A(t) \rangle_{>27}$), and second it avoids a
discontinuity in the model compaction formulation caused by a change from
the spinup A to the model mean. The second reason occurs because as the
modeled mean accumulation replaces the spinup mean, starting at the surface
and propagating downward with time, the replacement introduces an artificial

trend in the modeled surface *H*(*t*) of several cm a⁻¹, thereby obscuring or
falsely indicating an elevation trend of several cm a⁻¹. Proper demonstration
of this effect requires a time-dependent formulation in the FC model as dis-
cussed above.