A sufficiently accurate grid of crustal thickness throughout the oceanic realm is not yet avail- On the continents, we have estimated dynamic topography by scaling long-wavelength free-9 air gravity anomalies using a constant admittance of +50 mGal km −1 . It is important to excise 10 all regions where gravity anomalies are affected by active/ancient orogenic belts and by flexure 11 at subduction zones. We used a combination of gravity, topographic and geological maps to 12 manually excise these regions. The resultant exclusion polygons are shown in Figure S1b.

We perform joint nonlinear inversions of glacial isostatic adjustment (GIA) data, including the following: postglacial decay times in Canada and Scandinavia, the Fennoscandian relaxation spectrum (FRS), late‐Holocene differential sea level (DSL) highstands (based on recent compilations of Australian sea level histories), and the rate of change of the degree 2 zonal harmonic of the geopotential, J2. Resolving power analyses demonstrate the following: (1) the FRS constrains mean upper mantle viscosity to be ∼3 × 1020 Pa s, (2) postglacial decay time data require the average viscosity in the top ∼1500 km of the mantle to be 1021 Pa s, and (3) the J2 datum constrains mean lower mantle viscosity to be ∼5 × 1021 Pa s. To reconcile (2) and (3), viscosity must increase to 1022–1023 Pa s in the deep mantle. Our analysis highlights the importance of accurately correcting the J2 observation for modern glacier melting in order to robustly infer deep mantle viscosity. We also perform a large series of forward calculations to investigate the compatibility of the GIA data sets with a viscosity jump within the lower mantle, as suggested by geodynamic and seismic studies, and conclude that the GIA data may accommodate a sharp jump of 1–2 orders of magnitude in viscosity across a boundary placed in a depth range of 1000–1700 km but does not require such a feature. Finally, we find that no 1‐D viscosity profile appears capable of simultaneously reconciling the DSL highstand data and suggest that this discord is likely due to laterally heterogeneous mantle viscosity, an issue we explore in a companion study.

Earth's body tide-also known as the solid Earth tide, the displacement of the solid Earth's surface caused by gravitational forces from the Moon and the Sun-is sensitive to the density of the two Large Low Shear Velocity Provinces (LLSVPs) beneath Africa and the Pacific. These massive regions extend approximately 1,000 kilometres upward from the base of the mantle and their buoyancy remains actively debated within the geophysical community. Here we use tidal tomography to constrain Earth's deep-mantle buoyancy derived from Global Positioning System (GPS)-based measurements of semi-diurnal body tide deformation. Using a probabilistic approach, we show that across the bottom two-thirds of the two LLSVPs the mean density is about 0.5 per cent higher than the average mantle density across this depth range (that is, its mean buoyancy is minus 0.5 per cent), although this anomaly may be concentrated towards the very base of the mantle. We conclude that the buoyancy of these structures is dominated by the enrichment of high-density chemical components, probably related to subducted oceanic plates or primordial material associated with Earth's formation. Because the dynamics of the mantle is driven by density variations, our result has important dynamical implications for the stability of the LLSVPs and the long-term evolution of the Earth system.

We present a method for calculating the derivatives of measurements of glacial isostatic adjustment (GIA) with respect to the viscosity structure of the Earth and the ice sheet history. These derivatives, or kernels, quantify the linearised sensitivity of measurements to the underlying model parameters. The adjoint method is used to enable efficient calculation of theoretically exact sensitivity kernels within laterally heterogeneous earth models that can have a range of linear or non-linear viscoelastic rheologies. We first present a new approach to calculate GIA in the time domain, which, in contrast to the more usual formulation in the Laplace domain, is well suited to continuously varying earth models and to the use of the adjoint method. Benchmarking results show excellent agreement between our formulation and previous methods. We illustrate the potential applications of the kernels calculated in this way through a range of numerical calculations relative to a spherically symmetric background model. The complex spatial patterns of the sensitivities are not intuitive, and this is the first time that such effects are quantified in an efficient and accurate manner.

No abstract