# Dense mantle flows periodically spaced below ocean basins

Isabelle Panet, Marianne Greff-Lefftz, Barbara Romanowicz

<figure><img src="image_1.png" /></figure>

## To cite this version:

Isabelle Panet, Marianne Greff-Lefftz, Barbara Romanowicz. Dense mantle flows periodically spaced below ocean basins. Earth and Planetary Science Letters, 2022, 594, 10.1016/j.epsl.2022.117745. insu-03776402

HAL Id: insu-03776402

https://insu.hal.science/insu-03776402v1

Submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dense mantle flows periodically spaced below ocean basins

Isabelle Panet<sup>a,b,*</sup>, Marianne Greff-Lefftz<sup>a</sup>, Barbara Romanowicz<sup>a,c,d</sup>

<sup>a</sup>Université Paris Cité, Institut de physique du globe de Paris, CNRS, IGN, F-75005 Paris, France

<sup>b</sup>ENSG-Géomatique, IGN, F-77455 Marne-la-Valle, France

<sup>c</sup>Collège de France, Paris

<sup>d</sup>Berkeley Seismological Laboratory, University of California, Berkeley

**Abstract**

Understanding mantle flow is key to elucidate how deep Earth dynamics relate to tectonics at the global scale. The convective mass transport is reflected in lateral variations of the gravity field, seismic velocities, as well as deformations of the Earth's surface. Yet, upper to mid-mantle dynamics have been difficult to constrain at the medium scales of thousands of km. Here, we analyze the second-order horizontal derivatives of seafloor topography and of the gravity potential over the Pacific and Northern Indian ocean basins, and provide evidence for periodic undulations of 1600-2000 km wavelength in both signals, elongated along the direction of absolute plate motion. We investigate potential crustal and lithospheric sources and show that at least part of this signal must originate below the lithosphere, with alignments of sub-lithospheric upper mantle mass excess below seafloor lows. Furthermore, we find that these alignments coincide geographically over wide areas with similarly periodic slow seismic velocity fingers located at upper mantle depths. These two fields may thus record an intermediate scale of mantle convection below ocean basins, which cannot be explained by purely thermal convection and requires instead lateral variations in composition in the upper mantle. Elucidating the nature of the detected mass excess sources coincident with the slow seismic velocities calls for a joint dynamical modelling of all observations in a thermo-chemical context.

**Keywords:** Mantle mass structure, Gravity gradients, Seafloor topography, Seismic tomography

2010 MSC: 00-01, 99-00

<footnote>*Corresponding author</footnote>
Introduction

Imaging mantle structure at high spatial resolution is a major goal in geophysics, as we are still striving to constrain the planform of convection and elucidate whether the mantle convects in one or several layers that could involve motions at various spatial scales. Lateral variations in the gravity field, surface deformations and seismic velocities can provide complementary constraints on the patterns of mantle flow thanks to their different radial and lateral sensitivities to mantle mass transport and elastic structure. At the largest scales, global mass redistributions in the Earth's mantle associated with plate subduction is informed by the coincidence of long-wavelength geoid lows with circum-Pacific fast seismic velocities in the lower mantle, first documented decades ago [1] and attributed to the graveyard of slabs [2]. While more difficult to resolve, various dynamic features have been proposed to exist at short to medium scales. They include thin and vigorous plumes rising vertically from the bottom of the mantle to the surface to explain hotspot volcanism [3], or secondary scale convection driven by moving plates in the upper mantle [4, 5]. Deciphering these different features in the geophysical observables and reconciling them with seismic tomographic images of the mantle remains a challenge.

At intermediate lateral scales, enigmatic signals have been found in the geodetic and seismological data. In two dimensions along a profile between Tonga and Hawaii roughly orthogonal to the Pacific absolute plate motion (APM), a regional 2D seismic tomographic study, combined with geoid and seafloor topography data, showed a series of ~1500-km wavelength seafloor swells and geoid highs correlated with similarly periodic high shear velocities extending to depths in excess of 200 km [6]. As pointed out by the authors, this relationship is puzzling from the dynamic point of view: a thermal interpretation of these fast seismic velocities, representing cold anomalies and therefore downwelling flow, predicts seafloor lows where the highs are observed. Recently, global 3D upper mantle seismic shear velocity tomography at improved lateral resolution [7, 8] revealed quasi-periodic, ~1800-2000 km wavelength horizontally elongated slow velocity "fingers" aligned with present-day APM, ubiquitous below all ocean basins. At the same time, a directional geoid fabric at the same wavelength and along a similar direction was detected in GRACE satellite gravity data in the Pacific Ocean [9], beyond the 100's km scales of the previously known geoid undulations derived from satellite altimetry over ocean basins [10, 11, 12, 13, 14]. A dynamically unexpected relation with the slow seismic signals was again observed [8]. Today, the processes at the origin of these combined geodetic and seismic observations remain unresolved, and the characterization of the
signals in the geodetic data needs to be pursued. Their predominance in the gravity data has been demonstrated
in [9], but their spatial mapping still needs to be performed together with investigations of a possible seafloor
topographic expression.

Here, we continue and extend these previous studies by building an exhaustive set of geodetic observations in
order to constrain the processes at the origin of these elongated structures, and in particular, the depth at which
they originate. We search for the spatial signature of medium-scale, APM-oriented geoid and seafloor topography
signals over broad oceanic domains and assess their consistency with the slow velocity fingers of the SEMum2
upper mantle tomographic model [8]. We analyze whether the obtained signals can originate from the crust or
the lithosphere by comparing their pattern and amplitude with those of potential crustal or lithospheric sources.
We conclude that at least part of the observed anomalies must have a dynamical origin within the convecting
mantle.

# 1. Data and Methods

We focus our analysis primarily on the Pacific ocean, since the slow velocity fingers in the SEMum2 upper
mantle tomographic model [8] are most prominent below this fast-moving plate. We also consider the fast
Australian-Indian plate in the northern Indian ocean and the slow-moving Antarctic plate. In these areas, we
analyze global geoid and seafloor topography datasets as described in SI Section A. The geoid is based on GRACE
and GOCE satellite gravity data, and the seafloor topography is based on a global bathymetric grid, corrected
for the isostatic contribution of the load due to sediments. Rather than remove a priori models of shallow layers
(e.g. [15]), we thus seek to extract the target signal in direction and scale directly from the geodetic data, after
correcting only for sediment loading. In these datasets, we search for medium-scale undulations using an analysis
of second-order horizontal gradients.

## 1.1. *Geopotential and seafloor topography second-order gradients*

The gravity field is the first-order gradient of the gravity potential. The so-called gravity gradients, used in
this study, are the second-order gradients of this potential. In Cartesian coordinates, these gradients represent the
spatial derivatives of the gravity field in the three different directions of space, transformed through appropriate
rotations into spherical coordinates [16]. Because of the second-order differentiation, the spatial resolution of the
gradients above the source is higher compared to that of the geoid. Furthermore, these directional derivatives
efficiently bring out elongated features orthogonal to the direction of differentiation (SI Section B). They have
been used either individually or combined, to investigate lithosphere and deeper mantle structure at high spatial
resolution, e.g. [17, 18, 19, 20, 21]. The first and second-order gradients of the seafloor topography can be defined
in a similar fashion (SI Section B). Here, we refer to the first-order gradients of the seafloor topography as the
seafloor slopes, and to the second-order gradients as the seafloor slope gradients. While the first-order spatial
gradients record the orientation of the gravity vector and that of the seafloor topographic slope, the second-order
gradients highlight their rapid changes in the vicinity of localized mass sources.

At a given geographical location, we build a multi-scale gradient analysis by computing the horizontal second-
order gradients in the local spherical frame from the wavelet-filtered geoid and seafloor topography, for wavelet
scales from 500 to 3600-km. The methodology is presented in detail in [16]. This way, we obtain gravity
gradients and seafloor slope gradients at different spatial scales, which provide us with an enhanced description
of the geometry of the geoid and seafloor topography signals at each spatial scale. In addition, we progressively
rotate the spherical frame around the radial axis in order to search for directions that show significant strength.
For instance, when the *e*θ basis vector of the spherical frame (Suppl. Fig. B1) is orthogonal to a plate motion
direction, the θθ gradients emphasize signals orthogonal to eθ, thus parallel to the considered plate motion
direction. Note finally that, in order to increase the signal-to-noise ratio in real data analyses, we may also stack
the gradient maps for a range of frame rotation angles.

## 1.2. Scale-orientation diagrams

At each spatial point in the considered areas, we compute the horizontal gravity and seafloor slope gradients
for each scale and orientation of the spherical frame (note that the θθ and φφ gradients are simply related by
a 90° rotation around the radial axis). At a given location, their amplitude variations as a function of scale
and orientation provides a pointwise bi-dimensional spectrum [16]. Its maxima indicate a well-resolved elongated
signal at the considered location, and provide the corresponding characteristic scale and orientation. For all points
spanning a chosen region, we then count the number of occurrences of each characteristic scale and orientation
identified in the pointwise spectra, to derive regional scale-orientation diagrams and marginal distributions of
scales for a selection of characteristic orientations. These regional spectra provide a diagnostic of the different
contributions to the total gravity and seafloor topography variations in a chosen region as a function of their scale and orientation. They are a useful tool to investigate whether the gravity field and the seafloor topography signals exhibit a preferred orientation over ocean basins, and find the corresponding spatial scale.

## 2. Prominent signals at medium spatial scales

### 2.1. Predominant directions from scale-orientation diagrams

We compute the gravity gradients and seafloor slope gradients as explained above, and first build regional scale-orientation diagrams in the Central and South Central Pacific, near Antartica and in the Northern Indian ocean, for a 0.25° spacing of the grid points. Over wide areas in all investigated ocean basins, we find that the distribution of characteristic orientations systematically peaks around the present-day (in the Pacific ocean) or recent (in the Indian Ocean) APM direction [22] for intermediate scales (Fig. 1ab; SI Section C). South of India, the orientation of the medium-scale signals is close to that of the recent plate motion (present-day to 30 Myrs ago, label 1 in Suppl. Fig. C3). East of the Ninety East Ridge they can be tilted by up to ~15° with respect to the APM (label 1 in Suppl. Fig. C3). Note that the existence of these ~1000-km scale signals (i.e. ~2000 km wavelength) is furthermore marked by a distinct peak, often prominent, in the distribution of scales obtained within transverse profiles orthogonal to the present-day APM, for characteristic orientations selected within a range around the present-day APM (Figs. 1c, 2af; SI Section D). Suppl. Table D1 summarizes the characteristic scales of the gravity and seafloor topography signals found on these profiles. Signals from other known structures are also identified in these diagrams, such as those aligned with the end-of-Cretaceous (~70 Myrs ago) APM direction [42] in the case of the Line Islands chain, or following the ~N80°E oriented fracture zones, which delimit lithosphere sections of different age and thickness.

### 2.2. Maps of the APM-oriented signals

We then reconstruct the APM-oriented signals identified in the previous diagrams in the space domain, by computing the gravity gradients and seafloor slope gradients for the corresponding scales and spherical frame orientations. Accordingly, the obtained maps of the GRACE gravity gradients and seafloor slope gradients reveal a series of medium-scale undulations aligned over thousands of kilometers with the present-day APM direction in the investigated basins (Figs. 1de, 2bcgh; SI Section E). As could be expected from the scale-orientation
diagrams, these elongated anomalies vanish when rotating the spherical frame away from the APM-orientation (Movies E4-E5). Over the Central Pacific and the Northern Indian oceans, this directional fabric appears consistently in both the total gravity field and the seafloor relief, with a ~3-4 larger amplitude in the older Pacific plate, in the northern Indian ocean and near Antarctica as compared to the younger Pacific plate. It can also be recognized directly in the geoid and the seafloor topography by band-pass filtering of these datasets at intermediate wavelengths as done in [6] (SI Section F), even in the presence of other superimposed signals. Suppl. Table E1 summarizes the amplitudes of the gradient undulations obtained by visual inspection of the maps and cross-sections, at the 1100-km scale and at their characteristic scale. Finally, the consistency between the gravity gradient and the seafloor slope gradient undulations is confirmed by a correlation analysis over the corresponding ocean basins (SI Section G). Over the Central Pacific ocean, the correlation between these two fields amounts to 0.67, with similar values over the younger or the older side of the Pacific plate. Over the Northern Indian ocean, the correlation is even higher and reaches 0.77.

Before investigating the depth range in which these features may originate, let us discuss some simple synthetic examples, that inform us on the source structure of these gradient anomalies. Their characteristic scale corresponds to the half-wavelength of sinusoidal sources [16]. Note that the sign of the diagonal horizontal ($\theta\theta$ or $\phi\phi$) gradients computed along a direction orthogonal to the main axis of a given mass source is opposite to that of its geoid signal for both oscillating and localized sources [16], whereas the radial (rr) gravity gradients keep the same sign as the geoid. Thus, the direct Newtonian attraction of an elongated *mass excess* generates a *geoid high* and a $\theta\theta$ or $\phi\phi$ *gradient low* above the source, as illustrated in SI Section B for a 600-km wide, 2400-km long, 10 kg.m$^{-3}$ rectangular mantle mass excess source located between 110 and 210-km depth within a rigid Earth, oriented along the azimuth $N120^\circ E$. Fig. 3a shows the gravity gradient signal for a cross-section orthogonal to the source major axis. In the same way as for gravity, horizontal *seafloor slopes gradients highs* reflect *seafloor topography lows*, as illustrated in a similar cross-section in Fig. 3b for a 530m high, 600-km wide Gaussian-shaped seafloor depression following the same azimuth.

Fig. 3c illustrates the gravity gradient signals of a vertically equilibrated system of two opposite mass sources, the seafloor topography low (Fig. 3b) and the mantle source of mass excess (Fig. 3a). The sign of the total
gravity gradient signal indicates the predominant contribution to the gravity field, i.e. the seafloor topography in this example. Subtracting the topographic contribution from the observed total gravity, we obtain Bouguer gravity gradients which provide information on the support of this topographic relief. Here, large negative anomalies in the Bouguer gradients reflect internal sources of mass excess below the topographic lows, as in Fig. 3a. In the framework of isostasy, the support of the topography could be provided, for instance, by variations in crustal thickness. In the case of a dynamic support in a deformable Earth, a negative Bouguer gravity gradient anomaly indicates that, excluding the surface topography, the combined gravitational effect of all the internal mass anomalies supporting this topography (the moving source, the induced internal interface deflections due to phase changes and viscous flow and the induced density variations in a compressible Earth), is equivalent to that of a source with mass excess at depth.

Such lines of dense sources clearly show up in the Bouguer gravity gradient maps built from our observations over the Central Pacific and the Northern Indian Ocean (Fig. 4ab; note the flipped colorbars in the Bouguer maps). These maps feature large negative residual anomalies that coincide with maxima in the seafloor slope gradients, and therefore correspond to alignments of dense mass sources beneath the seafloor lows.

## 2.3. *Comparison with upper mantle seismic tomography*

As part of our investigation of the potential sources of these geodetic signals, we tested their correspondence with the seismic structure of the upper mantle below all investigated ocean basins. We observe that the APM orientation and the nearly 1000-km half-period of the geodetic anomalies coincide with those of the slow seismic velocity fingers, best detected in the 60-350 km depth range of the SEMum2 [8] upper mantle tomographic model (Fig. 2di, 5b), but extending to depths in excess of 660 km in some areas (e.g. Fig 5d). Comparing the locations of the SEMum2 slow velocity anomalies with those of the positive and the negative gradient anomalies, we find that in many areas, the seafloor topographic lows (reflected in the positive seafloor slope gradients) are located above the slow seismic fingers, and the seafloor highs, between the fingers. This is emphasized in profiles taken in directions orthogonal to these fingers, which also illustrate the periodicity of the signals and, in addition, the correlations with the gravity gradients (Figs. 2ej, 5c). In these profiles, the correspondence of the negative Bouguer gradients with the slow seismic fingers indicates that mantle mass excess laterally coincides with the slow seismic velocities. Generalizing the 2D observation along the Tonga-Hawaii profile [6], we detect this surprising
correlation west of the East Pacific Rise, south of India and in the western part of the northern Indian ocean (Fig. 2e). On the Antarctic plate, the continental borders and ridges approximately aligned with the APM limit the detectability of such patterns especially in the seafloor slope gradients, still two positive total gravity gradient undulations overlie velocity lows between the ridge and the continent (Fig. 2g).

170

## 3. Investigation of a crustal origin for the APM-oriented medium-scale signals

We first tested a crustal origin for the APM-oriented mass excess sources beneath the seafloor lows, focusing on the sediment correction and on the Moho topography as a major crustal contribution due to the significant crust/mantle density contrast.

### 3.1. Sediment correction

175

We estimated the seafloor slope gradient signals of the isostatic sediment correction applied to the bathymetric grid (SI Section A), and also analyzed the correction obtained from another global model of sediment thickness [30], which is used in CRUST1 [27] and LITHO1.0 [28]. Over the Central Pacific, the amplitude of the sediment signal does not exceed ~ 10% of the discussed seafloor signals, except in an equatorial band north of the Polynesian Superswell, where it reaches ~ 50% of the seafloor signals. The spatial pattern of these sediment signals reflects older, fracture zone directions rather than recent ones. Over the Northern Indian ocean, the sediment signal does not exceed ~ 15% of the seafloor signals except off the Eastern coast of India above the equator, where it reaches ~ 50% of the seafloor signals. Except there, it does not exhibit an APM-oriented structure as observed in the geodetic data. Thus, errors from this contribution are not expected to significantly change our results.

180

### 3.2. Isostatic crust models

185

We computed the Moho depth and the gravity gradients predicted for a crustal isostatic model, with contributions from seafloor topography, Moho depth and LAB (lithosphere-asthenosphere boundary) depth. We consider two different cases: 1. crustal isostasy without a lithosphere, and 2. isostasy including a lithospheric layer. The principle of the calculation of the isostatic Moho depth in each of these two cases is presented in SI Section H.

190
We first considered the amplitude of the isostatic Moho topography. In case 1, the amplitude of the iso-static Moho variations is proportional to the amplitude of the seafloor topographic undulations, which can be derived from the observed seafloor slope gradients using forward modelling. In the Northern Indian Ocean, the 70m/(100km)² seafloor slope gradient signal at 1100-km scale thus corresponds to a 700m amplitude, 1000km half-period sinusoidal topography, and an isostatic Moho undulation of amplitude ~ 3.8km (peak-to-peak 7.6km). The peak-to-peak amplitude, which reflects the relative variations in crustal thickness, is well above the 3-5km level of error of the Moho model by [24] in the Indian Ocean and commensurate with the average oceanic crustal thickness, and yet these large undulations over broad areas are hardly found in the investigated seismological Moho models. This rules out a purely isostatic Moho support of the observed seafloor topographies in the North-ern Indian Ocean. On the younger side of the Central Pacific, the obtained amplitude of the isostatic Moho depth undulations is smaller (1.1 km, corresponding to a 2.2 km peak-to-peak amplitude), and we investigate the gravity gradients of the crustal isostatic model.

We compute the gravity gradients of the crustal isostatic model for the Pacific ocean by summing the con-tributions from the seafloor topography, the Moho depth and in a second calculation (corresponding to case 2), the LAB depth. The results are presented in Fig. 6. In both cases, the pattern of undulations has weak amplitudes (± 100 mEötvös), compared to the observed undulations in the GRACE data (± 400 mEötvös, e.g., Fig. 1d) over most of the Pacific plate, because the Moho mass source compensating the seafloor topography variations is close to the seafloor source. This conclusion is also supported when computing the amplitude ratio of the topographic gravity gradients versus the GRACE gravity gradients (SI Section G2). This ratio can be viewed as an adaptation to the second-order gradients of the geoid-to-topography ratio, which informs on the compensation depth of the topography and can be compared to modelled values (e.g. [25]). Over the Central Pacific ocean basin, the observed ratio (6.7) is lower than predicted from the isostatic crust models (10.7 to 11.2, with or without a cooling lithosphere). Thus the observed GRACE signal cannot be explained fully by an isostatic support of the seafloor undulations by crustal thickness variations. In addition, the gravity gradients scale-orientation diagram for the isostatic crust model does not show any structure along the APM orientation (see Fig. 6ab for both cases 1 and 2). Instead, the spatial structure of this crustal model is dominated by signals along the fracture zone orientations. This difference in spatial structure of the crustal gravity gradients
compared to the seafloor and the Moho undulations (Suppl. Fig. H1) probably results from the near-cancellation of these two contributions in the gravity gradients, as the distance between these two opposite mass sources is small as compared to the scale of the signals, combined with spatial variations in the attenuation of the gravity signal of the Moho source as a function of its depth below the Earth's surface. This depth depends in particular on the spatial variations in seafloor depth. Accordingly, we found a lower correlation between the predicted crustal gravity gradients and the seafloor topographic signal than the observed correlation (Section 2.2).

These conclusions are consistent with results by [15, 26], who generally omit the isostatic crustal correction and estimate dynamic topographies over the oceans after plate cooling corrections and excision of the areas of anomalous crustal thickness. We calculated the seafloor slope gradients derived from their model of dynamic topography and found it in broad agreement with our results based on seafloor topographic data uncorrected for crustal isostasy and lithospheric plate cooling (SI Section J). Additionally, we confirmed the limited impact of the cooling plate correction on the discussed signals, and their robustness to corrections for anomalous crustal thickening around seamounts and plateaus (SI Section I).

### 3.3. Seismological Moho models

We compared the predictions obtained for the isostatic Moho undulations with the 2nd-order gradients of the Moho depths in two different global seismological crustal models: the model of Szwillus et al. [24], obtained by geostatistical interpolation of seismological data from the USGS database, which comprises mainly active seismic measurements, and the LITHO1.0 model [28], in which the CRUST1.0 crustal model of constant crustal thickness below the seafloor topography is modified in order to fit surface wave dispersion data. We note that the Moho depth model of Szwillus et al. [24] (Suppl. Fig. H1c) is weaker, with amplitudes several times smaller than expected from isostatic models of the oceanic crust (Suppl. Fig. H1a,b). Its gravity signal is far too small to compensate that of the seafloor reliefs as observed; the same holds for the Northern Indian ocean (not shown). The LITHO1.0 model has no significant structuration along the APM-orientation corresponding to the observed one, as shown by scale-orientation diagrams over the Pacific and the Northern Indian oceans (Fig. 7). In these diagrams, we do not find any structure corresponding to the GRACE gravity gradients and the seafloor slope gradient signals, marked by the black ellipses. Note that we did not consider CRUST1.0 [27], which is inappropriate for our study, as it has a constant thickness of the crust in the oceans, and its Moho topography
is the same as the seafloor topography.

### 3.4. Lateral variations in crustal structure above a hot mantle?

We finally investigated how crustal density structure might change in the presence of a hot underlying mantle. The slower than average shear velocity "fingers" in the seismic model SEMum2 suggest the presence of a positive thermal anomaly mostly below the lithosphere, which is also consistent with the absence of volcanism in the seafloor topography lows. In order to explain a mass excess in the crust above the corresponding thermal upwellings, either crustal thinning or a dense intrusion in the lower crust is required. Instead, observations in the vicinity of many oceanic islands [29] show crustal thickening by an underplated lower crustal layer formed as a result of partial melting and crystallization of the hot oceanic upper mantle, resulting in a mass default instead of the needed mass excess.

## 4. Comparison with plate ages and spreading directions

To assess potential lithospheric sources, we first compared the ages of the plates and their former spreading direction with the location, orientation and amplitude of the geodetic signals, as reflected in the negative Bouguer gravity gradients and associated alignments of mantle mass excess. Fig. 4ac shows such a comparison for the Australian-Indian plate in the Indian ocean. The difference between the direction of the gravity/seafloor relief signals and that of the fossil spreading is clear: the Bouguer anomalies follow the recent (0-30 Myr) north-easterly plate motion (Fig. 4a) over areas with different crustal ages (25-110 Myrs), with no clear relation between their amplitude and the crustal age (Fig. 4c). For ages above 40 Myr, the plate spreading direction at the time of formation was northerly, as shown by the Chagos-Maldives-Laccadives and the Ninety East Ridge orientations in Fig. 4c. This indicates that the source of the Bouguer gravity gradients is posterior to the formation of the plate. The same conclusion holds for the Pacific Ocean as shown in Fig. 4bd, where significant differences can be observed between the former plate spreading direction (as indicated by the fracture zones) and direction of the Bouguer gravity gradients. Note that this comparison also rules out intra-lithospheric density variations from hydration near fracture zones, which are oriented in the SW-NE direction. For plates younger than 100 Myr, there is again no obvious systematic link between the plate age and the amplitude of the undulations. For instance, comparable amplitudes of the Bouguer signal, between -400 and locally down to -600 mEötvös, are observed for plate ages from 10 to 100 Myr in the Central Pacific (undulations east of Hawaii, east of the
Marquesas chain and north of the Louisville chain, away from the hotspot). This suggests a relatively limited sensitivity of the observed gradients to lithosphere thickness variations due to conductive plate cooling. This signal is indeed filtered out in our analysis, as it predominates at longer wavelengths and in a different direction. All these comparisons indicate that the source of the observed intermediate-scale undulations must be recent and independent of plate formation processes.

We note however, that the amplitude of the Bouguer gravity gradients over the Pacific plate is larger for plate ages larger than 100-120 Myr, which is also the case for the seafloor slope gradients (e.g. Fig. 1e). This reflects a dichotomy of the signals, where larger seismic velocity anomalies tend to be associated with geodetic signals of smaller amplitude, as observed over the younger side of the Pacific plate, whereas the opposite holds over the older side of this plate.

# 5. Investigation of present-day lithospheric sources

## 5.1. Passive upwellings from lithospheric extension

Seafloor lows above a hot mantle could be associated with passive upwellings due to extensional stresses in the lithosphere, as is observed in rift zones, where asthenospheric material penetrates into a thinned lithosphere [31]. However, instead of the detected mass excess, this thinning would likely lead to mass default at the base of the lithosphere as it is replaced by a hotter asthenosphere, less dense than the stretched lower lithosphere [32, 33] especially for increasing plate age [34]. Furthermore, volcanism would be expected in the seafloor troughs instead of the highs [35]. We finally notice that the periodicity of the undulations is also difficult to explain: the wavelengths of crustal and lithospheric stretching under extensional stresses do not exceed hundreds of kilometers, much less than the 2000 km observed here [36]. Thus, this hypothesis appears unrealistic.

## 5.2. Failure of the model with lithospheric underplating by a dense basaltic layer

As mentioned above, the slow seismic velocity fingers suggest that positive thermal anomalies are present in the mantle and could a-priori be interpreted as hot and light upwellings. In this context, variations in lithosphere-asthenosphere boundary (LAB) topography and density of the lower lithosphere could be expected as a result of decompression melting at the top of these upwellings. As the hot mantle partially melts in the upper mantle low
velocity zone and the melt percolates towards the surface [37], a dense basaltic layer would accumulate at the base of the lithosphere. At the same time, a light depleted residuum would be left behind in the sub-lithospheric mantle. We thus use a two-layer source system to represent a dense basaltic layer accumulated at the base of the lithosphere (eventually forming its lower part), and a layer of light depleted residuum which is left behind in the asthenosphere after partial melting. The layer depths are set to 50-60 km for the underplated layer, and 60-160 km for the depleted residuum. These values are respectively representative of the deeper lithosphere, and of the expected depth range for peridotite partial melting (e.g. [37]). As the mass is conserved when the basaltic constituents are segregated into the melt, each of these two layers carries an opposite mass anomaly, +δM in the basaltic layer, -δM in the residuum. This system of sources is depicted in the inset of Fig. 8.

We apply a dynamic modelling approach and consider two viscosity profiles, with or without a sublithospheric low viscosity zone (LVZ) (see SI Section K). We estimate the expected mass δM of lithospheric underplating in two ways: (1) by fitting the observed Bouguer gravity gradient anomalies and (2) by fitting the observed seafloor slope gradients, assuming they represent dynamic topography. From the mass δM we infer the density anomaly δρR of the depleted residuum. In each case, the two estimates should plot on the diagonal in Fig. 8. Without a sublithospheric LVZ, the results show that the two measurements are incompatible: they differ by a factor of ~24, which furthermore involves unrealistic values of the densities to fit the Bouguer gravity gradient anomalies. With a sublithospheric LVZ, the two estimates are close to the diagonal; however, they correspond to a very large amount of melting. To quantify this, we convert the residuum density anomaly into a percentage of melting and depletion using the reference values in [38]. Because the melt depletion effect on the density depends on depth, we use two reference values: the 60-160 km depth-averaged depletion density anomaly, and the maximum depletion density anomaly at 130 km depth. For 100% depletion, the 60-160 km average anomaly (resp. the maximum 130-km depth anomaly) is Δρ = 117kg.m⁻³ (resp. Δρ = 165kg.m⁻³). With these values, fitting the observed Bouguer gravity gradients requires nearly 50% (resp. 35%) depletion of the whole 60-160 km depth layer, assuming an average (resp. maximum 130-km depth) depletion density anomaly. Such extensive melting is not consistent with seismic velocity anomalies of only 2-3% (Figs. 2, 5). It is associated with significant amounts of underplating, corresponding to the mass of a 10-15 km thick layer with a 250-400 kg.m⁻³ basaltic melt density anomaly [37].
Note that, in our modelling, we did not include the light mass anomalies that correspond to the presumed upwellings associated with the slow seismic velocities in the upper mantle. Further tests show that, when including such thermal anomalies, the value of $\delta\rho_R$ fitting the Bouguer gravity gradients shows a much larger increase than that fitting the seafloor slope gradients, making the two observations incompatible. This is due to the fact that, in a dynamic modeling framework, the masses beneath a LVZ contribute weakly to the surface dynamic topography. We have also considered the case where the melt is partially separated from the depleted residuum and find that the amount of underplated mass increases as compared to the case where the melt is fully separated from the residuum. Finally, we investigated the case where the depleted residuum accumulates in the lower part of a thicker lithosphere, between 60 and 100-km depth. Considering the viscosity profiles of Suppl. Fig. K1, modified to include a 100-km thick lithosphere, we find that the densities $\delta\rho_R$ fitting the Bouguer gravity gradients increase by a factor 4.2 to 4.6 compared to Fig. 8, with and without a sublithospheric LVZ respectively. They remain incompatible with those fitting the dynamic topography by a factor -6 to -10, due to the small amount of dynamic topographic high generated by this system of sources, poorly sensitive to the viscosity profile. Including hot thermal anomalies in the underlying upper mantle further increases these values as well as the incompatibility of the two measurements.

## Discussion: a dynamical origin within the mantle

The above analysis confirms the failure of crustal and lithospheric sources to explain the pattern and amplitude of the gravity and seafloor anomalies, considered jointly, calling for a deeper mantle process. This conclusion is supported by the correlation of the geodetic signals with the seismic structure of the upper mantle extending through the transition zone at the same scale and orientation below all investigated ocean basins, which also suggests a dynamical origin within the convecting mantle below the lithosphere. Also, in the presence of hot and light mantle mass anomalies within the slow seismic fingers, the amplitudes of the crustal and lithospheric contributions needed to explain the detected mass excess would significantly increase as compared to the estimates presented here, making these shallow sources even more unlikely.

Between the smaller scales of the previously known geoid undulations [10, 11, 12, 13, 14] and the long-
wavelength mantle structure, the observed 800–1000-km scale (1600–2000-km wavelength) pattern of gravity and seafloor topography undulations may thus record an intermediate scale of mantle convection below ocean basins. This intermediate scale is particularly clear over the Pacific ocean where it is manifested as a distinct peak in the distribution of scales, between a peak around a scale of 3200-km in the region of the Pacific LLSVP, which may contain the root of broad mantle plumes imaged in recent whole mantle seismic tomography [39, 40], and another peak around ~ 300-km corresponding to the smaller-scale undulations (see Fig. 1c). In addition, the detection of these elongated intermediate-scale signals over both fast-moving plates and the slow-moving Antarctic plate suggests that the corresponding mantle flow is likely not entirely driven by plate motions, as would be the case in the "Richter roll" model [4]. Instead, it probably originates in the deep mantle and may drive the slow motions of the Antarctic plate; it may interact with, or even partially control the motion of a fast-moving plate such as the Pacific plate [41]. Indeed, the Richter roll model would be able to explain the geometry of the observations; however, for a plate velocity of 2.0 cm/yr (larger than the 1.2–1.4 cm/yr velocity of the Antarctic plate) it would take more than 200 Myrs for the convective structures to become aligned as longitudinal rolls [4]. This timescale is much larger than the 10’s Myr over which major changes in the motions of this plate have occurred [42].

We investigated whether purely thermal convection could explain both the magnitudes and the signs of the signals, considering Rayleigh-Benard convection with a ~1:1 cell aspect ratio, a process a-priori able to explain the nearly parallel, 2000-km wavelength periodic alignments in the data and the depth extent of the slow velocity fingers, from the asthenosphere down through the transition zone [8, 39], and, in some areas, down to 1000-km depth (Fig. 5d). We thus considered an alignment of horizontally elongated Rayleigh-Benard cells between the base of an extended transition zone and the base of the sub-lithospheric low velocity zone, where each upwelling or downwelling limb of the convective system is represented by a mass anomaly of width 1000-km extending from 1000-km to 170-km depth, elongated in the APM direction [4]. In a purely thermal interpretation of the seismic velocities, the low velocity fingers would coincide with light, hot upwellings, and the detected mass excess would be entirely due to the corresponding upward deflections of the 660-km discontinuity, attributed to the ringwoodite phase transition which is known to have a negative Clapeyron slope [43]. However, this thermal convection model does not work (SI Section L), as previously argued for the Tonga-to-Hawaii corridor [6]. Even with a large negative 660-km Clapeyron slope of -2.8 10⁶MPa/K, we find that the upward deflections of the 660-km
discontinuity by far do not cancel the negative gravity effect of the light upwellings, which are amplified by the
corresponding downward deflection of the 410-km discontinuity. For example, below the Central Pacific ocean,
the predicted Bouguer gravity gradient signal above such upwellings, assuming a ±100K temperature anomaly in
each upwelling or downwelling limb (corresponding to a 200K lateral variation), exceeds 2400 mEötvös, compared
to the observed ~-400 mEötvös Bouguer anomaly (Suppl. Fig. and Table L1). The 410 and 660-km deflections
respectively reach -8 and +6km, the latter contributing only -400 mEötvös (Suppl. Table L1). Thus, even a
factor 2 enhancement of the Clapeyron slope for the 660-km discontinuity is insufficient (Suppl. Table L1): more
than 45km displacement of this interface is needed in order to explain the observed gravity gradients, which
translates into an increase of the Clapeyron slope by a factor of 8. Such a deflection would carry a negative
buoyancy more than 5 times larger than the positive buoyancy of the upwelling limb of the roll in the lower
mantle and would likely impede any such upward flow. The same holds for the case of deeper and cooler rolls.
Another dynamic issue with the thermal convection model is that, in all cases, it predicts seafloor highs above
the actively ascending flows [44, 6]. Switching these highs into lows as large as observed requires to substantially
reduce the upward vertical velocities when reaching shallower depths, which cannot be achieved even by intro-
ducing a sub-lithospheric channel with unrealistically low viscosities, 10¹⁷ Pa.s over broad areas of the Pacific
(Suppl. Table L1, L2).

Conclusions

From an analysis of gravity gradients and seafloor slope gradients, we have found evidence for intermediate
scale geodetic signals elongated in the direction of APM over broad areas of the Pacific and northern Indian
Oceans. We have shown that the source of these signals must reside in the convecting part of the mantle
and examined their relationship with the shear wave velocity "fingers" present in recent seismic tomographic
images in the upper mantle, down to transition zone depths. This relationship is puzzling, as seafloor lows and
sublithospheric mantle mass excess are found to coincide geographically with the slow seismic velocity fingers and
cannot be explained by a purely thermal convection model, suggesting instead, the presence of lateral variations
in composition in the upper mantle. Reconciling the different sets of observations and their regional variations
in amplitude, and understanding the sources of the mantle density anomalies now calls for a joint dynamical
410 modelling in a thermo-chemical context.

Acknowledgements

This work was supported by CNES. It is based on observations from the GRACE satellite mission as well
as observations from Cryosat-2, Jason-1 and Envisat. BR acknowledges support from NSF grant EAR-1800324.
We thank Bernhard Steinberger and an anonymous reviewer for their constructive reviews, which contributed to
improve our manuscript. The authors used the GMT software [45] for plotting.

415

References

[1] Dziewonski, A.M., Hager, B.H. and O'Connell, R.J. (1977). Large-scale heterogeneities in the lower mantle, *J. Geophys. Res.*, 82(2), 239–255.

420

[2] Richards, M.A., Engebretson, D.C. (1992). Large-scale mantle convection and the history of subduction, *Nature*, 355, 437440.

[3] Morgan, W. J. (1971). Convection plumes in the lower mantle, *Nature*, 230, 4243.

[4] Richter, F. and Parsons, B. (1975). On the interaction of two scales of convection in the mantle, *J. Geophys. Res.*, 80, 25292541.

425

[5] Ballmer, M.D., Ito, G., van Hunen, J. and Tackley, P. (2010). Small-scale sublithospheric convection recon-
ciles geochemistry and geochronology of 'Superplume' volcanism in the western and south Pacific, *Earth and
*Planetary Science Letters*, 290, 224-232.

[6] Katzman R., Zhao L., and Jordan T.H. (1998). High-resolution, two-dimensional vertical tomography of the central Pacific mantle using ScS reverberations and frequency-dependent travel times, *J. Geophys. Res.*, 103 (B8), 17933-17971.

430

[7] Lekic, V. and Romanowicz, B. (2011). Inferring uppermantle structure by full waveform tomography with the spectral element method, *Geophysical Journal International*, 185(2), 799-831.
[8] French, S.W., Lekic, V., Romanowicz, B. (2013). Waveform Tomography Reveals Channeled Flow at the Base of the Oceanic Asthenosphere, *Science*, 342, 227-230.

[9] Hayn, M., I. Panet, M. Diament, M. Holschneider, M. Mandea, and A. Davaille (2012). Wavelet based directional analysis of the gravity field: evidence for large-scale geoid undulations, *Geophysical Journal International*, 189(3), 1430-1456.

[10] Haxby, W.F. and Weissel, J.K. (1986). Evidence for small-scale convection from Seasat altimeter data, *J. Geophys. Res.*, 91, 3507-3520.

[11] Cazenave, A., Houry, S., Lago, B. and Dominh, K. (1992). Geosat-derived geoid anomalies at medium wavelength, *J. Geophys. Res.*, 97, 70817096.

[12] Cazenave A., Parsons B. and Calcagno P. (1995). Geoid lineations of 1000 km wavelength over the Central Pacific, *Geophysical Research Letters*, 22(2), 97–100.

[13] Maia, M. and Diament, M. (1991). An analysis of the altimetric geoid in various wavebands in the Central Pacific Ocean: constraints on the origin of intraplate features, *Tectonophysics*, 190, 133–153.

[14] Wessel, P., Bercovici, D. and Kroenke, L. (1994). The possible reflection of mantle discontinuities in Pacific geoid and bathymetry, *Geophys. Res. Lett.*, 21(18), 19431946.

[15] Hoggard, M.J., White, N. and Al-Attar, D. (2016). Global dynamic topography observations reveal limited influence of large-scale mantle flow, *Nature Geoscience*, 9, 456-463.

[16] Panet, I. (2018). An analysis of gravitational gradients in rotated frames and their relation to oriented mass sources, *Journal of Geophysical Research*, 123, 11062–11090.

[17] McKenzie, D., Yi, W. and Rummel, R. (2014). Estimates of $T_e$ from GOCE data, *Earth and Planetary Science Letters*, 399, 116-127.

[18] Panet, I., Pajot-Métivier, G., Greff-Lefftz, M., Métivier, L., Diament, M. and Mandea, M. (2014). Mapping the mass distribution of Earth’s mantle using satellite-derived gravity gradients, *Nature Geoscience*, 7, 131135.

[19] Greff-Lefftz, M., Métivier, L., Panet, I., Caron, L., Pajot-Métivier, G. and Bouman, J. (2016). Joint analysis
of GOCE gravity gradients data, of gravitational potential and of gravity with seismological and geodynamic
observations to infer mantle properties, *Geophys. J. Int.*, 205(1), 257-283.

[20] Ebbing, J., Haas, P., Ferracioli, F., Pappa, F., Szwillus, W. and Boumann, J. (2018). Earth tectonics as seen by GOCE - Enhanced satellite gravity gradient imaging, *Scientific Reports*, 8:16356.

460

[21] Afonso, J.C., Salajegheh, F., Szwillus, W., Ebbing, J. and Gaina, C. (2019). A global reference model of the lithosphere and upper mantle from joint inversion and analysis of multiple data sets, *Geophysical Journal International*, 217, 16021628.

[22] Altamimi, Z., Métivier, L. and Collilieux, X. (2012). ITRF2008 plate motion model, *J. Geophys. Res.*, 117, B07402.

465

[23] Clouard, V. and Gerbault, M. (2008). Reply to Break-up spots: Could the Pacific open as a consequence of plate kinematics? Comment by R. Pilger, *Earth and Planetary Science Letters*, 275, 196199.

[24] Szwillus, W., Afonso, J.C., Ebbing, J. and Mooney, W. (2019). Global crustal thickness and velocity structure from geostatistical analysis of seismic data, *J. Geophys. Res.*, 124, 1626-1652.

[25] Haxby, W.F. and Turcotte, D.L. (1978). On Isostatic Geoid Anomalies, *J. Geophys. Res.*, 83, B11, 5473-5478.

470

[26] Hoggard, M. J., Winterbourne, J., Czarnota, K. and N. White (2017). Oceanic residual depth measurements, the plate cooling model, and global dynamic topography, *J. Geophys. Res. Solid Earth*, 122, 2328-2372.

[27] Laske, G., Masters, T. G., Ma, Z. and Pasyanos, M. E. (2013). Update on CRUST1.0. A 1-degree global model of Earths crust, *Geophysical Research Abstracts*, 15.

[28] Pasyanos, M. E., Masters, T. G., Laske, G. and Ma, Z. (2014). LITHO1.0: An updated crust and lithospheric model of the Earth, *J. Geophys. Res.*, 119, 21532173.

475

[29] Leahy, G.M. and Park, J. (2005). Hunting for oceanic island Moho, *Geophys. J. Int.*, 160, 1020-1026.

[30] Laske, G. and G. Masters (1997). A global digital map of sediment thickness, *EOS Trans. AGU*, 78(46), Fall Meet. Suppl. F48397.

[31] Sengör, A.M. and Burke, K. (1978). Relative timing of rifting and volcanism on Earth and its tectonic implications, *Geophysical Research Letters*, 5(6), 419-421.

480
[32] McKenzie, D. (1978). Some remarks on the development of sedimentary basins, *Earth and Planetary Science Letters*, 40, 25-32.

[33] Sclater, J. and Celerier, B (1987). Extensional models for the formation of sedimentary basins and continental margins, *Norsk Geologisk Tidsskrift*, 67, 253-267.

485

[34] Oxburgh, E.R. and Parmentier, E.M. (1977). Compositional and density stratification in oceanic lithosphere - causes and consequences, *J. geol. Soc. Lond.*, 133, 343-355.

[35] Sandwell, D. T., Winterer, E.L., Mammerickx, J., Duncan, R.A., Lynch, M. A., Levitt, D. A. and Johnson, C. L. (1995). Evidence for diffuse extension of the Pacific plate from Pukapuka ridges and cross-grain gravity lineations, *J. Geophys. Res.*, 100, 15087-15099.

490

[36] Ricard, Y. and Froidevaux, C. (1986). Stretching instabilities and lithospheric boudinage, *J. Geophys. Res.*, 91, B8, 8314-8324.

[37] Sakamaki T., et al. (2013). Ponded melt at the boundary between the lithosphere and asthenosphere, *Nature Geoscience*, 6, 1041–1044.

[38] Schutt, D.L. and Lesher, C.E. (2006). Effects of melt depletion on the density and seismic velocity of garnet and spinel lherzolite, *Journal of Geophysical Research*, 111, B05401.

495

[39] French, S.W. and Romanowicz, B. (2015). Broad plumes rooted at the base of the Earth's mantle beneath major hotspots, *Nature*, 525, 95–99.

[40] Lei W., Ruan Y., Bozda E., Peter D., Lefebvre M., Komatitsch D., Tromp J., Hill J., Podhorszki N., Pugmire D. (2020). Global adjoint tomography-model GLAD-M25, *Geophys. J. Int.*, 223, 121.

500

[41] Coltice, N., Gérault, M. and Ulvrova, M. (2017). A mantle convection perspective on global tectonics, HYPERLINK "https://www.sciencedirect.com/science/journal/00128252", *Earth Science Reviews*, 165, 120-150.

[42] Seton, M. et al. (2012). Global continental and ocean basin reconstructions since 200 Ma, *Earth-Science Reviews*, 113, 212–270.

505
[43] Stacey, F.D. and Davis, P.M. (2008). Physics of the Earth, Fourth Edition, Cambridge University Press, Cambridge, UK.

[44] Richards, M.A. and Hager, B.H. (1984). Geoid anomalies in a dynamic Earth, *J. Geophys. Res.*, 89, 5987-6002.

510

[45] Wessel, P. and Smith, W.H.F. (1995). New version of the generic mapping tool, *EOS Trans. Am. Geophys. Un.*, 76(33), 329.

[46] Bird, P. (2003). An updated digital model of plate boundaries, *Geochem. Geophys. Geosyst.*, 4, 1027–1078.

[47] Steinberger, B. (2000). Plumes in a convecting mantle: Models and observations for individual hotspots, *J. Geophys. Res.*, 105, 11127–11152.

515

[48] Conrad, C.P. and Lithgow-Bertelloni, C. (2006). Influence of continental roots and asthenosphere on plate-mantle coupling, *Geophys. Res. Lett.*, 33, L05312.

[49] Müller, R. D., Sdrolias, M., Gaina, C. and Roest, W.R. (2008). Age, spreading rates, and spreading asymmetry of the world's ocean crust, *Geochemistry, Geophysics, Geosystems*, 9(4), Q04006.
## Figures

**Figure 1: Medium-scale periodic patterns in the gravity field and the seafloor relief for the Pacific plate.** Bottom: scale-orientation diagrams for the Pacific plate, showing the number of occurrences of the characteristic orientations and scales in the GRACE gravity gradients (a) and the seafloor slope gradients (b), in the area delimited by the black box in the panel (c) map. The azimuth clockwise from north is given in abscissa. The numbers in panels (a) and (b) indicate: orientations of elongated signals sub-parallel to the directions of fracture zones (1), present-day APM (2) (also marked by black ellipses) and end of Cretaceous APM (3) [42]. Black vertical bar: direction of the present-day APM. Panel c: Distribution of the characteristic scales of undulations elongated in the direction of azimuths $N110°-140°$ in the GRACE gravity (blue) and the seafloor slope (pink) gradients, sampled within transverse profiles across the undulations shown in the inset map. Top: GRACE gravity gradients (d) and seafloor slope gradients (e) around the Pacific APM direction (1100km scale, $θθ$ gradients, 20-50° frame rotation). Black dashed lines locate the SEMum2 [8] slow velocity fingers on all maps. Orange lines: plate boundaries [46]; black arrows: present-day APM direction [42]; green dots with red circles: hotspots [47].

**Figure 2: Medium-scale periodic patterns in the upper mantle shear velocities, seafloor slope gradients, and GRACE gravity gradients for the Australian-Indian and Antarctic plates.** (a) and (f) same as Fig. 1c for the Indian (a) and Antarctic (f) plates (azimuths $N30°-60°E$ for both plates). (b,c,g,h) same as Fig. 1de for the Indian (b-c) and Antarctic (g-h) plates: (b) and (g) GRACE gravity gradients; (c) and (h) seafloor slope gradients around the respective APM directions (Indian ocean: 1100km scale $φφ$ gradients, 30-60° frame rotation; Antarctic plate: 800km scale $φφ$ gradients, 30-60° frame rotation). (d) and (i) Same as Fig. 5b for the Indian and Antarctic plates, respectively. Note that black dashed lines indicate the location of SEMum2 slow velocity fingers on all maps. For each region, variations along the profiles marked by thick black lines in the maps (d,i) are shown in panels (e,j) for all observables, as in Fig. 5. For the Antarctic plate, the Bouguer gradients are not shown because they are strongly affected by the Antarctic continent boundaries, and a short additional cross-section is shown for SEMum2, along the green dashed line in the map (i).

**Figure 3: Seafloor slopes gradients and gravity gradients associated with an equilibrated system**
of two sources of opposite sign: a seafloor topographic load and its internal support. Left: 600-km scale, $\theta\theta$ gravity gradients of a 600-km wide, 2400-km long, 10 kg.m$^{-3}$ rectangular mantle mass excess source located between 110 and 210-km depth within a rigid Earth. In all panels, the spherical frame is rotated around the radial axis such that the $e_\theta$ axis is orthogonal to the source major axis. Middle: 600-km scale, $\theta\theta$ seafloor slope gradients (orange curve) and gravity gradients (blue curve) of a 530m high, 600-km wide and 2400-km long Gaussian-shaped seafloor depression following the same azimuth. Seafloor slopes (resp. gravity) gradients maxima reflect seafloor topography (resp. geoid) lows. Right: total 600-km scale, $\theta\theta$ gravity gradients of the system of these two sources carrying mass anomalies of opposite signs, the seafloor low and the internal mass excess.

**Figure 4: Bouguer gravity gradient signals compared with plate age and spreading direction.** Left panels (a,b): Maps of 1100-km scale Bouguer gravity gradients derived from GRACE and seafloor data over the Pacific and Australian-Indian plates along the respective APM directions (Top left, Indian Ocean: $\phi\phi$ gradients, 30-60° clockwise frame rotation; Bottom left, Pacific plate: $\theta\theta$ gradients, 20-50° clockwise frame rotation). Note the flipped colorbar of the Bouguer gravity gradients. Purple (resp. green) arrows: present-day (resp. 30 Myr old, for the Indian ocean) APM direction [42]; green dots with red circles: hotspots [47]. Right panels (c,d): map showing the ages of the oceanic plates [49]. Black lines show the locations of the negative Bouguer gradient anomalies; the age of the plates at the beginning and the end of each line is indicated. The orange lines in the bottom panel indicate the fracture zones of the Pacific plate: Mendocino (a), Murray (b), Molokai (c), Clarion (d), Clipperton (e), Galapagos (f), Marquesas (g).

**Figure 5: Correlation of the periodic gravity field and seafloor relief signals with SEMum2 (8) and SEMUCB-WM1 (15) slow velocity anomalies in the upper mantle and transition zone.** (a): SEMum2 [8] slow velocity anomalies at 260km depth (in yellow to red colors), superimposed with the positive GRACE gravity gradients (dark blue contours) and seafloor slope gradients (pink contours) signals in three areas (East of Hawaii, Polynesia, Louisville) delimited by the black boxes in panel (b). The contours of the negative Bouguer gravity gradients (light blue) are also included for Polynesia. In this Figure, all gradients are built using a 10-40° clockwise frame rotation, well-suited for the regions close to the East Pacific Rise. (b) SEMum2 slow
velocity anomalies as in (a) across the whole Pacific Basin with location of areas shown in (a) and of cross section shown in (c) and (d) (thick black line with white and purple dots). Green dots with red circles are hotspot locations. (c): Variations along the profile indicated in panel b, for all 4 observables (shear velocity anomalies at 260 km depth). Note that the signs of the shear velocity anomalies and Bouguer gradients have been reversed. Green rectangles indicate regions of lower than average shear velocity and pink rectangles indicate regions of seafloor lows. (d): depth cross-section in SEMUCB-WM1 seismic model, showing relative shear velocity anomalies with respect to the global mean, along the same profile.

**Figure 6: Spatial structure of an isostatic crust model for the Central Pacific.** Top: scale-orientation diagrams for the Central Pacific, showing the number of occurrences of the characteristic orientations and scales in the gravity gradients of an isostatic crust model including (a) or not (b) an age-dependent lithospheric layer [48]. The considered area is the same as in the diagrams of Fig. 1ab. The black vertical bar indicates the direction of the present-day APM and the black ellipses show the orientation of the undulated signals sub-parallel to the APM in the geodetic data, as in Fig. 1ab. Bottom: maps of the 1100-km scale, $\theta\theta$ gravity gradients of these isostatic crustal signals with (c) or without (d) the age-dependent lithospheric layer [48] in the Central Pacific; 20–50° stacked frame rotation. Compared with the observed ones (Fig. 1d), these predicted amplitudes are 4 times smaller.

**Figure 7: Spatial structure of a seismological Moho model for the Pacific and the Northern Indian ocean.** Scale-orientation diagrams for the Central and South Central Pacific (panels a, b) and the Northern Indian Ocean (panels c, d), for the second-order gradients of the Moho depth from the LITHO1.0 seismological model [28]. The considered areas are indicated on the insert maps. The APM orientation is indicated on each diagram by the solid blue line. The blue dashed line marks the Indian plate motion orientation 30 Myrs ago. Black ellipses indicate the observed GRACE and seafloor signals, not found in this Moho model.

**Figure 8: Density anomaly of the light depleted residuum in a model of lithospheric underplating.** Density anomaly $\delta\rho_R$ of the light depleted residuum calculated for two viscosity profiles: with a low viscosity channel LVZ (red circle) and without a LVZ (blue circle). The densities fitting the Bouguer gravity gradients
605

(resp. the surface dynamic topography gradients) are shown in the vertical axis (resp. the horizontal axis). The inset shows a zoom comparing the obtained values with those associated with 50 and 100 % depletion in the 60–160 km layer [38]. Lower bounds: average 60-160km depth value ; upper bounds: 130 km depth maximum value. The two-layer sources system is shown in the right panel, where $\delta M$ indicates the total underplated / depleted mass.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 1</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 2</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 3</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 4</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 5</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 6</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 7</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 8</figcaption></figure>