Postglacial rebound and Earth’s viscosity structure from GRACE

[ 1 ] The Gravity Recovery and Climate Experiment (GRACE) satellite mission was launched in March 2002 and has an expected 5-year lifetime. One potential application of GRACE measurements of time-variable gravity will be to isolate the postglacial rebound signal, which can then be used to estimate the Earth’s viscosity structure. In this paper we present a sensitivity analysis of simulated GRACE data, designed to assess the accuracy with which those data can be used to recover a simple model of Earth viscosity. We find that without combining with any other data type, but ignoring complications caused by uncertainties in the global ice loading history, GRACE data alone would allow us to determine the viscosity of a uniform lower mantle layer and an upper mantle/transition zone layer to within ±30–40% and to estimate lithospheric thickness to within ±15–20%. GRACE will have a harder time differentiating between the separate viscosities of the transition zone and upper mantle, but accuracies of within a factor of 2 might still be achievable for those parameters. Errors in the ice loading history could significantly degrade these viscosity estimates, particularly for the transition zone and upper mantle. The accuracy of recovery of the true Earth viscosity will depend in part on how well the model parameterization used for the grid search can represent the true Earth structure. However, combining GRACE data with data from other more traditional measurements of postglacial rebound has the potential of dramatically improving viscosity estimates throughout the Earth, particularly in the lower mantle.


Introduction
[2] Geodetic observations are playing an increasingly important role in constraining postglacial rebound (PGR) of the solid Earth. The geological evidence, especially observations of relative sea level variations during the late Pleistocene and Holocene (i.e., over the last 20,000 years or so), has proven exceptionally useful for learning about upper mantle viscosity and lithospheric thickness. These observations, though, are relatively insensitive to deepmantle viscosity. Certain types of geodetic observations, on the other hand, are primarily sensitive to lower mantle viscosity. These include observations of the large negative free-air gravity and geoid anomalies centered over the original ice sheets, of temporal variations of the Earth's gravity field, of positions of points on the surface, and of changes in the Earth's rotation (for recent reviews, see Peltier [1998] and Wahr and Davis [2002]). Furthermore, geodetic observables and relative sea level variations are sensitive to the history and spatial distribution of the ice sheets in different ways. Combining those two data sets can help resolve the trade-off between viscosity and ice sheet models.
[3] The ongoing PGR process of redistribution of mass deep within the Earth causes temporal variations in the gravity field. Satellites are far better at mapping gravity fluctuations over large (hundreds to thousands of kilometers) spatial scales than are surface gravimeters. Satellite measurements of the time-varying gravity field are relatively insensitive to localized non-PGR gravity effects and to the free-air effect of vertical crustal motion seen in surface gravity measurements. Measurements of long-wavelength, PGR-induced gravity fluctuations are useful because they are sensitive to viscosity deep within the Earth.
[4] One of the most serious obstacles to satellite observations of PGR-induced gravity fluctuations is that these signals are contaminated by other secularly varying mass anomalies. Particularly troublesome for existing measurements have been the possible effects of present-day changes in mass of the Antarctic and Greenland ice sheets. Results from satellite laser ranging (SLR), especially from LAGEOS, have been used to determine secular changes in the longestwavelength zonal (i.e., independent of longitude) gravity coefficients. So far, published SLR results have not provided secular coefficients for any nonzonal terms, which carry information about longitude-dependent mass variability. Although some attempt can be made to use the zonal coefficients to separate the PGR and polar ice contributions, it is not possible to use only the zonal terms to meaningfully separate the Canadian or Scandinavian PGR signals from each other, or from the PGR or present-day signals over Greenland.
[5] Satellite measurements of time-variable gravity will become much less ambiguous using data from the Gravity Recovery and Climate Experiment (GRACE). The GRACE mission is under the joint sponsorship of NASA and the Deutsches Zentrum fü r Luft-und Raumfahrt (DLR). Launched in March 2002 and 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 about 200 -250 km apart that range to each other using microwave phase measurements. Onboard GPS receivers determine the position of each spacecraft in a geocentric reference frame. Onboard accelerometers are used to detect nongravitational accelerations so that their effects can be removed from the satellite-to-satellite distance measurements. The residual range rates are used to map the gravity field, orders of magnitude more accurately, and to considerably higher resolution than by any other satellite. Temporal variations in gravity are determined down to scales of about 200 km every 30 days. These can be used to study problems in a number of disciplines, from monitoring changes in water, snow, and ice on land, to determining changes in seafloor pressure, to studying PGR within the solid Earth. Comprehensive descriptions of these and other applications are given by Dickey et al. [1997] and Wahr et al. [1998].
[6] In this paper we simulate the recovery of the PGR signal using GRACE data alone. We are interested in assessing how accurately a simple approximation of the Earth's viscosity structure can be recovered using GRACE data, particularly given that the GRACE secular coefficient rates will contain many different secular mass flux signals in addition to PGR. Toward this end, we include the effects of GRACE errors and of other possible secular signals, such as those from hydrology, ocean circulation, and present-day ice mass trends, as well as the possible effects of errors in the late Pleistocene and Holocene ice history. We do not address here the inaccuracies that might result in attempting to recover more complex (e.g., three-dimensional) viscosity structures; nor do we consider the more general and favorable situation in which GRACE measurements are combined with other types of PGR observations to improve the viscosity estimates.

Data Simulation
[7] The secular rate of change in the geoid can be expanded in a spherical harmonic representation as [e.g., Kaula, 1966] where a is the Earth's mean radius, q and f are colatitude and east longitude, _ C lm and _ S lm are the secular rates of change of the dimensionless Stokes' coefficients C lm and S lm , andP lm are normalized associated Legendre functions [e.g., Heiskanen and Moritz, 1967;Chao and Gross, 1987]. The spatial scale (half wavelength) of any (l, m) term in this expansion is approximately equal to 20,000/l km.
[8] GRACE will deliver measurements of C lm and S lm , up to a maximum degree and order (i.e., l and m) of 100, every month. Secular terms can be fit to these monthly values to determine _ C lm and _ S lm . It is these secular rates of change that will be used to recover the PGR signal.
[9] To evaluate the accuracy with which GRACE can recover this signal and to explore the implications for estimating the Earth's viscosity profile, we construct 5 years of synthetic, monthly Stokes' coefficients C lm and S lm , which include simulated GRACE measurement errors as well as plausible gravitational effects of PGR, present-day ice mass imbalances in Antarctica and Greenland, redistribution of water mass in the ocean and in its storage on continents other than Antarctica and Greenland, and errors in the corrections applied for changes in atmospheric mass. The secular rate of change for each dimensionless coefficient is estimated by least squares fitting to the synthetic monthly values.

Simulation of Postglacial Rebound
[10] Our synthetic GRACE data include contributions from PGR, which we calculate using a dynamical model. That model requires knowledge of the Earth's viscoelastic structure and of the spatial and temporal distribution of global ice during the late Pleistocene and Holocene.
[11] For our estimate of the deglaciation history we generally use the global ICE-3G Pleistocene ice model of Tushingham and Peltier [1991], with an additional 90 kyr linear glaciation phase added at the beginning. Though in one experiment, we use ICE-4G [Peltier, 1994] instead.
[12] We adopt a viscosity profile consisting of three uniform viscosity layers: the upper mantle (between the base of the lithosphere and 400 km depth), the transition zone (between 400 and 670 km depth), and the lower mantle (between the core-mantle boundary and 670 km depth). The overlying lithosphere and underlying fluid core are assumed to be elastic. We simulate results for many viscosity values and lithospheric thicknesses. For our ''reference'' model, we assume a lithospheric thickness of 100 km, an upper mantle viscosity of 1.0 Â 10 21 Pa s, a transition zone viscosity of 1.0 Â 10 21 Pa s, and a lower mantle viscosity of 1.0 Â 10 22 Pa s.
[13] For the elastic structure, we assume these five spherical layers are incompressible and homogeneous, with densities and shear wave velocities as given in Table 1. Although Table 1 implies a lithospheric thickness of 100 km, we will vary that thickness in the experiments below. By adopting this simplified elastic structure we are able to compute secular Stokes' coefficients for a wide range of viscosity profiles and for a large number of spherical harmonic degrees, with a minimum of effort.
[14] To make sure this simplified elastic model leads to reasonable results, we have compared our predicted values of _ C lm and _ S lm for a few viscosity profiles, with those computed by Han and Wahr [1995] for the same viscosity profiles. Han and Wahr use a compressible Earth with elastic values taken from the Preliminary reference Earth model (PREM) [Dziewonski and Anderson, 1981]. We have found agreement to within 10% or better for angular degrees less than about 25, with better than 5% agreement for degrees less than 20. We will use our model results only to simulate the extraction of a PGR signal from noisy GRACE data. We assume this level of discrepancy from PREM does not notably affect our conclusions.

Satellite Measurement Errors
[15] We include estimates of satellite measurement errors in the synthetic data. We use preliminary accuracy estimates for the Stokes' coefficients provided by B. Thomas and M. Watkins at the Jet Propulsion Laboratory. These errors are described in more detail by Wahr et al. [1998]. We include them in our synthetic monthly secular Stokes' coefficients by generating Gaussian random numbers with variances s 2 consistent with these preliminary error estimates. We assume the errors are uncorrelated from one month to the next, in which case we expect the 1-sigma RMS error on the trend of n monthly measurements to be 12s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 12=½ðn À 1Þnðn þ 1Þ p mm/yr. Here, the factor inside the square root comes from the covariance matrix of the fit; see, e.g., section 15.6 of Press et al. [1992].
[16] The GRACE errors for the secular components can be directly compared with the PGR signal, for a preliminary estimate of the sensitivity of GRACE to PGR. Figure 1a compares the degree amplitudes of the secular geoid variation, defined as a , predicted from our reference PGR model (the pluses), with the corresponding degree amplitudes of the secular GRACE errors (the solid black line). Note that the GRACE errors are smaller than the PGR signal for degrees less than about l = 40. Since there are 2l + 1 values of m for every l, and since the GRACE errors are expected to be roughly independent of m, there are likely to be on the order of 40 2 secular Stokes' coefficients where the PGR signal exceeds the GRACE measurement errors. In fact, Figure 1b shows the number of Stokes' coefficients, for each l, where the expected secular errors in the GRACE measurements are smaller than the PGR contributions predicted using the reference model. The dotted line in this figure shows the total number of Stokes' coefficients (i.e., 2l + 1) for each l. In all, there are between 2100 and 2200 Stoke's coefficients where the predicted PGR signal is larger than the expected GRACE measurement errors, as opposed to maybe five coefficients, all zonal, which have been detected using SLR.

Other Contributions to GRACE
[17] The GRACE gravity measurements will also include contributions from variable mass in the ocean and from the storage of water and snow over continental regions. We simulate the effects of the ocean on the geoid by using output from a variant of the Parallel Ocean Program, a general ocean circulation model developed at Los Alamos National Laboratory [Dukowicz and Smith, 1994] which is forced by both winds and pressure. The model output is used to construct monthly Stokes' coefficients (M. Molenaar and F. Bryan, personal communication, 1998), as described by Wahr et al. [1998], which are then added to the synthetic GRACE data. The inverse barometer response to pressure variability is removed from the oceanic solution prior to constructing the Stokes' coefficients. This removal is equivalent to adding, to the full oceanic solution, the time-variable mass signal from the atmosphere above the oceans.
[18] 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. These are used to construct monthly averaged Stokes' coefficients, which are then included as part of the synthetic GRACE data. The model's predictions of water (i.e., snow) storage on Antarctica and Greenland are not included, since we will model ice mass balance from these regions separately.
[19] Variations in the distribution of atmospheric mass also contribute to the time-variable geoid measured by GRACE. The contributions from the atmosphere over the oceans are included in the simulation through the removal of the inverse barometer component of the ocean simulation. Over land, the atmospheric mass variability can be modeled using global, gridded pressure fields, such as those routinely produced by the European Center for Medium Range Weather Forecasts (ECMWF) and the U.S. National Centers for Environmental Prediction (NCEP). The atmospheric contributions estimated from these pressure fields over land will be routinely removed from the GRACE measurements. Pressure corrections will not be made over the ocean. The pressure fields over land are not perfect, and errors in those fields could degrade the monthly geoid estimates [Velicogna et al., 2001]. To model this degradation, we estimate the errors in the monthly averaged pressure fields over land as dP ¼ ðP ECMWF À P NCEP Þ= ffiffi ffi 2 p , where P ECMWF and P NCEP are the monthly averaged surface pressure fields predicted by the ECMWF and NCEP, respectively. We use dP to estimate monthly Stokes' coefficients and include those coefficients in our synthetic GRACE data. Velicogna et al. [2001] note that the difference between ECMWF and NCEP analyzed pressure fields underestimates the true pressure error because both models assimilate the same observations. However, the error in the secular Stokes' coefficient rates introduced by imperfect correction for atmospheric mass is negligible relative to other sources of error. The upper mantle radius listed here implies a lithospheric thickness of 100 km. Other values are also used (see text). The Earth is assumed to be incompressible.
[20] We also include the effects of present-day changes in polar ice sheet mass. The present-day mass imbalance of the ice sheets is not well known, but recent estimates from mass imbalance studies [Church et al., 2001] suggest that Antarctic ice is melting at an average rate of 31 ± 32 mm/yr (corresponding to a rate of global sea level rise of 1.04 ± 1.06 mm/yr), and Greenland is melting at an average rate of 22 ± 27 mm/yr (corresponding to a sea level rise of 0.12 ± 0.15 mm/yr.
[21] We added these best estimates of the total ice mass balance for Antarctica and Greenland, distributed uniformly over each land mass, to estimate their contributions to the secular geoid change. We also added a uniform global, nonsteric sea level rise of 1 mm/yr. In this study we did not impose conservation of total water on the planet.

Method
[22] We examine the generalized prediction error between models and simulated data to assess how well we might expect to recover the ''true'' viscosity profile from GRACE estimates of the PGR signal. We expect the noise in the Figure 1. (a) Degree amplitudes of the secular geoid signal for the reference PGR model (pluses); the GRACE measurement errors (solid black line); the present-day mass imbalance for Greenland (triangles) and for Antarctica (diamonds) using the best estimates from mass imbalance studies summarized by Church et al. [2001] and assuming that the ice changes are distributed uniformly over each ice sheet; a 1 mm/yr nonsteric sea level change (light gray line); the ocean currents (dark gray dashed line); the continental hydrology (small solid circles) and for the entire synthetic GRACE data (large light gray circles). The latter includes all of the other signals except PGR. (b) The number of Stokes' coefficients at a given value of l, for which the amplitude of the predicted PGR signal is larger than the expected secular GRACE measurement error. The PGR signal is predicted using the reference model described in section 2.1. The dotted line shows the total number of Stoke's coefficients as a function of l. Altogether, there are a total of 2166 Stokes' coefficients where the predicted PGR signal exceeds the expected GRACE measurement errors.
GRACE PGR estimates to fall somewhere between the expected GRACE instrumental error (in the event that it is possible to remove the effects of other secular signals) and the GRACE instrumental error plus all other secular signals (if it turns out that we cannot adequately remove them). In section 4 we will examine the recovery of the viscosities and lithospheric thickness for both these end-members, as well as for a case where we first fit and remove the largest contaminating signals. Here, we describe the method we use to determine the viscosity profile from the GRACE data in each of these cases.
[23] Inverse solution for parameters of a model can benefit significantly from a priori knowledge of the expected measurement error by minimizing the generalized prediction error R 2 = e T We, where e is a vector containing the misfit between a model and the observations and W is a diagonal tensor containing the inverse of the variance for each observation [e.g., Menke, 1984]. In the case of GRACE data, we have an a priori estimate of variance that is wavelength-dependent. Measurements of the geoid rate expressed in the spatial domain will contain all the errors at all wavelengths. Thus there is some advantage of expressing the norm of a PGR model as a generalized prediction error in the spherical harmonic domain, instead.
[24] We simulate solutions for the Earth's viscosity structure from GRACE data by constructing a synthetic secular signal similar to that which will be measured by GRACE, as described in section 2. We subtract a series of PGR models, one at a time, for a range of Earth viscosity structures and lithospheric thicknesses, varying two parameters at a time in a grid search algorithm. The PGR model that minimizes the generalized prediction error for the GRACE data is taken to be the preferred model. For our simulated data set, the generalized prediction error can be written in a normalized form as where s(l) corresponds to the estimated GRACE measurement error, _ C lm ref and _ S lm ref are the Stokes' coefficient rates for the reference (''true'' Earth) model of PGR, _ C lm err and _ S lm err are the coefficient rates of the GRACE errors (and other secular signals, if included), and _ C lm mod and _ S lm mod are the PGR coefficient rates predicted by the Earth model we are comparing with. The normalization of the L 2 norm results in a weighted average of the root mean square error between the model and the simulated data, rather than the c 2 parameter that would result if the weighting were not normalized.
[25] Using the norm in equation (2) and after finding R min 2 , i.e., the generalized prediction error for the model that best fits the measurements, one can estimate confidence intervals via the likelihood ratio method [Beck and Arnold, 1977]. If the errors are jointly normal, zero mean, and uncorrelated, then the confidence region having probability a of containing the solution corresponds to the entire volume of the model parameter space for which where M is the number of model parameters, n is the number of measurements, and F À1 is the inverse of the F cumulative distribution function.
[26] Our errors can include the non-PGR secular signals described in section 2. The secular Stokes' coefficients caused by these signals will be neither normally distributed or uncorrelated with the PGR model, and so the confidence levels computed using equation (3) are not strictly appropriate. We will conclude below that the most accurate results are obtained by first fitting and removing as much of the spatially coherent non-PGR signals as possible. This will reduce the nonnormal errors and the correlation of error with the PGR signal, and so will increase the relevance of equation (3).

Results
[27] To examine the potential of using GRACE data to discriminate between different plausible viscosity structures, we divide our analysis into three parts. First, we examine how accurately we can recover the ''true'' viscosity structure when GRACE measurement errors are the only source of error in the secular Stokes' coefficients (i.e., under the assumption that other secular signals can somehow be removed from the Stokes' coefficients either prior to or during the recovery of the PGR signal). Second, we examine the recovery of the Earth's viscosity structure when the other secular signals contribute to the Stokes' coefficients (corresponding to the scenario in which these signals cannot be removed). Third, we consider a case where the signals from present-day mass variability in Antarctica, Greenland, and globally averaged sea level are fit and removed before comparing the data with the PGR models.
[28] We consider radially symmetric viscosity structures with uniform viscosities in the lower mantle, the transition zone, and the upper mantle. The viscosity values for the reference model (i.e., the simulated ''true'' Earth model) are described in section 2.1. We compare the simulated GRACE data with models in which we vary pairs of parameters (e.g., lower mantle viscosity and upper mantle viscosity, lower mantle viscosity and lithospheric thickness, etc.) in a grid search fashion, while the remaining parameters were kept fixed to the reference values. We include one scenario in which the upper mantle and transition zone are assumed to have the same viscosity, which is varied simultaneously with the lower mantle viscosity. In each case, we calculate the generalized prediction error as defined in equation (2) for 441 models (i.e., a 21 Â 21 grid of the parameter space). Confidence intervals were estimated using equation (3), with M = 2.

Simulation With GRACE Errors Only
[29] Figure 2 depicts the square root of the generalized prediction error calculated using equation (2) for different pairs of parameters, for the case where GRACE measurement errors are the only non-PGR contribution to the synthetic GRACE data. In each panel the viscosity values that produce the minimum prediction error (the center of the cross) are almost exactly coincident with the reference model (the center of the diamond). The solid contours, which enclose such small areas on the scale of Figure 2 that they appear as dots at the centers of the crosses, represent the 65% and 95% confidence contours. The implication is that if GRACE measurement errors were the only error source, if the ice history is known accurately, and if our general model of an Earth with a layered, Newtonian viscosity profile were correct, then that viscosity profile could be determined with high accuracy using the GRACE data alone.
[30] The viscosity solution in these cases is well constrained both because the number of measurements n is large (we used all degrees and orders up to l = 70, giving n = 5037), and because the secularly varying geoid is highly sensitive to the viscosity values. Evidently, even a small change in viscosity causes enough of a secular change in the geoid to be statistically significant. The confidence intervals are approximately circular (which is impossible to verify The center of the diamond corresponds to the true answer, and the center of the cross corresponds to the best fit. The 65% and 95% confidence intervals are included. In Figures 2a -2c those confidence intervals are so tightly bound at the center of the cross that they appear as a dot. Note that viscosity is shown on a logarithmic scale, whereas lithospheric thickness is linear. from Figure 2), which indicates that the two parameters are about equally well determined and that there is not much cross correlation between them.

Simulation Including Other Signals
[31] GRACE instrumental errors will not be the only error source for the GRACE estimate of PGR. The other gravity signals affecting GRACE will have secular components which contaminate the PGR solution. Efforts can be made to reduce the effects of those contaminating signals. For example, in section 4.3 we will describe a method in which some of those signals are fit and removed from the data prior to the recovery of the PGR signal. However, the contaminating secular terms cannot be removed entirely, so that some or (in the worst case) all of them will remain in the GRACE solution. Consequently we examine here the worst-case scenario in which all other secular signals are present as errors in the GRACE data, and so contaminate the PGR solution.
[32] Figure 1a provides a preliminary indication of how serious the problem could be. In that figure, the solid circles are the degree amplitudes for the total secular error in the synthetic GRACE gravity data: the GRACE measurement errors plus the secular terms from all the non-PGR signals described in section 2.3. Note that this total error is far larger than the GRACE measurement errors alone, and is as large as or larger than the PGR degree amplitudes at degrees larger than about 10. Much of this secular error is caused by present-day ice mass changes in Antarctica and Greenland, which are also included in Figure 1a, calculated using the preferred estimates from mass imbalance studies summarized by Church et al. [2001]. Continental hydrology and nonsteric sea level change (also shown in Figure 1a) contribute significantly to the secular GRACE signal as well. However, the potential for hydrology and sea level signals to contaminate estimates of solid Earth properties will be less than implied by the degree amplitudes, because their spatial distribution is quite different from that of PGR. The present-day Antarctic and Greenland ice sheet signals, which have a spatial distribution similar to that of PGR, are likely to be of the same order as the PGR signal and larger than the GRACE measurement errors out to degrees of l = 30-40.
[33] These results support the concern that the effects of present-day polar mass imbalance could contaminate the GRACE PGR estimates. To assess how serious a problem this could be, we repeat the analysis described in section 4.1, except that here we include all the secular signals in the simulated GRACE measurements.
[34] Results for different pairs of parameters are shown in Figure 3. In each case the confidence intervals are larger than those obtained when only GRACE measurement errors are included (Figure 2). Furthermore, the preferred values from the GRACE data analysis do not agree with the ''correct'' values used to construct the simulation. For three of the four panels, the correct value lies outside the 95% confidence interval.
[35] This discrepancy arises because one of the underlying assumptions of the method we use to estimate confidence intervals (equation (3)) is that the signal and error are uncorrelated. That is not true here. In particular, the effects of present-day melting of Antarctic and Greenland ice are correlated with the PGR model, since that model includes the Earth's viscoelastic response to melting of Antarctic and Greenland ice during the late Pleistocene and Holocene. The gravity signal caused by the viscoelastic redistribution of rock underneath those ice sheets is highly correlated with the gravity signal caused by present-day changes in the ice sheets themselves.

Fitting and Removing Other Secular Signals
[36] The results in section 4.2 show that the presence of non-PGR secular terms that are correlated with the PGR signal can degrade the recovery of the Earth's viscosity profile. We expect that the effects of this correlation could be reduced by removing as much of that correlated signal as possible from the GRACE data.
[37] To do this, we expand the definition of our model. Instead of modeling the simulated GRACE data as the secular geoid coefficient rates for PGR corresponding to a particular Earth viscosity, we model the coefficient rates as PGR plus present-day secular trends in Antarctic ice mass, Greenland ice mass, and sea level. We first calculate three sets of secular geoid coefficient rates, corresponding to (1) a 1 cm/yr loss of ice mass uniformly distributed over Antarctica ( _ Y lm Ant ); (2) a 1 cm/yr ice mass loss uniformly distributed over Greenland ( _ Y lm Grn ), and (3) a 1 mm/yr water mass increase uniformly distributed over the oceans ( _ Y lm Slv ). Here, _ Y lm is used as a shorthand for ( _ C lm , _ S lm ). For each PGR model in the grid search, we first subtract the PGR model from the simulated GRACE data. Then, using linear weighted least squares, we solve for the constant coefficients a Ant , a Grn , and a Slv that minimize an expanded definition of the generalized prediction error in which we substitute into equation (2): Fitting and removing the ice mass and sea level signals in this manner should not only reduce the correlation between the PGR and the residual secular terms in the gravity field: It should also remove the largest nonnormally distributed portion of the GRACE errors, and so would increase the relevance of equation (3) for estimating confidence limits. Although other contributions (particularly from continental hydrology and ocean currents) are as large as or larger than the ice mass and sea level change signals, we do not attempt to fit these because they are not highly correlated to the PGR signal and because they can not be approximated by a simple uniform change over some geographical region.
[38] We have applied this procedure to our synthetic data. Results are shown in Figure 4 for a grid search over pairs of parameters. Note that now the preferred and correct values agree to well within the confidence intervals. Those intervals are still much larger than those shown in Figure 2 for the case where only GRACE measurement errors are included. In fact, the area of the confidence intervals is not much different than that shown in Figure 3, where the contaminating secular signals have not been fit and removed.
[39] In general, the 95% confidence limits on the lower mantle viscosity are between about 6 Â 10 21 and 14 Â 10 21 Pa s. Since the value used to create the synthetic GRACE data was 10 Â 10 21 Ps s, we infer that GRACE is capable of inferring the lower mantle viscosity to about ±40%. Similarly, for the combined upper mantle/transition zone viscosity shown in Figure 4a, the 95% confidence interval extends between 0.75 Â 10 21 and 1.3 Â 10 21 Pa s, which agrees with the correct value (1.0 Â 10 21 Pa s) to within ±30%. The results shown in Figure 4d suggest that the lithospheric thickness can be recovered to within ±15 -20% with 95% confidence (the 95% confidence interval extends between 80 and 115 km).
[40] Figure 4c shows that GRACE would have a harder time discriminating between separate viscosities in the upper mantle and transition zone. For each of those parameters, the 95% confidence interval includes numbers that are between one half and twice the correct value.
[41] So far, we have assumed there are no errors in the ice history model used to calculate the PGR signals. Although this assumption is certainly false, its effects are difficult to determine since there is no obvious way to estimate the errors in an ice model. To obtain what we hope is an upper bound on these effects, we consider a second ice history model: ICE-4G from Peltier [1994]. ICE-4G is an improved model that uses new constraints which were not available when ICE-3G was developed. [42] Our approach is to generate the synthetic GRACE data using a PGR signal calculated from ICE-3G; but to compute the PGR signals used in the grid search by using ICE-4G. The synthetic GRACE data include all of the non-PGR secular signals. The results, shown in Figure 5, show that the differences between the ice models can have significant effects on the recovered viscosity. The preferred upper mantle/transition zone viscosity is about 1.6 Â 10 21 Pa s, or 1.6 times the correct value. The preferred lower mantle mantle viscosity is (6 Â 10 21 Pa s), 40% less than the correct value. It is of course possible that this analysis underestimates the effects of true uncertainty in the ice model, as ICE-3G and ICE-4G were derived by the same group and all of the data used to constrain ICE-3G was also incorporated into ICE-4G. However, this is perhaps the best estimate we can make for potential uncertainties.

Discussion and Conclusions
[43] The comparisons of simulated GRACE measurements and Earth models presented in this paper suggest that by using GRACE data alone, it should be possible to obtain useful information about the Earth's viscosity structure. It had been suggested previously that a time variable satellite gravity mission could provide useful information on lower mantle viscosity [Dickey et al., 1997]. The results described above confirm that suggestion, indicating that GRACE could conceivably determine the lower mantle viscosity, with 95% confidence, to ±40%. However, what  Figure 3, but where the signals from present-day changes in Antarctic and Greenland ice and from a global sea level rise have been fit and removed from the synthetic GRACE data prior to recovering the viscosity profile.
we found surprising was that GRACE could also determine the combined viscosity of the upper mantle/transition zone, and to even slightly better accuracy: ±30%. We found that we could recover the lithospheric thickness to within ±15 -20%, which was also much better than we would have expected. It will be more difficult to use GRACE to solve for separate viscosities in the upper mantle and transition zone. Those parameters could probably only be determined to within a factor of 2. Note that the accuracy of recovery of the true Earth viscosity will depend in part on how well the model parameterization used for the grid search can represent the true Earth structure.
[44] The accurate recovery of the model parameters is made possible in part by weighting of the coefficients according to the inverse variance of GRACE measurement errors in defining the generalized prediction error (equation (2)). We also performed a grid search in which the minimum was defined for an unweighted-residual error norm. Without error weighting, best fit solutions for lower mantle viscosity differed by a factor of 5 from the input reference model, and best fit viscosities of other layers differed by even more. Clearly, the inverse variance weighting is crucial to achieve the performance described here. These simulations assumed that the covariance matrix of GRACE secular coefficients is perfectly diagonal. Full covariance matrices to degree 120 have been estimated at JPL and Texas from numerical simulations of the GRACE gravity recovery (M. Watkins, personal communication, 1999). However, these have not been used in the science analyses of GRACE applications due to the numerical cumbersomeness of the extremely large files and due to the fact that the assumption of diagonality is generally quite reasonable, especially for the low to middle degree and order terms. Since the correlation properties are also affected by the parameterization used for the gravity solution, the effect of the full covariance on science products is a topic for additional study when the actual GRACE gravity products become available. However, the analysis described in this paper can easily be generalized to use the full covariance matrix by substituting the matrix inverse for W in the expression for the generalized prediction error R 2 .
[45] The results these conclusions are based on all use the same relatively simple reference model. To examine whether there is significant dependence on the choice of reference model, we performed the calculations described above using different reference models, and found that the solution error space can be somewhat sensitive to the ''true Earth'' viscosity structure. For example, using a reference model with a lower mantle viscosity of 5 Â 10 23 and an upper mantle and transition zone viscosity of 6 Â 10 21 , and including GRACE measurement errors but no non-PGR secular signals, causes an additional local minimum in the grid search over lower mantle versus coupled upper mantle/ transition zone viscosities ( Figure 6). Nevertheless, the grid search recovers the correct answer and the confidence intervals remain small, much as before.
[46] The results described at the end of section 4.3 suggest that a limiting factor in recovering the Earth's viscosity structure using GRACE data may be uncertainties in the ice loading history, rather than GRACE measurement errors or contamination from non-PGR secular signals. The effects of errors in the ice model can conceivably be reduced by parameterizing that model and using the data to solve for those parameters and the viscosity profile simultaneously [e.g., Lambeck et al., 1998; L. Tarasov and W. R. Peltier, Greenland glacial history and local geodynamic consequences, submitted to Geophysical Journal International, 2001].
[47] It is unlikely that the GRACE data alone will have enough resolving power to permit such a comprehensive inversion. The problem would be much better constrained if other PGR data types were combined with the GRACE measurements. The upper mantle/transition zone recovery, although surprisingly good, is probably not as accurate as results that can already be obtained using geological measurements of relative sea level change and geodetic observations of present-day uplift rates. What the GRACE data can provide is sensitivity to lower mantle viscosity that is more accurate and less ambiguous in its interpretation than any other PGR data type. It is the combination of PGR data types, along with the sort of information normally used to help constrain the ice models (e.g., positions of moraines, beach heights, coral based data) that should provide the most dramatic improvement in estimates of the Earth's viscosity profile.