1. Introduction

The principal processes affecting the mass balance and dynamics of the ice sheets are illustrated in Figure 1, including interactions with the atmosphere and ocean. The principal ice-mass input is from atmospheric precipitation in the form of snowfall, with subsequent losses from sublimation, snow redistribution by surface winds, and drift removal at the margins. In Greenland, surface melting in the ablation zone causes water runoff from the grounded ice and acceleration of the ice flow as meltwater propagates to the ice-sheet base to enhance basal sliding. In contrast, surface melting on the grounded ice of Antarctica is very small, and subject to refreezing in the firn, with negligible mass loss. Interaction with the ocean occurs at the undersides of the floating ice shelves and glacier tongues, with ocean melting of the basal ice or accretion of ice from basal freezing. Consequent changes in the thickness of floating ice affect the rate of ice flow from the grounded ice.

Fig. 1. The principal processes affecting the mass balance and dynamics of the ice sheets are ice mass input from snowfall with losses from sublimation and drifting. Surface melting on the grounded ice of Antarctica is very small, and subject to refreezing in the firn. Interaction with the ocean occurs at the undersides of the floating ice shelves and glacier tongues, and consequent changes in thickness affect the rate of ice flow from the grounded ice.

The annual mass input to the Antarctic ice sheet (AIS) from snow accumulation is ∼2000 Gt a⁻¹, based on field data (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999; Reference Giovinetto and ZwallyGiovinetto and Zwally, 2000) and meteorological data (Reference Lenaerts, Van den Broeke, Van de Berg, Van Meijgaard and Kuipers MunnekeLenaerts and others, 2012). By 2011, estimates of the mass balance (input–output) of Antarctica had improved from roughly 0 ± 400 Gt a⁻¹ in the mid-1990s to a range of approximately +50 to −250 Gt a⁻¹ for the period 1992–2009 (Reference Zwally and GiovinettoZwally and Giovinetto, 2011). The principal techniques have been based on satellite technology: (1) elevation and volume change from radar altimetry (RA) and laser altimetry (LA), (2) mass change from gravimetry (GR), and (3) the input–output method (IOM) using ice-velocity sensing for mass output estimates. Subsequent evaluation (Reference ShepherdShepherd and others, 2012) eliminated some of the larger estimates of Antarctic mass loss and gave a mean estimate of −72 ± 43 Gt a⁻¹ from the four techniques (RA, LA, GR and IOM) for the intercomparison period October 2003 to December 2008. Despite the progress, significant uncertainties remained about rates of mass change for various reasons (Reference HannaHanna and others, 2013), especially for East Antarctica (EA). Despite residual uncertainties, the general pattern of gains and losses has been a net mass loss from coastal regions of West Antarctica (WA) and the Antarctic Peninsula (AP), acceleration of those losses during the last 20 years, and a net mass gain in EA (Reference Zwally and GiovinettoZwally and Giovinetto, 2011; Reference ShepherdShepherd and others, 2012; Reference HannaHanna and others, 2013).

Variations in the surface mass balance (SMB) and vertical ice velocity of the grounded ice sheet both cause changes in net mass balance (dM/dt) and surface elevations (dH/dt). Over Antarctica, the rate of total mass change per unit area is principally the sum of surface balance, dM _a/dt, and dynamic-driven mass change, dM _d/dt:

(1)

The dynamic dM _d/dt is defined by differences between the vertical ice flux from the surface and the long-term (>decades) average accumulation rate 〈A〉. dM _a/dt is driven by short-term (<decades) variations in the accumulation rate. In our analysis, any net changes in ice thickness caused by melting or freezing at the base of grounded ice are neglected and are therefore essentially included in the dynamic term. The dM _d/dt may be caused by either dynamic variations in ice velocity on various timescales or by persistent long-term changes in 〈A〉 with a consequent long-term dynamic response.

The essential characteristic of both dynamic thickening (dH _d/dt > 0 and dM _d/dt > 0) (Reference Li and ZwallyZwally and others, 2011) and dynamic thinning (dH _d/dt < 0 and dM _d/dt < 0 ) (Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009; Reference Filament and RémyFilament and Rémy, 2012) is a difference between the vertical ice flux from the surface and 〈A〉. Information on the timing of changes in dynamic thinning or thickening is only provided if observed dH _d/dt differ between measurement periods. Interpretation of changes in dH _d/dt on longer timescales (>decadal) requires additional information such as variations in accumulation rates from ice cores.

In this paper, we determine the values of the mass terms in Eqn (1) in 50 km × 50 km gridcells over the AIS using satellite-altimeter measurements of dH/dt along with meteorological data on accumulation variations. We use radar-altimeter measurements of dH/dt by European Remote-sensing Satellites 1 and 2 (ERS-1/-2) for the period 1992–2001 and satellite laser-altimeter measurements by the Ice, Cloud and land Elevation Satellite (ICESat) for 2003–08. Our procedures (Reference Li and ZwallyLi and Zwally, 2011; Reference ZwallyZwally and others, 2011) for deriving dM/dt from measured dH/dt apply height corrections for factors that do not change the ice mass, specifically vertical motion of the underlying bedrock (dB/dt) and changes in the vertical velocity (V _fc) of the surface due to variations in the rate of firn compaction (FC). The dM _a/dt during elevation measurement periods are calculated directly from the accumulation variations, δA(t), provided by meteorological reanalysis data relative to their 27 year average (1982–2008).

We analyze our derived dM/dt, and its dM _a/dt and dM _d/dt components, in 27 drainage systems (DS) (http://icesat4.gsfc.nasa.gov/cryo_data/ant_grn_drainage_systems.php), and the major regions of WA, EA and the AP as shown in Figure 2. Grounded ice on adjacent islands and ice rises within ice shelves is included in the respective DS. WA is divided into WA1 (consisting of three dynamically active coastal DS that include Pine Island, Thwaites and Smith glaciers and the coastal DS of Marie Byrd Land) and WA2 (consisting of three inland DS flowing into the Ross and Filchner Ice Shelves and a small coastal DS at the base of the AP). EA is also divided into the westernmost EA1 consisting of DS2–11 and the easternmost EA2 consisting of DS12–17, in order to illustrate differences in their meteorological variations. We discuss the glaciological and climatological significance of the derived rates of mass change, including significant temporal changes between the 1992–2001 and 2003–08 measurement periods and the absence of significant temporal changes in certain DS and regions.

Fig. 2. Antarctic drainage systems and regions. The Antarctic Peninsula (AP) with DS24–27. West Antarctica (WA) is divided into WA1 (Pine Island Glacier DS22, Thwaites and Smith Glaciers DS21 and the coastal DS20) and WA2 (coastal DS23 and inland DS1, DS18 and DS19). East Antarctica (EA) is divided into EA1 (DS2–11) and EA2 (DS12–17).

2. Altimeter Data, Elevation Change Solutions and dH/dT Maps

For 1992–2001, we use previous gridded dH/dt (Reference ZwallyZwally and others, 2005) that were calculated from time series of elevation differences at orbital crossovers of ERS-1 and ERS-2 radar altimeter data (http://icesat4.gsfc.nasa.gov/) from mid-April 1992 to mid-April 2001. The instrument, range retracking, ERS-1/ERS-2 bias and other corrections, as well as the use of DUT DGM-E04 orbits which have a radial orbit precision of 5–6 cm (Reference Scharroo and VisserScharroo and Visser, 1998), are described in Zwally and others (Reference Zwally2005; http://icesat4.gsfc. nasa.gov/). Principal changes in our analysis for the ERS period are the separation of accumulation-driven and dynamic-driven changes and the dB/dt correction that uses a more recent model on bedrock motion (Section 3).

For 2003–08, we use ICESat data release 634 (http://nsidc.org/data/icesat/index.html) for 15 laser campaigns (each with data for 33 to ∼36 days) from October 2003 to October/mid-December 2008. The surface elevations in release 634 data products were changed from release 633 by application of the Gaussian–centroid (G-C) range correction (see Appendix), which also significantly increased the need for laser campaign bias corrections to obtain accurate elevation changes (Reference ZwallyZwally, 2013). We apply ICESat laser campaign bias corrections derived from laser measurements of the sea surface height (SSH) over open water and thin ice in leads and polynyas within the Antarctic and Arctic sea-ice packs, with adjustments for SSH variations measured by Envisat radar altimetry concurrently with ICESat measurements (see Appendix). We calculate dh/dt at equally spaced (172 m) reference points on a set of ICESat reference tracks for points where sufficient data (i.e. N ≥4 satellite passes) are available, solving simultaneously by linear least-squares minimization for dh/dt, cross-track surface slope (α) and elevation (h ^trk(t ₀)) on the reference track at time t ₀ (Reference ZwallyZwally and others, 2011). Data on parallel repeat passes are first interpolated to the locations on the tracks perpendicular to each reference point. Solutions at each reference point also include data at three points before and after the reference point. We account for the along-track slope by first calculating it separately from a fit of the data from all tracks to a surface curved in the along-track direction. To extend the analysis beyond the ICESat campaigns through 2007 that have cloud flags in the data records previously used for editing, we developed multi-stage data-editing criteria that gave essentially the same dM/dt results (i.e. −173 Gt a⁻¹ versus −171 Gt a⁻¹ for 2003–07 for Greenland in Reference ZwallyZwally and others, 2011). In the first stage, we use all ICESat data for which valid surface elevations are indicated and iteratively exclude 3σ outliers from the curved surface in the least-squares solution. The maps of dh/dt and dh/dt over grounded ice are presented in Figure 3. The second stage excludes dh/dt outliers in the averaging of dH/dt into gridcell maps of dH/dt.

Fig. 3. Along-track solutions at 172 m spacing in Antarctica from ICESat data. (a) dh/dt, and (b) σ _dh/dt showing σ _dh/dt is mostly <0.05 m a⁻¹, with larger values at margins where dh/dt tend to be larger.

The calculated along-track profiles of dh/dt and Δh(t_i ) over Smith Glacier in DS21 in Figure 4 illustrate the quality of the solutions in a region of large thinning near the coast in WA. The maximum rate of surface lowering on track 1300 over Smith Glacier is 9 m a⁻¹, resulting in a 48 m lowering in 5.5 years. Figure 5 illustrates the variations of dh/dt in a region of small changes over Vostok Subglacial Lake in EA, where dh/dt varies from approximately +0.0 ± 1 cm a⁻¹ to +3.0 ± 1 cm a⁻¹ with km-scale oscillations and significant spatial variability of average values over the lake. Other features of the dh/dt variability are low dh/dt in the trough on the southwestern end of track 0330, in contrast to the average value there on track 1312, and the above average on track 0077. Near the small ridge on the northeastern side, dh/dt is approximately +2 cm a⁻¹ on track 0330, +1.5 cm a⁻¹ on track 1312, and +0.5 cm a⁻¹ on track 0077.

Fig. 4. Along-track profiles in West Antarctica from ICESat data. (a) Profiles of derived h, dh/dt and σ _dh/dt from along-track solutions at 172 m spacing along track 1300 across Smith Glacier in DS21. Maximum dh/dt is −8 m a⁻¹ near the center of the glacier. (b) Surface elevation (Δh(t _i)) profiles (ten-point smoothing) relative to the first profile from October 2003 to April 2009. Maximum surface lowering is 48 m in 5.5 years. (c) Along-track dh/dt for the 33 day repeat-cycle tracks plotted on ICESat elevation map.

Fig. 5. Along-track profiles over Vostok Subglacial Lake from ICESat data. (a) Profiles of derived dh/dt, σ _dh/dt and h from along-track solutions at 172 m spacing along six tracks. The dh/dt varies from <0 cm a⁻¹ near the northern end of track 0071 to >3.0 cm a⁻¹ on the southern end of track 0077. The σ _dh/dt are mostly 1.0 cm a⁻¹. The vertical blue lines indicate crossings of the lake boundary (Reference Tikku, Bell, Studinger and ClarkeTikku and others, 2004). (b) dh/dt for the 33 day repeat-cycle tracks plotted on ICESat elevation map showing outline of the lake.

The spatial variability of dh/dt may be caused by spatial variability of several cm in local values of A(t) from variations in snowdrift, precipitation and sublimation, which also affect the temporal variability of A and elevations at specific locations. Therefore, comparisons between dh/dt from altimetry and surface-based GPS estimates (see Appendix) are highly dependent on the specific locations. A characteristic feature illustrated in Figure 5 is the continuity of dh/dt values on the grounded ice surrounding the lake and the values on the floating ice on the lake. Another interesting feature is the 20 km oscillation in dh/dt between approximately +5 and −2 cm a⁻¹ on the southern end of track 0190, indicating downslope migration of the snow dunes shown on the elevation profiles.

Average values (dH/dt) of the dh/dt in each 50 km gridcell are calculated as previously (Reference ZwallyZwally and others, 2011), and the resulting dH/dt maps for 1992–2001 and 2003–08 are presented in Figure 6. We use kriging (optimal interpolation) to fill the dH/dt gridcells without data south of the 86°S coverage of ICESat, and for some cells at the ice-sheet margins for both ERS and ICESat. For 1992–2001, we use the ICESat dH/dt to fill in the area south of the 81.5°S coverage of ERS. The average dH/dt of 1.96 cm a⁻¹ south of 81.5°S is therefore the same for both ERS and ICESat, but their respective dM/dt values, of 56.0 and 45.9 Gt a⁻¹, differ because dM _a/dt, dH ^a/dt and dC _A/dt (Section 3) differ between periods. The results of using the ICESat dH/dt south of 81.5°S in the ERS dH/dt grid compare with the alternative method of kriging, which gave an average dH/dt of 2.35 cm a⁻¹ and dM/dt of 64.1 Gt a⁻¹. However, the distribution among the DS south of 81.5°S and the division between WA and EA are better using the ICESat dH/dt.

Fig. 6. Maps of dH/dt (a) for 1992–2001 from ERS-1 and -2 data (dH/dt from ICESat used south of 81.5°E) and (b) for 2003–08 from ICESat data.

In area-weighted summations for DS and regions, cells partially on grounded ice are weighted accordingly. Cells on DS boundaries are partially assigned in proportion to their coverage in each DS. For ICESat data, the partial-cell coverage of the along-track data at the ice margin is determined by an ice-edge boundary with 1 km resolution. For ERS, dH/dt values in the partial cells (as derived from crossover analysis) are weighted by the area fraction on grounded ice, also determined by the 1 km boundary. For ERS, 152.4 out of 4810.9 cells (counting partial cells as fractions) at the ice-sheet margins are kriged and their average dH/dt is −4.9 cm a⁻¹ compared to +0.50 cm a⁻¹ for all cells. For ICESat, 1.6 out of 4810.9 cells at the ice-sheet margins are kriged and their average dH/dt is −20.3 cm a⁻¹ compared to +0.76 cm a⁻¹ for all cells. Previously, cells for ERS were taken as either totally on-ice or off-ice (Reference ZwallyZwally and others, 2005).

We examine the relative calibration of the ERS radar and ICESat laser altimetry over Vostok Subglacial Lake in central EA, which is a dynamically stable region with low accumulation and small vertical velocities. In this region, short-term dynamic changes that might affect the derived H(t) should be very small. Vostok Subglacial Lake also has technical advantages, because of its flatness and consequent minimization of effects of surface slopes on both the radar and laser altimeter measurements. We use ERS crossover analysis and ICESat along-track solutions within a box (78.6–76.19°S, 101.3–107.0°E) of approximately 63 km × 260 km over Vostok Subglacial Lake.

The resulting H(t) time series from ERS, both before and after the backscatter correction for the variable penetration depth of the radar signal in the firn, and the H(t) time series from ICESat, along with their linear trends are shown in Figure 7. As shown in Table 1, the respective trends for ERS and ICESat are very similar, at 2.03 and 2.02 cm a⁻¹. After correction for dB/dt and dC _T/dt (defined in Section 3), the respective dI/dt are 2.01 and 2.08 cm a⁻¹. After further subtracting the accumulation-driven change , dH _d/dt are 2.18 and 2.05 cm a⁻¹. The small differences (≤0.13 cm a⁻¹) in these height parameters are a good indication of the close relative calibration of ERS and ICESat, for which we estimate respective systematic calibration errors of 0.20 and 0.15 cm a⁻¹. In addition, over the whole of EA (Table 2), the average dH/dt from ERS and ICESat are 1.11 and 1.30 cm a⁻¹, the dI/dt are 1.19 and 1.44 cm a⁻¹ and the dH _d/dt are 1.58 and 1.59 cm a⁻¹, indicating minimal short-term dynamic change as discussed in Section 5.

Fig. 7. H(t) time series on Vostok Subglacial Lake. From ERS 1992–2003 with trend of +2.03 cm a⁻¹ after backscatter correction (red) and +2.18 cm a⁻¹ before backscatter correction (black). From ICESat 2003–08 with trend of +2.02 cm a⁻¹ (blue). The backscatter correction significantly reduces the amplitude of the seasonal variability in the ERS signal.

Table 1. Components of elevation change (cm a⁻¹) over Vostok Subglacial Lake

Table 2. Components of surface elevation change (cm a⁻¹), dH/dt, dC _T/dt, dI/dt, , dH _d/dt, dB/dt by region

3. Firn Compaction Modeling, Bedrock Motion, and Calculation of Mass Changes

Calculation of dM/dt mass changes from measured dH/dt requires correction for the dH/dt changes that do not involve changes in ice mass, specifically changes in the vertical velocity (V _fc) of the surface due to variations in the rate of FC and motion of the underlying bedrock, dB/dt. Defined as perturbations from steady state, the vertical components of dH/dt are

(2)

where dH ^a/dt is the direct height change driven by accumulation variations (δA(t) = A(t) − 〈A〉), dC _A/dt and dC _T/dt are changes in V _fc(t) driven by δA(t) and temperature variations (T(t)), dH _d/dt = −(V _ice(t) −〈V _ice〉) is the dynamic height change determined by physical changes in velocity relative to the long-term (>decades) average velocity (〈V _ice〉=〈A〉/ρ _i) at the firn/ice transition, where ρ _i is the relative density of ice, and dB/dt is taken to be time-independent over decades.

The dynamic-driven height change is then

(3)

where dH/dt is the measured height, and dC _T/dt and dC _A/dt are calculated with the FC model driven by δA(t) from meteorological reanalysis data and by T(t) from satellite measurements, as described below. The dH ^a/dt is also calculated from δA(t) using a near-surface relative density of 0.3. The dB/dt is from a recent model (Reference Ivins, James, Wahr, Schrama, Landerer and SimonIvins and others, 2013) of glacial isostatic adjustment. Combining terms in Eqn (3) gives

(4)

where dI/dt = dH/dt − dC _T/dt − dB/dt is an effective rate of ice thickness change taking into account the bedrock motion and the FC changes driven by temperature. combines direct height change and FC change which are both driven by accumulation variations.

The density associated with the dynamic-driven height changes is ρ _i and the dynamic mass change for Eqn (1) is therefore

(5)

The accumulation-driven mass change is

(6)

where τ is the time period of the measurements. To clarify the timescales involved, dM _a/dt is determined only by δA(t) during τ, and not by prior anomalies or variations in V _fc. Secondly, variations in V _fc, which depend on the time history of T(t) and δA(t) before and during the measurements, only affect dC _T/dt and dC _A/dt. The values of dC _T/dt and dC _A/dt affect the apportionment of dI/dt = dH/dt − dC _T/dt − dB/dt between and dH _d/dt according to Eqn (4). Effects of δA(t) and T(t) variations before 1982 on the V _fc during 1992–2008 are shown to be limited to minimal values by the response times of the firn to A(t) and T(t) perturbations (Reference Li and ZwallyLi and Zwally, 2015).

Our FC model is initiated with steady-state density profiles calculated for each gridcell (approximately 50 km 50 km) for thousands of years to convergence using estimates of long-term accumulation rate 〈A〉 and mean annual temperature 〈T〉. We use 〈A〉 modified from a compilation and interpolation of field data (Reference Giovinetto and ZwallyGiovinetto and Zwally, 2000) and 〈T〉 as the average for 1982–84 of measured T(t) from satellite Advanced Very High Resolution Radiometer (AVHRR) measurements (Reference Li, Zwally and ComisoLi and others, 2007). From 1982 onward, the model is driven by monthly values of the measured T(t) and δA(t) from meteorological data. We use δA(t) = A(t) − 〈A(t)〉₂₇, where A(t) is precipitation minus sublimation (P − S) from the European Centre for Medium-Range Weather Forecasts’ ERA-Interim atmospheric-model reanalysis (Reference DeeDee and others, 2011) and 〈A(t)〉₂₇ is the 27 year average (1982–2008).

Most important for our purposes is the accuracy xof the temporal variability of δA(t). A comparison (Reference Bromwich and NicolasBromwich and Nicolas, 2011) of five reanalyses concluded that the ERA-Interim data provide the most realistic depiction of interannual variability and trends in Antarctic P − S during 1989–2009. Comparison of the time series of A(t) from the ERA-Interim data with the SMB from a regional atmospheric climate model (RACMO) (Reference Lenaerts, Van den Broeke, Van de Berg, Van Meijgaard and Kuipers MunnekeLenaerts and others, 2012) shows an overall correlation of 0.79 over Antarctic grounded ice for 1979–2010, indicating similar overall temporal variability. However, comparison of the distributions of the dM _a/dt we obtain using ERA-Interim (Section 4; Fig. 10, further below) with the distributions we obtain using RACMO shows significant differences in the spatial distribution of the temporal variability, particularly in coastal regions. We also find that both datasets give dM _a/dt losses in EA for both measurement periods (−6 and −33 Gt a⁻¹ for RACMO and −11 and −11 Gt a⁻¹ for ERA-Interim). However, the large difference between periods suggested by the RACMO data gives a corresponding large difference in the dynamic thickening in EA, which we believe is not realistic.

Further support for our use of ERA-Interim is provided by a detailed analysis (Reference MedleyMedley 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 and (1) four reanalyses (including ERA-Interim and RACMO) and (2) ice cores. The snow radar flight lines covered ∼90 000 km² of the eastern part of the 217 404 km² of our DS21 extending slightly into DS22. The temporal correlation for ERA-Interim from 1980 through 2009 was 0.93 compared to only 0.68 for RACMO, 0.91 and 0.92 for the other two reanalyses and 0.80 for the ice cores. Furthermore, the time series of the snow radar data (fig. 3 of Reference MedleyMedley and others, 2013) shows a negative anomaly for 1992–2001 of −0.008 m w.e. a⁻¹ versus a positive anomaly for 2003–08 of +0.024 m w.e. a⁻¹. Assuming those anomalies extended over the entire DS21 gives net values of −1.6 and +5.1 Gt a⁻¹ for the two periods, highly consistent with our respective dM _a/dt of −1 ± 1 and +4 ± 2 Gt a⁻¹ for DS21 using the ERA-Interim data. In contrast, our values for DS21 using RACMO are less satisfactory at −1 and +12 Gt a⁻¹.

Regional values of the dH/dt, dC _T/dt, dI/dt, and dH _d/dt components of elevation change are shown in Table 2. Mean values over the AIS for the two measurement periods are dH/dt = 0.50 and 0.76 cm a⁻¹; dC _T/dt = −0.30 and −0.27 cm a⁻¹; dI/dt = 0.71 and 0.95 cm a⁻¹; and 0.27 cm a⁻¹; and dH _d/dt = 1.20 and 0.69 cm a⁻¹. For EA, the accumulation-driven of −0.40 and −0.16 cm a⁻¹ are small relative to the dynamic dH _d/dt = 1.58 and 1.59 cm a⁻¹. The are also small relative to dH _d/dt in WA1 and the AP where dynamic losses are large.

In EA, where dH/dt are 1.18 and 1.30 cm a⁻¹ for the two periods, dC _T/dt temperature-induced surface lowerings of −0.13 and −0.19 cm a⁻¹ represent ∼10% corrections to the measured dH/dt. Applying those dC _T/dt corrections to the measured dH/dt adjusts the derived dM/dt for EA by +12.1 and +17.6 Gt a⁻¹ for the two periods as shown in Table 3. For WA, the respective temperature-induced corrections of 19.3 and 9.2 Gt a⁻¹ are comparable to EA, because the temperature increases in WA are larger although their area is only 18% as large. For the AIS, the temperature-induced corrections are 33.8 and 30.4 Gt a⁻¹, which emphasizes the importance of FC corrections for obtaining accurate mass-balance estimates during periods of climate warming (or cooling). As previously shown for Greenland (Reference Zwally and GiovinettoZwally and others, 2011), a 4.1 cm a⁻¹ surface lowering induced by higher temperatures during 2003–07 might incorrectly be interpreted as a 54 Gt a⁻¹ additional mass loss when the dC _T/dt compaction is not taken into account.

Table 3. Mass effects (Gt a⁻¹) of dC _T/dt and dB/dt corrections on mass change estimates

4. Firn Densities (σ _a) Associated with Accumulation-Driven dM _a/dt and Pseudo-Densities Relating dH/dT to dM/dT

We also calculate the density, ρ _a, associated with dM _a/dt and the corresponding height changes using

(7)

Although we do not use ρ _a with to calculate dM _a/dt as we did previously (Reference ZwallyZwally and others, 2011), the values of ρ _a are of interest because they represent the density of the firn added (or not added) as a result of the sum of the δA(t) anomalies over time . We use the regional average ρ _a in our estimation of errors, because the values of ρ _a affect how systematic errors in dM _a/dt cause corresponding systematic errors in dM _d/dt and dM/dt though their coupling in Eqns (1) and (4).

The ρ _a represent firn distributed over a range of depths depending on the temporal distribution of the δA(t) and how the anomalies propagate into the firn. Therefore, ρ _a do not represent the density of a particular firn layer at a specific depth. The distribution of ρ _a calculated for all gridcells over the AIS for the ERS and ICESat time periods is shown in the histograms of Figure 8. Realistic values of ρ _a cannot be calculated for all gridcells, because singularities occur when approaches zero causing unrealistically large positive or negative values of v_a. For the ERS and ICESat periods, 76% and 86% respectively of the ρ _a values lie between 0.2 and 0.92; the outliers are excluded from calculations of regional averages. The most probable ρ _a for the AIS is ∼0.33 in both periods. The average regional values of ρ _a range from 0.37 to 0.61, with an average of 0.39, over the AIS (Table 4). The average ρ _a is lowest at 0.37 in the colder EA where the compaction is slower, and is larger in the warmer areas at 0.46 in WA and 0.60 in the AP.

Fig. 8. Distribution of firn densities, ρ _a, associated with dM _a/dt accumulation-driven changes for AIS. In 1992–2001 (red) and 2003–08 (black), the average ρ _a are 0.39 excluding outliers from singularities in Eqn (7).

Table 4. Accumulation density (ρ _a), average density (ρ _avg) and pseudo-densities (ρ _pseudoH and (ρ _pseudoI ) by region

In contrast to our calculation of dM _d/dt = ρ _i(dH _d/dt) using Eqn (5) with the well-defined ρ _i, use of to calculate the accumulation-driven mass change can only provide a rough approximation to dM _a/dt, because of the difficulties of calculating ρ _a. Furthermore, the local and regional values of dM _a/dt and dM _d/dt are highly variable and often of opposite signs, factors which preclude a priori selection of single or multiple densities to calculate dM/dt from measured dH/dt as has been done by several investigators. Examples are the selection of a single density (e.g. Reference Davis, Li, McConnell, Frey and HannaDavis and others, 2005; Reference ZwallyZwally and others, 2005) or binary densities depending on locations and assumptions about the dynamic or accumulation-driven character of the dH/dt (e.g. Reference SørensenSørensen and others, 2011; Reference ShepherdShepherd and others, 2012; Reference McMillanMcMillan and others, 2014). To further illustrate this issue, we consider several alternate definitions of densities. Reference Li and ZwallyLi and Zwally (2011) defined an average density as , which has values within the range of firn/ice densities (e.g. 0.80 and 0.86 in EA and 0.75 and 0.79 for the AIS; Table 4). However, there is no corresponding single dH/dt to use with ρ _avg to obtain the correct dM/dt. We also calculate ρ _pseudoH dM/dt/(dh/dt × Area) and ρ _pseudoI dM/dt/(dI/dt × Area) that give the correct dM/dt using our dM/dt, dH/dt and dI/dt from Tables 2 and 5. The values of the pseudo-densities vary widely as shown in Table 4 (e.g. maximal regional values of ρ _pseudoH = 0.35 in 1992–2001 and 2.6 in 2003–08 in WA and 7.1 in 2003–08 in EA2, which are caused by local values of and dH _d/dt that are in opposite directions). More reasonable values of ρ _pseudoH in the range 0.86–1.8 are calculated separately for the WA1 and WA2 subregions of WA and for the whole of EA where are small relative to dH _d/dt. The results for ρ _pseudoI are similar to those for ρ _pseudoH .

Table 5. Rates of total (dM/dt), dynamic-driven (dM _d/dt) and accumulation-driven (dM _a/dt) mass changes (Gt a⁻¹) and surface mass balance (SMB) (Gt a⁻¹) by drainage system (DS) and regions

In order to estimate dM/dt from either dH/dt or dI/dt, neither a single density nor binary densities dependent on location can properly account for the typical mixtures of accumulation-driven and dynamic-driven mass changes. The same problem applies to both the large areas of regions and sub-regions and to the smaller sizes of gridcells, as shown by the mixtures of accumulation-driven and dynamic-driven mass changes (often with opposite signs) at most locations in the maps of dM _a/dt and dM _d/dt (Figs 10 and 11). In particular, while the choice for relative ice density of 0.917 for the dynamically active portions of DS22, DS21, DS20, DS18 and DS13 (fig. 1 inset in Reference McMillanMcMillan and others, 2014) may seem appropriate, those areas also have significant dM _a/dt for which McMillan and others snow density of 0.35 would be more appropriate. Furthermore, the spatial distribution of dM _a/dt in those dynamically active areas changes with time as shown in Figure 10. More seriously, McMillan and others’ assignment of snow density to the rest of the AIS, in particular the very large areas of EA and WA2, with mostly dynamic thickening (Fig. 11) and small or negative dM _a/dt (Fig. 10), causes underestimates of the mass gains.

5. Discussion of Results

Maps of the derived dM/dt, dM _a/dt and dM _d/dt for the two periods are presented in Figures 9–11, values by DS and defined regions are shown in Table 5, and values only by DS in Figure 12. For the whole of the AIS, mass gains from snow accumulation exceeded losses from ice discharge by 112 ± 61 and 82 ± 25 Gt a⁻¹ respectively during the 1992–2001 and 2003–08 measurement periods. The overall positive balance is due to large net gains of +136 ± 50 and +136 ± 28 Gt a⁻¹ in the two periods in EA, plus smaller net gains of +44 ± 14 and +72 ± 9 Gt a⁻¹ in WA2. Those net gains together exceed the total losses of −60 ± 12 and –97 ± 6 Gt a⁻¹ in WA1 and the losses of −9 ± 10 and –29 ± 2 Gt a⁻¹ in the AP.

Fig. 9. Maps of total mass changes, dM/dt: (a) during 1992–2001 with overall positive balance of +112 ± 61 Gt a⁻¹ (dH/dt from ICESat used south of 81.5°E); (b) during 2003–08 with +82 ± 25 Gt a⁻¹. The distribution of dM _d/dt in Figure 11 is very similar to dM/dt. In WA1, the mass-loss rate increased from 60 Gt a⁻¹ to 97 Gt a⁻¹ mainly due to a 51 Gt a⁻¹ increase in dynamic thinning (note inland spreading of mass loss in 2003–08). In WA2 and EA together, the total mass gains of 180 and 208 Gt a⁻¹ in the two periods are caused mainly by persistent dynamic thickening (deficiency of ice flow relative to long-term 〈A〉 ) of 202 and 211 Gt a⁻¹, which is a residual of the dynamic response to a marked increase in precipitation at the beginning of the Holocene.

Fig. 10. Maps of the accumulation-driven mass changes, dM _a/dt, during (a) 1992–2001 and (b) 2003–08. dM _a/dt are generally smaller than dM _d/dt (note 8× larger scale than in Figs 9 and 11). In WA1, the dM _a/dt increase of 14 Gt a⁻¹ in the 2003–08 period partially offset the 51 Gt a⁻¹ increase in dynamic thinning. The increase in dM _a/dt during 2003–08 also extended over DS23 and DS1, causing a 20 Gt a⁻¹ increase in WA2, with negligible net changes in DS18 and DS19. In the AP (DS24–27), the dM _a/dt variation is slightly positive during 1992–2001 and slightly negative during 2003–08. In EA, the variability of dM _a/dt ranges from –8 Gt a⁻¹ in DS12 during 1992–2001 and −9 Gt a⁻¹ in DS14 during 2003–08 to +6 Gt a⁻¹ in DS3 during 2003–08, but is unchanged at –11 Gt a⁻¹ between periods over the whole of EA. (dM _a/dt are calculated over the whole area.)

Fig. 11. Maps of the dynamic-driven mass changes, dM _d/dt, (a) for 1992–2001 (dH/dt from ICESat used south of 81.5°E) and (b) for 2003–08. In WA1, the net dynamic-loss rate increased from 55 Gt a⁻¹ to 106 Gt a⁻¹. Dynamic thickening (excess of long-term accumulation over ice flux) occurred over WA2 and EA in both periods. Dynamic thickening (27 Gt a⁻¹ in 2003–08) is strongest in DS18 in an area inland from Kamb Ice Stream that stagnated 150 years ago. Similarly, the dynamic thinning in Eastern DS17 and Western DS18 is inland of Mercer and Whillans Ice Streams, which restarted flowing 400 years ago.

Fig. 12. Rates of mass change and SMB by DS. (a) Total mass change and SMB, showing net mass losses in the AP and WA1 regions and gains in WA2 and EA. (b) Accumulation-driven mass change, showing the negative anomaly in WA1 and WA2 during 1992–2001 and the positive anomaly during 2003–08. (c) Dynamic-driven mass change, showing the increase in loss from dynamic thinning in the AP and WA1 and the gains from long-term dynamic thickening in WA2 and EA.

The indicated change of −29 ± 66 Gt a⁻¹ between measurement periods in the overall AIS mass balance is not significant relative to the estimated uncertainties. In the EA region, dM/dt gains of 136 Gt a⁻¹ in both periods also indicate no significant change (δ = 0 ± 57 Gt a⁻¹), which is also true separately for the EA1 and EA2 subregions (δ = −7 ± 40 Gta⁻¹ and +7 ± 33 Gt a⁻¹). In contrast, regional changes in mass gains and mass losses in WA and the AP are significant. In WA2, the increase in mass gain from 44 ± 14 Gt a⁻¹ to 72 ± 9 Gt a⁻¹ (δ = +28 ± 16 Gt a⁻¹) is statistically significant and is largely due to a 20 ± 8 Gt a⁻¹ increase in dM _a/dt from greater snowfall. Increases in mass loss in WA1 and the AP are also significant, and consistent with observations of accelerated mass loss in those regions (discussed below). Specifically in WA1, the mass loss increase is from −60 ± 12 Gt a⁻¹ to −97 ± 6 Gt a⁻¹ (δ = –37 ± 13 Gt a⁻¹), and in the AP the mass loss increase is from −9 ± 10 Gt a⁻¹ to −29 ± 2 Gt a⁻¹ (δ = −20 ± 10 Gt a⁻¹).

The dM _a/dt maps (Fig. 10) illustrate the large spatial variability of accumulation in each period and the large temporal variability between periods. Regionally, dM _a/dt increased by 14 ± 5 Gt a⁻¹ in WA1 and by 20 ± 8 Gt a⁻¹ in WA2 and decreased by 5 ± 2 Gt a⁻¹ in the AP (Table 5). Although dM _a/dt over EA is unchanged at −11 ± 6 Gt a⁻¹ in both periods, it increased by 21 ± 8 Gt a⁻¹ in EA1 and decreased by 21 ± 13 Gt a⁻¹ in EA2, indicating a large regional shift in accumulation anomalies. In general, dM _a/dt are small compared to dM _d/dt. Overall, dM _a/dt as a fraction of dM/dt are −21% and +6% for the two periods, whereas dM _d/dt are much larger fractions at +121% and +95%.

In EA, the large dM/dt gains of 136 Gt a⁻¹ in the two periods consist of dM _d/dt gains of 147 Gt a⁻¹ which are much larger than the small dM _a/dt losses of 11 ± 6 and 11 ± 6 Gt a⁻¹. The small losses from accumulation anomalies relative to the 27 year average accumulation show that the elevation increases and mass gain in EA are not caused by increasing snowfall during the 17 year measurement period. This result contradicts previous conclusions that growth in EA during the first period was driven by contemporaneous increases in snowfall (Reference Davis, Li, McConnell, Frey and HannaDavis and others, 2005), a conclusion questioned by Reference MonaghanMonaghan and others (2006), or by accumulation variability and depth of the firn layer (Reference HelsenHelsen and others, 2008).

Instead, the dM/dt gains in EA are caused by a large dynamic thickening of 147 Gt a⁻¹ due to a deficiency of ice flow relative to the long-term 〈A〉 (or equivalently an excess of 〈A〉 over the ice flow). Dynamic thickening extends over most of EA, and none of the DS have significant net dynamic thinning (i.e. <−6 ± 6 Gt a⁻¹). DS13, which is close to dynamic balance in both periods due to near-coastal thinning on Totten Glacier, differs from the adjacent DS12 and DS14 which exhibit dynamic gains of 22 and 24 Gt a⁻¹ and coastal thickening.

A dominant characteristic of the dynamic thickening in EA is its persistence, with no significant change between periods. This thickening is likely a result of the marked increase in accumulation that began in the early Holocene, ∼10 ka ago, with a 67–266% increase from the Last Glacial Maximum as derived from six ice cores (Reference SiegertSiegert, 2003). A residual thickening from the Holocene increase in accumulation is consistent with the characteristic slow response time of the ice flow to accumulation changes in EA. Comparing the 147 Gt a⁻¹ dynamic gain with the SMB gives a current residual dynamic imbalance of ∼13%.

Using 100% as a rough estimate of the average accumulation increase in the early Holocene implies that the dynamic adjustment over the last 10 ka is ∼87% complete on average over EA. In the interior areas of the ice sheet, dynamic response times are much longer than in the higher-accumulation and thinner areas closer to the coastal margins. For example, in the interior of EA where the accumulation is low (e.g. 0.03 m a⁻¹) and the ice is thick (e.g. 2500 m), the ice thickness increase over 10 ka is only ∼300 m. The corresponding fractional increase in driving stress over that time is only ∼12%, implying the flow takes a very long time to adjust to the change.

The WA2 portion of WA also exhibits a significant persistent dynamic thickening of ∼60 Gt a⁻¹. Together, the combined dM _d/dt of 202 and 211 Gt a⁻¹ shows no significant short-term change. However, the timescale for some of the dynamic thickening in WA2 is shorter than the 10 ka thickening in EA based on other observations. Much of the 27 ± 2 Gt a⁻¹ dynamic thickening in DS18 is inland of Kamb Ice Stream flowing into the Ross Ice Shelf, and is likely a result of ice-stream stagnation there ∼150 years ago (Reference Joughin, Tulaczyk, Bindschadler and PriceJoughin and others, 2002). Similarly, the dynamic thinning in the eastern part of DS17 and western part of DS18 is inland of Mercer and Whillans Ice Streams, which restarted flowing 400 years ago (Reference Hulbe and FahnestockHulbe and Fahnestock, 2007).

In contrast, WA1 shows a significant short-term increase in dynamic thinning from −55 ± 14 Gt a⁻¹ to −106 ± 11 Gt a⁻¹ (δ = −51 ± 17 Gt a⁻¹). DS22, with ice discharge through Pine Island Glacier, has the largest increase in mass loss, from −8 ± 14 Gt a⁻¹ to −35 ± 11 Gt a⁻¹ (δ = −26 ± 9 Gt a⁻¹), consistent with recent glacier acceleration (Reference Rignot, Vaughan, Schmeltz, Dupont and MacAyealRignot and others, 2002; Reference Wingham, Wallis and ShepherdWingham and others, 2009). The increased thinning in DS21 and DS22 extends >200 km inland from the front to the ice divide (Fig. 11). DS21, which discharges through Thwaites and Smith Glaciers, shows the largest dynamic thinning of any DS in Antarctica, increasing from −39 ± 6 Gt a⁻¹ to –56 ± 5 Gt a⁻¹ (δ = −16 ± 8 Gt a⁻¹). DS20 shows a smaller increase in dynamic thinning of −8 ± 3 Gt a⁻¹, which may be associated with the thinning of the Dotson and Getz Ice Shelves (Reference ZwallyZwally and others, 2005; Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012; Reference Paolo, Fricker and PadmanPaolo and others, 2015).

Much of the 51 ± 17 Gt a⁻¹ increase in dynamic loss in WA1 is offset by a contemporaneous increase in snowfall of 34 ± 12 Gt a⁻¹ in the whole of WA. The dM _a/dt patterns in DS21 and DS22 of WA1 and DS23 and DS1 of WA2 differ between periods (Fig. 10), with a negative anomaly of 16 Gt a⁻¹ during 1992–2001 and a positive anomaly of 20 Gt a⁻¹ in 2003–08. For the whole of WA, the dM _a/dt increase of 34 ± 12 Gt a⁻¹ offset most of the −43 ± 37 Gt a⁻¹ increase in dynamic thinning. In EA, snowfall anomalies subtracted 11 Gt a⁻¹ from the total dM/dt in both periods. Over the whole AIS, accumulation anomalies reduced the balance by 24 Gt a⁻¹ during 1992–2001 and added 5 Gt a⁻¹ during 2003–08.

Between periods, the dynamic thinning in the AP increased by 15 ± 12 Gt a⁻¹, which is consistent with (1) an increased loss (Reference Rott, Müller, Nagler and FloricioiuRott and others, 2011) of 4.3 Gt a⁻¹ from glacier acceleration and thinning in the Larsen B basin following ice-shelf disintegration in 2002, (2) a loss (Reference Shuman, Berthier and ScambosShuman and others, 2011) of 11.2 Gt a⁻¹ in the Larsen A and B basins, (3) widespread acceleration of tidewater glaciers (Reference Pritchard and VaughanPritchard and Vaughan, 2007), and (4) extensive retreat of marine glacier fronts (Reference Cook, Fox, Vaughan and FerrignoCook and others, 2005).

As noted in Section 1, the evaluation of Reference ShepherdShepherd and others (2012) eliminated some larger estimates of Antarctic mass loss and gave a reconciled mean of −72 ± 43 Gt a⁻¹ for the technique intercomparison period October 2003 to December 2008 (fig. 4 of Reference ShepherdShepherd and others, 2012). For WA, Shepherd and others’ reconciled mean of −67 ± 21 Gt a⁻¹ contained the means and most of the ranges of the RA, LA, GR and IOM results from several research groups. Also, their reconciled mean of −28 ± 10 Gt a⁻¹ for the AP contained the means and most of the ranges of the LA, GR and IOM results. In contrast, for EA the reconciled mean of +24 ± 36 Gt a⁻¹ contained the means and most of the ranges of only the RA and GR results. The EA mean from IOM was more negative at −30 ± 76 Gt a⁻¹ and the mean from LA was more positive at +109 ± 57 Gt a⁻¹. Variations among the LA estimates were due to different methods of estimating dM/dt from dH/dt, differing choices for the ICESat inter-campaign biases, differing methods of dh/dt solutions, and different data-editing procedures, which have been improved for our EA mass gain of +136 ± 28 Gt a⁻¹ in Table 5.

Possible causes for the lower RA estimate of +22 ± 39 Gt a⁻¹ for EA are inadequate corrections to the Envisat data for variable radar penetration depth (as discussed in the Appendix) and the use of a low density of 0.35 to estimate dM/dt from dH/dt assuming the changes were due to snowfall anomalies (as discussed in Section 4). A likely cause for the lower GR estimate is the sensitivity of the GR estimates to the glacial isostatic adjustment (GIA) correction, as discussed in the Appendix where we note that a −1.6 mm a⁻¹ change in the modeled dB/dt would bring the GR and our dM/dt into agreement at approximately +150 Gt a⁻¹. The additional ice loading from a dynamic thickening of 1.59 cm a⁻¹ over EA (Table 2) for 10 ka implies an additional bedrock depression of 27 m continuing at a rate of 2.65 mm a⁻¹ assuming full long-term isostatic adjustment. Therefore, the −1.6 mm a⁻¹ needed to bring the gravimetry and altimetry dM/dt estimates into agreement is only 60% of the full isostatic adjustment rate, and therefore within the range of what can be expected if the ice loading implied by the long-term dynamic thickening is accounted for in the GIA models.

A principal issue regarding the IOM concerns the extrapolation procedure used to estimate the ice discharge from non-observed areas, for which ice velocity measurements were not available, as shown in Reference Zwally and GiovinettoZwally and Giovinetto (2011) regarding the IOM results of Reference RignotRignot and others (2008, Reference Rignot, Velicogna, Van den Broeke, Monaghan and Lenaerts2011). The estimate in Reference RignotRignot and others (2008) for EA was near-zero at −4 ± 61 Gt a⁻¹, which was based on an input estimate of 1131 Gt a⁻¹, an observed output of 784 Gt a⁻¹ and an extrapolated non-observed output of 349 Gt a⁻¹. Therefore, a disproportionate 31% of Rignot and others’ total output from EA came from only 15% of the area that was non-observed. If the non-observed output had been scaled in proportion to the area as stated, then the non-observed output would be only 173 Gt a⁻¹ (i.e. 15% from 15% of the area) and the net IOM balance would be +174 Gt a⁻¹. Or if the non-observed slower-moving areas were only 70% as effective at discharging ice as the faster-moving observed areas, then the modified non-observed output would be +121 Gt a⁻¹ and the net IOM balance estimate for EA would be +226 Gt a⁻¹ (Reference Zwally and GiovinettoZwally and Giovinetto, 2011).

6. Conclusions

During the period 1992–2001, the Antarctic mass gain from snow accumulation exceeded the mass loss from ice discharge by 112 ± 61 Gt a⁻¹. During 2003–08, the gain exceeded the loss by a similar 82 ± 25 Gt a⁻¹, which is 4% of the SMB and equivalent to 0.23 mm a⁻¹ sea-level depletion. The mass-balance distribution is mostly positive in EA and WA2, and mostly negative in the AP and WA1. Although the overall balance of the AIS shows no significant change over the 17 years, significant increases in the dynamic losses occurred from the AP and several coastal DS in WA; however, those increases in dynamic losses were partially compensated by accumulation increases in WA. In EA, no change occurred in the short-term contributions of short-term accumulation, which are slightly negative by 11 Gt a⁻¹ in both periods compared to the 27 year average accumulation rate. Therefore, the positive mass balance in EA clearly has not been caused by contemporaneous increases in snowfall.

Major characteristics of the Antarctic imbalance are the large long-term dynamic thickening of 147 Gt a⁻¹ in EA and 59 Gt a⁻¹ in WA2, and the short-term increase of dynamic thinning from −55 Gt a⁻¹ to −106 Gt a⁻¹ in WA1 and from −9 Gt a⁻¹ to −29 Gt a⁻¹ in the AP. While continued increases in dynamic thinning are expected in response to thinning of adjacent ice shelves in WA1 and the AP (Reference Cook, Fox, Vaughan and FerrignoCook and others, 2005; Reference Pritchard and VaughanPritchard and Vaughan, 2007; Reference Rott, Müller, Nagler and FloricioiuRott and others, 2011; Reference Shuman, Berthier and ScambosShuman and others, 2011; Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012), the long-term dynamic thickening in EA and WA2 should provide a significant buffer against such continued increases in mass loss. If dynamic thinning continues to increase at the same rate of 4 Gt a⁻² with no offset from further increases in snowfall, the positive balance of the AIS will decrease from the recent 82 Gt a⁻¹ to zero in ∼20 years. However, compensating increases in snowfall with climate warming may also be expected (Reference Gregory and HuybrechtsGregory and Huybrechts, 2006; Reference Winkelmann, Levermann, Martin and FrielerWinkelmann and others, 2012).