A method for separating Antarctic postglacial rebound and ice mass balance using future ICESat Geoscience Laser Altimeter System, Gravity Recovery and Climate Experiment, and GPS satellite data

[ 1 ] Measurements of ice elevation from the Geoscience Laser Altimeter System (GLAS) aboard the Ice, Cloud, and Land Elevation Satellite can be combined with time-variable geoid measurements from the Gravity Recovery and Climate Experiment (GRACE) satellite mission to learn about ongoing changes in polar ice mass and viscoelastic rebound of the lithosphere under the ice sheet. We estimate the accuracy in recovering the spatially varying ice mass trend and postglacial rebound signals for Antarctica, from combining 5 years of simulated GRACE and GLAS data. We obtain root-mean square accuracies of 5.3 and 19.9 mm yr (cid:1) 1 for postglacial rebound and ice mass trend, respectively, when smoothed over 250 km scales. The largest source of error in the combined signals is the effect of the unknown time-variable accumulation on the density of the ice column. To estimate this contribution and so obtain better estimates of ice mass trend and postglacial rebound, we add Global Positioning System (GPS) measurements of vertical velocities as additional constraints. Using an empirical relation between the errors in postglacial rebound and ice mass trend that result from the unknown density variation within the ice column, we are able to solve for all three unknowns in the problem: ice mass trend, postglacial rebound, and the snow compaction trend. The addition of a plausible distribution of GPS measurements reduces the errors in estimates of postglacial rebound and ice mass trend to 3.4 and 15.9 mm yr (cid:1) 1 , respectively. GRACE, Antarctica, postglacial rebound, ice mass imbalance, altimeter, GPS Citation:


Introduction
[2] Current best estimates from tide gauge data suggest that sea level rose about 1.5 ± 0.5 mm yr À1 during the last century [Church et al., 2001], but the causes of global sea level rise (i.e., steric expansion versus mass influx, and the possible sources of additional water) are not well determined. The Antarctic ice sheet is the largest reservoir of fresh water on Earth, and estimates of the Antarctic contribution to sea level rise during this century are in the range 1.04 ± 1.06 mm yr À1 [Church et al., 2001]. However, that estimate assumes the ice shelves are in balance. If that assumption is incorrect then the true uncertainty is underestimated. Improved knowledge of the Antarctic contribution would help to better understand the nature of sea level rise and would improve the reliability of climate models used to study the effects of atmospheric greenhouse gases. This in turn would lead to better estimates of possible changes in global temperature and sea level in the coming century.
[3] An important step toward estimating present-day Antarctic and Greenland ice mass balance will be the 2002 launch of NASA's Ice, Cloud, and Land Elevation Satellite (ICESat), which will carry the Geoscience Laser Altimeter System (GLAS) and will have a mission lifetime of 3 -5 years. To study the polar ice sheets, a laser pulse generated by the altimeter will reflect off the snow/ice surface and return to the satellite. Measurements of the round-trip distance, combined with Global Positioning System (GPS) measurements of the geocentric position of the spacecraft, will map changes in the surface elevation of the polar ice sheets at regular time intervals. The exact repeat period of the ground tracks will be 183 days, though there will be a 25day near-repeat subcycle in which measurement locations shift by 15 km at the equator (and much less over the poles) from those made 25 days previously. Changes in the ice sheet elevation will be determined from crossover differences. Using GLAS data, it will be possible to estimate the rate of change in Antarctic ice mass over the mission lifetime between the ice sheet margins and 88°S latitude.
[4] GLAS measurements of ice sheet elevation will reflect more than just the ice mass change. Uplift and subsidence due to postglacial rebound (PGR) and variable compaction of snow will also influence the height change measured by GLAS. Wahr et al. [2000] showed that it is possible to reduce the PGR contribution to the GLAS estimate of Antarctic ice mass balance by combining GLAS measurements with time-variable gravity from the Gravity Recovery and Climate Experiment (GRACE) dedicatedgravity satellite mission. The GRACE mission is under the joint sponsorship of NASA and the Deutsches Zentrum für Luft-und Raumfahrt (DLR). Scheduled for launch in March 2002 with a nominal 5-year lifetime, GRACE consists of two satellites in low-Earth orbit (an initial altitude in the range of 450 -500 km) and a few hundred kilometers apart that range to each other using microwave phase measurements. Onboard GPS receivers will allow positioning of each spacecraft in a geocentric reference frame, and onboard accelerometers will detect nongravitational accelerations so their effects can be removed from the satellite-to-satellite distance measurements. The residuals will be used to map the gravity field. The geoid will be determined orders of magnitude more accurately, and to considerably higher resolution, by GRACE than by any existing satellite.
[5] In this paper we use a technique analogous to that of Wahr et al. [2000] but with several differences. First, we attempt to recover the spatial variability of both the ice mass trend (IMT) and PGR signals, whereas Wahr et al. [2000] estimated only the average of IMT over the entire ice sheet and made no attempt to assess the accuracy of the recovered PGR signal. Second, we simulate a spatially varying IMT from Bentley and Giovinetto [1991] instead of assuming an average trend applied uniformly throughout Antarctica. Third, although we initially solve for PGR and IMT by combining only GLAS and GRACE data, we later add an additional constraint to the simulation using GPS measurements of vertical velocity to solve for PGR, IMT, and the temporal variation of density within the ice column.

Synthetic Data
[6] To simulate the recovery of Antarctic postglacial rebound and ice mass trend from GRACE, GLAS, and GPS, we construct 5 years of synthetic satellite data as the sum of geophysical signals and measurement errors. Almost all of the contributions to these signals are identical to those described by Wahr et al. [2000], and a more detailed discussion of the simulated data can be found therein. The new signal contributions are those found in the GPS data and the spatially varying ice mass trend estimated by Bentley and Giovinetto [1991]. We briefly summarize the contributions to simulated signals here.
[7] Monthly ice sheet elevations in the simulated GLAS measurements consist of contributions from snow accumulation and horizontal ice flow on Antarctica, from the Earth's elastic response to loading by these mass fields, from Earth's viscoelastic response to past loading (i.e., PGR), and from GLAS measurement errors.
[8] For GRACE, we simulate 5 years of monthly measurements of the geoid: the equipotential surface corresponding to mean sea level over the oceans. The geoid can be expanded in a spherical harmonic representation as [Kaula, 1966] where a is the Earth's mean radius, q and f are colatitude and east longitude, C lm and S lm are dimensionless coefficients, and theP lm are normalized associated Legendre functions [e.g., Chao and Gross, 1987]. GRACE will deliver values of C lm and S lm , up to a maximum degree and order (i.e., l and m) of 100, every month. The simulated data include GRACE measurement errors, as well as the gravitational effects of snow accumulation and ice flow on Antarctica, of the elastic response to loading, of PGR, of redistribution of water mass in the ocean and on landmasses other than Antarctica, and of errors in correction for atmospheric pressure.
[9] The simulated GPS data consist of 5 years of daily height coordinates. The coordinates include a constant vertical velocity contributed by PGR, the Earth's elastic response to accumulation and ice flow on Antarctica, and GPS measurement errors estimated to be 1.5 cm root-meansquare (RMS) for daily coordinates.

Snow and Ice Simulation
[10] Contributions of Antarctic snow accumulation to our simulated GLAS, GRACE, and GPS data are derived from monthly precipitation in the National Center for Atmospheric Research (NCAR) Climate System Model (CSM-1) [see, e.g., Boville and Gent, 1998;also G. Bonan, personal communication, 1997]. This model couples component models of atmospheric, oceanic, sea ice, and land surface processes into a 300-year integration which is not meant to represent any particular 300-year time span. Our assumption is that the amplitudes and the spatial and temporal characteristics of Antarctic precipitation are reasonably well represented by the model output (see Wahr et al. [2000] for further details). To isolate the time-variable accumulation for inclusion into our synthetic satellite data, we calculate at every grid point a detrended time series of the accumulated mass. When we refer to the accumulation fields in the remainder of this paper, we mean these residual, detrended fields.
[11] We assume that horizontal ice flow has a constant rate [see, e.g., Oerlemans, 1981], and we use an estimate of the present-day mass balance of different Antarctic drainage basins that is based on accumulation data and drainage basin analyses [Bentley and Giovinetto, 1991]. The IMT estimate includes both the century scale accumulation and the mass lost to horizontal outflow.
[12] To incorporate the snow and ice effects into the GRACE simulated geoid, we sum the CSM-1 estimates of time-variable accumulation and the Bentley and Giovinetto [1991] model of the long-term mean, and integrate over latitude and longitude to generate a time series of the geoid coefficients C lm and S lm (see equations (11) and (12) of Wahr et al. [1998]). To simulate changes in GLAS ice heights corresponding to the change in mass fields, we use the density properties of snow and ice modified by the timeand accumulation-dependent compaction model of Wingham [2000].
[13] We estimate the elastic vertical displacement U of the solid Earth caused by the mass load, using standard methods (see, e.g., equations (40) and (42) of Mitrovica et al. [1994]). We use elastic load Love numbers, h l , computed for the preliminary reference Earth model PREM [Dziewonski and Anderson, 1981] and provided by D. Han (personal communication, 1995). The elastic response in our synthetic GRACE geoid coefficients is represented by the load Love number k l in equation (12) of Wahr et al. [1998].

Simulation of Postglacial Rebound
[14] The synthetic GLAS, GRACE, and GPS data used in this paper each include PGR as a signal component. The Earth model used to estimate PGR assumes a 120 km thick elastic lithosphere, an upper mantle viscosity (between the base of the lithosphere and 670-km depth) of 1.0 Â 10 21 Pa s, and a lower mantle viscosity of 1.0 Â 10 22 Pa s. The Earth's core is assumed to be inviscid. We estimate the PGR vertical velocity by convolving viscoelastic load Love numbers, computed as described by Han and Wahr [1995], with estimates of the Antarctic deglaciation history. We use the ICE-3G Pleistocene ice model of Tushingham and Peltier [1991], with the addition of a preceding 90-kyr glaciation phase to build the ice sheet. We include not only the Antarctic ice sheet but all the other ICE-3G ice sheets as well (e.g., Laurentia, Fennoscandia, Greenland).

Satellite Measurement Errors
[15] We include estimates of satellite measurement errors in the synthetic data for both GLAS and GRACE. For GRACE we use preliminary accuracy estimates for the geoid coefficients provided by B. Thomas and M. Watkins at the Jet Propulsion Laboratory. These errors are described in slightly more detail by Wahr et al. [1998]. We assume the errors are uncorrelated from one month to the next.
[16] We estimate GLAS measurement errors using an indirect approach described by Wahr et al. [2000]. We construct monthly, gridded elevation fields over regions of varying area with errors designed to be consistent with estimates of 6-monthly crossover error by Choe [1997]. We assume that there is no correlation between errors from one month to the next, or from one grid element to another. We also include an additional constant-rate GLAS error of RMS 0.1 mm yr À1 (C. K. Shum, personal communication, 1998), which represents the possible effects of orbital drift errors during the 5-year mission.

Other Contributions to GRACE
[17] The GRACE gravity measurements will also include contributions from variable mass in the ocean, the atmosphere, and the storage of water and snow over continental regions outside of Antarctica. We simulate the effects of the ocean by including monthly geoid coefficients derived from a variant of the Parallel Ocean Program general ocean circulation model developed at Los Alamos National Laboratory [Dukowicz and Smith, 1994]. We include the effects of continental water storage by using output from a land-surface water and energy balance model coupled to a high-resolution climate model at the Geophysical Fluid Dynamics Laboratory in Princeton (C. Milly and K. Dunne, personal communication, 1999). The model generates daily, gridded estimates of soil moisture, snowpack, surface water, and groundwater. The model's predictions of water (i.e., snow) storage on Antarctica are not included, since we are using output from the NCAR climate model for that.
[18] Variations in the distribution of atmospheric mass contribute to the time-variable geoid measured by GRACE. Over land, the contribution of atmospheric mass to the geoid will be removed using global, gridded pressure fields routinely produced by the European Center for Medium Range Weather Forecasts (ECMWF). Errors in the ECMWF pressure fields could degrade the Antarctic results [Velicogna et al., 2001]. To model this degradation, we estimate the errors in monthly averaged pressure fields over land as where P ECMWF and P NCEP are the monthly averaged surface pressure fields predicted by the ECMWF and the National Center for Environmental Prediction (NCEP), respectively. We use dP to estimate monthly geoid coefficients and include those coefficients in our synthetic GRACE data. Pressure corrections to GRACE data will not be made over the ocean. The contributions from the atmosphere over the oceans are included in the simulation via removal of the inverse barometer component of the ocean simulation.
[19] Wahr et al. [2000] show that none of these additional contributions to GRACE (i.e., from the ocean, from continental water storage outside of Antarctica, or from atmospheric pressure errors) has significant impact on the uncertainty of the rate of Antarctic mass change and the PGR signal. This is because these contributions either originate from regions far removed from Antarctica, or have only a small secular component, or both. The effects of the GLAS and GRACE measurement errors are more important but are still considerably smaller than the effects of time-varying density in the ice column resulting from a variable accumulation rate.

Estimating PGR and IMT From GLAS and GRACE Data
[20] To estimate IMT and PGR from the simulated data, we use an iterative method in which we try to avoid dependence on any a priori assumptions about Earth viscosity and ice loading history. The iterative approach is analogous to that described by Wahr et al. [2000]. Initially, we assume that GLAS is sensitive only to the ice thickness change and not to the PGR uplift, and we determine the ice thickness change using only the GLAS elevation data. We refer to this initial estimate as the zeroth iteration.
[21] We compute the secular rate of change in the geoid caused by the zeroth order GLAS mass balance estimate and remove that geoid signal from the GRACE data. The residual secular geoid change is then interpreted as being due entirely to PGR. We use the geoid residual to estimate the corresponding crustal uplift due to PGR, using the method of Wahr et al. [2000], in which we let The explanation for equations (3) and (4) is that the change in the geoid caused by PGR is mostly due to mass anomalies associated with vertical motion of the surface [see . Wahr et al. [2000] showed that the crustal uplift rate _ U PGR calculated from a PGR model differs negligibly from the approximate uplift computed using equation (5) for a wide variety of viscosity profiles and ice models.
[22] We estimate PGR uplift from the residual (GRACE minus the IMT estimate of geoid) using equation (5). Subtracting this PGR estimate from GLAS then results in a better estimate of ice balance. We refer to this estimate as the first iteration. We then repeat the process: compute the geoid contributions from this new GLAS ice mass estimate and remove them from GRACE, interpret the secular component of the new GRACE residuals as the effect of PGR, calculate the PGR uplift and remove it from GLAS, and use the GLAS residuals to construct a further improved ice balance estimate. We iterate this procedure eight times, after which the improvement is negligible.
[23] The resulting PGR and IMT estimates have large errors at short wavelengths because the _ C lm PGR and _ S lm PGR values inferred from the GRACE measurements become increasingly inaccurate as l gets large. To reduce those short-wavelength errors, we smooth the PGR and IMT results by constructing Gaussian averages of those fields. The Gaussian average of the recovered PGR field is [Wahr et al., 1998] where _ U PGR GLAS=GRACE is the unsmoothed field, and the averaging function is where a is the angular distance between (q, f) and (q 0 , f 0 ) and with R the distance along the Earth's surface at which W has dropped to 1/2 its value at a = 0 [see Jekeli, 1981, equation (59)]. We apply this smoothing process to the recovered PGR and IMT fields, using R = 250 km. Thus the end products of our estimation algorithm can best be described as the PGR and IMT fields smoothed over 250 km scales.
To determine the accuracy of these smoothed, recovered fields, we also apply this smoothing process to the input PGR and IMT fields used to construct our simulated data, and we assess the difference between those smoothed input fields and the smoothed recovered fields.

Error Caused by Variations in Accumulation Rate
[24] ICESat and GRACE will have lifetimes of about 5 years. Interannual and interdecadal variations in accumulation rate cause the mass trend on 5-year timescales to differ from the century-scale trend. Part of this difference occurs because the normal variability of climate includes nonsecular components with periods greater than 5 years. For GRACE, which is directly sensitive to mass change, the difference between the 5-year and century-scale trends is the only contribution of the time-variable accumulation to the error in estimating the century-scale trend. For GLAS, which is sensitive to ice sheet thickness rather than to mass, there is the additional problem that the relation between thickness and mass is not simple multiplication by r i but is complicated by variable compaction of the snow-ice column.
[25] Because the density profile in the upper layers of the snow-ice column depends on prior accumulation rate, the assumption of a constant ice density r i introduces an error in the estimate of IMT from GLAS data. It is unlikely that the variable accumulation at interannual and interdecadal periods can be reconstructed accurately enough to directly calculate variable compaction, although GLAS measurements will help to constrain it for the short term, and microwave remote sensing of scattering from ice lenses and pipes [e.g., Cogley et al., 2001] will contribute additional information. Nevertheless we expect that, as a first approximation, researchers will convert GLAS estimates of ice thickness rate of change to mass rates by simply multiplying by the maximum density of compact ice, r i = 917 kg/m 3 , since this is the density appropriate for the century-scale mass variability. This will introduce an additional error into the GLAS mass balance estimates that does not affect GRACE [Wahr et al., 2000].
[26] Wahr et al. [2000] consider two separate effects of time-varying accumulation as two different error types. They refer to the difference between a 5-year mass trend and the century-scale mass trend as the ''undersampling error,'' and the error in the GLAS mass balance estimate caused by approximating the time and accumulation dependence of snow compaction with _ h = _ m/r i as the ''compaction error.'' Note that the compaction error causes an error in the GLAS estimate of the 5-year trend, and not just the century-scale trend. We only concern ourselves with how accurately we will be able to recover the 5-year trend in this paper. Wahr et al. [2000] conclude that the compaction error in estimates of the 5-year ice mass trend using GLAS alone is likely to be ±4.5 mm yr À1 (water thickness equivalent) when averaged over the entire Antarctic ice sheet, which is equivalent to an error of about ±0.15 mm yr À1 in global sea level rise. However, without combining with GRACE data, the GLAS IMT estimates would also include errors due to unmodeled PGR crustal uplift, which maps directly into a change in ice sheet elevation.
[27] The addition of the GRACE data permit the separation of the PGR and IMT signals, and so reduce the contribution of the unmodeled PGR signal to the IMT error. Wahr et al. [2000] conclude that the PGR contribution to the IMT error averaged over the ice sheet is on the order of 5 mm yr À1 RMS. Wahr et al. [2000] also find that the GRACE/GLAS iteration described here would remove all of the PGR contribution to the IMT error if there was no time-variable accumulation. However, time-variable compaction violates the assumption of r = r i used to estimate mass from GLAS heights, and the mass inconsistency leads to an error in the final estimate of PGR contribution to the surface heights. Thus the GRACE residuals used to infer the PGR signal suffer indirect effects of the compaction error. Wahr et al. [2000] find that the final PGR contribution to the IMT error, averaged over the entire ice sheet, is approximately 0.31 times the IMT compaction error. The value of the proportionality factor is an artifact of the method used to combine the GLAS and GRACE data. Specifically, it arises from combining equation (15) of Wahr et al. [1998] to estimate the geoid from the GLAS surface mass (i.e., IMT) estimate, with equation (5) above, relating the GRACE PGR geoid to an inferred uplift rate. The proportionality factor between PGR and IMT after any single iteration is about equal to d = 0.24. After N iterations, the proportionality factor = d (1 À d N )/(1 À d) % 0.31 when N is large.

PGR and IMT From GLAS, GRACE, and GPS Data
[28] In the simulations performed by Wahr et al. [2000], the limiting error in the IMT recovery arises principally because only two observables (GLAS ice sheet elevations plus the GRACE geoid) are used to determine three unknowns (IMT, PGR, and the time-variable density of the ice column). To solve for all three unknowns, we must add an additional constraint. We consider here the assimilation of continuous GPS point measurements of vertical velocity, simulated as described in section 2. Although the simulated data include the Earth's elastic response to loading, the velocity is dominated by PGR. We estimate PGR and IMT using the iterative algorithm described in section 3. Then, following the last iteration, the PGR error is estimated from the difference between GPS vertical velocities (slope of the daily height coordinates) and the PGR estimate _ U PGR derived by applying equation (5) to the secular GRACE geoid minus the geoid from the estimated IMT. That PGR error is used to estimate the ice compaction trend (see below), and a corresponding correction to the estimate of IMT.
[29] We expect that the GPS measurements available to estimate the PGR correction will be sparse and irregularly distributed, so that assimilating them into the data constraints will require some care. We first interpolate the gridded _ U PGR GLAS=GRACE to each GPS location using two-dimensional optimal interpolation. The PGR correction at a GPS site is simply Á _ U PGR compaction = _ U PGR GPS À _ U PGR GLAS=GRACE . Then, we interpolate the PGR correction back to the grid points using an optimal interpolation. However, we expect the GPS sites to be sparse and widely distributed, so we multiply the interpolated correction at each grid point by a Gaussian function, W(a), given by equation (7), where R in equation (8) is 500 km and a is the angular distance to the nearest GPS site (see section 7 for more details about the choice of R). By thus downweighting the correction, we effectively localize the estimate of PGR error near the GPS sites where it is measured. As we will address further in the discussion section, the accuracy of the estimate of PGR error calculated using this method depends fundamentally on the correlation length scale of the time-varying compaction effects that cause the error.
[30] Local correction of the GLAS/GRACE PGR estimate with GPS vertical velocities improves the PGR estimate, but subtracting the corrected PGR estimate from GLAS heights does not dramatically improve the estimate of IMT. However, the empirical relation between the IMT compaction error and the compaction error in PGR Á _ U PGR compaction = À0.31Á _ h PGR compaction (see the discussion at the end of section 4), allows us to use the GPS-derived PGR errors to also estimate the compaction error. In practice, we use the iterative routine that combines GRACE and GLAS (see description in section 3) and after the last iteration we use the GPS to correct the GLAS/GRACE estimate of PGR by removing the difference between GPS and GLAS/ GRACE PGR vertical velocities (Á _ U PGR compaction ) estimated as described above. Then we correct the IMT by removing the PGR error and the compaction error estimated as the PGR error divided by 0.31. In this manner we effectively solve for all three unknowns in the Antarctic mass balance (i.e., IMT, PGR, and the time variation of density within ice columns). Once we have obtained the PGR and IMT estimates, we smooth those fields using the 250-km Gaussian averaging function defined in equation (7).

Results
[31] In this section we present results of the simulations described in sections 3 and 5, using 5 years of simulated GLAS, GRACE, and GPS data constructed as described in section 2. In the figures, we will show results for a single 5year time span from among the 219 years of CSM-1 accumulation fields that we used. In Table 1, we will compare results from that 5-year time span as well as from multiple independent 5-year periods.

Ice Mass Trend and Postglacial Rebound Estimation From GRACE and GLAS Data
[32] Figure 1a shows one example of an input 5-year, smoothed IMT field, consisting of the century-scale IMT from Bentley and Giovinetto [1991] plus the 5-year trend from one period of the Antarctic accumulation fields. Figure  1b shows the smoothed IMT field recovered using only GLAS and GRACE observations, and Figure 1c depicts the difference between these two smoothed fields. The error in recovering IMT is dominated by the compaction error. If we run the simulation without including time-variable accumulation in the simulated data, the difference between the VELICOGNA AND WAHR: GLAS, GRACE, AND GPS IN ANTARCTICA recovered and input smoothed IMT fields is small, with an RMS over the entire ice sheet of 2.1 mm yr À1 (Table 1). The residual has RMS 18.4 mm yr À1 when variable accumulation is included for this particular example 5-year period. That corresponds to $68% of the RMS of the smoothed IMT signal (the IMT RMS = 27 mm yr À1 ), and so we are able to recover the variance of the input signal with an error of 18.4 2 /27 2 = 46%.
[33] Figure 2 depicts the recovery of the smoothed PGR signal from GLAS and GRACE data. The difference between the smoothed input ( Figure 2a) and smoothed recovered (Figure 2b) signals is given in Figure 2c. The largest errors in the PGR recovery, as for the IMT recovery, result from the compaction error. When the time-variable accumulation is not included in the simulated data, the residual between the input and recovered PGR signals is relatively small (1.3 mm yr À1 ; see Table 1). With timevariable accumulation included, the RMS error in the recovered PGR during this particular example 5-year period is 5.5 mm yr À1 . For comparison, the smoothed PGR signal depicted in Figure 2a has an RMS of 8.3 mm yr À1 , so that the variance in that signal has been recovered with an error of 5.5 2 /8.3 2 = 43%.
[34] The example calculations shown in Figures 1 and 2 are representative of errors that might be expected from this approach. We ran this simulation for all of the independent 5-year periods. The RMS error of the spatially varying signals for individual 5-year periods ranges from 8.4 to 31.5 mm yr À1 for IMT and from 2.9 to 8.5 mm yr À1 for PGR ( Table 1). The RMS average of the RMS errors is 19.9 mm yr À1 for IMT and 5.3 mm yr À1 for PGR. Not surprisingly, 5year periods with larger errors correspond to periods with unusually large accumulation rates.

Ice Mass Trend and Postglacial Rebound Estimation From GRACE, GLAS, and GPS Data
[35] To better estimate PGR and IMT, we include GPS measurements of vertical velocity in the simulation. Initially, we assume continuous GPS measurements are made at 45 existing and proposed sites, including previous campaign sites along the Antarctic coast [Dietrich, 2001] and existing and proposed continuous sites along the Transantarctic Mountains and Marie Byrd Land in West Antarctica [Raymond et al., 1998]. Thirty-six of the 45 sites used are on the Antarctic mainland, while the remainder are on islands at high southern latitudes.
[36] Figure 3 shows the recovery of the smoothed IMT after assimilation of the GPS vertical velocities, for the same 5-year period examined using only GLAS and GRACE in Figure 1. Figure 4 depicts the smoothed PGR signal recovery given this GPS site distribution for that 5-year period. By incorporating the GPS data, we obtain a smoothed PGR estimate with just 2.8 mm yr À1 of RMS error (Table 1), a 49% improvement relative to that recovered from GRACE and GLAS data alone. The corresponding estimate of the smoothed IMT has 12.0 mm yr À1 of RMS error, a 35% improvement. The similarity of the error reductions for IMT and PGR indicates that the empirical scaling used to relate PGR error to IMT compaction error is justified. Using this distribution of 45 GPS sites, we calculated estimates for all of the independent 5-year periods of the CSM-1 model, and the RMS average of the RMS errors was 3.4 mm yr À1 of PGR error and 15.9 mm yr À1 error in IMT. The PGR error ranged from 1.4 to 5.4 mm yr À1 and the error in the IMT estimate ranged from 5.8 to 26.9 mm yr À1 .
[37] Because the largest error source in the PGR and IMT estimates is time-variable accumulation, we also examine how the estimates of IMT and PGR are affected by using just a few GPS sites located where the compaction error is largest in Figure 1. The compaction error is largest where the accumulation rates are largest. The results after locating the GPS measurements near the largest accumulation centers are shown in Figure 5. With a configuration of only 10 GPS sites, the RMS error in the smoothed PGR estimate for the same example 5-year period shown in Figure 2 is 2.3 mm yr À1 . This is a reduction of 58% from the GLAS/GRACE estimate. The RMS error of the smoothed IMT estimate is 10.5 mm yr À1 , a similar 43% reduction. When calculated for all of the independent 5-year periods, the RMS average of the RMS errors was 2.6 mm yr À1 of PGR error and 13.2 mm yr À1 error in IMT. The PGR error ranged from 1.5 to 3.6 mm yr À1 and the error in the IMT estimate ranged from 6.1 to 19.1 mm yr À1 . However, even if we ignore the logistical difficulties of siting continuous GPS instruments in the Antarctic interior, the locations where the accumulations are largest are among the least likely places to find exposed bedrock. Moreover, there is no reason to expect that the precipitation rate always will be highest at the locations indicated in Figure 1c. The loci of the maximum compaction errors will change with time, and we cannot predict where the compaction error will be largest during the GRACE/GLAS missions.
[38] To obtain an extreme lower bound for the recovery errors, we also examined the accuracy of signals recovered 10 GPS sites 2.3 2.6 1.5 -3.6 10.5 13.2 6.1 -19.1 50-km-spaced GPS 1.5 --7.9 -in the unlikely limit of very dense GPS coverage. By putting GPS sites at a regular 50 km spacing everywhere south of À60°, we retrieved the smoothed PGR with RMS error of just 1.5 mm yr À1 and the smoothed IMT with 7.9 mm yr À1 for the example 5-year period shown in Figures 1 -5 (Table 1). We suspected that a portion of the remaining error was a result of the uncorrected elastic load response contained in the GPS velocities. To test this we ran the simulation again but did not include the elastic load response in the GPS heights. The errors in the signal recovery were reduced to 0.1 mm yr À1 (PGR) and 7.7 mm yr À1 (IMT). We thus infer that the portion of the PGR RMS error caused by uncorrected secular elastic loading in the GPS signal is % ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1:5 2 À 0:1 2 p mm yr À1 = 1.5 mm yr À1 , and the corresponding portion of the IMT error is % ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7:9 2 À 7:7 2 p mm yr À1 = 1.8 mm yr À1 . We expect that the effects of the uncorrected elastic errors would be similar for any distribution of GPS sites but that they could be

Discussion and Conclusions
[39] Our simulations show that by combining GLAS and GRACE data it should be possible to recover the spatially varying Antarctic PGR and IMT signals, smoothed over 250 km scales, to accuracies of about 5 and 20 mm yr À1 , respectively. These errors are due almost entirely to the compaction error associated with variability in the accumu-lation. When GPS measurements are combined with the GLAS and GRACE data, the accuracy of the recovered fields is improved. This improvement depends on the distribution of the GPS receivers. When we assume there are continuous data from 45 existing and proposed GPS sites, the RMS errors in the smoothed PGR and IMT fields drop to about 3 and 16 mm yr À1 , respectively.
[40] When GPS data are added to the simulations, compaction error in the IMT is reduced by exploiting the relationship observed by Wahr et al. [2000], that compac-  tion errors in PGR are 0.31 times the compaction errors in IMT. However, differences between the GPS vertical velocity and the GRACE/GLAS estimate of PGR that arise from sources other than the compaction trend will cause errors in the correction of IMT, and the correction scales these differences as approximately IMT = (1 + 1/.31) PGR = 4.2 PGR . One example of a difference between the GPS vertical velocity and the GRACE/GLAS estimate of PGR that does not result from the compaction trend arises because of the smoothing which must be applied to GRACE geoid estimates. The PGR signal used in this paper, when smoothed using a Gaussian averaging function with 250 km radius, differs from the unsmoothed PGR signal by $1.4 mm yr À1 RMS. We are able to attenuate the effect of this difference somewhat by applying Gaussian smoothing to the estimate of the PGR correction, Á _ U PGR compaction , before using it to correct the IMT. Another difference between the GPS and GLAS/GRACE estimate of PGR arises from the approximation used to convert geoid to PGR uplift in equation (5). This approximation produces an RMS error of 0.9 mm yr À1 in PGR for the signal modeled in this paper.
[41] Errors in the estimate of IMT described here are substantially larger than errors in the estimate of total IMT for all Antarctica, of ±4.5 mm yr À1 cited by Wahr et al. [2000], primarily because we have considered the spatially varying IMT in this analysis. Errors scale inversely proportional to the square root of the area considered, and improvement in the RMS estimate of the spatially varying IMT that results from including GPS should translate to improvements in the estimate of the total IMT for Antarctica.
[42] All of the error estimates in this paper depend heavily on the assumed time-variable accumulation field and GPS receiver distribution. It is difficult to assess the accuracy of the CSM-1 accumulation fields [Wahr et al., 2000]. We expect that Antarctic GPS receivers will be sparsely distributed, with the 45 sites representing a bestcase scenario. Moreover we expect these sites will be much more sparsely sampled in time than the 5 years of daily coordinates we simulated. Among existing sites in the Transantarctic Mountains, the coordinate time series from MCM4 (the IGS site at McMurdo station) is quasi-continuous, but remote stations COAT and MTCX have just 1 -5 months of data per year as a result of power and other system failures. Also many of the sites we incorporated are campaign sites, which may have only a few days to weeks of data per year. Vertical velocities estimated from the continuous GPS sites currently have errors of between $0.4 and 5.9 mm yr À1 , as compared to 0.3 mm yr À1 RMS errors in recovery of PGR vertical velocity from the simulated GPS time series. Errors in GPS vertical velocity estimates of more than about 4.4 mm yr À1 would introduce errors into the correction of IMT that exceed the RMS errors in GLAS/GRACE recovery without GPS. However, estimates of vertical velocity at continuous sites will improve with more observation, and spatial averaging will reduce the error further in the more densely instrumented regions (e.g., the Transantarctic Mountains and Marie Byrd Land).
[43] Another caveat in applying this method relates to the true scales of correlation of the signals we are solving for. The GLAS, GRACE, and GPS data sets sense changes at different length scales, ranging from sparsely distributed point velocity measurements (GPS) to dense altimetric heights with 70-m footprints (GLAS) to surface mass density integrated over scales >200 km (GRACE). Consequently the accuracy of the estimates of PGR, IMT, and the correction for the time-varying compaction trend depends on the length scales for which these signals are self-correlated. If all three signals self-correlate at scales consistent with the minimum resolution scale of GRACE and the scale of distribution of the GPS measurements, then we expect the results to be quite good. If on the other hand one of the signals is extremely variable at scales of a few hundred kilometers, we expect a poor result. [44] As evident from Figure 2, the PGR signal used in our simulation does have long characteristic length scales. This partially reflects the fact that the Antarctic component of the ICE-3G deglaciation model used to generate this PGR signal is dominated by these same long scales. However, even if the true deglaciation pattern had significant power at short wavelengths, viscoelastic rebound would be dominated by longer wavelengths because the stress induced by short-wavelength loading tends to be concentrated near the Earth's surface, and so within the elastic lithosphere. Thus it is not subject to viscoelastic relaxation.
[45] To estimate the minimum plausible correlation scale for PGR, we create an ice load with equal power at all wavelengths, distributed over a disc with angular radius 20°( the approximate angular radius of Antarctica). We assume this ice load, shown in Figure 6a, melted at a constant rate between 10,000 and 4000 years ago. The load has zero mean when averaged over all of Antarctica, and the spatial RMS of the ice thickness is arbitrarily chosen to be 1 km. We convolve this deglaciating load with viscoelastic Green's functions, to predict the present-day vertical velocities shown in Figures 6b and 6c. We assume the Earth has the viscoelastic profile described in section 2.2, but with a lithospheric thickness of either 120 km Figure 6b) or 80 km (Figure 6c).
[46] We calculate the RMS of the difference between the vertical velocities shown in Figures 6b and 6c, and a Gaussian average (equation (6)) of those same vertical velocities, using various averaging radii. The ratio of this RMS difference to the total RMS of the vertical velocity is shown in Figure 6d as a function of averaging radius. A small RMS ratio at an averaging radius of R implies there is little power in the velocity field at scales smaller than R. Note from Figure 6d that an 80 km lithosphere (dashed line) leads to predictions with larger RMS ratios, implying more short-wavelength power, than does a 120 km lithosphere (solid line). Both models, though, have RMS ratios <1/e for radii of up to about 250 km: 290 km for the 120 km lithosphere and 230 km for the 80 km lithosphere. This suggests that no matter what scales are present in the deglaciation pattern, the PGR signal is significantly selfcorrelated out to lengths of 250 km. Note that this method underestimates the true PGR correlation scale because, unlike the ice model shown in Figure 6a, the real deglaciating Antarctic ice is unlikely to have as much power at short wavelengths as at long wavelengths.  [47] The self-correlation length scale of the compaction error depends on the corresponding scale of the time-variable accumulation. Wingham et al. [1998] note that 5-year accumulation trends measured by the ERS-1 and ERS-2 satellite radar altimeters are correlated out to a length scale of $400 km. Using the same CSM-1 accumulation fields used to generate our simulated GLAS/GRACE/GPS data, we compare the compaction trends with Gaussian averages of those same trends using a 250 km averaging radius. We find that the Gaussian-averaged compaction trends are $0.5 that of the unsmoothed, suggesting that indeed the compaction predicted from the CSM-1 accumulation is correlated on these length scales.
[48] We also examine the effects of changing the radius of the Gaussian cap (equation (7)) that was used to downweight the estimated compaction trend correction far from the nearest GPS site. We expect that if GPS corrections are interpolated to distances where the compaction trend is decorrelated from that of the GPS sites, the overall error in the estimates of PGR and IMT will actually increase. On the other hand, if the correction is dampened too much at correlated distances, the GPS correction will not be optimal. Hence a minimum in the RMS error of the estimates can be expected for some radius of the Gaussian dampening, which will be related to the correlation length scale of the compaction trend. Using the 45 continuous and campaign sites (e.g., Figure 4), we ran the simulation for 20 different Gaussian radii and all of the different 5-year periods, and found that a minimum in the RMS error of recovery occurs at a radius of about 500 km (Figure 7). This minimum reflects correlation properties of the CSM-1 accumulation model as opposed to real data, and the spatial correlation properties of the model may not be representative. However, the result is consistent with what we would expect to see given the $400 km length scale of decorrelation for accumulation observed by Wingham et al. [1998], and the >700 km length scale of decorrelation of the PGR signal ( Figure 6). This was the reason we adopted 500 km for the radius of dampening in the results presented here.