We develop a new formalism for the prediction of secular variations in the gravitational potential field of a spherically symmetric, self-gravitating, (Maxwell) viscoelastic planetary model subjected to an arbitary surface load which may include a gravitationally self-consistent ocean loading component. The theory is applied to generate the most accurate predictions to date, of the present-day secular variations in the zonal harmonics of the geopotential (the so-called J˙l for degree l) arising as a consequence of the late Pleistocene glacial cycles. In this respect, we use the very recent ICE-3G reconstruction of the last late Pleistocene deglaciation event (Tushingham and Peltier, 1991). A comparison of these predictions with those generated using simplified disk models of the ice sheets, which have been used in all previous studies of the J˙l harmonics (l>2), (indicates that the disk model approximation introduces unacceptably large errors at all spherical harmonic degrees except perhaps l=2. Predictions have also been made using a eustatic loading approximation (also used in previous studies) in place of a gravitationally self-consistent ocean loading component, and we have found that the resulting discrepancy is largest at degrees 2, 8 and 10. In the case of J˙2 the magnitude of the error incurred using the eustatic approximation can be as large as order 10--15% of the predicted value. We have attributed this discrepancy to the present day net flux of water away from the equatorial regions arising from the remanant present-day adjustment associated with the late Pleistocene glacial cycles. The effect represents a heretofore unrecognized contribution to the J˙l harmonics, or alternatively the nontidal acceleration of Earth's axial rate of rotation. In terms of the later, the maximum anormaly in the length of day is approximately 1.7 μs/yr.

We also consider the sensitivity of the J˙l data to variations in the radial mantle viscosity profile by using a suite of forward calculations and an examination of Fr¿chet kernls. The theory required for the computation of those kernls is described herein. We find that the radial variation in sensitivity can be a strong function of the viscosity model used in the calculations. For models with a uniform upper mantle viscosity (vum) of 1021 Pa s, forward predictions of the J˙l harmonics exhibit a pronounced peak when a wide enough range of lower mantle viscosities (vLM) are considered (we denote the vLM value at this peak as vˆlLM). At the lowest degrees (l≤4), Fr¿chet kernels computed for a series of increasing vLM values (1021 Pa s≤vLM<1024 Pa s) indicate a migration of the dominant sensitivity of the J˙l data to variations in viscosity from regions below approximately 1200 km depth (for vLM≤vˆlLM) to regions above this depth in the lower mantle (for vLM≥vˆlLM). The sensitivity of the J˙l data to variations in the viscosity profile in the shallowest parts of the lower mantle, for the case vLM≤vˆlLM, is also reflected in a set of forward calculations described herein. As an example, J˙l predictions made using Earth models in which the viscosity above 1200 km depth is constrained to be 1021 Pa s, do not exhibit the multiple solutions characteristic of the vUM=1021 Pa s calculations.

The same is true of Earth models in which the upper mantle viscosity is weakened an order of magnitude to 1020 Pa s. The theory described herein is also applied to commute the J˙l signal (l≤10) arising from the retreat of small ice sheets and glaciers described by Meier (1984) and also from any potential variations in the mass of the Antarctic and Greenland ice sheets. The present day J˙l signal due to the late Pleistocene glacial cycles dominates the signal from Meier's sources at all degrees except l=3. In contrast, the J˙l