Click
Here

JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 114, B06406, doi:10.1029/2008JB006161, 2009

for

Full
ArticleGlacial isostatic adjustment and nonstationary
signals observed by GRACEPaul Tregoning,1 Guillaume Ramillien,2,3 Herbert McQueen,1 and Dan Zwartz1
Received 19 October 2008; revised 27 March 2009; accepted 28 April 2009; published 23 June 2009.

[1] Changes in hydrologic surface loads, glacier mass balance, and glacial isostatic

adjustment (GIA) have been observed using data from the Gravity Recovery and Climate
Experiment (GRACE) mission. In some cases, the estimates have been made bycalculating a combination of the linear rate of change of the time series and periodic
seasonal variations of GRACE estimates, yet the geophysical phenomena are often notstationary in nature or are dominated by other nonstationary signals. We investigate thevariation in linear rate estimates that arise when selecting different time intervals of
GRACE solutions and show that more accurate estimates of stationary signals such as GIA
can be obtained after the removal of model-based hydrologic effects. We focus on North
America, where numerical hydrological models exist, and East Antarctica, where suchmodels are not readily available. The root mean square of vertical velocities in North
America are reduced by 20% in a comparison of GRACE- and GPS-derived uplift rates
when the GRACE products are corrected for hydrological effects using the GLDAS
model. The correlation between the rate estimates of the two techniques increases from0.58 to 0.73. While acknowledging that the GLDAS model does not model all aspects of
the hydrological cycle, it is sufficiently accurate to demonstrate the importance ofaccounting for hydrological effects before estimating linear trends from GRACE signals.
We also show from a comparison of predicted GIA models and observed GPS uplift
rates that the positive anomaly seen in Enderby Land, East Antarctica, is not a stationary
signal related to GIA.
Citation: Tregoning, P., G. Ramillien, H. McQueen, and D. Zwartz (2009), Glacial isostatic adjustment and nonstationary signals
observed by GRACE, J. Geophys. Res., 114, B06406, doi:10.1029/2008JB006161.

1. Introduction[2] One aim of the Gravity Recovery and Climate Experiment (GRACE) mission [Tapley et al., 2004] is to identify
the gravity signal associated with glacial isostatic adjustment (GIA) of continental regions that were ice covered
during the Last Glacial Maximum. Satellite altimetry and
GRACE space missions detect not only the present-day
mass change but also the remnant GIA signal. It is therefore
not possible to estimate the component of present-day massbalance changes without first removing the GIA signal.
[3] On decadal timescales, GIA causes essentially linearchanges in the geopotential over the affected regions. Thus
attempts have been made to extract linear variations fromexisting GRACE solutions in order to estimate GIA signals[e.g., Tamisiea et al., 2007]. In particular, a positive gravity
rate anomaly seen in Enderby Land, East Antarctica [Chen1
Research School of Earth Sciences, Australian National University,
Canberra, ACT, Australia.
2
Laboratoire d’Etudes en Géophysique et Océanographie, Centre
Nationale de la Recherche Scientifique, Toulouse, France.
3
Now at Dynamique Terrestre et Planétaire, Centre Nationale de laRecherche Scientifique, Toulouse, France.

Copyright 2009 by the American Geophysical Union.
0148-0227/09/2008JB006161$09.00

et al., 2006; Ramillien et al., 2006; Lemoine et al., 2007]
has been interpreted as either an unmodeled GIA uplift or
recent snow accumulation [Chen et al., 2006] or, more
recently, as reflecting errors in the GIA model or being
related to snow accumulation [Chen et al., 2008]. Despite
this uncertainty in geophysical origin of the signals, linear
rates are estimated from GRACE time series [e.g., Chen etal., 2008].
[4] Hydrologic signals are typically cyclic in nature, withdominantly near-annual periods [Schmidt et al., 2008].
Hydrologic processes cause significant variations in land
surface heights and both good [e.g., Davis et al., 2004] andpoor [e.g., van Dam et al., 2007] agreements have beenfound between estimates from GRACE and the Global
Positioning System (GPS). However, interannual trends
associated with, for example, droughts can cause multiyear
increases/decreases and departures from simple periodic
variations. If sufficiently large, such nonstationary variations (that is, variations that do not repeat in a regular,
predictable pattern) in surface mass will be present inmonthly GRACE solutions in addition to long-term, stable
trends such as GIA. Global hydrological models such as
GLDAS [Rodell et al., 2004] have been used to modelsurface and soil moisture signals on broad scales [e.g., Syedet al., 2008] and to mitigate the effects on GRACE
estimates of GIA [Tamisiea et al., 2007].

B06406

1 of 10

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

B06406[5] The aim of this paper is to demonstrate the nonstationary nature of many of the observed GRACE signals andto caution against the process of simply estimating a linear
trend from monthly GRACE fields. We focus on two regions
where the GIA signal is not well constrained by ice models(East Antarctica and Laurentia) and show that accurate
secular trends can be estimated from GRACE after removing
the dominant hydrologic signals, thereby permitting the longterm GIA signals to be identified. We also show from a
combination of observed GPS uplift rates and predicted geoid
rates from numerical models that the positive anomaly seen inEnderby Land is most likely not related to GIA.
[6] This paper is divided into three main parts. Firstly, we
describe the GPS, GLDAS and GRACE data sets andrelevant components of the analysis and utilization of thetechnique products. Next, we look in detail at the observed
GIA signals in North America, a region covered by the
GLDAS hydrology model, to demonstrate that a large
portion of the observed GRACE signals are nonstationary
and to highlight the importance of correcting for hydrologic
effects before estimating linear rates from GRACE. Then,being mindful of the fact that hydrologic effects can affectsignificantly linear rate estimates from GRACE, we investigate the unexplained positive anomaly in Enderby Land,
East Antarctica. Given that there is no available hydrologic
model (such as GLDAS) that covers the Antarctic continent,
we invoke different analytical techniques to distinguishbetween GIA and hydrological causes for the observed
GRACE signals and conclude that the latter is more plausible.

2. Vertical Deformation From GRACE
Observations[7] The observed GRACE anomalies are a combinationof components related to both elastic and viscoelastic
effects, and the Stokes coefficients of the temporal sphericalharmonic fields are the sum of the two effects.
dC nm ðt ÞdS nm ðt Þ¼

dC enm ðt Þ
dS enm ðt Þ
þdC vnm ðt ÞdS vnm ðtÞð1Þ

where n, m are the degree and order, dCnm, dSnm are the
Stokes coefficients of the GRACE anomaly fields at time t,
and the superscripts e and v refer to the elastic and viscoelastic components, respectively. At monthly timescales,anomalies caused by surface loads can be considered togenerate only elastic deformation. The elastic component of
the anomalies observed by GRACE is [Wahr et al., 1998]dC enm ðt Þ
dS enm ðtÞ0

3r 1 þ k n¼ w
rav 2n þ 1^ nm ðt ÞdC
d S^nm ðt Þð2Þ

where rw is the density of fresh water, rav is the averagedensity of the Earth, k0n are elastic Love loading numbers^ nm, d S^nm are dimensionless[Pagiatakis, 1990] and d C
Stokes coefficients that represent the surface load anomalies
at time t. We can use the coefficients dCenm, dSenm to calculatevertical elastic deformation [Davis et al., 2004].
N X
n
X

0

hn
1
þ
kn0
ð3Þ
n¼0 m¼0 ee Pnm ðcosqÞ dCnm
ðt Þcosml þ dSnmðt Þsinml

dU e ðq; l; t Þ ¼ R

B06406

where R is the mean radius of the Earth (6371 km), Pnm are
the fully normalized Legendre functions and h0n, k0n areelastic Love loading numbers [Pagiatakis, 1990] computedfor the Preliminary Reference Earth Model (PREM)[Dziewonski and Anderson, 1981].
[8] Long-term viscoelastic deformation associated with
GIA can be approximated by [Wahr et al., 2000]
dU v ðq; l; t Þ ¼ R

N Xn
X
2n þ 1

2
ð4Þ
 v

v
ðt Þsinml
 Pnm ðcosqÞ dCnm ðt Þcosml þ dSnm
n¼0 m¼0

where, in this case, dCvnm, dSvnm are the coefficients representing the viscoelastic deformation components ofthe gravity changes observed by GRACE at time t. Thus in
the elastic case, we multiply by h0n/(1 + k0n) whereas in theviscoelastic case we multiply by (2n + 1)/2.
[9] It is not possible from the GRACE spherical harmonic
fields alone to separate the elastic and viscoelastic components of the total Stokes coefficients. The surface load
anomalies vary in an unpredictable manner with the changes
in the hydrological cycles, hence are nonstationary innature, whereas the viscoelastic signals associated with
GIA will be essentially linear on timescales of years to
decades.

3. Data
3.1. GPS[10] It has been demonstrated that observations with the
Global Positioning System (GPS) can detect GIA uplift
patterns [e.g., Milne et al., 2001; Lidberg et al., 2007]. We
analysed data from 80 globally distributed GPS sites
with the GAMIT/GLOBK software [Herring et al., 2008]
in a global analysis following, for example [Tregoning
and Watson, 2009]. We used the VMF1 mapping function[Boehm and Schuh, 2004; Boehm et al., 2006] with a priori
zenith hydrostatic delays from ray-tracing through the
ECMWF global weather model and we applied nontidal
atmospheric pressure loading at the observation level. We
transformed daily free-network solutions onto the terrestrialreference frame using 30 globally distributed stabilizationsites. We used coordinates for the stabilization sites fromthe ITRF2005, although we first removed a 1.8 mm/yr
Z-translation velocity, being the difference between the
ITRF2000 and ITRF2005 reference frames [Altamimi et al.,2007]. This provides a reference frame with the internalconsistency of the ITRF2005 (Z. Altamimi, personal communication, 2007) without the apparently geophysically implausible translation rate along the Z-axis [Argus, 2007].
[11] We estimated the GPS uplift rate at the Antarctic IGS
sites at Mawson, Davis and our remote, solar-powered siteat Richardson Lake in Enderby Land, East Antarctica aspart of our global analysis (site locations are indicated on
Figure 3). The latter site, affixed firmly to exposed bedrock,
was installed in January 2007 specifically to estimate thepresent-day uplift rate in light of the positive gravity
anomaly rate that had been identified in GRACE solutions[Chen et al., 2006; Ramillien et al., 2006]. To date, 153 daysof data have been transmitted from the site by satellite communications, spanning 1.2 years.

2 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS[12] The most comprehensive estimate of vertical velocities in the North American region, from a combination of
permanent and campaign-based observations, has been
undertaken by Sella et al. [2007]. Their published velocities
agree to within 1 mm/yr of our estimates at common sites(not shown); thus, we adopt their velocity estimates forthe assessment of GRACE-derived GIA in North America(see section 5.1).
3.2. GRACE[13] We used the monthly GRACE solutions from theGroupe de Recherche de Géodesie Spatiale as described by
Lemoine et al. [2007]. These solutions use a 10-day sliding
window approach to analyze the data in 30-day segments,
with the first 10 days weighted by 0.25, the middle 10 days
weighted by 0.5 and the last 10 days weighted by 0.25. An
empirical constraint approach is used when estimating thespherical harmonic coefficients which negates the need toapply subsequent spatial filtering techniques. These solutions have been used previously to estimate, for example
mass balance in Antarctica and Greenland [Ramillien et al.,2006], co-seismic deformation [Panet et al., 2007], hydrologic [e.g., Ramillien et al., 2004] and oceanic signals [e.g.,
Lombard et al., 2007; Tregoning et al., 2008]. We calculated Stokes coefficient anomalies by subtracting the meanvalue over the time series 2002.6 to 2007.6 to derive valuesof dCnm, dSnm of equation (1).
3.3. Hydrology[14] We used the GLDAS global hydrologic model[Rodell et al., 2004] to generate estimates of changes incontinental water mass, incorporating four soil moisturelayers, accumulated snow quantity and total canopy storage.
Previous studies have shown agreement between GRACE
and GLDAS estimates of hydrologic signals [Syed et al.,2008], although the GLDAS model does not include Antarctica nor valid modeling of Greenland [Rodell et al., 2004].
Using the NOAH10 3-hourly values from 2002 to 2008, wecomputed the mean value at each 1° space grid node, thensubtracted it from each 3-hourly field to generate anomaly
values at each node (we set the anomaly values to zero overGreenland). We then generated monthly anomaly estimatesat the epochs of the GRACE solutions, using the same
weighting procedures invoked when the GRACE solutions
are generated (see previous section). We converted theseestimates into spherical harmonic coefficients, which are the^ nm, d S^nm coefficients of the right-hand side of equationdC(2), then converted them into elastic Stokes coefficients,
dCenm, dSenm, using equation (2).
[15] The hydrology signals will affect the GRACE estimates of gravity change through the potential of thechanging surface loads themselves and through the elastic
deformation of the surface caused by the loads [e.g., Wahr etal., 1998]. The latter effect will also affect the GPS estimateof station coordinates: since the GPS antennae are fixed
to the surface of the Earth, they will detect any elastic
deformation that might occur. We can thus utilise GPS
coordinate estimate anomalies to validate the GLDAS modelover the Laurentide region. The local deformation recorded
by a GPS site will only be well modeled by the GLDAS
model if the dominant signals are long wavelength; how-

B06406

ever, since the same is true of the relation between GRACE
and hydrologic signals, it is an interesting comparison tomake.
[16] We convolved the hydrology anomalies (in terms ofequivalent water height) with elastic Greens functions[Farrell, 1972] using Load Love numbers from thePREM model and generated deformation anomalies using
equations (2) and (3) (Figure 1). These agree well in phasewith the GPS height anomalies at sites located around
Churchill Bay in Canada, a region of significant GIA. The
amplitude of the deformation from the GLDAS model
underestimates the movement detected by GPS.
[17] We estimated the scale factor per site that would need
to be applied to the GLDAS elastic deformations in order tomatch the GPS height anomalies and found that it varied
considerably from site to site with a range from 1.1 to 2.7(Figure 2). The differences may be caused by local hydrological phenomena that affect individual GPS sites but are
not captured by the 1° space GLDAS model or deficiencies
in the GLDAS model (which does not include all hydrological processes, e.g., groundwater), or perhaps errors in
the ocean tide loading model (we used the FES2004 model)that alias to low-frequency signals. Lacking any better
alternative, we used the GLDAS to model the hydrologic
signals and assume that it captures the majority of the longwavelength components in Canada. While this is not a
perfect approach, the spatial variation in the scale factor
estimates and the scarcity of permanent GPS sites in thisregion of Canada implies that we cannot verify whether
modifying the GLDAS model through some spatially averaged scaling process would lead to a better result.

4. Linear Rate Estimates[18] We first generated GRACE time series of verticaldeformation, without removing the hydrologic component,
to derive least squares estimates of linear rates and amplitudes/phases of annual periodic signals. This is the approachthat has been used to study geoid rate signals [e.g., Tamisiea
et al., 2007; Davis et al., 2008] as well as the GIA of
Antarctica [e.g., Chen et al., 2006, 2008]. Figure 3a shows
the resulting linear trends obtained using the whole timeseries. Some particular significant features are immediately
apparent, for example the mass loss occurring in West
Antarctica [e.g., Rignot et al., 2008] and southern Greenland [e.g., Luthcke et al., 2006]. Other signals are either
present when perhaps not expected (e.g., the positive signal
in East Antarctica), are not clearly present or are not of the
expected magnitude (e.g., the GIA signals in North America
and Fennoscandia).
[19] We next considered rate estimates using subsets ofthe total GRACE data to investigate the spatial stability ofthe observed signals. We first separated the time series into
two, utilizing 3 years of data for each estimate. Figures 3band 3c show the rate field for 2002.6 – 2005.6 and 2004.6 –2007.6, respectively. A comparison of these fields showsthat some of the signals seen in Figure 3a are always
present, indicating the stationary nature of the signal. The
mass loss in West Antarctica decreases slightly in magnitude
while there is a 50% increase in mass loss in southeast
Greenland. The mass loss in southern Alaska appears to

3 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

B06406

Figure 1. Comparison of residual height variations estimated from GPS (pink) and calculated byconvolving GLDAS hydrologic anomalies into elastic vertical deformation (blue) at GPS sites Baker
Lake (BAKE), Churchill Bay (CHUR), and Scherrerville (SCH2), all in eastern Canada (see Figure 2).

have decelerated, while mass loss may have commenced innorthwestern Greenland to the south of Thule. Many ofthese features are either not visible or are significantly
damped in the rate estimates from the whole time series.
These are the geophysical interpretations that could be
made, but we are not suggesting that these are validinterpretations since seasonal hydrological signals still remain in these GRACE results.
[20] To separate the stationary, viscoelastic effects fromthe nonstationary, hydrologic effects, we subtracted fromthe monthly GRACE spherical harmonic anomaly coefficients (dCnm, dSnm) the monthly GLDAS coefficients (i.e.,the dCenm(t), dSenm(t) terms of equation (1)), making theassumption that the hydrological signals captured by the
GLDAS model will account for all elastic-related surface
mass variations detected by GRACE. The remaining signals
should be related mainly to the nonelastic processes of GIA,
but will also include any unmodeled hydrological or oceanic
elastic signals. We can then generate trend estimates fromtime series of values for each 10-day epoch on a global grid

4 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

B06406gauge observations [e.g., Snay et al., 2007] as well as thefree-air gravity signal from GRACE observations [Tamisiea
et al., 2007]. The predicted uplift signal over Laurentia
computed from GRACE solutions (2002.6 – 2007.6) whenignoring the possible presence of elastic effects is shown in
Figure 5a. Three positive uplift regions are visible, alongwith a zone of subsidence around (E250° space N59°). In

Figure 2. Scale factor estimates at GPS sites in North
America between observed GPS height anomalies andelastic deformation calculations using water load anomalies
from the GLDAS model.

using equation (4), thus deriving viscoelastic vertical deformation as observed by GRACE. This approach is more
rigorous than that of Tamisiea et al. [2007] who removed
the effects of hydrological changes when computing free-air
gravity trends over northern Canada by subtracting from the
GRACE trends a trend estimated from GLDAS data (onedoes not obtain the same velocity estimate from (1) estimating a velocity from raw GPS heights, separately fromGLDAS vertical elastic deformation computations anddifferencing the velocities; and (2) correcting the GPS
heights for the GLDAS-derived elastic deformation, thenusing the corrected time series to estimate the velocity).
Figure 4 shows time series of the original GRACE solutions, the GLDAS elastic deformation predictions and
GRACE minus GLDAS at Flin Flon in southern Canada(54.7°N, 258.0°E). By propagating the formal uncertainties
of the GRGS monthly GRACE fields, we calculated that theuncertainty of each monthly estimate of vertical position
anomaly is around 15– 25 mm (Figure 4), leading to aformal uncertainty of the rate estimates of 1 mm/yr. We
consider this to be an optimistic estimate of the uncertainty
since it does not include any estimate of error in the
GLDAS fields. We acknowledge also that any errors in
the GLDAS representation of the hydrological effects will
still remain in our uplift rate estimates; however, there is not
currently a superior approach available to that which wehave used here.

5. Discussion
5.1. Laurentia[21] GIA models for North America predict significant
vertical uplift rates, with the magnitude depending on the
choice of ice history and viscosity models [e.g., Peltier,2004]. Numerous studies have observed uplift from GPSdata [e.g., Calais et al., 2006; Sella et al., 2007] and tide

Figure 3. Linear trend estimates (in terms of geoid rate)using GRACE solutions spanning (a) 2002.6 to 2007.6,(b) 2002.6 to 2005.6, and (c) 2004.6 to 2007.6. GPS sites in
Antarctica are indicated as D for Davis, M for Mawson, andR for Richardson Lake.

5 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

Figure 4. (a) Time series of surface height (in terms ofvertical viscoelastic deformation) at 54.7°, 258.0°E, using
GRACE-only solutions. (b) Time series of GLDAS
anomalies, expressed in terms of equivalent water height.
(c) As for Figure 4a but with the nonstationary hydrologicalsignals removed using the GLDAS model.

general, the pattern does not correspond well with thatpredicted from GIA models [e.g., Peltier, 2004].
[22] Interannual hydrologic variations that occur in North
America create nonstationary signals in the GRACE observations, which contaminate the purely stationary GIA
signals. Accounting first for the hydrological signals in
the GRACE spherical harmonic anomalies, then estimating
the linear trends yields a significantly different result(Figure 5b). The positive zone around (54°N, 260°E) is
reduced in magnitude while the zone of subsidence at 250°E

B06406is eliminated completely. Neither of these signals feature in
GIA models for Laurentia and are not related to the GIA
processes. The negative signal below 50°N at around 260°E
is more likely to be remaining hydrologic signal not captured by GLDAS than to be related to a peripheral bulge
signal of GIA.
[23] The uplift pattern of the ICE-5G (VM2) model has
three notable centers of maximum uplift, associated with thethree primary ice dome complexes in that model: one overKeewatin near Yellowknife, one in southeast Hudson Bay,and the third in the Foxe Basin to the west of Baffin Island[Peltier, 2004, Figure 21]. Our uplift estimates show clearly
the first two regions but, because the hydrology of Greenland is poorly represented in the GLDAS model [Rodell etal., 2004], we do not show on Figure 5 nor discuss here thethird center.
[24] The pattern of uplift west of Hudson Bay has amaximum of 12 mm/yr and is slightly elongated in an E-W
direction. There are a number of interesting features of thiscenter: firstly, there is a clear separation between thismaximum and that of the center to the southeast of Hudson
Bay. Secondly, the extension of the region toward Yellowknife is supportive of substantial ice loss in this region.
Available shoreline data in the Canadian region cannot
provide constraints on the likely thickness of the ice
shield in this region at the Last Glacial Maximum (LGM)(K. Lambeck, personal communication, 2008), nor are there
any GPS sites in the vicinity from which to estimate directly
the ongoing uplift. Thus this estimate from GRACE solutions provides new spatial observations of present-day
signals. The shape of the pattern is in general agreementwith the hydrology-corrected free-air gravity anomaly mapof Tamisiea et al. [2007]. A clear positive region appears
southeast of Hudson Bay, corresponding with the locationof a dome in the ICE-5G (VM2) model [Peltier, 2004,
Figure 21] but with a larger amplitude (17 mm/yr compared with 14 mm/yr). This zone of uplift is also present
in the models of Lambeck et al. [2008].
[25] Predictions of uplift rate from GIA models are
dependent on both the ice history and viscosity models
used; therefore, the usefulness of a comparison of GRACEderived and GIA model-derived uplift rates is limited by theaccuracy of each of the models, where the true accuracy ofeach is unknown. We compare our GRACE-derived velocities with GPS-derived uplift rates of Sella et al. [2007]
covering the region. The latter estimates are independent of
ice model and viscosity and so provide a better alternative
for assessing the accuracy of the two GRACE models. The
RMS differences between our GRACE velocity fields andthe GPS uplift rates are 3.7 mm/yr and 3.0 mm/yr for
GRACE-only and GRACE-GLDAS, respectively, showinga 20% improvement in agreement with the GPS velocities
once the hydrologic signals are removed. There is also an
improvement in the correlation between the GPS andGRACE velocities (0.58 for GRACE-only, 0.73 for
GRACE-GLDAS) (Figures 5b and 5d). A similar patternof uplift rates has recently been found by van der Wal et al.
[2008] in a study of GIA in North America.
5.2. Enderby Land, Antarctica[26] The above investigation demonstrates that not accounting for nonstationary hydrological effects can signif-

6 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

B06406

Figure 5. (a) Uplift rates from GRACE for the Laurentide region including hydrologic signals and(c) with hydrologic signals removed. Vertical velocities estimated from GPS observations by Sella et al.
[2007] are shown using the same color scheme. (b) Comparison of GRACE-derived and GPS uplift
velocities including hydrologic signals and (d) with the GLDAS model removed from each monthly field.
The locations of Yellowknife (Y) and Flin Flon (F) are indicated (red squares) in Figure 5a.
icantly bias the linear rate estimates of GIA. Despite any
possible shortcomings of the GLDAS model, the fact thatsuch models exist at all for North America make suchstudies possible. On the other hand, no such models are
publicly available that cover the Antarctic continent so thesame approach cannot be used to obtain the GIA estimatesfrom GRACE over Antarctica. In this section we invoke adifferent approach to study positive anomaly signals whichhave been detected in Enderby Land, East Antarctica. We
compare viscoelastic modelling and a comparison with
GPS-derived uplift rates to assess whether GIA or interannual snow accumulation are the most likely cause of thesignals observed by GRACE.
[27] Chen et al. [2008] estimated a rate of change of 33.7 ±0.65 mm/yr of equivalent water height in the center of thepositive anomaly feature in Enderby Land, confirming theearlier detection of a positive anomaly in this region [Chenet al., 2006; Ramillien et al., 2006]. Figure 6 shows (in red)a time series at this location (point D in Figure 7a) from the
GRGS GRACE solutions. This is comparable to the timeseries shown by Chen et al. [2008, Figure 5b] except that
we express the gravity field changes as vertical deformation
rather than surface load in terms of equivalent water height.
Estimating the uplift rate from the last three rather than thefirst three years of the total time series changes the rateestimate from 19.3 ± 2.3 mm/yr to 2.5 ± 2.1 mm/yr, with anaverage rate over the entire series of 10.3 ± 1.1 mm/yr(Figure 6). The rates are not constant, indicating that thetime series is dominated by a signal(s) other than linear
GIA. Possible interpretations of the time series are that no
significant GIA is occurring (the variations in the time series
are related to variations in accumulated snow/ice) or that

7 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

B06406matches the GPS-observed surface uplift rates. This arises
because one cannot have a geoid rate signal without having
an associated uplift of the continental surface.

6. Conclusions[31] Nonstationary hydrologic signals distort estimatesof long-term GIA from GRACE data if a linear rate isestimated through uncorrected GRACE fields. Removing
the signals using a global hydrologic model enables theextraction of stationary, linear signals, which are moreconsistent with ground-based estimates of vertical uplift

Figure 6. Time series of gravity variations (expressed asvertical deformation using equation (4)) at point D in
Figure 7a. Uplift rates (in mm/yr) shown are estimated fromthe first three years (blue), the last three years (orange), andthe entire time series (red).
destructive interference has occurred in the latter part of thetime series between an underlying uplift signal and significant hydrologic mass loss.
[28] Because we can’t remove the hydrology signals, we
instead simulated an ice load over the region using aviscoelastic earth model with 65 km lithospheric thickness,upper mantle viscosity of 4  1020 Pa s and lower mantle
viscosity of 1022 Pa s [Lambeck et al., 1998] to generate ageoid rate anomaly comparable to that seen in the GRACE
solutions.
[29] The choice of a larger lithospheric thickness would
reduce the high-frequency component of the modeled GIA
signal and may affect the conclusions; however, the value
used is typical of values used in other studies [e.g., Ivins
and James, 2005]. There is not a unique solution to thisforward modelling process: a model with more ice loss
earlier will generate a similar present-day signal to a modelwith less ice loss more recently. Therefore we attempted tocover the range of possibilities by generating simulated
geoid rate patterns for two models, one with ice loss of a
1200 m cylinder occurring linearly between 11 ka and 3.5 ka(the time interval of Antarctic ice melt in the ICE-5G model[Peltier, 2004]) (Figure 7b) and the other with 600 m iceloss occurring between 3.5 ka and 1 ka (Figure 7c). We
stress that there are not currently any geophysical observations to support the presence of such ice quantities in thepast; this is simply an exercise to develop a GIA model thatcan reproduce the positive geoid rate anomaly seen by
GRACE in Enderby Land.
[30] Figures 7d and 7e show the present-day verticaldeformation associated with each of the models, withpredicted uplift rates of 10 mm/yr at RICH and MAW1.
The uplift rates at these sites estimated from GPS observations
are near zero and statistically significantly different from thepredicted rates at the 99% confidence level (Figure 8). While,
admittedly, the time series at RICH spans only 1.2 years, theheight estimates clearly do not align with the predicted rate
(red line in Figure 8c). In other words, it is not likely thata GIA model can generate a positive gravity anomaly thatmatches the Enderby Land feature and, simultaneously,

Figure 7. (a) GRACE geoid rate over Enderby Land, East
Antarctica. The cross labeled D (S68°, E54°) indicates theapproximate location used by Chen et al. [2008] to computethe linear trend in Figure 6. Geoid rates associated with(b) 1200-m ice loss from 11 ka to 3.5 ka and (c) 600-m iceloss from 3.5 ka to 1 ka. (d) Present-day vertical velocities
associated with Figure 7b and (e) associated with Figure 7c.
The locations of the nearest GPS sites are shown, whereas
the white outlines show the spatial extents of the regions
where extra ice has been melted.

8 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

B06406logistics personnel who contributed to the fieldwork. We thank K. Lambeck
and T. Purcell for the use of their GIA modeling software. The GPS data
were computed on the Terrawulf II computational facility at the Research
School of Earth Sciences, a facility supported through the AuScope
initiative. AuScope Ltd. is funded under the National Collaborative
Research Infrastructure Strategy (NCRIS), an Australian Commonwealth
Government Programme. We thank the associate editor and three anonymous reviewers whose helpful review comments improved this manuscript.

References

Figure 8. GPS height time series at (a) Mawson (MAW1),(b) Davis (DAV1), and (c) Richardson Lake (RICH) in
Antarctica. The yellow lines are the modeled values of a rate
and annual term, the black lines show the rate only. The redline in Figure 8c represents the 10 mm/yr uplift rate
predicted by the GIA models that approximate the observed
geoid signal seen by GRACE.
rates. The GRACE-derived uplift pattern in the Laurentide
region provides spatial information on the GIA process inregions that lie between GPS sites and also where theavailable geomorphological and shoreline observations cannot constrain the ice thicknesses during the last glacial period.
Such GRACE estimates provide important additional constraints in inversions of geological and geophysical data toconstruct a more accurate ice model for North America.
[32] The positive anomaly in Enderby Land identified inother studies [Chen et al., 2006; Ramillien et al., 2006;
Chen et al., 2008] is not a GIA signal that has been
overlooked in recent GIA models for Antarctica (ICE-5G[Peltier, 2004], IJ05 [Ivins and James, 2005]): GPS uplift
rates estimated both in Enderby Land and at Mawson andDavis are incompatible with the predicted uplift rates
associated with ice loss models implied by the geoid rate
signal seen in the region. Thus we attribute the cause of thepositive anomaly in Enderby Land to the accumulation ofsnow over the 2002 – 2005 period. The recent GRACE
solutions indicate a relatively constant mass since 2005.
[33] Acknowledgments. We thank the GRGS GRACE team formaking their solutions freely available and the IGS for the global GPS
data. The GPS fieldwork undertaken in Enderby Land was supported by theAustralian Antarctic Division and we thank the many field, station, andAltamimi, Z., X. Collilieux, J. Legrand, B. Garayt, and C. Boucher (2007),ITRF2005: A new release of the International Terrestrial Reference
Frame based on time series of station positions and Earth orientation
parameters, J. Geophys. Res., 112, B09401, doi:10.1029/2007JB004949.
Argus, D. (2007), Defining the translational velocity of the reference frameof Earth, Geophys. J. Int., 169, 830 – 838.
Boehm, J., and H. Schuh (2004), Vienna mapping functions in VLBI analyses, Geophys. Res. Lett., 31(1), L01603, doi:10.1029/2003GL018984.
Boehm, J., B. Werl, and H. Schuh (2006), Troposphere mapping functions
for GPS and very long baseline interferometry from European Centre for
Medium-Range Weather Forecasts operational analysis data, J. Geophys.
Res., 111, B02406, doi:10.1029/2005JB003629.
Calais, E., J. Y. Han, C. DeMets, and J. M. Nocquet (2006), Deformation of
the North American plate interior from a decade of continuous GPS measurements, J. Geophys. Res., 111, B06402, doi:10.1029/2005JB004253.
Chen, J. L., C. R. Wilson, D. D. Blankenship, and B. D. Tapley (2006),
Antarctic mass rates from GRACE, Geophys. Res. Lett., 33, L11502,
doi:10.1029/2006GL026369.
Chen, J. L., C. R. Wilson, B. D. Tapley, D. Blankenship, and D. Young(2008), Antarctic regional ice loss rates from GRACE, Earth Planet. Sci.
Lett., 266, 140 – 148.
Davis, J. L., P. Elósegui, J. X. Mitrovica, and M. E. Tamisiea (2004),
Climate-driven deformation of the solid Earth from GRACE and GPS,
Geophys. Res. Lett., 31, L24605, doi:10.1029/2004GL021435.
Davis, J. L., M. E. Tamisiea, P. Elósegui, J. X. Mitrovica, and E. M. Hill
(2008), A statistical filtering approach for Gravity Recovery and Climate
Experiment (GRACE) gravity data, J. Geophys. Res., 113, B04410,
doi:10.1029/2007JB005043.
Dziewonski, A. M., and D. L. Anderson (1981), Preliminary reference
Earth model, Phys. Earth Planet. Inter., 25, 297 – 356.
Farrell, W. E. (1972), Deformation of the Earth by surface loads, Rev.
Geophys., 10, 761 – 797.
Herring, T. A., R. W. King, and S. McClusky (2008), Introduction to
GAMIT/GLOBK, Mass. Inst. of Technol., Cambridge, Mass.
Ivins, E., and T. James (2005), Antarctic glacial isostatic adjustment: A newassessment, Antarct. Sci., 17, 541 – 553, doi:10.1017/S0954102005002968.
Lambeck, K., C. Smither, and P. Johnston (1998), Sea-level change, glacial
rebound and mantle viscosity for northern Europe, Geophys. J. Int., 134,
102 – 144.
Lambeck, K., T. Purcell, and J. Zhao (2008), Sea levels and ice sheets
during the last glacial cycle: New results from glacial rebound modelling,
paper presented at Geological Society William Smith 2008 Meeting,Observations and Causes of Sea-Level Changes on Millennial to DecadalTimescales, London, U. K., 1 – 2 Sept.
Lemoine, J.-M., S. Bruinsma, S. Loyer, R. Biancale, J.-C. Marty,
F. Perosanz, and G. Balmino (2007), Temporal gravity field models inferred from GRACE data, Adv. Space Res., 39, 1620 – 1629.
Lidberg, M., J. M. Johansson, H.-G. Scherneck, and J. L. Davis (2007), An
improved and extended GPS-derived 3D velocity field of the glacialisostatic adjustment (GIA) in Fennoscandia, J. Geod., 81, 213 – 230.
Lombard, A., D. Garcia, G. Ramillien, A. Cazenave, R. Biancale, J.-M.
Lemoine, F. Flechtner, R. Schmidt, and M. Iishi (2007), Estimation ofsteric sea level variations from combined GRACE and Jason-1 data,
Earth Planet. Sci. Lett., 254, 194 – 202.
Luthcke, S. B., H. J. Zwally, W. Abdalati, D. D. Rowlands, R. D. Ray,R. S. Nerem, F. G. Lemoine, J. J. McCarthy, and D. S. Chinn (2006),
Recent Greenland ice mass loss by drainage system from satellite gravity
observations, Science, 314, 1286 – 1289, doi:10.1126/science.1130776.
Milne, G. A., J. L. Davis, J. X. Mitrovica, H.-G. Scherneck, J. M. Johansson,
M. Vermeer, and H. Koivula (2001), Space-geodetic constraints on glacialisostatic adjustment in Fennoscandia, Science, 291, 2381 – 2385.
Pagiatakis, S. D. (1990), The response of a realistic earth to ocean tideloading, Geophys. J. Int., 103, 541 – 560.
Panet, I., V. Mikhailov, M. Diament, F. Pollitz, G. King, O. de Viron,
M. Holschneider, R. Biancale, and J.-M. Lemoine (2007), Coseismic
and post-seismic signatures of the Sumatra 2004 December and 2005
March earthquakes in GRACE satellite gravity, Geophys. J. Int., 171,
177 – 190.

9 of 10

B06406

TREGONING ET AL.: NONSTATIONARY GRACE SIGNALS

Peltier, W. R. (2004), Global glacial isostasy and the surface of the ice-age
Earth: The ICE-5G (VM2) model and GRACE, Annu. Rev. Earth Planet.
Sci., 32, 111 – 149.
Ramillien, G., A. Cazenave, and O. Brunau (2004), Global time variationsof hydrological signals from GRACE satellite gravimetry, Geophys. J. Int.,
158(3), 813 – 826.
Ramillien, G., A. Lombard, A. Cazenave, E. R. Ivins, M. Llubes,
F. Remya, and R. Biancale (2006), Interannual variations of the massbalance of the Antarctica and Greenland ice sheets from GRACE, Global
Planet. Change, 53, 198 – 208.
Rignot, E., J. L. Bamber, M. R. van den Broeke, C. Davis, Y. Li, W. Jan van
de Berg, and E. van Meijgaard (2008), Recent Antarctic ice mass loss
from radar interferometry and regional climate modelling, Nat. Geosci.,
1, 106 – 110, doi:10.1038/ngeo102.
Rodell, M., et al. (2004), The global land data assimilation system, Bull.
Am. Meteorol. Soc., 85, 381 – 394.
Schmidt, R., S. Petrovic, A. Güntner, F. Barthelmes, J. Wünsch, andJ. Kusche (2008), Periodic components of water storage changes fromGRACE and global hydrology models, J. Geophys. Res., 113, B08419,
doi:10.1029/2007JB005363.
Sella, G. F., S. Stein, T. H. Dixon, M. Craymer, T. S. James,
S. Mazzotti, and R. K. Dokka (2007), Observation of glacial isostatic
adjustment in ‘‘stable’’ North America with GPS, Geophys. Res. Lett., 34,
L02306, doi:10.1029/2006GL027081.
Snay, R., M. Cline, W. Dillinger, R. Foote, S. Hilla, W. Kass, J. Ray,
J. Rohde, G. Sella, and T. Soler (2007), Using global positioning systemderived crustal velocities to estimate rates of absolute sea level change
from North American tide gauge records, J. Geophys. Res., 112, B04409,
doi:10.1029/2006JB004606.
Syed, T. H., J. S. Famiglietti, M. Rodell, J. Chen, and C. R. Wilson (2008),
Analysis of terrestrial water storage changes from GRACE and GLDAS,
Water Resour. Res., 44, W02433, doi:10.1029/2006WR005779.
Tamisiea, M. E., J. X. Mitrovica, and J. L. Davis (2007), GRACE gravity
data constrain ancient ice geometries and continental dynamics overLaurentia, Science, 316, 881, doi:10.1126/science.1137157.

B06406

Tapley, B. D., S. Bettadpur, M. Watkins, and C. Reigber (2004), The
gravity recovery and climate experiment: Mission overview and early
results, Geophys. Res. Lett., 31, L09607, doi:10.1029/2004GL019920.
Tregoning, P., and C. Watson (2009), Atmospheric effects and spurious
signals in GPS analyses, J. Geophys. Res., doi:10.1029/2009JB006344,
in press.
Tregoning, P., G. Ramillien, and K. Lambeck (2008), GRACE estimates ofsea surface height anomalies in the Gulf of Carpentaria, Australia, Earth
Planet. Sci. Lett., 271, 241 – 244, doi:10.1016/j.epsl.2008.04.018.
van Dam, T., J. Wahr, and D. Lavallée (2007), A comparison of annual
vertical crustal displacements from GPS and Gravity Recovery and Climate Experiment (GRACE) over Europe, J. Geophys. Res., 112, B03404,
doi:10.1029/2006JB004335.
van der Wal, W., P. Wu, M. G. Sideris, and C. K. Shum (2008), Use of
GRACE determined secular gravity rates for glacial isostatic adjustment
studies in North America, J. Geodyn., 46, 144 – 154.
Wahr, J., M. Molenaar, and F. Bryan (1998), Time-variability of the Earth’s
gravity field: Hydrological and oceanic effects and their possible detection using GRACE, J. Geophys. Res., 103, 30,205 – 30,230.
Wahr, J., D. Wingham, and C. Bentley (2000), A method of combining
ICESat and GRACE satellite data to constrain Antarctic mass balance,
J. Geophys. Res., 105, 16,279 – 16,294.
H. McQueen, P. Tregoning, and D. Zwartz, Research School of Earth
Sciences, Australian National University, OHB-A, Building 61A, Mills
Road, Canberra, ACT 0200, Australia. (paul.tregoning@anu.edu.au)
G. Ramillien, Dynamique Terrestre et Planétaire, UMR 5562, Observatoire, Midi-Pyrénées, Centre Nationale de la Recherche Scientifique,
14 avenue Edouard Berlin, F-31400 Toulouse, France. (ramillie@legos.obsmip.fr)

10 of 10