## Abstract

A comprehensive analysis of satellite datasets has estimated that the ice sheets of Greenland, West Antarctica, the Antarctic Peninsula, and East Antarctica experienced a net mass loss of −100 ± 92 Gt yr^{−1} over the period 1992–2000 and −298 ± 58 Gt yr^{−1} over the period 2000–11, representing an increase of −198 ± 109 Gt yr^{−1} between the two epochs. The authors demonstrate that the time rate of change of the degree-four zonal harmonic of Earth's gravitational potential provides an independent check on these mass balances that is less sensitive to uncertainties that have contaminated previous analyses of the degree-2 zonal harmonic [e.g., due to ongoing glacial isostatic adjustment (GIA), solid Earth body tides, and core–mantle coupling]. For the period 2000–11, the signal implied by the ice sheet mass flux cited above is (3.8 ± 0.6) × 10^{−11} yr^{−1}, whereas the change in the harmonic across the two epochs is (2.3 ± 1.1) × 10^{−11} yr^{−1}. In comparison, using satellite laser ranging (SLR) data, the authors estimate a GIA-corrected value of (3.8 ± 0.6) × 10^{−11} yr^{−1} for the epoch 2000–11 and a change across the two epochs of (5.3 ± 1.6) × 10^{−11} yr^{−1}. The authors conclude that the former supports recent estimates of melting over the last decade, whereas the latter suggests either that estimated melt rates for the earlier epoch were too high or that the uncertainty associated with the SLR-based inference of during the earlier epoch is underestimated.

## 1. Introduction

A recent, combined reanalysis of satellite-based regional altimetry, interferometry, and gravimetry datasets (Shepherd et al. 2012, hereafter S2012) has estimated the mass balance of ice sheets covering Greenland (GIS), the Antarctic Peninsula (APIS), West Antarctica (WAIS), and East Antarctica (EAIS) over the period 1992–2011 as well as during various epochs within this two-decade time window (Table 1). The uncertainties ascribed to these estimates account for the limited and unique time span associated with each of the three space geodetic techniques; observation errors; and, through modeling, uncertainties in corrections for ice sheet surface mass balance and glacial isostatic adjustment (GIA). Over the period 2000–11, for example, the total estimated mass flux from these ice sheets, −298 ± 58 Gt yr^{−1} (1-*σ* uncertainty), equates to a eustatic sea level (ESL) rise of 0.83 ± 0.16 mm yr^{−1}. [We adopt the term “eustatic” as the geographically uniform shift in sea level over the oceans that would yield a volume equal to the meltwater addition. Gregory et al. (2013) have suggested the term “barystatic” for this quantity.] In this paper, we demonstrate that this estimate may be independently tested by invoking the constraint provided by the time rate of change of the zonal harmonic of Earth's gravitational potential at spherical harmonic degree 4 (or, alternatively, the trend in the Stokes coefficient; Hofmann-Wellenhof and Moritz 2006).

Estimates of secular trends in the low-degree zonal harmonics of the geopotential based on satellite laser ranging (SLR) have been available for over a quarter of a century, initially at spherical harmonic degree 2 (Yoder et al. 1983; Rubincam 1984) and later at higher degrees (Cheng et al. 1989). The first analyses of the harmonic interpreted it as being dominated by ongoing GIA (Yoder et al. 1983; Rubincam 1984; Wu and Peltier 1984; Peltier 1983; Yuen and Sabadini 1985; Ivins et al. 1993; Vermeersen et al. 1998) and commonly used the harmonic as a basis for inferring the average viscosity of the lower mantle. Motivated, at least in part, by the growing appreciation of the impact of modern climate change on Earth's gravitational field, other, generally later, analyses of these harmonics highlighted the potential signal from recent changes in the mass balance of polar ice sheets and mountain glaciers (Sabadini et al. 1988; Mitrovica and Peltier 1993; James and Ivins 1997; Trupin and Panfili 1997). In addition, coupling between the fluid outer core and the overlying mantle has also been identified as a potentially important contributor to the harmonic (Fang et al. 1996).

A significant change in the secular trend of the coefficient in the mid-1990s, as estimated from SLR data (Cox and Chao 2002), has led to renewed interest in the geophysical interpretation of this harmonic. Explanations for the observation have included a possible error in the standard geodetic correction for the 18.6-yr body tide signal (Benjamin et al. 2006), a change in the dynamics of the El Niño–Southern Oscillation (Cheng and Tapley 2004; Marcus et al. 2009) and increased melting of either mountain glaciers (Marcus et al. 2009; Dickey et al. 2002) or polar ice sheets (Nerem and Wahr 2011).

The spherical harmonic basis function at degree 2 and order zero has the same sign in both high northern and southern latitudes. As such, the datum provides an integrated measure of the combined (recent) mass balance of the GIS and the Antarctic ice sheet (AIS; Ivins et al. 1993; Mitrovica and Peltier 1993). In contrast, basis functions associated with odd harmonics have a different sign near the north and south poles, and they thus provide a measure of the difference in mass balance of the two polar ice sheets. Because the mass balance of the AIS and GIS are of comparable magnitudes, this difference will be sensitive to their relative sizes. However, while the integrative nature of the harmonic is an advantage over the odd harmonics for inferring total polar ice sheet mass changes, there are other sources of uncertainty that decrease with higher spherical harmonic degrees (e.g., 18.6-yr body tide, core–mantle coupling, etc.). Consequently, the optimal spectral component to examine for estimating net polar mass change is an even harmonic that is of a low enough degree that the basis function integrates over the region of interest (in this case the polar ice sheets) but is of high enough degree to reduce the contribution from nonmelt signals.

In this regard, the harmonic, which has been overlooked in recent geodetic analyses of climate signals, provides a potentially more rigorous constraint on integrated polar ice sheet mass balance than . The signal due to recent mass loss from either the GIS or AIS should only be ~10% smaller than the analogous signal in the harmonic (Ivins et al. 1993; Mitrovica and Peltier 1993), and hence it should be significant and observable (as we demonstrate below). Furthermore, the datum avoids many of the various complications that arise in analyzing the harmonic. For example, the signal associated with the 18.6-yr body tide and core–mantle coupling will be about three orders of magnitude smaller than the analogous signal in .

The and harmonics will both have a significant contribution from ongoing GIA and any effort to estimate recent mass flux from polar ice sheets based upon these harmonics must address the uncertainty associated with this signal. In the analysis described below, we implement two strategies for this purpose. The first is to consider observed changes in the rates of the harmonics: for example, from the epoch 1992–2000 relative to the epoch 2000–11 (Table 1, column 4). This will remove any sensitivity to the GIA process because the change in the GIA-induced rate will be negligible over the two-decade time window of the present analysis. A second strategy for dealing with the contamination of the harmonic due to GIA, which we turn to next, is forward modeling of the GIA process.

## 2. Analysis and results

### a. GIA predictions

The results in Fig. 1 show the predicted GIA-induced signal in and as a function of the adopted lower-mantle viscosity of the earth model. The calculations are based on spherically symmetric, Maxwell viscoelastic Earth models (Peltier 1974) with density and elastic structure given by the Preliminary Reference Earth Model (PREM) seismic model (Dziewonski and Anderson 1981). The lithospheric thickness and upper-mantle viscosity are set to 70 km and 5 × 10^{20} Pa s, respectively. Our predictions of the long-wavelength zonal harmonics are insensitive to the former, and in the calculations below we explicitly investigate the sensitivity to the latter. The predictions adopt the ICE-5G model for the evolution of ice thickness over the last glacial cycle (Peltier 2004), and they incorporate a gravitationally self-consistent sea level theory (Kendall et al. 2005). This theory outputs the present-day rate of change of relative sea level and sea surface height; the spherical harmonic coefficients of the latter are proportional to the rates of change of the Stokes coefficients (Mitrovica and Peltier 1993).

In addition to the total GIA-induced perturbation in each harmonic (solid lines in Fig. 1), we decompose the perturbation into contributions from the GIS, AIS, and all other ice sheets included in the ICE-5G inventory (e.g., Laurentide, Fennoscandian, etc.). The results indicate that, although ice sheets other than the AIS and GIS dominate the GIA-induced signal in the harmonic, the signal associated with AIS evolution over the last glacial cycle dominates the GIA prediction of the harmonic. Moreover, for lower-mantle viscosities greater than ~3 × 10^{21} Pa s, which is a hard lower bound for ice age–based inferences (Peltier 2004; Nakada and Lambeck 1989; Lambeck et al. 1998; Mitrovica and Forte 2004), the GIA-induced signal driven by ice sheets other than the AIS is relatively insensitive to uncertainties in the radial profile of mantle viscosity. We have verified that this insensitivity extends to predictions based on Earth models in which the upper-mantle viscosity is varied from 3 × 10^{20} to 10 × 10^{20} Pa s (Fig. 2) and an independent ice history (Fleming and Lambeck 2004). These predictions are also insensitive to the adopted lithospheric thickness.

The decomposition in Fig. 1 highlights another important advantage of using the observed harmonic, rather than the harmonic, as a basis for estimating recent polar ice sheet mass flux. In particular, any uncertainty in the GIA correction to is largely isolated to the contribution to the GIA signal associated with the AIS. In this regard, GIA models for Antarctica over the last glacial cycle have recently been improved, relative to earlier models such as ICE-5G/VM2 (Peltier 2004), by combining numerical ice sheet models with geological observations of ice extent, local relative sea level histories, and GPS-derived uplift rates (Whitehouse et al. 2012). In adopting the harmonic to estimate recent polar mass flux, a second strategy for dealing with the GIA signal is to use the correction generated by this new, improved class of ice-age AIS models.

We begin by focusing on the period 2000–11 and compute a GIA correction based on a combination of our GIA modeling and the W12a Antarctic GIA model (Whitehouse et al. 2012). The results in Fig. 2 indicate that the GIA contribution to from sources other than the Antarctic falls within the (conservative) range (−0.7 ± 0.2) × 10^{−11} yr^{−1}. Moreover, using the upper- and lower-bound estimates of the W12a model (Whitehouse et al. 2012) yields an Antarctic GIA contribution to of (−0.96 ± 0.14) × 10^{−11} yr^{−1}. Thus, we estimate the total GIA contribution to be (−1.7 ± 0.2) × 10^{−11} yr^{−1}. We analyzed the SLR-derived times series for over the time periods 1992–2000 and 2000–11 (Fig. 3). Correcting the observed trend in during each of these epochs using the total GIA signal yields residual trends of (−1.6 ± 1.5) × 10^{−11} yr^{−1} and (3.7 ± 0.6) × 10^{−11} yr^{−1}_{,} respectively.

We note that the Antarctic ice volume change in the W12a model (Whitehouse et al. 2012) is significantly smaller than in the ICE-5G history. The difference amounts to approximately 10 m of equivalent eustatic sea level rise, or ~10% of the total non-Antarctic ice volume in ICE-5G. This reduction would have to be compensated by an increase in the excess volume of Northern Hemisphere ice relative to ICE-5G in order to maintain a fit to far-field relative sea level records (Austermann et al. 2013). This suggests that the amplitude of our computed non-Antarctic GIA contribution (Fig. 2) will be biased low by ~10%. Accounting for this effect yields GIA-corrected trends in of (−1.5 ± 1.5) × 10^{−11} yr^{−1} and (3.8 ± 0.6) × 10^{−11} yr^{−1} for the epochs 1992–2000 and 2000–11, respectively.

### b. Polar mass change contribution to

We next turn to the estimate of polar ice sheet mass balance based on the comprehensive S2012 analysis of satellite-based measurements (Table 1). What signal in does this level of modern melting imply? To answer this, we adopted a special, elastic case of our sea level software to compute the signal associated with uniform melting over 1) EAIS, 2) WAIS, 3) APIS, and 4) GIS (Table 2). The calculations are normalized so that they represent perturbations to the trend per millimeter per year of ESL rise. Using these results, we converted the total mass flux estimates in Table 1 into a signal (last row of Table 1). For the epoch 2000–11, we obtain the range (2.9 ± 0.6) × 10^{−11} yr^{−1}. A final, significant contribution to the trend is associated with melting of mountain glaciers. A tabulation of mountain glacier mass balance by Jacob et al. (2012), based on Gravity Recovery and Climate Experiment (GRACE) gravity data, covers the period 2003–10. Using this database, we computed a trend of 0.9 × 10^{−11} yr^{−1}. If we assume that this rate is applicable to the 2000–11 epoch adopted in Table 1, we can augment the above estimate [(2.9 ± 0.6) × 10^{−11} yr^{−1}] to include the mountain glacier signal; in particular, we arrive at an estimate of (3.8 ± 0.6) × 10^{−11} yr^{−1}. The consistency between this estimate and our GIA-corrected SLR-derived trend, (3.8 ± 0.6) × 10^{−11} yr^{−1}, provides independent support for recent estimates of polar ice sheet (S2012) and mountain glacier (Jacob et al. 2012) mass flux over the last decade.

Alternatively, in inferring polar ice sheet mass balance, we can avoid entirely the contaminating effect of GIA by considering the change in the trend of the harmonic from 1992–2000 to 2000–11 (i.e., the first strategy for dealing with GIA discussed above). From Table 1, the change in polar ice mass flux across these time periods inferred in the S2012 study maps into a change in the signal of (1.9 ± 1.1) × 10^{−11} yr^{−1}. Taking into account a change in the mountain glacier signal using the tabulations of Kaser et al. (2006) for the first epoch and (as above) Jacob et al. (2012) for the second epoch raises this value to (2.3 ± 1.1) × 10^{−11} yr^{−1}. As noted above, our direct, SLR-derived estimate of the change in trend is (5.3 ± 1.6) × 10^{−11} yr^{−1}. This analysis suggests either that estimates of melt rates during the earlier epoch are too high and/or that the uncertainty in the SLR-derived change in the signal (Fig. 3) has been underestimated. In regard to the latter, we note that the observed trend across the earlier epoch [(−3.3 ± 1.5) × 10^{−11} yr^{−1} as cited in the caption of Fig. 3] is sensitive to the start time of the epoch.

## 3. Conclusions

We have demonstrated that estimates of net polar ice sheet mass balance may be independently tested by invoking the long-neglected constraint associated with the rate of change of the zonal harmonic degree-4 of Earth's potential. The harmonic is relatively insensitive to the uncertainties associated with signals from the 18.6-yr Earth tide and core–mantle coupling that must be accounted for in analyses of the harmonic. Analyses of the harmonic must, however, address the potentially significant signal from GIA, either by considering a change in trend across two epochs, rather than the absolute value of the harmonic, or through numerical modeling. In regard to the latter, we have shown that the dominant GIA contribution to the is associated with the Antarctic ice sheet, and thus the uncertainty in the GIA correction can be reduced significantly by adopting recent, well-constrained models of Antarctic GIA. In any event, our application of this approach indicates that the trend supports recent inferences of increased melting from polar ice sheets and mountain glaciers over the last decade. This demonstrates that analyses of the harmonic can play an important role in advancing our understanding of ice sheet stability in a progressively warming world.

## Acknowledgments

We thank Felix Landerer and two anonymous reviewers for helpful comments and suggestions and Minkang Cheng for providing the time series of SLR zonal harmonics used in the present study. The study was funded by Harvard University, Princeton University, the Canadian Institute for Advanced Research, the Natural Sciences and Engineering Research Council of Canada, and NSF Grant EAR-1014606.

## REFERENCES

_{2}anomaly

*Physical Geodesy.*2nd ed. Springer-Verlag Wien, 403 pp.

_{2}coefficient from LAGEOS and non-tidal acceleration of earth rotation

_{2}observation