Antarctic contribution to sea level rise observed by GRACE with improved GIA correction

(1)

Antarctic contribution to sea level rise observed by GRACE

with improved GIA correction

Erik R. Ivins,1Thomas S. James,2,3John Wahr,4 Ernst J. O. Schrama,5 Felix W. Landerer,1and Karen M. Simon3,2

Received 31 August 2012; revised 25 April 2013; accepted 3 May 2013.

[1] Antarctic volume changes during the past 21 thousand years are smaller than

previously thought, and here we construct an ice sheet history that drives a forward model prediction of the glacial isostatic adjustment (GIA) gravity signal. The new model, in turn, should give predictions that are constrained with recent uplift data. The impact of the GIA signal on a Gravity Recovery and Climate Experiment (GRACE) Antarctic mass balance estimate depends on the specific GRACE analysis method used. For the method described in this paper, the GIA contribution to the apparent surface mass change is re-evaluated to be +55 ˙ 13 Gt/yr by considering a revised ice history model and a parameter search for vertical motion predictions that best fit the GPS observations at 18 high-quality stations. Although the GIA model spans a range of possible Earth rheological structure values, the data are not yet sufficient for solving for a preferred value of upper and lower mantle viscosity nor for a preferred lithospheric thickness. GRACE monthly solutions from the Center for Space Research Release 04 (CSR-RL04) release time series from January 2003 to the beginning of January 2012, uncorrected for GIA, yield an ice mass rate of +2.9˙ 29 Gt/yr. The new GIA correction increases the solved-for ice mass imbalance of Antarctica to –57 ˙ 34 Gt/yr. The revised GIA correction is smaller than past GRACE estimates by about 50 to 90 Gt/yr. The new upper bound to the sea level rise from the Antarctic ice sheet, averaged over the time span 2003.0–2012.0, is about 0.16˙ 0.09 mm/yr.

Citation: Ivins, E. R., T. S. James, J. Wahr, E. J. O. Schrama, F. W. Landerer, and K. M. Simon (2013), Antarctic contribution to sea level rise observed by GRACE with improved GIA correction, J. Geophys. Res. Solid Earth, 118, doi:10.1002/jgrb.50208.

1. Introduction

[2] Mantle creep deformation allows any large-scale

surface-load-induced deviation from gravitational hydro-static equilibrium in the Earth to be re-established over time by slow viscous flow [Kaula, 1980]. The most ubiquitous and widespread mantle flow of this type during the Pleis-tocene is caused by large exchanges of water mass between the oceans and the great continental ice sheets [Milne et al., 2009]. The slow viscous mantle flow causes global changes in the shape and the Earth’s gravitational field during the

Additional supporting information may be found in the online version of this article.

1_{Jet Propulsion Lab, California Institute of Technology, Pasadena,} California, USA.

2_{Geological Survey of Canada, Sidney, British Columbia, Canada.} 3_{School of Earth and Ocean Sciences, University of Victoria, Victoria,} British Columbia, Canada.

4_{Department of Physics and Cooperative Institute for Research in} Environmental Sciences, University of Colorado, Boulder, Colorado, USA. 5_{Faculty of Aerospace Engineering, Delft University of Technology,} Delft, Netherlands.

Corresponding author: E. R. Ivins, Jet Propulsion Lab, California Insti-tute of Technology, MS 300-233, 4800 Oak Grove Dr., Pasadena, CA 91109-8099, USA. (erik.r.ivins@jpl.nasa.gov)

present day. The phenomenon is known as glacial isostatic adjustment, or GIA. One of the main goals of modern geodesy is to measure global trends in fresh water transport that are separable from other geophysical signals, includ-ing GIA [Wahr et al., 1998; Wu et al., 2010; Kusche et al., 2012]. The largest potential source of sustained transport is found in the 57 m of equivalent global mean sea level rise stored in the continental Antarctic ice sheet. Geodetic sciences have a major role to play in constraining models of past and present-day climate change.

[3] The idea that the geodetic measurement of elastic

crustal motions using Global Positioning System (GPS) sta-tions on rock outcrops adjacent to the great ice sheets could be sensitive to interannual variability of the net water mass balance above the grounding line was shown by Wahr et al. [1995] and Conrad and Hager [1995]. That concept is now in use for monitoring the mass balance of the Greenland ice sheet [Bevis et al., 2012]. Largely because of ice sheet reconstructions, the vertical motions predicted due to man-tle viscous ﬂow response to glacial unloading since the Last Glacial Maximum (LGM) at 21 ka are thought to be many times larger than the elastic signature associated with changes in the Antarctic ice sheet (AIS) [James and Ivins,

1995]. Some evidence bearing this out came from both14_C

dates of retreat in the Ross Sea [Conway et al., 1999] and

(2)

Figure 1. Antarctic ice velocity map [Rignot et al., 2011] with paleoconstraint data locations

(yellow-numbered regions are discussed in the text). “IC” refers to deep ice core sites and “TAM” to the Transantarctic Mountain Ranges. The basic boundaries of the numbered regions are deﬁned by the veloc-ity reconstruction using radar data and extend into the ice shelf environment. (Adapted from Rignot et al. [2011]).

lowered by some 700 m in the Sarnoff Mountains (Figure 1) during the last 10,000 years of glacial evolution [Stone et al., 2003]. Such evidence gave credence to the idea that regional GIA might generate uplift at a rate of 20 mm/yr, larger than observed in the Gulf of Bothnia (Scandinavia) or Hudson Bay (Canada) [Lidberg et al., 2010; Mazzotti et al., 2011]. In anticipation of two simultaneously operating space missions designed to recover ice sheet mass balance, one measuring gravity changes (Gravity Recovery and Climate Experiment or GRACE) and the other measuring ice height changes (Ice, Cloud, and Land Elevation Satellite or ICESat), Wahr et al. [2000] devised a scheme to use both measurements to decou-ple the two geophysical signals and recover ice mass balance of Antarctica. The paucity of data to support GIA models is a significant source of uncertainty in determination of Antarc-tic ice mass balance from the monthly time-varying GRACE gravity field solutions provided by the various analysis cen-ters [Velicogna and Wahr, 2006]. Lack of a robust sampling by ICESat and the difficulty of obtaining accurate trend data on Antarctic bedrock using GPS have meant that this source of uncertainty has plagued the recovery of an accurate mass balance for almost a decade. The uncertainty is quite large, 50–150 Gt/yr, as reported by Gunter et al. [2009], and could rival the amplitude of the mass balance signal itself.

[4] This paper revisits the secular gravity signal retrieved

for Antarctica via the Center for Space Research (CSR),

University of Texas at Austin, Release 04 (RL04) for time-varying, de-aliased, spherical harmonic (SH) ﬁelds from January 2003 to January 2012. Using a new Antarctic ice retreat model, a viscoelastic GIA response model examines a large variation of mantle viscosity parameters using 240 computations of Earth responses, retrieving an improved GIA correction that is compatible with the GPS data. There-fore, we provide a more reﬁned estimate of Antarctic ice mass balance averaged over this period than has been avail-able to date. Both our past analyses of elastic responses to ice loading/unloading [Ivins et al., 2005] for a present-day glacial mass balance study produced by Rignot and

Thomas [2002] and a more recent analysis by King et al.

[2012] indicate that elastically induced vertical crustal motions can compete with those of GIA models only within the coastal region adjacent to the Amundsen Sea and in the northernmost Antarctic Peninsula.

2. Emergence of New Data

[5] Two new data sets have become available that allow

the GIA correction uncertainty to be lowered. First, a spa-tially and temporally improved GPS data set now exists [Bevis et al., 2009; Thomas et al., 2011], with many

sta-tions providing vertical trends with 1 uncertainty less than

(3)

core data have become available that constrain past ice sheet reconstructions. Critical to the improved constraints are new and widespread cosmogenic dating [Bentley et al., 2010] and mapping of the offshore extent of the AIS at global LGM during global at sea level lowstand, some 19–26 ka [Livingston et al., 2012]. These glacial-geological data are supported by atmospheric and ice ﬂow modeling of deep ice core data throughout the Antarctic interior [Mackintosh

et al., 2011; Buiron et al., 2011]. Whitehouse et al. [2012a]

have assembled all the pertinent glacial-geological data and used a shallow ice numerical code to examine a full spectrum of possible ice change scenarios. In addition, the relatively small number of relative sea level data available at coastal sites were modeled by Whitehouse et al. [2012b] using a comprehensive global GIA model from Mitrovica and Milne [2003]. Whitehouse et al. [2012b] used variations of the ice model to compare to recently published GPS data [Thomas

et al., 2011] and determined an optimum ice model: W12a.

[6] In this paper we revisit the IJ05 model [Ivins and

James, 2005] and develop a revision, IJ05_R2, which

incor-porates largely the same new geochronological constraint data as given in Table 1 of Whitehouse et al. [2012a]. However, the approach to interpreting the ice thickness and timing data assumes both the “thickest” and “youngest” interpretation of LGM through Holocene evolution, in some contrast to the Whitehouse et al. [2012a] approach. (A full list of the differences in approach by Whitehouse et al. [2012a] and that undertaken here is given in the online sup-porting material S4.2 of Shepherd et al. [2012].) The new ice model, IJ05_R2, is forced to have no ice changes after 1 ka. What is demonstrated in this paper is that new GRACE GIA corrections can be generated by comparing predicted uplift from forward modeling, with the revised ice model, to ﬁeld GPS uplift data. A 9 year time series of monthly GRACE solutions (2003.0 to 2012.0) is then used to re-evaluate Antarctic mass loss averaged over this period.

2.1. Data: Glacial Geology

[7] The model IJ05 presented a reconstruction of the ice

load since the LGM based on Antarctic glacial geology and present-day mass balance estimates. By comparison to ear-lier models, such as the maximum model proposed in the CLIMAP (Climate: Long-range Investigation, Mapping, and Prediction) project of the late 1970s with total meltwater

volume 1.0 107_km3

(similar to the ICE-3G model of

Tushingham and Peltier [1991]), the ice load model IJ05

incorporated substantial improvements in Antarctic recon-struction by folding in regional constraints available from ice cores, paleogrounding line positions on the continental shelf, moraine positions, and rock exposure dating. The net effect was to lower the total meltwater volume expulsion to near0.4 107km3.

[8] In Figure 2 we show the AIS volume changes as a

function of time since LGM (21 ka) for a number of

con-temporary models including M11 [Mackintosh et al., 2011], W12a [Whitehouse et al., 2012b], and PDcPaSEb [Briggs, 2012], which incorporate local geological constraints, and the ICE-5G (version 1.2) [Peltier, 2004] that uses limited regional constraints and ice dome histories inferred from far-ﬁeld sea level models and records. The most recent new model using the regional glacial chronology constraints is by

Briggs and Tarasov [2013] and uses a sophisticated scoring

0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 time (ka) e. s. l. ( m) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 11 12 13 14 15 16 17 18 19 20 total ICE-5G Briggs (2012) PDcPaSEb

total IJ05_R2 IJ05_R2 (East Ant.) total W12a

total M11 M11 (East Ant.)

(total)

Figure 2. Ice history of IJ05_R2, ICE-5G, W12a, M11, and

PDcPaSEb from Briggs [2012] in terms of volume changes as equivalent mean sea level changes (e.s.l.).

system to evaluate the performance of numerical simulations versus the data and error statistics. What is important to note from Figure 2 is that all of the more recent models tend to be smaller in total volume of meltwater by a factor of more than 2 in comparison to ICE-5G. The revised thickness history in ice model IJ05_R2 assumes a “maximal” interpretation of the paleo-ice thickness and area extent records, as is clearly seen in comparisons to other models in Figure 2. The only exception where more volume change occurs in one of the new models is in the comparison of IJ05_R2 to the period 21 to 9 ka for model M11. Thus, our interpretation will lead to larger predictions of overall uplift and therefore provide an the upper bound on the newly derived GIA corrections for GRACE determination of mass balance trends.

2.2. Reconstruction of Ice Loading

[9] Over the last decade, there has been a

substan-tial increase in the number and distribution of constraints on ice load history. Important sites where new data are located include Mac. Robertson Land, Talos Dome, Dome C, Dome A, and Dome Fuji (Figure 1). Interpretation of this new data set suggests that the volume of ice at the LGM was considerably smaller than previously thought [Bentley et al., 2010; Mackintosh et al., 2011] and the IJ05 model was mod-iﬁed accordingly. The following three sections provide an overview of the new information included in the revision. A view of the grid caps that are used is given in section A and shown in Figure A1.

[10] A major breakthrough in understanding ice

thick-ness changes in East Antarctica comes from detailed snow accumulation modeling of ice core annual layers. For the Talos Dome region, for example, estimates are that the max-imum elevation changes occurred at 12.7–17 ka, and are limited to about 170 m of ice height differential (IHD)

(4)

[Buiron et al., 2011]. These constraints, and the new expo-sure dates on erratics in coastal mountains of East Antarctica and those in the central TAM (Figure 1), are important to the revision.

2.2.1. Dronning Maud Land (15ıW) to Mac. Robertson Land (75ıE)

[11] Mackintosh et al. [2007, 2011] presented new data

and model for the region 75ıE to 55ıE constraining

past ice thickness from the continental shelf to at least 500–600 km inland to be lower than that in the IJ05 model.

To the west of the Napier Mountains (66ı_{S, 56}ı_{E), both}

coastal emergence data and 10_Be, 26_{Al exposure dating of}

erratics indicate that past ice thicknesses were at least 350 m thicker than at present and that substantial retreat occurred between about 11 and 7 ka [Yamane et al., 2011]. This new evidence is one of the few regions indicating that more local mass transport from continent to ocean should exist than in IJ05. Further to the west in Dronning (Queen) Maud

Land in the Heimefrontfjella Mountain Range (74ı_{S, 14}ı_W;

see Figure 1), radiocarbon-dated snow petral nest materi-als restrict the proximal ice thickness increase at LGM to lower elevations [Thor and Low, 2011]. These data require the regional IJ05 model to be reduced in thickness by 200–300 m.

2.2.2. Mount Manthe (100ıW) to the Wisconsin Range (135ıW)

[12] New evidence in West Antarctica is found in the

central Transantarctic Mountain Ranges (TAM), close to the geographic South Pole. Ackert et al. [2007, 2011] demonstrated that near the Ohio Mountain Ranges, LGM ice was about 125 m higher than at present, and this increased thickness persisted until about 9–10 ka. Todd

et al. [2010] dated moraines in the Wisconsin Ranges

(Figure 1) and found that ice reached a maximum height of roughly 250 m higher than at present between about 17 and 14 ka. These new constraints lower the regional IHD values in the new model (IJ05_R2) by about 70–200 m. At Siple Dome, there has been about 350 m of thinning between 15 and 14 ka [Price et al., 2007]. Near the coast, new data obtained at Mount Manthe and adjacent mountain outcrops indicate that deglaciation was underway by about 14 ka [Johnson et al., 2008]. This information only slightly changes the IHD magnitudes of the IJ05 model, reducing past thickness by 50–75 m. More importantly, the new data from Johnson et al. [2008] show that the region was actively changing up to at least 2.2 ka. In sum, IJ05_R2 lowers IHD values adjacent to the Ross Sea in the TAM by roughly 30–180 m. We should also note that the vast inte-rior of East Antarctica adjacent to the TAM remains rather poorly constrained.

2.2.3. Weddell Sea (15ıW) to Bellingshausen Sea (90ıW)

[13] A considerable amount of new information has been

collected in the mountain ranges and nunataks adjacent to glaciers that ﬂow to the Weddell Sea. For example, Hein

et al. [2011] ﬁnd that in the Shackleton Range, exposure

dating reveals a thickening of no more than 165 m at the site where erratics were sampled. We therefore reduce the IJ05 thicknesses at LGM by 250–350 m in the region. New data from the Ellsworth Mountain Ranges [Bentley

et al., 2010] play a pivotal role in lowering the

esti-mate of past ice volume in the eastern Weddell region.

Table 1. Grounded Volume of Deglaciating Ice (105km3₎a

Region Time (kyr B.P.)

15 11.5 7.6 6.8 3.2 2.2 Wilkes B. 3.14 2.81 2.22 1.56 0.0 0.0 (4.62) (4.14) (3.26) (2.30) (0.0) (0.0) Lambert 0.99 0.82 0.76 0.65 0.12 0.01 (1.45) (1.21) (1.12) (0.96) (0.18) (0.13) DML 4.15 3.33 2.24 0.75 0.0 0.0 (6.08) (4.90) (3.29) (1.11) (0.00) (0.00) Weddell 8.13 5.45 4.36 2.76 0.73 0.06 (11.95) (8.01) (6.41) (4.06) (1.08) (0.75) Ant. Pen. 3.51 2.61 1.66 1.78 0.31 0.03 (5.16) (3.84) (2.45) (2.62) (0.46) (0.33) Amundsen 1.93 1.08 0.87 0.72 0.20 0.03 (2.59) (1.37) (1.10) (0.90) (0.25) (0.18) Ross Sec. 4.60 3.60 2.84 1.76 0.49 0.05 (6.76) (5.29) (4.18) (2.59) (0.72) (0.51) EAP 2.62 5.23 5.64 3.01 0.0 0.0 (3.85) (3.78) (3.26) (2.30) (0.0) (0.10) Total 29.05 24.93 20.59 12.99 1.84 0.27 (42.97) (32.54) (23.91) (15.42) (2.78) (1.90)

a_{Values in parentheses are for the original IJ05 model. DML}₌_Dronning

Maud Land, EAP=East Antarctic Plateau.

The erratics analyzed by Bentley et al. [2010] suggest that regional ice was thicker by 230–480 m at 15 ka, and conse-quently, the IJ05 model thicknesses are reduced in IJ05_R2 by 250–380 m.

[14] Analysis of sedimentary core data in the Weddell Sea

and Eltanin Bay by Hillenbrand et al. [2010, 2012] reveals that regional deglaciation via ice streams ﬂowing from the mountains of the southernmost Antarctic Peninsula, north of

Ellsworth Land and south of the Behrendt Mountains (74ıS,

74ı_{W), was continuous from about 11 to 4 ka and that retreat}

on the inner shelf of the southern Weddell Sea was under-way as late as 8 ka. While this puts no constraint on past ice thicknesses, it supports a storm- and accumulation-driven climatic environment of the past ice sheet, wherein post-Antarctic cold reversal Paciﬁc westerlies allowed thicker ice than at present to be maintained well after 13 ka [Johnson

et al., 2008; Livingston et al., 2012]. However, in the

inte-rior, just west of the Rutford ice stream and the Ellsworth Mountain Ranges, the mass changes are minimal from about 7 ka to present [Ross et al., 2011].

[15] The differences between IJ05 and IJ05_R2 are shown

in the maps of Figures 3a–3c. HereıIHD is the difference

in IHD between the version of the history model of Ivins

and James [2005] and the new model introduced here. The

yellow and red regions of Figures 3a–3c show where the new model has ice heights that are increased in value rela-tive to IJ05, whereas the green to blue colors show where the new model paleo-ice heights are decreased. Figure 3d shows the IHD values for the revision at 15 ka, a time near the maximum Antarctic LGM expansion in the model. In Figure 3, the spherical harmonic representation of caps (also see Figure A1) is truncated at degree and order 256 and forms the forcing function for the spherical harmonic solution of the differential equation system used for GIA predictions. Using the eight regions deﬁned in Ivins and

James [2005], Table 1 compares the ice volumes to be lost

to sea level rise at six time slices between the old (values in parantheses) and new model.

(5)

-30˚

Figure 3. (a–c) Differences in model ice heights (ıIHD) between the original IJ05 [Ivins and James, 2005] and the revised model presented here (IJ05_R2) for three time slices. (d) New model IJ05_R2

heights at 15 ka. Here “IJ05_R0” refers to the original IJ05. The quantity is the AIS total collapse

volume represented as its global sea level rise equivalent (SLE) (Figure 3d). In Figures 3a–3c,ıis the

incremental change in that value at the various time slices in a one-to-one comparison between R0 and R2 models. We note that the use of “IJ05_R0” implies an intermediate model “R1,” which has been used in some unpublished work but not reported in the literature.

3. Ice to Solid Earth Surface Loading

[16] Many of the technical aspects of modeling the

pertur-bation of the solid Earth away from gravitational equilibrium by surface ice and ocean loading have been described in

James and Ivins [1998] and Ivins and James [2005]. The

main assumptions that underlie the ice history reconstruction IJ05_R2 are summarized in this section.

[17] A collection of 455 spherical caps of varying radii

(see section A and Figure A1) forms the continental and shelf loading grid. The height of each cap is the estimated IHD, with a partial isostatic adjustment being achieved by computing one glacial cycle of 100 ka. Over a surface that is currently covered by ocean, or by an ice shelf, loads are set to zero if, at the time of maximum glacial extent, the grounded ice mass did not exceed a value sufﬁcient to have produced

an incrementally larger radial stress load at the ice-bedrock interface than at present. All caps shown in Figure A1 have a nonzero ice load at some time in the history used for the GIA model, though only that grounded portion of the past load above freeboard of present-day sea level is incorporated. A series of tests performed by Simon et al. [2010] give con-ﬁdence to the accuracy of this approximation. In this way, much of the past mass lost from regions that are currently marine-based do not contribute directly to a change in past global mean sea level.

[18] Geochronological data and past grounding line

posi-tion data, where available, are applied directly to each grid at 9 times after the LGM. The time slices in IJ05_R2 may be seen in Figure 2. Linear extrapolations are applied between time slices. Further linear interpolations for each

(6)

disk are used to decimate the model into 1 ka increments for GIA runs. The load is set to zero at 1 ka. The model lacks both ﬂow dynamics and explicit accumulation ﬂux and hence is quasi-static. Temporal variability in accumulation is accounted for when possible.

[19] In the interior of East Antarctica, and far from ice

core sample locations, there are no constraints on IHD val-ues. Such regions include portions of numbered basins 8–16 and 18–19 (Figure 1) where there are vast areas having no rock outcrops. Here we assume ice domes, and not depres-sions, that reach ice ﬂow equilibrium values. However, no massive ice domes are assumed (i.e., that approach 1 km in excess ice thickness) as these cannot be supported by modeled accumulation records from ice core data [Parrenin

et al., 2007, Neumann et al., 2008, Verleyen et al., 2011, White et al., 2011, Siddall et al., 2011, Buiron et al., 2011].

[20] The model reconstruction is biased to large IHD

values at each of the 10 time slices. The bias allows the model to be used as a way of upper bounding the volume among competing ice forward models, while obeying the constraints provided by the IHD data. In Mac. Robertson Land (Figure 1), many of the IHD constraints are located within 100 km of one another, and each may be discordant by a factor of 2, or more. During the next decade, availabil-ity of new sites, where either additional rock exposure dating can be recovered or new deep ice core data are extracted, holds the promise of additional reﬁnements for the load history.

3.1. Data: GPS

[21] A great deal of effort has gone into obtaining

uplift data on bedrock in Antarctica for constraining GIA [Tregoning et al., 2000; Raymond et al., 2004; Dietrich and

Rülke, 2008; Wilson et al., 2008]. Because of the time span

of these records, and because of various improvements in geodetic reference frames, hardware standards, processing strategies, and time series analyses, GPS vertical motion trends are now retrieved from the data for analysis of GIA [Bevis et al., 2009; Ivins et al., 2011; Argus et al., 2011]. The most comprehensive analysis of these data for Antarctica is the report by Thomas et al. [2011] using 39 individual station records spanning the period 1995–2010. Here we use 18 of

these trend data and their 1errors, eliminating records from

the Antarctic Peninsula north of 72ıS, due to the difﬁculty of

dealing with both large elastic and transitional viscoelastic [e.g., Ivins et al., 2000; Simms et al., 2012] signals present there. We average together the values from stations located within 100 km of one another and eliminate some stations with reported errors greater than the amplitude of the signal.

4. GIA Modeling of the Solid Earth Deformation 4.1. Model Computation Architecture

[22] A spherical code that includes ocean loading,

self-gravitation, and a layered mantle has been used for the computations in the revised predictions from the load model. A simple incompressible Earth model with four Maxwell viscoelastic layers above the core-mantle boundary (CMB) are retained, along with a ﬂuid core. The elastic param-eters and density have been chosen to approximate the PREM model [Dziewonski and Anderson, 1981] and are constrained to ensure that buoyancy forces at mantle

inter-Table 2. Earth Structure Parameters for GIA Model

h (km) UM(Pa s) LM(Pa s) (2)LM(km) (1)LM(km)

65, 115 1–10 1020 _0.8–_{80 10}21 ₅₇₀₁ ₄₅₀₁

h(GPa) UM(GPa) (2)LM(GPa) (1)LM(GPa) sg (m/s2)

57.8 74.9 240.0 293.4 9.799

h(kg/m3) UM(kg/m3) (2)LM(kg/m3) (1)LM(kg/m3) cmbg (m/s2)

2890 3900 4701 4988 10.675

faces and the lithosphere are reasonably close to those of PREM. The mathematical structure of the numerical code is in the time domain, using the system of partial differ-ential equations in the spherical harmonic (SH) domain as derived from mass conservation, momentum balance, Poisson equation for the gravitational potential, and linear Maxwell viscoelastic constitutive equations. Each layer is assumed to be homogeneous, a feature that generally ren-ders solutions to the radial-dependency analytic [Vermeersen

et al., 1996]. A matrix system of ordinary differential

equations deﬁned by mantle and mantle-core interface boundary conditions is then solved as an eigenvalue-eigenvector problem as described by Ivins et al. [1993]. Benchmark comparisons to the code used by Simon et al. [2010] demonstrated that the required accuracy for uplift and geoid computations reported here is achieved by this eigenvalue-eigenvector routine. In this paper, the model described here is used to evaluate both the GIA correction for GRACE and the crustal motion uplift.

[23] The fundamental parameters of the study for the GIA

intercomparison are given in Table 2 wherein (i)

LM and (i)

LMare the elastic rigidity and density within theith layer

of the lower mantle and(i)

LMare the radial position to the

top boundary of these layers. Subcripts ‘h’ and ‘UM’

indi-cate lithosphere and upper mantle, respectively. For each model calculation there is a predictive set of quantities:

Pur( , , r, t)andˆP0( , , r, t), wherePuris the radial

compo-nent of the bedrock velocity andˆP0_{the time derivative of the}

external potential ﬁeld caused by the rock motions

through-out the Earth, including at the surface, r = Re. Predicted

values at present day (’tpd’) ofPur( , , Re, tpd)(at co-latitude

and longitude ) may be compared to GPS monitoring of

vertical bedrock motions andˆP0₍_{, , R}

e+ a0, tpd) to

space-craft accelerations at altitudea0. The use of such solutions

with GRACE data is straightforward: Each of the analysis

center releases are in terms of SH coefﬁcients,CN`mj, from

which a secular trend can be determined, PNCS

`mj and can be

directly compared to the expansion coefﬁcients,GIAPNC

`mj, that

are formed from each GIA computation [Ivins et al., 2011]. For completeness, we write out the expansion of the time-varying perturbed gravitational potential about a mean (or “static”) value at the surface of the Earth:

ˆ0( , , Re, t) = GM Re L X `=1 ` X m=0 2 X j=1 N C`mj(t) NY`mj( , ), (1)

where the subscripts `mj on CN are the degree and order,

andsin orcospart of SH coefﬁcient set, respectively. The

(7)

Figure 4. Uplift map of prediction with IJ05_R2. The solid

Earth model assumes lithospheric thicknessh = 115 km,

upper mantle viscosityUM= 21020Pa s, and lower mantle

viscosityLM= 41021Pa s. Station locations used for

eval-uation of GIA model are from Thomas et al. [2011]. Stations north of the southernmost Antarctic Peninsula were elimi-nated from the analysis of model performance here. See text for explanation of the color symbols and note that there is no difference in classiﬁcation of green stars and blue diamonds, as the former are due to the background color contour near the zero value.

harmonics common to geodesy is employed [Tapley et al., 1988; Torge, 1989; Ivins et al., 1993]. For GRACE appli-cations, these coefﬁcients are also directly related to the apparent surface mass redistribution by eliminating the elas-tic solid Earth deformation-related part of the potential by

using the multiplicative factor,1/(1+k0

`), wherek

0

`is the load

Love number [Wahr et al., 1998].

[24] For the purposes of improving our understanding of

the GIA correction for Antarctica, we need only compute

GIAPNC

`mjfor determiningGIAˆP0( , , Re, t)and then calculate

the apparent mass by integrating over a mask defined by boundaries of the target region on the surface. This region represents the space over which the GRACE satellites resolve the gravity change signal. The area overlaps onto the part of the continental shelf having no grounded ice by 100 km, since any “real” or “leaked” signal there is effec-tively indistinguishable from signal that originates from the nearby grounded ice sheet. No tapered mask is employed. The method has been well studied by Swenson and Wahr [2007], and tests of the mask and spherical harmonic fit with degree and order truncation of 60, with Gaussian filter radius

G = 300 km, yield the important scale factor. All of the

GIA computations for the 240 mantle-parameter study are each sampled using this mask and rescaling. Although many

other lithospheric thicknesses were studied, two thickness values were selected as representative, and then eight values of the upper mantle viscosity and 15 values of lower mantle viscosity were computed in the parameter sweep.

4.2. Integration of the Vertical Motion Data With Mantle Parameters

[25] In Figures 4 and 5, two uplift predictions are given

in map view for differing values of lithospheric thickness, one more consistent with the West Antarctica rift

sys-tem (h = 65 km) and the other with the East Antarctic

craton (h = 115 km). The East Antarctic lithosphere is

topped by a thick crust, having a very deep seismic Moho boundary at 40–45 km depth [Hansen et al., 2009]. The mantle portion of the lithosphere beneath this cratonic crust is much thicker. This value of effective elastic thickness is a reasonable guess, yet is rather poorly constrained.

Estimates of lateral variations in seismic Q from core

reﬂections of shear waves (ScS) also indicate a thick Pre-cambrian cratonic lithosphere for East Antarctica [Souriau

et al., 2012], possibly thicker than in the present study.

A new model study of paleo sea level change during the Oligocene growth of the Antarctic ice sheet with compre-hensive viscoelastic loading and self-gravitation favors a

lithospheric thickness in East Antarctica of h 100 km

[Stocchi et al., 2013].

[26] Among the broad range of model parameters for

upper and lower mantle viscosity, forh = 65 km, values

ofUM = 2 1020 andLM = 1.5 1021 Pa s provided

a reasonable set of predictions of the uplift data at the 18 stations selected from Thomas et al. [2011]. For the thicker

lithosphere at h = 115 km, UM = 2 1020 Pa s and

LM = 4 1021Pa s also provided a reasonable set of

pre-Figure 5. Uplift map of prediction with IJ05_R2. Solid

Earth model assumesh = 65km,UM = 2 1020Pa s, and

(8)

Table 3. GPS Station Uplift (mm/yr) With 3000 Days of Observationa

Station MAW1 DUM1 DAV1 SYOG CAS1 VESL

GPS_Pu

r 0.06 ˙ 0.39 –0.79 ˙ 0.46 –0.94 ˙ 0.49 2.6 ˙ 0.36 1.18 ˙ 0.43 1.06 ˙ 0.45a

Obs. period 2000–2010 1998–2010 2000–2010 1995–2010 2000–2010 1999–2010

GIA IJ05 R2b _0.76 _0.60 _0.73 _0.71 _0.22 _0.03

GIA IJ05 R2c _0.50 _1.18 _0.57 _0.51 _0.49 _0.33

a_{Errors reported here are from formal 1}_{estimation. The data analysis is from Thomas et al. [2011].} b_{Fit assumes}_{h = 115}_km,

LM= 4 1021Pa s, andUM= 2 1020Pa s. c_{Fit assumes}_{h = 65}_km,

LM= 1.5 1021Pa s, andUM= 2 1020Pa s.

dictions. The upper mantle structure parameters are slightly lower than those reported for Fennoscandia by Zhao et al. [2012] who determined that a comprehensive regional GPS

data set there is best ﬁt usingh = 93to 110 km andUM= 3.4

to5.0 1020 _{Pa s. The main point of the GIA model}

solu-tions here is that they recover reasonably justiﬁed values by a ﬁrst-order comparison to Fennoscandia. However, it is gen-erally cautioned that the geodetic data for Antarctica lack the spatial density necessary to form more formal solutions for mantle viscosity and lithospheric thickness.

[27] To quantify the misﬁt of the models, we perform

two evaluations. One assessment simply uses the standard

normed2_{assessment which we deﬁne as}

2 ₌ 1 N N X ı=1 GPS_Pu rı – GIAPurı Q ı 2 , (2)

with N = 18. Here the superscripts GIA and GPS

indi-cate model and observation, respectively. However, theQıis

not just the formalized 1 errors reported by Thomas et al.

[2011] but includes an addition of 1.05 mm/yr as an assess-ment of an error level affecting all stations due to elastic deformation or other uncalibrated errors and biases, but still downweights the noisiest uplift values that rely on less than 160 days of observations. Part of the justiﬁcation for uncer-tainties slightly higher than 1 mm/yr, for example, comes

from the total ambiguity of reference frame issues ( ˙0.4

mm/yr) [Collilieux et al., 2011]. The additional uncertainty of 1.05 mm/yr added to the formal errors reported in Thomas

et al. [2011] is not entirely arbitrary. We selected reference

frame drift uncertainty from Collilieux et al. [2011], wherein an 8 year long drift of global GPS and VLBI-based frames is reported at 0.4 mm/yr (see Table 4 therein), doubled the average formal RMS errors of the highest-quality stations (Table 3) to 0.44 mm/yr, and employed the average elastic corrections of Thomas et al. [2011] for those same stations to characterize an additional of 0.21mm/yr. The total, 1.05 mm/yr, is then added as an additional estimate of error at

Table 4. GPS Station Uplift (mm/yr) With 1000 Days of Observationa

Station ABOA BELG SVEA MAIT AV

GPS_Pu

r 1.4 ˙ 0.84 2.97 ˙ 1.47 2.07 ˙ 1.95 0.12 ˙ 0.40a

Obs. period 2003–2010 1998–2005 2004–2008 1996–2009

GIA IJ05 R2b _0.23 _1.11 _0.48 _0.01

GIA IJ05 R2c _0.62 _1.19 _0.64 _0.11

a_{Errors reported here are from formal 1}_estimation. b_Assumes_{h = 115}_km,

LM= 4 1021Pa s, andUM= 2 1020Pa s. c_Assumes_{h = 65}_km,

LM= 1.5 1021Pa s, andUM= 2 1020Pa s.

each station employed in our analysis. It is important to note

that we have included this addition to Qı not because we

feel that Thomas et al. [2011] underreported their errors, but rather because it allows us to make assessments that are more sensitive to the larger uplifting regions of West Antarctica, while downweighting the inﬂuence of the very small, but more certain, values measured along the coast of East Antarctica. The second way of making a quantitative evaluation is simply to compute the quantity:

N ˛2 = 1 N N X ı=1 GPS_u rı – GIAurı 2 , (3)

which ignores errors altogether and which uniformly weights GPS and GIA model uplift rates.

[28] Figure 4 shows the station positions that are used

and a contoured prediction of GIA uplift. Tables 3–6 give

observations with 1errors and the GIA predictions for the

two cases of contoured uplift shown in Figures 4 and 5. The station positions are marked with symbols that indicate their location in the hierarchical Tables 3–6 with blue circles indicating the highest-quality continuous (cGPS) observing locations, while red circles and blue diamonds/green stars indicate those having greater than 1000 and 400 days of observing, respectively. Stations reporting signiﬁcant uplift, but having less than 160 days of observation, are shown with red diamonds. The station data from MBL1 AV (green dia-mond) played a special role, as they were tested for inﬂuence on the overall GIA correction. The impact of removing the station on our retrieval of the GIA correction for Antarctic mass loss is negligible.

5. GIA Model Correction for GRACE

Time-Varying Gravity

[29] Within the context of the recent Ice Mass Balance

Intercomparison Exercise (IMBIE) [Shepherd et al., 2012],

Table 5. GPS Station Uplift (mm/yr) With 400 Days of Discon-tinuous Observationa

Station TNB1 CRAR FOS1 BREN

GPS_Pu

r –0.23 ˙ 0.81 1.00 ˙ 0.65 2.14 ˙ 0.40a 3.85 ˙ 1.60

Obs. period 1999–2008 2002–2010 1995–2010 2006–2010

GIA IJ05 R2b _–0.22 _0.28 _5.11 _2.33

GIA IJ05 R2c _–0.31 _–0.01 _2.50 _2.45

a_{BREN has only 436 days of observations and the other three greater than}

1000 days.

b_Assumes_{h = 115}_km,

(9)

Table 6. GPS Station Uplift (mm/yr) With 160 Days of Discon-tinuous Observationa Station HAAG W05 AV W07 AV MLB1 AV GPS_Pu r 3.47 ˙ 0.71 4.86 ˙ 1.01 3.61 ˙ 1.58 3.28 ˙ 1.09 Obs. period 1996–2006 2002–2008 2002–2006 1998–2003 GIA IJ05 R2b _4.55 _4.05 _4.10 _1.62 GIA IJ05 R2c _3.89 _4.52 _5.05 _1.40

a_{HAAG, W07, W05, and MLB have 22, 56, 153, and 119 observation days,}

respectively.

b_Assumes_{h = 115}_km,

LM= 1.5 1021Pa s, andUM= 2 1020Pa s.

a number of GRACE analyses groups used several test patterns of Antarctic mass loss to determine the scaling factor needed to restore amplitude that is damped by both truncation and Gaussian smoothing [Velicogna and Wahr, 2006; Swenson and Wahr, 2007]. The scaling factor for

the algorithms used in this paper is 1.256 and the error in retrieval is estimated to be of order 2%.

5.1. Results of the Parameter Study

[30] In the parameter study, with GPS constraints, we

were able to determine a new set of optimum GIA

correc-tions to GRACE. The values of the square roots of˛N2_,2_,

and of the mass rate predicted by the GIA model, GIA_M_P_,

are shown in three frames, Figures 6a–6c, respectively, con-toured for all 120 model parameter computations, for both the h = 65 km andh = 115 km cases, respectively. The

minima in the mean of the deviation of model and data˛N

(Figures 6a and 6b) correspond to the map views of the predictions shown in Figures 4 and 5.

[31] The regions with contour minima in Figures 6a and

6b for the quantity ˛N near the white or light pink color

(˛ 0.2N ) reveal those parts of the upper and lower man-tle viscosity space best satisfying the GPS uplift data, and

21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 20.0 20.1 20.2 20.3 20.4 20.5 20.6 20.7 20.8 20.9 21.0 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 (b) (d) (f) 20.0 20.1 20.2 20.3 20.4 20.5 20.6 20.7 20.8 20.9 21.0 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8

mean model difference compared to uplift data

Log₁₀η_LM Log 10 ηUM Log 10 ηUM Log₁₀_ηLM Log₁₀η_LM Log₁₀_ηLM Log 10 ηUM Log₁₀η_LM Log₁₀η_LM 0.66 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8

χ2_{fit to uplift data.}

(a) (c) (e) 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 21.0 21.2 21.4 21.6 21.8 22.0 22.2 22.4 22.6 22.8 h = 65 km λ = 300 km_G

GIA correction in Gt/yr 20.0 20.1 20.2 20.3 20.4 20.5 20.6 20.7 20.8 20.9 21.0 20.0 20.1 20.2 20.3 20.4 20.5 20.6 20.7 20.8 20.9 21.0 h = 115 km 20.0 20.1 20.2 20.3 20.4 20.5 20.6 20.7 20.8 20.9 21.0 20.0 20.1 20.2 20.3 20.4 20.5 20.6 20.7 20.8 20.9 21.0

Figure 6. GIA model predictions using ice loadIJ05_R2 with (a, c, e)h = 65 km and (b, d, f)h =

115km. GIA solutions that have a Gaussian ﬁlter of radiusG= 300km, rescaling factor 1.256, and mask

are all applied. Figures 6a and 6b, 6c and 6d, and 6e and 6f correspond top˛N2_,_{, and}GIA_M_P_{, respectively,}

(10)

Table 7. GIA Correction Comparisona

GIA Model GIA_M_P _(Gt/yr)

ICE-5G VM2 (IC) –126

ICE-5G VM2 (C) –128

IJ05_R2ave.˛N115 km+65 km –53

IJ05_R2ave.˛Nfor all 1 –52

IJ05_R2min.˛N65 km –45

a_{ICE-5G assumes VM2 mantle structure.}

therefore deﬁning the range of permissible GIA corrections

to GRACE (Figures 6e and 6f). By selecting this˛N

crite-rion for “goodness-of-ﬁt,” each of the subset of “best” model predictions maintained two characteristics: The rebound amplitudes are each within 0.75 mm/yr of the values of

max-imum amplitude reported in Tables 3–6 and the2 _values

are less than approximately 0.8. The permissible range of

GIA corrections is from as low asGIA_M_P _{= 37.5 Gt/yr (}_{h = 65}

km case in Figure 6e) to as high asGIA_M_P _{= 67.5 Gt/yr (}_{h =}

115km case in Figure 6f). The limit of the permitted range

was tested by allowing different ranges of allowed˛N from

0.2 to 0.3 mm/yr. Such liberalization of the limit does not change our inference of GIA corrections for GRACE trends by more than 5–8 Gt/yr, and we account for this by adding this as part of our error estimate.

[32] Although now clearly in need of revision [e.g., Wang

et al., 2013], a ubiquitous GIA model used for correcting

GRACE ice mass balance trends has been ICE-5G VM2 [Peltier, 2004, 2009; Riva et al., 2009; Wu et al., 2010]. In Table 7, we compare ICE-5G GIA corrections [A et al., 2013] (computations provided by Geruo A) to those from

IJ05_R2. For ICE-5G, the Earth structure is identical to

the published VM2, ocean loading with solution of the sea level equation is included, and the modern rotational feed-back theory [Mitrovica et al., 2005] is employed with “C” and “IC” representing assumptions of a compressible and incompressible viscoelastic rheology. Three different

repre-sentations are in Table 7 for IJ05_R2 are referenced: (i)

averaging the best˛N ﬁt to GPS between theh = 65km and

115 km cases, (ii) all˛Nﬁts that also lie within 1of the data

in the2_{test, and (iii) the best}_˛_N _{ﬁt to the}_{h = 65}_{km case.}

The comparison shows that theIJ05_R2 model reduces the

GRACE GIA correction to roughly 42% of the ICE-5G VM2 values.

5.2. GRACE Analysis

[33] The GRACE monthly time series for Level 2 data

has now been extended into calendar year 2012. Here we use monthly releases from CSR Release 04 for January 2003 to January 2012 to estimate a 9 year trend in mass balance for all of Antarctica. The focus on the entire ice sheet avoids having to treat basin-to-basin leakage issues [Horwath and Dietrich, 2009] and more tedious rescaling formulations. Leakage from outside Antarctica is estimated

to be2.6 ˙ 2.0 Gt/yr using terrestrial water storage

simu-lated with the Global Land Data Assimilation (GLDAS-2) hydrology model and the ECCO-2 ocean model.

[34] The monthly SH coefﬁcients CNMc_`mj for Release 04,

corrected for elastic Earth deformation, are analyzed for trend using standard time series analysis and decomposed as follows using least squares [Bloomﬁeld, 2000]:

N CMc_`mj(t) = NCo_`mj+ PNCS_`mj t + K X =1 N C_`mj cosh2 ft + `mj i , (4)

where the ﬁrst term is a simple offset, NCo

`mj, the

coefﬁ-cient superscripted “S” is the “secular” rate, and the sum

on is formed from the decomposition of semiannual,

annual,S2, andK2tidal alias terms [Ray and Luthcke, 2006]

plus any additional known frequencies,f. Here the

super-script Mc indicates the monthly release coefﬁcients. (As

explained below, the CSR harmonic coefﬁcients are aug-mented slightly to recover the mass balance). For each SH

(`mj), these amplitudes,CN_`mj, and phases,_`mj, are solved

for and separated from the analysis for secular trends. A

principal goal of the time series analysis is to obtain PNCS

`mj

with as little corruption from tidal aliasing and interan-nual oscillations as possible. A pure trend component of the potential can then be deﬁned as

P ˆ0( , , Re) = GM Re L X `=1 ` X m=0 2 X j=1 PNCS `mj NY`mj( , ) . (5) WithGIAPNC

`mj from model calculation with theIJ05_R2 ice

model and PNCS

`mj from the analysis of the monthly ﬁelds,

we are now ready to compareGIA_M_P _with_M_P _{retrieved from}

GRACE. Assuming that all secular change originates from

surface mass, a surface density rate,P, may be related to the

secular coefﬁcients PNCS `mj: P ( , ) = N Re L X `=1 ` X m=0 2 X j=1 2` + 1 3 PNC S `mj NY`mj( , ) , (6)

withNthe mean density of the Earth [Ivins et al., 2011]. An

Antarctic grid mask,Aki, at positionsk, i, is then applied

over the globe, with the sum taken over each area

incre-ment, ıAki, to compute the mass change rate. The mask

extends seaward from the ice sheet grounding line by 150 km

to avoid under-sampling. The GIA equivalent rate, GIA_M_P_,

is computed simply by replacing PNCS

`mj withGIAPNC`mj in the

equation above and performing the same procedure. To recover the correct mass gain/loss for Antarctica from the GRACE SH ﬁelds, the scaling must be used and this includes GIA predictions when they are applied to the gravity data.

5.2.1. Additional Nuances of the GRACE Time Series

[35] The GRACE satellites are insensitive to

center-of-mass, center-of-ﬁgure shifts, and so contain no` = 1

spheri-cal harmonics. This lowest odd-parity harmonic, however, is important to closing the long wavelength mass budget, espe-cially for a clean separation of north-south hemispheric mass variability. A method developed by Swenson et al. [2008] uses the GRACE ﬁelds for isolating time-varying surface

loads, then reconstructs the` = 1mode from that. We

aug-ment the CSR GRACE ﬁelds by including monthly values

of those terms provided by Sean Swenson. Adding the` = 1

terms in this way approximately accounts for changes that are occurring both far from Antarctica, and on the continent

as well. The` = 1 mode has also been discussed recently

by Wu et al. [2012] and Rietbroek et al. [2012]. Additional parts of the low-order ﬁeld are also required to be restored. Antarctica is especially sensitive to the lowest zonal

(11)

CSR Jan 2003 to Jan 2012 trends for GRACE Release 04 in mm per yr WHE 0˚ 30˚ 60˚ 90˚ 120˚ 150˚ 180˚ 210˚ 240˚ 270˚ 330˚ -160 -150 -140 -130 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40

Trend WHE mm/yr

300˚ 300˚

Figure 7. GRACE trend withG = 300km, rescaling, annual, semiannual,K2,S2 tidal aliasing

fre-quencies removed, and all degree 1–4 restorations as discussed in the text. Note the three (burnt carmine)

centers of positive mass balance with contours at+30 mm/yr water height equivalent (WHE). Symbol

markings as in Figure 4.

harmonics, because of its area extent and position at the South Pole [James and Ivins, 1997].

5.2.2. Partial Restoration of the Low-Order Field

[36] The coefﬁcients for degree 1 used are from an

anal-ysis by Swenson et al. [2008] extended into 2012. We per-formed a time series analysis on this, similar to equation (4), employing the relations given in Swenson et al. [2008] for

conversion from frame trend (PX, PY, PZ) to SH coefﬁcients

(yielding PNCS 111= –0.27 10 –11_yr–1_, PNCS 112= –0.36 10 –12_yr–1_, PNCS 100 = –1.99 10

–11_yr–1_{). The} _{X, Y, Z}_{time series is}

care-fully corrected for leakage from non-Antarctic sources and iteratively freed from any GIA model dependence. (This is the same procedure employed for the GRACE intercom-parative study used in Shepherd et al. [2012]). In addition, we employ information from satellite laser ranging (SLR) studies [Cheng and Tapley, 2004; Cheng et al., 2013] for

properly restoring the low-order zonals PNCS

200, PNC S 300, PNC

S 400, as

recommended by Bettadpur [2007], and in GRACE Tech-nical Note 05 (update as of 12 June 2012) for the required

N

C200 time series. We also compared results by replacing

these with the trends of Maier et al. [2012], and this had

approximately a˙1.5 Gt/yr effect on estimated mass

bal-ance. Differences between our Antarctic solutions for degree 1 and the global values published by Wu et al. [2010] (over

a shorter a time span) are a˙2.5 Gt/yr effect. However, the

overall effect is not small. If we naively replace our Antarctic

time series for PNCS

100, PNCS111, and PNCS112over 2003.0–2012.0 with

the global analysis of Rietbroek et al. [2012] for 2003.0–

2009.0 (determining a value of PNCS

100 = –3.35 10–11yr–1,

or 65% larger), this would lower the absolute value of the estimated mass imbalance by 18.8 Gt/yr. Larger negative zonal degree 1 coefﬁcients predict larger Antarctic mass gain, in comparison to those of smaller negative ampli-tude. In fact, this is true of all satellite-determined very low

odd-degree zonalPNCS

l00[James and Ivins, 1997].

[37] The Antarctic mass balance sensitivity to the secular

trend of the degree 2 tesseral coefﬁcients (PNCS

211,PNC S

212) is

con-siderably smaller, probably smaller than about 1 Gt/yr for competing choices for independently determined rates. This can be veriﬁed by using different rate values selected from

Cheng et al. [2011] for the time span given (2002–2010), the

one that is most similar to 2003–2012. The results reported

here assume PNCS

211, PNCS212 that follow from the

documenta-tion of Bettadpur [2007] as updated on the website, ftp:// podaac-ftp.jpl.nasa.gov/allData/grace/L2/CSR/RL04/docs/.

5.2.3. GRACE Mass Balance Without GIA Correction

[38] Using the mask, rescaling, and a time series

analy-sis with no search for an ENSO-like interannual term (2.6 year period), and no search for tidal terms with aliasing periods longer than 4 years, we determined that the mass balance for Antarctica over the time span 2003.0–2012.0 is

+2.9 ˙ 29Gt/yr. The formal 1errors coming from the diag-onal terms of the covariance matrix are smaller than the total 29 Gt/yr reported here as error. Added are errors due to de-aliasing model dependency, differences among least squares solutions for trend, rescaling, and low-order ﬁeld reconstruc-tions, and these collectively amount to 13 Gt/yr. Figure 7 shows the pure GRACE trend in water height equivalent (WHE) that is recovered. Figure 8 shows two frames of map view GIA correction, one on Figure 8a is a WHE

(12)

predic--4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 180˚ 150˚W 120˚W 90˚W 60˚W 30˚W 0˚ 30˚E 60˚E 90˚E 120˚E 150˚E 180˚ 150˚W 120˚W 90˚W 60˚W 30˚W 0˚ 30˚E 60˚E 90˚E 120˚E 150˚E

WHE rate in mm/yr for IJ05_R2 GIA prediction

(a) (b)

mm/yr

Figure 8. GIA trend predictions in WHE with G = 300 km and rescaling 1.256, identical to

Figure 7. (a)h = 115km, with stiffer mantleUM = 8 1020Pa s andLM = 1 1022. (b)h = 115km,

with softer mantleUM = 4 1020Pa s andLM = 1.5 1021. Stiffer mantle predicts larger total

cor-rection (GIA_{M = 88}_P _{Gt/yr), but is unconstrained by GPS data in Figure 8a. Figure 8b shows a prediction}

with smaller total correction (GIA_{M = 48}_P _{Gt/yr), but is constrained by GPS data. An appropriate map of}

ice mass change with GIA correction (in terms of WHE) is to subtract map of Figure 8b from the map of Figure 7. The reader should note that the rescaling constant was tested for a number of conﬁgurations of gain/loss across Antarctica, though not for any smaller standalone region or drainage basin, such as the Antarctic Peninsula.

tive map that has no GPS constraint, whereas the case on Figure 8b is selected from among those that are compatible with the GPS data. Note that strong differences exist in both positive and negative parts of the color contours. Even in the case of larger GIA (Figure 8a), the largest corrections are generally less than half of the pure GRACE trends.

5.2.4. The Potential Role of New GPS Data and Improved GRACE Dealiasing Models

[39] No reliable GPS secular trend data are available

along the coastal rock outcrops of the Amundsen Sea sector (Figure 1, region 5) where dramatic ice mass loss is observed by GRACE (Figure 7). Using two campaigns (2006 and 2010) at three sites with continuous observations spanning from 4 to 16 days, Groh et al. [2012] determined the uplift to exceed that predicted by the elastic response to regional ice loss and estimated the possible increase in GIA correction to be about 30 Gt/yr, in contrast to either ICE-5G or the orig-inal IJ05 model and making the GRACE inferred ice loss more negative. While we have not incorporated such data into our analysis, the study by Groh et al. [2012] demon-strates the importance of retrieving more bedrock trend data in West Antarctica. The GRACE analysis centers have now released version 05 of the spherical harmonic coefﬁcients, and a preliminary estimate of the Antarctic trend for the time span 2003–2012 is that ice imbalance is more negative by roughly 20 Gt/yr than with version 04, used in this paper. The difference likely owes to improvements in the dealiasing models used to produce the Level 2 data product. Analyses with the improved ocean models for version 05 were con-sidered in the analysis of Shepherd et al. [2012]. Using the average of the IJ05_R2 and W12a GIA corrections, the later

study reported–81 ˙ 33Gt/yr for Antarctica during January

2003 to December 2010, which is generally consistent with results reported in this paper.

6. Conclusions

[40] Measurement of continent-wide GIA is needed to

interpret satellite-based trends for the grounded ice mass change of the Antarctic ice sheet (AIS). This is especially true for trends determined from the Gravity Recovery and Climate Experiment (GRACE) satellite mission. Three data sets have matured to the point where they can be used to shrink the range of possible GIA models for Antarctica: The glacial-geological record has expanded to include

expo-sure ages using10_Be,26_{Al measurements that constrain past}

thickness of the ice sheet, modeled ice core records now better constrain the temporal variation in past rates of snow accumulation, and Global Positioning System (GPS) verti-cal rate trends from across the continent are now available. Here we have derived a new GIA model predicting both the new GPS data and a new GIA correction for GRACE trend.

[41] The revised Antarctic mass balance is lowered by

60–120 Gt/yr relative to the past analyses of GRACE time series. The exact value of the reduction in the GIA depends on the details of the GRACE analysis method. Chen et

al. [2009], for example, used the CSR RL04 time series

and GIA model IJ05 [Ivins and James, 2005] and found

Antarctica to be changing by–190 Gt/yr during 2002.25–

2009.0 and noted that the loss rate would increase to–250

Gt/yr by using the ICE-5G GIA correction. For a shorter time series (2002.6–2007.0) and using CSR RL04, Peltier [2009] determined positive mass balance for Antarctica

(13)

85 85 85 85 85 85 85 85 42 30 6 6 11 12 149 149 129 85 280 187 85 17 6 6 6 11 17 17 17 149 258 237 211 235 187 ₁₈₇ 410 365 267 17 17 11 6 6 6 11 12 6 11 17 149 74 238 225 196 199 191 119 443 456 463 349 114 341 170 11 6 6 6 6 11 11 6 6 6 6 11 17 17 129 138 307 237 546 201 200 51 45 430 682 197 341 375357 121 341 204 11 11 6 6 6 17 136 136 17 6 6 6 6 6 17 17 3 367 174 172 363 221 204 166 25 59 54 157 512 497 478 478 444 546 478 248 180 17 11 6 6 6 6 16 100 212 138 17 6 6 6 6 6 6 17 17 17 17 93 93 93 264 337 303 204 208 206 235 244 238 324 576 446 450 333 299 478 470 444 375 150 17 6 6 6 17 10 207 207 109 109 212 43 6 6 6 6 6 6 6 7 17 17 17 17 17 93 179 254 375 409 426 422 337 144 119 68 68 142 238 444 344 378 344 93 179 303 176 269 ₄₄₃ 342 40 250 254 307 341 378 232 409 273 244 62 274 281 300 67 6 8 17 8 6 6 8 17 18 54 51 17 11 247 247 247 247 247 247 247 1 247 247 247 247 371 247 247 315 315 247 133 161 247 168 223 315 223 223 247 420 380 371 247371 247 189 303 409 312 408 408 78 281 267 349 170 68 17 28 28 51 33 34 30 34 95 170 341 307 223 371 223223223₃₄₉371223 349 349 274 274₂₇₂ 349 349 349 264267233 233 233 267 446 233 267 267 233 267 267267 1 267 267 336 274 1 1 267 267 267 301 267 567 542 437 336 349 571 272 272 64 375 331 297 307 341 314 272 248 341 375 306 371 349 371 267 1 1 1 341 1 255 255 1 1 378 341 1 497 480 493 459 562 497 478 463 470 341 206 225 ₂₂₅ 223 202 204 179 179 179230 179 179 179 179 179 179 110 110 145 110 110 110 110 316 93 341 403 375 403 403 230 196 230 307 452 314 341 170 204 136 314 136 102 44 230 230 230 44 44 44 44 23 23 44 IJ05_R2 Ice Height Differential (IHD) in meters at 21 ka

Figure A1. Cap locations and ice height differential to present-day forIJ05_R2 at 21 ka. The spherical caps shown are approximated in the plot with octagonal boundaries.

obtain –111 Gt/yr, implying an ICE-5G GIA correction of

133 Gt/yr. The upper bound on the GIA correction deter-mined in this paper is less than half this value. Both the new GPS uplift data and the wealth of new data available to con-strain past ice sheet volume and collapse history are key to coming to this conclusion. The new estimate for Antarctic

mass balance by GRACE during 2003.0–2012.0 is–57 ˙ 34

Gt/yr with the improvedIJ05_R2-GIA correction. Similar

conclusions concerning the lower Antarctic contribution to sea level rise during the GRACE era have also followed

from the implementation of the W12a model by King et

al. [2012] and by Sasgen et al. [2012], who used a hybrid

GPS-GIA model. The reduction in mass contribution to sea level rise from Antarctica is notable. Indeed, for future dis-cussions of closure scenarios for global water transport, it

will be important to incorporate thisŠ0.23 mm/yr

equiva-lent reduction in Antarctic GIA correction for GRACE, as it is an important component of the total sea level change bud-get that has been so well recorded over the past two decades [e.g., Masters et al., 2012].

(14)

Appendix A: Ice Loading Grid

[42] An example of the spherical grid caps that are used is

shown in Figure A1. The online supporting information pro-vides the complete grid model history and guide for imple-mentation as an ice history for geodetic and geophysical investigations.

[43] Acknowledgments. The research support for E.R.I. and F.W.L. comes from NASA’s Earth Science Division with grants from both the Cryosphere Program and the Earth Surface and Interior Focus Area as part of the GRACE Science Team effort. The work of E.R.I. and F.W.L. was performed at the Jet Propulsion Laboratory, California Institute of Tech-nology. J.W. was partially supported by NASA grants NNX08AF02G and NNXI0AR66G and by NASA’s Making Earth Science Data Records for Use in Research Environments (MEaSUREs) program. Robert Briggs is thanked for discussions and for providing data from his PhD Thesis. We thank Geruo A for calculating the ICE-5G VM2 GIA responses and three anonymous reviewers for thorough and helpful reviews. This is also a contribution of the Climate Change Geoscience Program of the Earth Sciences Sector (ESS) of Natural Resources Canada. T.S.J. and K.S. gratefully acknowl-edge support from the ArcticNet Networks of Centres of Excellence. The authors acknowledge Paul Wessel and the University of Hawaii for General Mapping Tools [Wessel and Smith, 1998]. This is ESS contribution 20120237.

References

A, G., J. Wahr, and S. Zhong (2013), Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: An application to glacial isostatic adjustment in Antarctica and Canada, Geophys. J. Int.,

192, 557–572, doi:10.1093/gji/ggs030.

Ackert, R. P., S. Mukhopadhyay, B. R. Parizek, and H. W. Borns (2007), Ice elevation near the West Antarctic Ice Sheet divide dur-ing the Last Glaciation, Geophys. Res. Lett., 34, L21506, doi:10.1029/ 2007GL031412.

Ackert, R. P., S. Mukhopadhyay, D. Pollard, R. M. DeConto, A. E. Putnam, and H. W. Borns (2011), West Antarctic Ice Sheet elevations in the Ohio Range: Geologic constraints and ice sheet modelling prior to the last highstand, Earth Planet. Sci. Lett., 307, 83–93.

Argus, D. F., G. Blewitt, W. R. Peltier, and C. Kreemer (2011), Rise of the Ellsworth mountains and parts of the East Antarctic coast observed with GPS, Geophys. Res. Lett., 38, L16303, doi:10.1029/2011GL048025. Bettadpur, S. (2007), GRACE 32-742, UTCSR Level-2 processing

stan-dards document for Level-2 Product Release 0004, 27 February 2007.

(Available at ftp://podaac.jpl.nasa.gov/allData/grace/docs).

Bentley, M. J., C. J. Fogwill, A. M. Le Brocq, A. M. Hubbard, D. E. Sugden, T. J. Dunai, and S. P. H. T. Freeman (2010), Deglacial history of the West Antarctic Ice Sheet in the Weddell Sea embayment: Constraints on past ice volume change, Geology, 38, 411–414.

Bevis, M., et al. (2009), Geodetic measurements of vertical crustal velocity in West Antarctica and the implications for ice mass balance, Geochem.

Geophys. Geosyst., 10, Q10005, doi:10.1029/2009GC002642.

Bevis, M., et al. (2012), Bedrock displacements in Greenland manifest ice mass variations, climate cycles and climate change, Proc. Nat. Acad. Sci.,

USA, 109, 11,944–11,948, doi:10.1073/pnas.1204664109.

Bloomﬁeld, P. (2000), Fourier Analysis of Time Series, pp. 261, Wiley and Sons, New York.

Briggs, R. (2012), Antarctic ice sheet evolution over the last glacial cycle: A data-constrained large-ensemble modelling approach, PhD Thesis, Dept. Physics and Phys. Oceanography, Memorial Univ. Newfoundland, Canada. pp. 391.

Briggs, R., and L. Tarasov (2013), How to evaluate model derived deglacia-tion chronologies: A case study using Antarctica, Quat. Sci. Rev., 63, 109–127, doi:10.1016/j.quascirev.2012.11.021.

Buiron, D., et al. (2011), TALDICE-1 age scale of the Talos Dome deep ice core, East Antarctica, Clim. Past, 7, 1–16.

Chen, J. L., C. R. Wilson, D. D. Blankenship, and B. D. Tapley (2009), Accelerated Antarctic ice loss from satellite gravity measurements,

Nature Geosci., 2, 859–862.

Cheng, M., and B. D. Tapley (2004), Variations in the Earth’s oblate-ness during the past 28 years, J. Geophys. Res., 109, B09402, doi:10.1029/2004JB003028.

Cheng, M., J. C. Ries, and B. D. Tapley (2011), Variations of the Earth’s ﬁgure axis from satellite laser ranging and GRACE, J. Geophys. Res.,

116, B01409, doi:10.1029/2010JB000850.

Cheng, M., B. D. Tapley, and J. C. Ries (2013), Deceleration in the Earth’s oblateness, J. Geophys. Res. Solid Earth, 118, 740–747, doi:10.1002/jgrb.50058.

Collilieux, X., L. Métivier, Z. Altamimi, T. van Dam, and R. J. Ray (2011), Quality assessment of GPS reprocessed terrestrial reference frame, GPS

Solutions, 15, 219–231.

Conrad, C. P., and B. H. Hager (1995), The elastic response of the Earth to interannual variations in Antarctic precipitation, Geophys. Res. Lett.,

22(23), 3183–3186, doi:10.1029/95GL03176.

Conway, H., B. L. Hall, G. L. Denton, A. M. Gades, and E. D. Waddington (1999), Past and future grounding-line retreat of the West Antarctic Ice Sheet, Science, 286, 280–283.

Dietrich, R., and A. Rülke (2008), A precise reference frame for Antarctica from SCAR GPS Campaign: Data and some geophysical implications, in Geodetic and Geophysical Observations in Antarctica, edited by A. Capra, and R. Dietrich, pp. 1–10, Springer, Berlin.

Dziewonski, A. M., and D. L. Anderson (1981), Preliminary reference Earth model, Phys. Earth Planet. Inter., 25, 297–356.

Groh, A., H. Ewert, M. Scheinert, M. Fritsche, A. Rülke, A. Richter, R. Rosenau, and R. Dietrich (2012), An investigation of glacial isostatic adjustment over the Amundsen Sea Sector, West Antarctica, Global

Planet. Change, 98-99, 45–53, doi:10.1016/j.gloplacha.2012.08.001.

Gunter, B., T. Urban, R. Riva, M. Helsen, R. Harpold, S. Poole, P. Nagel, B. Schutz, and B. Tapley (2009), A comparison of coincident GRACE and ICESat data over Antarctica, J. Geodesy, 83, 1051–1060.

Hansen, S. E., J. Juliá, A. A. Nyblade, M. L. Pyle, D. A. Wiens, and S. Anandakrishnan (2009), Using S wave receiver functions to estimate crustal structure beneath ice sheets: An application to the Transantarctic Mountains and East Antarctic craton, Geochem. Geophys. Geosyst., 10, Q08014, doi:10.1029/2009GC002576.

Hein, A. S., C. J. Fogwill, D. E. Sugden, and S. Xu (2011), Glacial/ interglacial ice-stream stability in the Weddell Sea embayment, Antarc-tica, Earth Planet. Sci. Lett., 307, 211–221.

Hillenbrand, C. D., R. D. Larter, J. A. Dowdeswell, W. Ehrmann, C. Ó. Cofaigh, S. Benetti, A. G. C. Graham, and H. Grobe (2010), The sed-imentary legacy of a palaeo-ice stream on the shelf of the southern Bellinghausen Sea: Clues to West Antarctic glacial history during the Late Quaternary, Quat. Sci. Rev., 29, 2741–2763.

Hillenbrand, C-D., M. Melles, G. Kuhn, and R. D. Larter (2012), Marine geological constraints for the grounding-line position of the Antarctic Ice Sheet on the southern Weddell Sea shelf at the Last Glacial Maximum,

Quat. Sci. Rev., 32, 25–47.

Horwath, M., and R. Dietrich (2009), Signal and error in mass change inferences from GRACE: The case of Antarctica, Geophys. J. Int., 177, 849–864.

Ivins, E. R., and T. S. James (2005), Antarctic glacial isostatic adjustment: A new assessment, Antarct. Sci., 17, 541–553.

Ivins, E. R., C. G. Sammis, and C. F. Yoder (1993), Deep mantle viscous structure with prior estimate and satellite constraint, J. Geophys. Res., 98 (B3), 4579–4609, doi:10.1029/92JB02728.

Ivins, E. R., C. A. Raymond, and T. S. James (2000), The inﬂuence of 5000 year-old and younger glacial mass variability on present-day crustal rebound in the Antarctic Peninsula, Earth Planets Space, 52, 1023–1029. Ivins, E. R., E. Rignot, X. Wu, T. S. James, and G. Casassa (2005), Ice mass balance and Antarctic gravity change: Satellite and terrestrial per-spectives, in Earth Observation With CHAMP: Results From 3 Years in

Orbit, edited by Ch. Reigber, H. Luhr, P. Schwintzer, and J. Wickert, pp.

3–11, Springer-Verlag, Berlin. ISBN: 3-540-22804-7.

Ivins, E. R., M. M. Watkins, D.-N. Yuan, R. Dietrich, G. Casassa, and A. Rülke (2011), On-land ice loss and glacial isostatic adjustment at the Drake Passage: 2003–2009, J. Geophys. Res., 116, B02403, doi:10.1029/2010JB007607.

James, T. S., and E. R. Ivins (1995), Present-day Antarctic ice mass changes and crustal motion, Geophys. Res. Lett., 22(8), 973–976, doi:10.1029/94GL02800.

James, T. S., and E. R. Ivins (1997), Global geodetic signatures of the Antarctic ice sheet, J. Geophys. Res., 102, 605–633.

James, T. S., and E. R. Ivins (1998), Predictions of Antarctic crustal motions driven by present-day ice sheet evolution and by isostatic memory of the Last Glacial Maximum, J. Geophys. Res., 103, 4993–5017.

Johnson, J. J., M. J. Bentley, and K. Gohl (2008), First exposure ages from the Amundsen Sea Embayment, West Antarctica: The Late Quaternary context for recent thinning of Pine Island, Smith and Pope Glaciers,

Geology, 36, 223–226.

Kaula, W. M. (1980), Problems in understanding vertical movements and Earth rheology, in Earth Rheology, Isostasy and Eustasy, edited by N-A. Mörner, pp. 577–589, J. Wiley, New York. ISBN: 0-4710-27593-X. King, M. A., R. J. Bingham, P. Moore, P. L.Whitehouse, M. J.Bentley, and

G. A. Milne (2012), Lower satellite-gravimetry estimates of Antarctic sea-level contribution, Nature, 491, 586–589, doi:10.1038/nature11621.