Continental Hotspots Tracks From an Analysis of GOCE Gravity Gradients Data

This study aims to detect residual thermal anomalies along the path of a hotspot beneath the continental lithosphere using gradiometry data. We calculate perturbations in gravity vectors and gradients by employing a simple thermal anomaly model. The results highlight that the horizontal components of Bouguer gravity gradients, sensitive to source structure and elongation, exhibit substantial magnitudes, exceeding current detection levels. To improve signals from sources of different sizes and orientations, we apply a filtering approach involving wavelet transforms and rotated local frames. We also assess the impact of lateral crustal variations on the gravity field by introducing random density anomalies and a seismic crustal model, testing the method's ability to distinguish between these anomalies and thermal sources. Using Gravity field and steady‐state Ocean Circulation Explorer satellite data, we generated scale‐orientation diagrams aiming to identify signals aligned with the direction of plate movement and with a spatial scale of a few hundred kilometers (scale of the plumes). Maps of filtered Bouguer gravity gradients aligned with hotspot trajectories are generated for all continental hot plumes, with the position and the age of the intraplate volcanism to support the interpretation. While we do not show a clear signal in regions that are too tectonically complex or with lateral variations in the crust or at the lithospheric boundary correlated with the hotspot trace, we found observed signals with expected sign and amplitude over Hoggar, Tibesti, Darfour and Cameroon tracks in Africa and over Iceland and Jan Mayen in Greenland. However, determining corresponding mass source depths remains a complex task.


Introduction
Hotspots are regions of thermal instability that arise in the mantle and lead to volcanic activity on the Earth's surface (Courtillot et al., 2003;J. P. Morgan, 1972).They are often associated with mantle plumes, which are column-like structures of hot rock ascending from the deep mantle.As the plume approaches the base of the lithosphere, the pressure and temperature cause the rock to melt and form magma.This magma then rises through the lithosphere, eventually erupting as volcanism on the surface (for a review, see Ballmer et al., 2015).
The stability of such plumes within the mantle is thought to depend on a variety of factors, including the composition of the plume material, the temperature gradient between the plume and surrounding mantle, and the dynamics of the convecting mantle itself (for an extensive review about mantle plumes, see Koppers et al., 2021).Recent research has suggested that plumes may be more stable than previously thought.Indeed some models suggest that they can remain stable over long periods of time and move primarily in a vertical direction within the mantle, with weak horizontal movement due to the strong convection currents (for a review see Arnould, 2018).This hypothesis has important implications for tectonic and geodynamic processes on Earth.As the lithospheric plates of the Earth move horizontally with respect to the mantle (with a horizontal velocity of a few centimeters per year), a hot plume that is fixed in the mantle will move with respect to the plate.This creates a chain of volcanoes that become progressively older and lower in elevation as we move away from the current position of the hotspot (J.P. Morgan, 1972;Richards et al., 1989).These features can be observed in the oceanic lithosphere, as alignments of islands and seamounts.The studies of hotspot tracks in the oceanic lithosphere have provided valuable information on the movement of tectonic plates and the characteristics of the mantle beneath them.They allowed researchers to construct reference frames linked to the mantle, in which the plumes would be fixed (Müller et al., 1993) or slightly variable (Doubrovine et al., 2012;O'Neill et al., 2005).
The continental lithosphere is thicker and more complex than the oceanic lithosphere, making hotspot traces more difficult to observe and interpret over continents.Previous research focused on the dynamics of mantle plumes and their interaction with the continental lithosphere has shown that mantle plumes have a significant impact on the lithosphere and vice versa (Agrusta et al., 2013;Moore et al., 1999;Shi et al., 2021).As the plume head reaches the base of the lithosphere, it spreads out laterally and becomes trapped beneath the lithosphere.This can lead to the formation of large, buoyant magma bodies that can cause uplift and deformation of the overlying lithosphere.Observed lithosphere uplifting together with mathematical modeling and geophysical data analyses could give us insight into the plume's activity (Burov & Gerya, 2014;Wang & Li, 2021).The interaction between the plume and the lithosphere can also change the chemical and thermal properties of the lithosphere, leading to the formation of different types of rocks and minerals.These rocks have distinctive chemical and isotopic signatures that can be used to trace their origin and the role of plumes in their formation (Allègre et al., 1987;Hart et al., 1992;Jackson et al., 2017).However, while the presence of the hotspot during its period of activity has left traces within the continental lithosphere that have been studied in the literature, the traces left by the passage of the plume under the lithosphere are much less well-known.One key feature of these tracks is the presence of volcanic and igneous rock formations different from the surrounding rocks.In addition, residual thermal anomalies (RTAs)-areas where the lithosphere has been heated by a hotspot in the past-can persist for millions of years whereas, since the plate has moved, the plume is no longer active there.Such a large amounts of heat can have a range of effects on the lithosphere.Several studies have been conducted to understand these features.One such study concerned the Yellowstone-Snake River Plain volcanic provinces (J.P. Morgan, 1972;L. A. Morgan et al., 2009) with the 10 Ma track of the hotspot (Pierce & Morgan, 2009) and the thermomechanical behavior linked to the plume tail imaged by Shervais and Hanan (2008).Another study by D. R. Davies et al. (2015) highlighted a long trace of 2,000 km in the eastern part of Australia that provides a record of volcanic activity between 33 and 9 million years ago.Based on a seismic study in Virginia (Chu et al., 2013), proposed that thermo-mechanical erosion of the lithosphere induced by a "hidden" hot plume that did not create surface volcanism could explain a slow anomaly of seismic velocity.These interactions between a fixed hot plume conduit and a moving depleted lithosphere have been modeled, taking into account both diffusion of heat and thermal erosion, using either a 1D model (G.F. Davies, 1994) or a 3D model (Yang & Leng, 2014).The plume can erode the lower part of the moving lithosphere, creating a corridor-like RTA associated with a low viscosity and low seismic velocities downstream in the direction of plate motion as well as surface heat flux anomalies (Heyn & Conrad, 2022).This thermal anomaly, which is approximately 300 km wide and 50 km height (D. R. Davies et al., 2015;Yang & Leng, 2014), could persist under the lithosphere for at least 100 million years and extend for thousands of kilometers along the past direction of plate motion.
Here, we aim to identify such RTA along hotspot trajectories using gradiometry data, which provide detailed information about the density and distribution of mass within the Earth's interior.From 2009 to 2013, the European Space Agency's Gravity field and steady-state Ocean Circulation Explorer (GOCE) satellite measured the six components of the gradiometric tensor, which represents the spatial derivatives of the Earth's gravity vector in three different directions.These data are highly sensitive to the directional structure of the gravity field and its sources and can highlight elongated geometries aligned with the orientations of the plate motions.We will use these data to create maps of gravity gradients at the surface.These maps will be filtered at a spatial scale of a few hundred kilometers, which is the scale of the hotspot, and oriented along the direction of the hotspot's trajectory.By analyzing these maps, we will search for traces of thermal anomalies associated with all continental hotspots, to gain a better understanding of their behavior and movement, with a particular focus on five regions: North Africa, Australia, Greenland, Europe and Virginia.

Modeling
The thermal anomaly resulting from the interaction between a hotspot conduit and a moving lithosphere was first studied by G. F. Davies (1994), considering both diffusion of heat and thermal erosion.The magnitude of the RTA depends on the size of the plume conduit, on the advection within the conduit, and on the plate velocity.A lower plate velocity results in a longer duration of the lithosphere heating by the hot plume, and a larger density anomaly.Accounting for the negative contribution (light) of the hot thermal anomaly and the positive contribution (heavy) of lithospheric thinning yields a density anomaly of ∼ 10 kg/m 3 for a lithospheric thinning of ∼50 km (D. R. Davies et al., 2015;Yang & Leng, 2014).
We modeled the RTA as a parallelepiped aligned along the equator, centered at [0°N, 210°E], with a height of 50 km and a bottom depth of 100 km.We assume a Gaussian density anomaly as defined in Wahr et al. (1998): the radius of the Gaussian (fixed to 280 km in latitude) is the distance at the Earth's surface for which the value of the Gaussian is equal to half the maximum value at the center ( 10 kg/m 3 here) as shown in Figure 1a.
To compute the associated perturbation of the gravitational potential and of the topography, we use classical geoid and topography kernels (Richards & Hager, 1984) assuming a radially stratified compressible Earth's model.The density is given by the PREM model (Dziewonski & Anderson, 1981) and we assume a viscous Newtonian Earth, with a simple viscosity profile.The viscosity is set to 2 × 10 22 Pa.s, 10 21 Pa.s and 4 × 10 22 Pa.s for the 100 km thick lithosphere, the underlying upper mantle and the lower mantle, respectively.The Moho is treated as a physical interface on which we assume that the momentum is continuous.The source is expanded in spherical harmonics up to the degree and order 200.Figures 1b and 1c show the associated geoid and topography at the surface.Above the source, the magnitudes of the signals are approximatively one m for the geoid and one hundred meters for the topography.Note that for a negative mass anomaly at the bottom of the lithosphere, we obtain a positive geoid because the positive contribution of the surface topography is larger than the negative direct effect of the sublithospheric mass source on the gravitational potential.The three components of the gravity vector are the first order spatial derivatives of the gravitational potential, whereas the nine components of the gravity gradient tensor are its second order spatial derivatives.Here, we compute the gravity gradient tensor in the local spherical basis (Casotto & Fantino, 2009).We note e → X the spherical vector along the meridian and directed toward the South pole, e → Y the Eastward spherical vector along the parallel, and e → Z the outward radial vector.We plot in Figure 1d the six components of this symmetric tensor.The order of magnitudes of the signals varies between 1 Eötvös (10 9 s 2 ) for the ZZ component and 1 Eötvös for XX.
This example illustrates the ability of gradients to highlight the geometry of the source.In particular, the component XX which corresponds to the second order derivative of the gravitational potential in the direction perpendicular to the equator clearly highlights the boundaries of the structure.As written above, for sources close to the surface, the signal associated with the mass redistribution induced by the viscous flow in the mantle is primarily created by surface topography.To isolate the sublithospheric source we are looking for, we will remove the topographic contribution from the surface gravity (Bouguer anomalies) and calculate Bouguer gravity gradients (Figure 1e) The amplitude of the Bouguer gravity gradient components is lower and of opposite sign than that of the gravity gradient components: 500 mEötvös for the vertical ZZ component and +500 mEötvös for the horizontal XX component.
In this simple calculation, we assumed that the lithosphere was a viscous layer.However, the thermal anomaly under investigation has a spatial scale smaller than the flexural wavelength of the lithosphere (approximately 400 km).Consequently, it would be necessary to consider the elastic strength of the lithosphere.To estimate this effect, we replicated the same calculation as presented above, but assuming that the lithosphere is an elastic layer.This results in a significantly reduced dynamic topography, to the extent that the geoid above the source becomes negative (this is the direct effect of masses beneath the lithosphere dominating the signal).Nevertheless, the Bouguer anomaly (sensitive to the source beneath the lithosphere) remains nearly identical, and thus, the order of magnitude of the associated Bouguer gravity gradient components will be around ±500 mEötvös.
In this specific scenario, the surface topography is the dynamic topography generated by the source at the base of the lithosphere and consequently, the interpretation of the Bouguer anomaly is easy.But in practice, the dynamic topography cannot be accurately measured: it cannot be easily separated from the topography associated with crustal thickness variations.In the real Earth, topography independent from the hotspot at the origin of the studied RTA may indeed also be present.This could be a topography from tectonic origin existing for much longer than the hot spot.Such a topography would be partially or isostatically compensated by a deformation of the Moho and thus have a crustal root.The two density anomalies (surface tectonic topography and its crustal root) create gravity gradients which differ with respect to the ones created by the RTA, making it possible to separate the sources (tectonic vs. hotspot).However, in some cases, the passage of the plume induces crustal density anomalies and volcanic chains that are too thin to be isostatically compensated and with a geometry very similar to that of the RTA, aligned in the direction of plate movement.Consequently, when a volcanic trace is associated with the plume, it may be challenging to distinguish the mass anomaly originating from the crust from that created by the RTA in the Bouguer anomaly maps.
To illustrate this, we present a synthetic example with three parallelepipeds, at different depths, of length L = 2,500 km aligned in the same direction (Figure 2a): • C1 represents the surface topography (height h = 400 m, radius R = 250 km, density ρ c = 2,700 kg/m 3 ) • C2 represents the crustal root (Radius R = 250 km, height w = ρ c ρ L ρ c h corresponding to isostatic equilibrium, located at the depth d c + w/2, with d c = 35 km the thickness of the crust and ρ L = 3,300 kg/m 3 the density of the lithosphere).The mass of C2 is the opposite of the mass of C1.
• C3 represents the RTA at the base of the lithosphere (see Figure 1) We plotted in Figure 2b the different topographies: the isostatic topography (C 1 ), the dynamic topography induced by the source C 3 , and the sum of these two contributions.The calculated horizontal gradient XX along a profile perpendicular to the direction of the considered sources is plotted in Figure 2c.We clearly see that it will be difficult to extract the signal created by the source C3 from the Bouguer anomaly (C2 + C3) when we investigate the signals along the path of a continental plume.Nevertheless, the amplitude of the expected signal is well above the data detectability threshold.
Another approach to address the separation of crustal mass anomalies from the sublithospheric RTA signal is to utilize independent data to correct the observed signal for lateral variations in crustal density, particularly for crustal roots.We calculated the gravity gradient predicted by models of crustal mass anomalies derived from seismology (models Crust1.0 (Laske et al., 2012) and LITHO1.0(Pasyanos et al., 2014)).These models do not assume mass balance (isostatic or not), while the gravimetric signal is highly sensitive to this assumption.This is why these models have produced large amplitude gravity anomalies that are challenging to interpret.This is a complex problem that is beyond the scope of this article (see for example Tenzer et al., 2015).To push the limits of exploration in the absence of a sufficiently realistic crust model, while still making the most of the information embedded in seismic crustal models, we will scrutinize the spatial scales and preferential orientations of their crustal density anomalies.This involves performing a multi-scale directional analysis of the gradients.Furthermore, random crustal density anomalies will be introduced to evaluate the signal-to-noise ratio.All these endeavors aim at fostering a discussion on the separation of the sources, be they RTA or crustal anomalies, impacting Bouguer gravity gradients.

Gravity Gradients in Rotated Frames, at Different Spatial Scales With Random Crustal Density Anomalies
To investigate the robustness of the detection of a hotspot related source in the presence of crustal noise, we added a random density anomaly in the upper 50 km of the Earth with an amplitude ranging from 0 to 50 kg.m 3 , sampled at intervals of one degree in latitude and longitude and 10 km in radius using the intrinsic fortran function random number: in each spherical volume element 1°× 1°× 10 km, we randomly generate a value between 0 and 1, which we multiply by the desired maximum density amplitude (10, 25 or 50 kg.m 3 ).To identify and separate these sources, we applied a filter at different spatial scales to the gravity potential, from which the gravity gradients are derived, using a wavelet transform (Holschneider et al., 2003) as described in Panet (2018aPanet ( , 2018b)).The characteristic scale of the thermal anomaly we are investigating is a few hundred kilometers.To demonstrate the noise-to-signal ratio as a function of filter width, we plot in Figure 3 the XX Bouguer gravity gradient component induced by the parallelepiped mass source C3 oriented along the azimuth 90°N superimposed with the crustal random density anomaly at amplitude of 10 kg/m 3 (left column), 25 kg/m 3 (middle column), and 50 kg/m 3 (right column), for different spatial scales.We can clearly see that the wavelet filtering allows us to separate the signal created by the plume on a larger spatial scale than the crustal density anomalies.However, when the amplitude of the crustal density is around 50 kg.m 3 , the signal becomes very noisy and difficult to interpret.To further separate these signals of different origins, we will also take into account the fact that they may have different orientations.In Figures 1 and 2, we underlined the ability of the second order derivative of the gravitational potential in the direction perpendicular to the parallelepiped to highlight the boundary of the structure.To In this rotated local frame, we compute the six components of the gravity gradient tensor (for more details on this method, see Panet, 2018aPanet, , 2018b)).As an example, we start with the same source as in Figure 1 but now the parallelepiped is aligned in a direction making an angle of 30°with respect to the North pole.We compute the horizontal gradient signal Y′Y′ (second order derivative of the potential in the direction e → Y′ ) for different clockwise rotations of the spherical frame by an angle A z ranging from 0°to 180°(Figure 4).We note that the horizontal gradient signal is strongly distorted as we rotate the frame away from the orientation of the source and that, for A z = 30°(i.e., along a direction perpendicular to the parallelepiped), it has the maximum amplitude and reflects the source orientation.
To determine if the hotspot signal can be detected by this analysis, we add a random crustal density anomaly to the parallelepiped RTA (as in Figure 3) and we examine different orientations and spatial scales.To present our results in a synthetic figure, we create a scale-orientation diagram (Panet et al., 2022) using the following steps: 1.At each point (colatitude θ, and longitude φ) in the studied area, we rotate the frame from 0°to 180°and vary the spatial scale from 150 to 500 km.Then, we record the maximal amplitude of the YY component (noted Geochemistry, Geophysics, Geosystems 10.1029/2023GC011200 YY Max ), along with the associated rotation angle α Max and spatial scale L Max .We obtain three functions: YY Max (θ, φ), α Max (θ, φ) and L Max (θ, φ). 2. We count the number of occurences of each rotation angle and spatial scale in α Max (θ, φ) and in L Max (θ, φ) over the studied region.If we denote the number of points in latitude as N lat and the number of points in longitude as N long , we can express this count as: where δ is the Dirac function.3. To compute the scale-orientation diagram Diag(α, L), we normalize this function by: We present the scale-orientation diagrams of the YY Bouguer gravity gradient component for the parallelepiped mass source C3 oriented along an azimuth of 90°N, with and without random crustal density anomalies of amplitudes 10 or 25 kg/m 3 (Figure 5).The dominant signal is observed at A z = 90°and L = 350 ∼ 400 km, for the cases without or with weak crustal anomalies.However, in the presence of stronger random crustal anomalies, the observed signal is dominated by the small spatial scales, and the large scale signal is too noisy to be detected.
In our search for a signal from the RTA along the trajectory of a hotspot, we will compute the scale-orientation diagram in each investigated region to estimate the direction and spatial scale at which the signal is maximum.We will investigate if there is a signal in the present-day Absolute Plate Motion direction (APM) on a spatial scale of a few hundred kilometers.We will place ourselves in a reference frame { e → X′ , e → Y′ , e → Z } where the vector e → X′ is aligned with the trajectory and pointing toward the present direction of the plate motion (APM), and e → Y′ is perpendicular as expected from a 90°counter-clockwise rotation of e → X′ .We will then calculate the two key components: the vertical component, noted ZZ, and the component perpendicular to the trajectory, noted Y′Y′ (denoted as YY in the figures for simplify), and we will search for a negative (blue) signal in ZZ and a positive (red) signal in YY component of the Bouguer gravity gradient along the path of the plume.
In each investigated region, to improve our interpretation, we will also plot the scale-orientation diagram associated with the gravitational potentials induced by the topography and by the crustal density anomalies.

Data
We use a global gravity field model based on GOCE satellite measurements and expanded in spherical harmonics up to degree/order 260 (Bruinsma et al., 2013) to compute the gravity gradients at scales greater than 200 km.In addition, we use the EGM2008 global gravity field model (Pavlis et al., 2008(Pavlis et al., , 2012)), based on satellite gravity and altimetry measurements and also including surface data sets, for scales smaller than 200 km.This latter gravity model is expanded in spherical harmonics up to degree/order 2190 and will be used only in regions densely covered in ground measurements.
A topography model is necessary for the calculation of the topography gravity gradients and of the Bouguer gravity gradients: we chose the DV-ELL-RET2012 model of the gravity field created by the Earth's topographic load (Claessens & Hirt, 2013).In this model, the density of the continental crust used for the topographic corrections is fixed at 2.67 and does not take into account any lateral variations (except from continental ice and water).However, in certain regions, the presence of sediments can significantly disturb the gravity signal.To solve this problem, we extracted from the model Crust1.0(Laske et al., 2012) the three layers of sediment located above altitude 0 and their density from which we subtracted the average value of 2.67 (used in RET2012).Then, by splitting the topography into 10 m thick radial slices and summing the contribution of each layer, we calculated the associated loading potential expanded into spherical harmonics as well as the associated gravity gradients.At the 250 km scale, the gravity gradient signal of this sediment correction can reach 300 mEötvös, therefore is quite significant in the interpretation of our maps.Finally, we have summed the coefficients of the spherical harmonics decomposition of this correction with those of RET2012 (see Text S1 in Supporting Information S1) to build our map of topography gravity gradients and of Bouguer gravity gradients.

Hotspot Trajectories and Intraplate Volcanism on Continental Lithosphere
To compute the trajectory of a hotspot, we need the Euler pole and angle of the plate in the geological past with respect to a mantle reference frame in which the hotspots are fixed or slightly moving.We used three mantle reference frames built thanks to the volcanoes alignments left on the oceanic lithosphere by the plumes: the first one assuming that hotspots are fixed in the mantle (Müller et al., 1993), the second one taking into account the slight advection of hotspots with a mantle circulation model (O'Neill et al., 2005), and the last one taking into account the advection of the hotspots and using both Indo-Atlantic and Pacific volcanic tracks with a refined plate circuit model linking the Pacific plate with the Indian Ocean bordering continent (Doubrovine et al., 2012).These three reference frames give three trajectories which are very close for recent periods (up to 40 Ma) but deviate significantly for older periods.
Same as for the oceanic lithosphere, the plumes may sometimes leave volcanic traces on the surface.We use a compilation of global surface volcanism (Greff-Lefftz & Besse, 2012) in which we extracted from three databases (the Georoc, the NGDB [see References] and the GPMDB v.4.6 database (McElhinny & Lock, 1996)) the locations and ages of volcanic edifices.In all our figures, intraplate volcanism will be represented by triangles of different colors depending on their age.

Hawaii Oceanic Hotspot
In this section, we validate our method using a major oceanic hotspot, Hawaii, as a proof of concept before delving into more intricate scenarios involving continental lithosphere.We computed the Bouguer scale-orientation diagram over a region between 20°-30°N and 178-208°E (Figure 6a) and found a signal around the APM direction ∼115°at the spatial-scale of 350 km (Figure 6b).We plotted in Figures 6c and 6d the vertical and the horizontal components-in the direction perpendicular to the APM-of the Bouguer gravity gradients: we clearly see a blue signal in the ZZ component as well as a red signal in the XX component with an amplitude about 3,000 mEötvös.

Continental Hotspots Tracks
For each investigated region, we computed the Bouguer, topographic and crustal gravity gradients at the surface in different rotated local frames and at various spatial scales.Then, we plotted the three corresponding scaleorientation diagrams and carefully examined the results around the APM direction to identify a signal with a spatial scale of a few hundred kilometers.The scale-orientation diagrams for the YY component of the topographic and crustal gravity gradients are presented in Text S2 in Supporting Information S1 for the four regions.In the main text, we present the scaleorientation diagram for the YY component of the Bouguer gravity gradient and maps for the radial ZZ and horizontal YY components of the Bouguer gravity gradients, in a local frame aligned with the current motion of the tectonic plate for the spatial scale with the highest amplitude.Figures for the other components are provided in Supporting Information S1.In these maps, we superimpose the expected trajectories of the hotspots according to the three mantle reference frames as well as the position and the age of the volcanism, to support the interpretation.The current positions of the hotspots are given in Table 1.We also indicate in this table the depth to which the current hot plume conduit is observed in the mantle by seismology (see Vicente de Gouveia (2017) for more details, and French and Romanowicz (2015), Montelli et al. (2006), andRitsema et al. (2011) for seismological studies).Knowing the rooting depth of the hotspots is crucial because the most stable and long-lived plumes (a prerequisite for our calculation of the reconstruction of the trajectory of hotspots in the lithosphere) are those anchored very deep in the mantle (Davaille et al., 2002;Jellinek & Manga, 2002).
In the next subsections, we will investigate five real cases in detail.These cases include the Northwestern part of Africa (with the Hoggar, Tibesti, Darfur and Cameroon hotspots), the Eastern part of Australia (Cosgrove), Greenland (Iceland and Jan Mayen), the central part of Europe (Eifel) and the Virginia (hidden hot plume).We have chosen these regions because, except for the central Europe, their tectonic structure is not too complex, allowing us to detect the signal we are looking for.

Hoggar, Tibesti, Darfur and Cameroon Hotspots in Africa
We initially focus on Northwestern Africa.We studied the trajectory of four hotspots (Hoggar, Tibesti, Darfur and Cameroon) located in a tectonically quiet region.
-The Hoggar hotspot is located in southern Algeria.It is associated with the formation of the Hoggar Mountains, which are composed of a variety of volcanic and sedimentary rocks that have been folded and faulted by tectonic activity.The hotspot is believed to have been active for at least 45 million years.The most recent volcanic activity in the area is thought to have occurred about 3 million years ago.The origin of this volcanism is still debated: plume or intraplate stress induced by the Africa-Europe collision (Liégeois et al., 2005).-The Tibesti hotspot is located in the central Sahara Desert, in northern Chad.The volcanic field it has created covers an area of about 31,000 km 2 (Vincent, 1970) and includes several volcanic mountains.The volcanic activity in this region started during Oligocene with major episodes from lower Miocene to Quaternary (Gourgaud & Vincent, 2004) and it is still on-going (Permenter & Oppenheimer, 2007).-The Darfur hotspot is located in western Sudan.It has been active for millions of years and has produced a number of volcanic features including a large volcanic field and several calderas.The most recent eruption occurred in the early Holocene, around 6,000 years ago (Franz et al., 1994).The hotspot has also caused the uplift of the Darfur Plateau, which is an elevated region of the Sahara desert that has an average elevation of around 900 m. -The Cameroon hotspot is located in the West African Rift system in Cameroon.It is associated with the formation of the volcanic line that extends from the Gulf of Guinea to Lake Chad.The hotspot is believed to have been active since at least 60 million years ago and is responsible for the formation of several volcanic mountains, including Mount Cameroon and Bioko Island.It is still considered to be active (Halliday et al., 1988).
We reconstructed the trajectories of these hotspots back to 90 Ma, using the three reference frames mentioned earlier.These trajectories are shown in Figure 7a, superimposed on the topography and intraplate volcanism.Despite the lack of clear consistency with the distribution of volcanic activity, the three obtained traces exhibit  remarkable similarity over the past 40 million years.For Darfur and Hoggar, we identified multiple volcanic episodes along their trajectories, indicating hotspot activity spanning tens of millions of years.However, the magmatic activity around Darfur follows a time sequence from N/E to S/W, which is opposite to the direction of the plume's track.In the case of the Cameroon chain, the ages of the volcanic events do not align with the hotspot tracks, suggesting a different origin, possibly of tectonic nature.In the vicinity of Tibesti, there is no volcanism older than 10-20 million years, indicating that this hotspot was not active prior to this period.
We estimated the average azimuth of the African plate motion between 40 Ma and the present, at a point at the center of the region and we found an azimuth of 35°.We subsequently generated the scale-orientation diagram for the YY component in this region (Figure 7b).We observed a noteworthy signal at a spatial scale of 300 km, with an azimuth that closely aligns with the direction of the APM.There is also a signal in the direction with an azimuth about 40°-90°away from absolute plate motion.This signal could be associated with the topography and/or the crust (see Text S2, Figures S2-1a and S2-2a and Text S3, Figure S3-b, middle column in Supporting Information S1), particularly in the NW/SE direction of the past episodes of extension over Africa, during the Mesozoic and Cenozoic (Guiraud et al., 2005).Additionally, there is a subtle peak at an azimuth of 125°, which appears to align with the direction of hotspots tracks older than 50 million years.It is worth noting that changes in plate motion direction could potentially complicate our analysis.
In this region, we calculate and plot, at different spatial scales, the isotropic vertical and directional horizontal components of the Bouguer gravity gradients in a local frame rotated by 35°E, that is, with e → X aligned with the present-day motion of the African plate (see Figure S3-a in Supporting Information S1).In Figures 7c and 7d, we plot the ZZ and YY components at the spatial-scale of 300 km.There is a significant signal in both components: we observed a blue negative signal for the ZZ component of the Bouguer gravity gradient and a red positive signal for the YY component over almost 20-40 millions years for the Hoggar, Darfur and Cameroon hotspots.The signal is more uncertain for Tibesti, but this may be due to the fact that this hotspot is not deeply anchored in the mantle (see Table 1) and, therefore, may be more affected by mantle-lithosphere interactions.A fifth similar signal in Figures 7c and 7d, located between the traces of Tibesti and Hoggar, is surprising because it does not correspond to a known hotspot.This anomaly could be associated with the Tibesti-Sirt high (plotted in black thick line) which is surrounded by two basins, the one in the West being the Murzig Basin with a trough axis aligned with the direction of the plate (Shalbak, 2015).We also clearly see, near the Cameroon volcanic chain, an anti-correlated anomaly related to the Benoue Trough (Braitenberg, 2015), schematized in the figure by a black thick line.
The scale-orientation diagram for the YY component of the topographic gradient shows a significant signal around the APM direction (Figure S2-1a in Supporting Information S1), and looking at all components of the topography gradients (Figure S3-b in Supporting Information S1), we found that we cannot separate the contribution of the RTA from that induced by the mass anomalies created by isostatic thickness variations of the crust (leading to negative anomalies under the mountains and positive ones under the basins).Even though the sources depths and widths are not the same, our filtering does not allow us to separate them.We are able to follow the surface hotspots tracks, but we are unable to interpret the gradiometric signal in terms of density.
To conclude this section on Africa, let us return to the sediments correction in the topographic gradients discussed in Part 3. We present different maps in Text S4 in Supporting Information S1 over West Africa and show that this correction has a smaller amplitude and is not correlated with the signals discussed here.

Cosgrove Hotspot Track on Australia
We move on to Australia, where there is a well-known northward drift of intraplate volcanism associated with the continental track of the hotspot linked to the Cosgrove volcanic field (D. R. Davies et al., 2015).This activity has left a record of volcanic activity between 33 and 9 Myr ago.The Cosgrove volcanism is 9 Ma old, and the current position of the associated plume (named East Australia or Cosgrove in the literature) can be calculated from the Australian plate motion relative to the mantle.It is located somewhere between King Island and Tasmania.
We have plotted the trajectory of this hotspot on the topography map (Figure 8a) for the three referentials described earlier.The three tracks are similar and cross Australia toward the north in less than 40 million years due to the fast velocity of the plate (around 5 cm/year).
The scale-orientation diagram of the YY component of the Bouguer gravity gradient is plotted in Figure 8b: it shows a dominant signal around A z = 20°, that is, in the direction of the Australian plate, whatever the spatialscale.This diagram is very similar to the scale-orientation diagram of the YY component of the topography gravity gradient plotted in Figure S2-1b in Supporting Information S1.
We then calculated the different gradients in a local frame rotated by 20°and at various spatial scales (Text S5 in Supporting Information S1).In Figures 8b and 8c, we plotted the vertical and horizontal Bouguer gravity gradients in the direction perpendicular to the trajectory (i.e., the quasi-WE direction) for a spatial scale of 300 km.The signal obtained from the horizontal component of the Bouguer gravity gradients does not show a red patch along the track and we do not find any corresponding blue patch in the vertical component, suggesting a complex origin.The signal reflected in Figure 8b may be perturbed by the N/S oceanic/continental boundary (see Text S5 in Supporting Information S1 for figures at different spatial scales) and/or by small scale convective instabilities induced by local variations in lithospheric thickness and the fast velocity of the plate (D. R. Davies & Rawlinson, 2014) or by bursts in slab flux associated with the past Pacific seafloor subduction zone (Mather et al., 2020).

Iceland Hotspot Track on Greenland
Iceland and Jan Mayen are two hotspots located in the North Atlantic Ocean.Iceland is a volcanic island located at the boundary between the North American and Eurasian tectonic plates, while Jan Mayen is a small volcanic island located about 600 km to the northeast of Iceland.Iceland is home to over 130 volcanoes, with ∼30 of them having erupted in the past 1,000 years.Jan Mayen, on the other hand, is a much smaller hotspot, with only one known volcano, Beerenberg.This volcano last erupted in 1985, and is thought to have been active for the past 5,000 years.
We reconstructed the trajectories of the Iceland and Jan Mayen plumes under Greenland since 80 million years ago.In Figure 9a, these trajectories are plotted on the topography for our three models of referentials, revealing an overall westward track with significant uncertainties in the plate motion reconstruction in this region.We also represented intraplate volcanism for different periods of time.At a latitude of 70°N, we note that Greenland experienced simultaneous volcanic episodes on its East and West coasts, between 70 and 50 Ma (cyan triangles).On the west coast, this volcanism is consistent with the passage of the Iceland plume 70 Myr ago.The volcanism of same age on the east coast may be explained by the rifting of the North Atlantic or by a complex interaction between the lithosphere and the plume (Steinberger et al., 2019).For period between 50-0 Ma, the westward trajectory of the Iceland plume seems consistent with the age of volcanism.
The trace of the Iceland hotspot over Greenland has been extensively studied in recent years to understand the history of the Iceland plume.Here is a brief summary of the main results according to the geophysical methods and tools used: -Seismic data: Tomographic studies (Knezevic Antonijevic & Lees, 2018;Kumar et al., 2005;Lebedev et al., 2018;Mordret, 2018;Pourpoint et al., 2018) have observed a thinning of the continental lithosphere, which could be associated with a thermal erosion created by the plume head but could also be interpreted as a tectonic boundary between cratonic blocks.-Geothermal flux data: Using ice-lithosphere coupled model and borehole measurements, thermal and compositional anomalies in the present-day Greenland lithosphere have been mapped and correlated with seismic velocity anomalies, which may be explained by the Iceland hotspot history (Rogozhina et al., 2016).A recent study (Alley et al., 2019) investigated the interaction between glacial-interglacial cycling and hot sublithosphere conditions left by the passage of the plume.-Geodetic data: Global Positioning System observations around Greenland of the vertical ground displacements created by the post-glacial rebound suggest a warmer sublithospheric mantle explained by the passage of the Iceland hotspot (Khan et al., 2016).-Magnetic data: An inversion of precise magnetic anomaly data has made possible to map the depth of the Curie isotherm and the geothermal heat flow (Martos et al., 2018), revealing a thicker crust and underplating by the Iceland plume, consistent with the measured Bouguer anomalies.
To the best of our knowledge, the gravity gradients have not been previously used to identify the plume track in this region.
In Figure 9b, we plotted the scale-orientation diagram of the YY component in this area.Due to the significant uncertainties surrounding the present-day direction of the plate, we denoted the APM with two vertical black lines that correspond to an azimuth ranging from 100°to 120°.Within this range of direction, we only observed a signal with a spatial scale of 200 km.
We then calculated the different gradients in a local 100°-120°rotated frame (in order to increase the signal-tonoise ratio, we have stacked the gradient maps for a range of frame rotation angles) at different spatial scales (as detailed in Text S6 in Supporting Information S1).In Figures 9c and 9d, we show the vertical and horizontal Bouguer gravity gradients for a frame rotation angle of A z = 100°-120°and a spatial scale of 200 km.We clearly observe in the YY component two coherent positive signals with an amplitude greater than 500 mEötvös, starting respectively at the present-day positions of Iceland and Jan Mayen (thick gray lines in Figure 10a).
The geometry of the topography is different from that of the plume track, which is promising for separating the sources and detecting the signal induced by the passage of the plume under the lithosphere.Notably, we found that the amplitude of the topographic gravity gradient perpendicular to the track is smaller than that aligned with the track, except in the northernmost and southernmost regions of Greenland (compare XX and YY components in Figure S6-b in Supporting Information S1).
In Figure 10a, we overlaid a possible trajectory for these two plumes (thick gray lines, defined as the maximum of the YY Bouguer gravity gradient) on the horizontal Bouguer gravity gradient for a spatial scale of 200 km in a local frame rotated by 100°-120°.Our Iceland track is within the range of tracks based on several geodynamic reconstructions (Rogozhina et al., 2016).Our interpretation of the Iceland hotspot trajectory from gravity gradient anomalies seems to differ from that obtained by calculating the thermal thickness of the lithosphere (e.g., Figure 9 of Artemieva ( 2019)) or from magnetic anomalies (Figure 1 of Martos et al. (2018)).To understand this difference, we compared maps of the YY horizontal component of the Bouguer gravity gradient at different spatial scales: 200 km in Figure 10a and 500 km in Figure 10b.We observe a large red patch (depicted by a thick dashed black line) in Figure 10b that seems consistent with previously published lateral variations in thermal thickness.If we associated this red patch with the passage of a single hotspot, it would follow a NW/SE trajectory, but it would not explain the volcanism of the west coast of Greenland, around the latitude of 70°N.However, the fine spatial resolution of the GOCE satellite allowed us to separate the trajectories of two hot plumes (Figure 10a): Iceland and Jan Mayen.In this case, the volcanism of the West coast is located on the trajectory of the Iceland plume.

Eifel Hotspot Tracks on Europa
The Eifel hotspot is a volcanic hotspot located in western Germany and active for the past 65 million years.There have been no known eruptions in the past 10,000 years.However, the region is still considered geologically active.
The European crust and lithosphere have been the focus of several seismological and gravimetric studies.These studies have highlighted the locations of the various tectonic units in the region (Kaban et al., 2018;Tessauro et al., 2008;Yegorova & Starostenko, 2002).This is a tectonically complex region, with the Alps in the south, the Pannonian Basin in the east and the German Polish Basin in the north.Geochemistry, Geophysics, Geosystems 10.1029/2023GC011200 We begin by examining the trajectory of the Eifel plume since 70 Myr in the three mantle reference frames plotted on the topography (Figure 11a).The path of the Eifel plume has roughly moved eastward with trajectories significantly different before 40 million years ago, depending on the chosen reconstruction.We have also superimposed the reconstructed positions of intraplate volcanism on the topography.We note that the eastward movement of the Eifel plume seems too slow to be correlated with the intraplate volcanism, suggesting another origin for the latter.In addition to recent volcanism around the present-day position of Eifel (zone labeled A in green dashed line in Figure 11a), there are two groups of peri-alpine tertiary alkaline volcanism correlated with the main faults in the peri-alpine rift zone: one between 10 and 20 Ma, near the Rhine graben (zone B), and another between 20 and 40 Ma, further east (Eger graben, zone C) (Dèzes et al., 2004).Separating the volcanism resulting from the passage of a hot spot from that linked to regional tectonics is a challenge that would require a detailed study of the chemistry of the lava, which is beyond the scope of this article.So, what information can we obtain from gravity gradiometry?
Recently, the gravity gradients measured by GOCE have been used to map the Moho discontinuity in Central Europe (Lenczuk et al., 2019).The gravity gradient maps obtained after various corrections were calculated at the altitude of the satellite (around 255 km) and highlight different geological structures on a large spatial scale.However, structures as fine as those associated with the trajectory of a hotspot did not appear in these maps.Therefore we computed gradient maps at the Earth's surface in different directions and at different spatial scales in order to search for a trace of the path of the Eifel hotspot.We created a scale-orientation diagram of the YY component on this region (Figure 11b).We notice a clear signal in the APM direction at the spatial-scale of 300-350 km.There also clear signals in the scale-orientation diagrams of the YY component of the topographic and crustal gravity gradients  in Supporting Information S1).
In Figures 11c and 11d, we plotted the vertical and horizontal components of the Bouguer gravity gradients in the direction of the present-day APM of Eurasian plate (A z = 50°) for a spatial scale of 300 km.The signal is strongly perturbed by the presence of the Alps, which created a large topographic gradient in APM direction.We do not see any red patch along the track nor any blue patch in the vertical component.
As this regions is densely covered in ground measurements, we used the gradients from the EGM2008 model, filtered at the spatial scale of 150 km to achieve a finer spatial resolution.We plotted the different gradients at various spatial scales (150, 250, and 350 km) (Text S7 in Supporting Information S1).We clearly see a correlation between the volcanism and the topographic gradients at the spatial scale of 150 km (Figure S7b, top in Supporting

Information S1
).Consequently, it seems difficult to find any trace of the hotspot path under the lithosphere in this region because of its geological complex and eventful tectonic.

Hidden Hotspot Tracks in Virginia
Let us finally revisit the previous seismic study conducted in Virginia by Chu et al. (2013), which revealed a thermomechanical erosion of the lithosphere by a "hidden" mantle plume.In that study, the authors identified a linear seismic anomaly in the lower lithosphere characterized by a reduction in P-wave velocity and high attenuation, interpreted as a hotspot track.
To determine the trajectory of this "hidden" plume, we calculated its present-day position from the 75 millionyear-old Elliot Kentucky kimberlite using the three mantle reference frames, as shown in Figure 12a: it is located to the north of the New England seamounts, in a region of the North Atlantic with no traces of volcanic activity.In addition, our calculations indicate that between 90 and 70 million years ago, the plume moved eastward with an azimuth of approximately 110°.We also plotted on Figure 12a the trajectory of the Bermuda hotspot which crossed the American plate about 1,000 km further south.
In the region around the position of the Elliot Kentucky kimberlite, the lithosphere is about 200 km thick, and the RTA model proposed by Chu et al. (2013) would be equivalent to a parallelepiped approximately 100 km thick, twice as thick as our example in Section 2. We computed a scale-orientation diagram for this synthetic anomaly and found a maximum for a 350 km spatial-scale.We obtained an amplitude of the horizontal component YY just above the parallelepiped of about 380 mëotvös for a density anomaly of 10 kg/m 3 .Using observed GOCE data, we calculated the ZZ and YY components of Bouguer gravity gradients in a local rotated frame with azimuth of 110°(Figures 12b and 12c) at a spatial scale of 350 km and superimposed the path of this hidden hotspot.We found no clear signal in the ZZ component, except for the signal induced by the Appalachian mountains in the east.However, in the gray box representing the seismic anomaly (see Figure 3 in Chu et al. ( 2013)), we observed a signal in the YY component with an amplitude of approximately 500 mëotvös, which could be associated with a RTA of approximately Δρ = 13 kg/m 3 under this part of the trajectory.In our synthetic model, such a Δρ will involve a dynamic topography just above the parallelepiped of 380 m which is precisely the value proposed by Chu et al. (2013) in their Figure S8.
Therefore, it appears that a thermal residual anomaly persists beneath the lithosphere and is detectable through gradiometry.

Other Hotspots Tracks Over the Continents
To conclude, we examined gravity signals from all the continental hotspots listed in Table 1 (SI) of Koppers et al. (2021) (see also Jackson et al., 2021): two in West America (Yellowstone and Raton in Figure S8 in Supporting Information S1), four in East Africa (from north to south: Afar, East Africa, Lake Victoria, and KenyaDome, in Figure S9 in Supporting Information S1), and two located near a boundary between a continent and an ocean (Erebus in Antarctica (Figure S10 in Supporting Information S1), and Hainan in Southeast Asia (Figure S11 in Supporting Information S1)).For each map, we present the 3 components (XX, YY, and ZZ) of the Bouguer gravity gradients and of the topographic gravity gradients, in the direction of present-day plate motion of the investigated region and at a 250 km spatial scale.
However, interpreting the results is difficult for two reasons: -First, tectonics may be very complex.For example, on the west coast of the United States, the presence of the Rocky Mountains and the Sierra Nevada prevents us from detecting a recent trace left by the hotspots Yellowstone and Raton.Similarly, in the eastern part of Africa, the trajectories of the Victoria, East-Africa and Afar hotspots are strongly disturbed by the opening of the Red Sea.Geochemistry, Geophysics, Geosystems 10.1029/2023GC011200 -Second, the signals may be strongly correlated with topography, or the thickness of the lithosphere may vary greatly, perturbing the gravity signal.For example, Erebus or Hainan hotspots may be affected by these factors.

Discussion
In this study, we investigated the traces in gradiometry data of a possible RTA along the trajectory of a hotspot under the continental lithosphere.This anomaly has been quantified by different models taking into account both conduction and thermal erosion at the base of the lithosphere: it is on average a few hundred kilometers wide and 50 km thick, with a density anomaly around 10 kg/m 3 .These parameters obviously depend on the plate velocity (the faster the plate, the shorter the duration of heating from below), on the temperature anomaly and on the velocity of advection in the plume conduit as well as its size.The intensity of the surface gravity perturbation also depends on the thickness of the plate: the thicker the continental plate, the more distant from the surface the source of mass.To be observed in the gradiometry data and in particular in a direction along the movement of the plate, the investigated region must not be tectonically too complex.For example, on the West Coast of the United States, the Bouguer gravity gradient anomaly associated with the Yellowstone or Raton hotspots since 40 Myr (with an amplitude of a few hundred mEötvös) is completely masked by the signal created by the crustal roots of the mountains (Text S10, Figure S10 in Supporting Information S1).
This demonstrates the importance of understanding the lateral variations in crustal thickness and density anomalies.If we were to calculate the gravitational potential and its derivatives induced by such a crustal model, we will obtain a signal with an amplitude one order of magnitude larger than the observed one.A portion of the crust is in isostatic equilibrium, leading to mass compensation and subsequently reducing the associated gravitational signal.A potential issue arises when starting with a model of crustal density anomalies and Moho depth (e.g., Crust1.0) and computing the radial displacement of the Moho depth to ensure isostatic equilibrium.But such a study goes far beyond the scope of this article.Correcting the observed gravimetric signals for the effects of the Earth's crust is a crucial challenge that could significantly improve our understanding of mantle dynamics.
We worked on the Northwestern part of Africa, a tectonically stable region, where we were able to highlight the trace of three of the four hotspots over a few tens of millions of years.The most significant scale is around 300 km.However due to a significant signal in the topographic gradient in this direction and at this spatial scale, it is difficult to interpret the observed amplitude in terms of density anomaly at the base of the lithosphere.
In eastern Australia, several studies have already highlighted a correlation between the South/North alignment of volcanoes according to their age and the northward motion of the Australian plate on the mantle.The plume associated with this volcanism is not currently observed: its last appearance, 9 Ma ago, had given rise to the Cosgrove volcanism.Its current position, reconstructed from an Australian plate motion model, would be close to Tasmania.It is difficult to find a signal in the second order derivative in the East/West direction of the gravitational field, even at small spatial scales, as the gravity anomaly is disturbed by the lateral East/West variations of the Lithosphere-Asthenosphere Boundary and of the crustal thickness, related to the proximity of the ridge, which, coupled with the fast velocity of the plate, can induce small scale convection instability (D. R. Davies & Rawlinson, 2014).
We continued our study through Greenland, a region intensively studied recently to better understand the melting of the ice sheet.The passage of the hotspot currently under Iceland is often invoked to explain heat flow anomalies at the base of the lithosphere, anomalies of viscosity and seismic velocity.The different models of Greenland plate motion on the mantle do not seem to be well constrained and the age of the volcanism suggests a complex coupling mechanism between both plume and tectonics.At small spatial scales (200 km), the second-order derivative of the gravitational potential perpendicular to the plate direction highlights two corridor-like parallel anomalies, not correlated with the topography.One of them corresponds to the trajectory of the Iceland plume, the other one to that of a hotspot located a little further North, named Jan Mayen, with an amplitude of a few hundred mEötvös.Investigations at larger spatial scales (500 km) would not have made it possible to separate these two distinct traces.
Next, we focused on the European plate and the Eifel hotspot.The intraplate volcanism that occurred between 0 and 40 million years ago displays a west-east trajectory that appears to be faster than the movement of the European plate in the mantle.We created vertical and horizontal Bouguer gravity gradients maps with a spatial resolution of 300 km, but due to the complex tectonic history of the region, identifying any evidence of the hotspot path beneath the lithosphere appears to be challenging.
Finally, we searched for the thermomechanical erosion of the lithosphere by a "hidden" plume in Virginia revealed by the seismic study of Chu et al. (2013).Our analysis suggests that a thermal residual anomaly could persists beneath the lithosphere and be detectable through gradiometry.
To end this study, gravity signals from all continental hotspots were examined, and maps presenting the Bouguer gravity gradients and topographic gravity gradients were created.However, the results are difficult to interpret due to the complex tectonics in some areas, such as the presence of the Rocky Mountains and the Sierra Nevada in the United States and the opening of the Red Sea in Africa.Additionally, for some hotspots, such as Erebus and Hainan, the signals may be affected by topography or variations in lithosphere thickness.
Better constraining the trajectory of hotspots in the continental lithosphere is an important challenge because it could significantly improve the design of a mantle-linked reference frame, which currently relies solely on oceanic hotspot traces.
In this study, we used a global approach to identify local anomalies.However, to improve the accuracy of the results, a more precise crustal density model should be used for the calculation of topographic and crustal corrections in certain regions.This is especially important in areas where volcanic structures align with the plume's path, as uncertainties in the density and thickness of crustal roots can produce gravity anomalies aligned with the trajectory, potentially leading to misinterpretation of the signal.Sediments can also contribute to such misinterpretations.

Conclusion
This study explores the potential of utilizing gravity gradients detected by the GOCE satellite to identify the RTA resulting from a plume track across continental lithosphere.We predicted the expected gravity gradients associated with a hotspot track and found that they exceed current detection levels.An analysis based on a scaleorientation diagram shows the feasibility of detecting such gradients within the background variations in the gravity field.However, if we find that the plume track should generate detectable gravity gradients, these may be overshadowed by gravity gradients arising from crustal density heterogeneities.Using a scale-orientation diagram, we attempted to isolate the gravity gradients attributable to plumes, applying this approach to all continental hotspots.We identified some plume tracks where we anticipated them, particularly in Africa and Greenland.Our findings indicate that gradiometric data have the potential to detect and track hotspots in the continental lithosphere for tens of millions of years, while determining the depth of the corresponding mass sources remains challenging.
The next phase of this research will apply our multi-scale directional analysis of the gradients to the highly precise World Digital Magnetic Anomaly Map on the continents (Lesur et al., 2016).The passage of a hotspot beneath the lithosphere disrupts the isotherms in the crust, especially those near the Curie point of crustal minerals.This Geochemistry, Geophysics, Geosystems 10.1029/2023GC011200 disruption has the potential to induce crustal magnetization, leading to the generation of an observable magnetic field at the Earth's surface.

Figure 1 .
Figure 1.(a) Simple model of a parallelepiped residual thermal anomaly; Induced perturbations in the geoid (b) and in the dynamic topography (c), in the six components of the gravity gradient tensor (d) and of the Bouguer gravity gradient tensor (e) obtained after removing the contribution of the (dynamic) surface topography.

Figure 2 .
Figure 2. Perturbations in the topography and in the XX component of the Bouguer gravity gradient tensor induced by the different sources.(a) simple model with two parallelepiped of constant density schematizing a tectonic topography (C 1 ) and its isostatic crustal root (C 2 ), with a width about 250 km, and the residual thermal anomaly (C 3 ) with a Gaussian density (Figure 1).The thick black line on the top of cylinder C 1 is the total topography (isostatic topography and dynamic topography as plotted in (b)) (b) Surface topography in meter: isostatic topography (gray line), dynamic topography induced by the cylinder C 3 (red line) and total topography (black line).(c) XX component of the Bouguer gravity gradients associated with the different sources (C2 in black line, C3 in red), and for the sum of the 2 sources (in green).

Figure 3 .
Figure 3. Analysis at different spatial scales of the XX Bouguer gravity gradients component induced by the parallelepiped mass source C3 oriented along the azimuth 90°N superimposed to a random density anomaly in a 50 km thick crust with an amplitude about 10 kg/m 3 (left column), 25 kg/m 3 (middle column), and 50 kg/m 3 (right column).From top to bottom: full synthetic signal, signal filtered at 300, 350, 400, and 450 km spatial scales.
enhance the gravity signals of sources oriented along orientations, at each point of the Earth surface, we rotate the local horizontal axes e → X and e → Y around the fixed radial axis e → Z by an angle noted A z .The obtained rotated local spherical frame is denoted { e →

Figure 4 .
Figure 4. Y′Y′ Bouguer gravity gradient component in rotated frames for a parallelepiped mass source oriented along the azimuth 30N.The angle of rotation A z of the spherical frame along its radial axis is indicated on each panel as well as the rotated frame.

Figure 5 .
Figure 5. Scale-Orientation diagrams of the YY Bouguer gravity gradient component for the parallelepiped mass source C3 oriented along the azimuth 90°N (top) superimposed to a random density anomaly in a 50 km thick crust with an amplitude about 10 kg/m 3 (middle), and 25 kg/m 3 (bottom) (d).

Figure 6 .
Figure 6.(a) Trajectory of Hawaii (in purple with dot every 10 Myr) for Doubrovine et al. (2012) reference frame superimposed to the topography (meter) and intraplate volcanism.(b) Scale-Orientation diagrams of the YY Bouguer gravity gradient component over Pacific (longitude between 178°and 208°and latitude between 20°and 30°N).ZZ isotropic vertical component (c) and YY directional horizontal component (d) of Bouguer gravity gradients computed at a 350 km spatial scale in a local reference frame rotated about 115°with respect to the vertical axis.Thin black lines: contour lines of the topography with a contour interval of 2,000 m.

Figure 7 .
Figure 7. (a) Trajectories of Cameroon, Hoggar, Tibesti, and Darfur (with dot every 10 Myr) for the three referentials described in the text (in white for O'Neillet al. (2005), in black forMüller et al. (1993) and in purple forDoubrovine et al. (2012)) superimposed to the topography (meter) and intraplate volcanism.The two thick black lines indicate the position of the Murzuq basin (to the north) and the Benoue basin (to the south) (b) Scale-Orientation diagrams of the YY Bouguer gravity gradient component over Africa (longitude between 7°and 35°and latitude between 5°and 35°N).ZZ isotropic vertical component (c) and YY directional horizontal component (d) of Bouguer gravity gradients computed at a 300 km spatial scale in a local reference frame rotated about 35°with respect to the vertical axis.Thin black lines: contour lines of the topography with a contour interval of 300 m.

Figure 8 .
Figure 8.(a) Trajectories of Cosgrove (with dot every 10 Myr) for the three referentials described in the text (in white for O'Neill et al. (2005), in black for Müller et al. (1993) and in purple for Doubrovine et al. (2012)) superimposed to the topography (meter) and intraplate volcanism.The Cosgrove volcanic field, 9 Ma old, used to calculate the present-day position of the plume is plotted as a green dot.(b) Scale-Orientation diagrams of the YY Bouguer gravity gradient component over East part of Australia (longitude between 145°and 150°and latitude between 40°and 21°S).ZZ isotropic vertical component (c) and YY directional horizontal component (d) of Bouguer gravity gradients computed at a 300 km spatial scale in a local reference frame rotated about 20°with respect to the vertical axis.Thin black lines: contour lines of the topography with a contour interval of 200 m.

Figure 9 .
Figure 9. (a) Trajectories of Iceland and Jan Mayen since 80 Ma (with dot every 10 Myr) for the three referentials described in the text (in white for O'Neill et al. (2005), in black for Müller et al. (1993) and in purple for Doubrovine et al. (2012)) superimposed to the topography (meter) and intraplate volcanism.The present-day position of these hotspots is plotted as a green dot.(b) Scale-Orientation diagrams of the YY Bouguer gravity gradient component over Greenland (longitude between 310°and 325°and latitude between 66°and 72°N and longitude between 305°and 335°and latitude between 72°and 80°N).ZZ isotropic vertical component (c) and YY directional horizontal component (d) of Bouguer gravity gradients computed at a 200 km spatial scale in a local reference frame rotated about 100°-120°with respect to the vertical axis.

Figure 10 .
Figure 10.YY component of the Bouguer gravity gradient: (a) at a spatial scale of 200 km for a 100°-120°rotation of the local spherical frame and (b) at a spatial scale of 500 km for the 100°-120°rotation of the local spherical frame.The two gray thick curves represent the possible paths of the two hotspots, Iceland and Jan Mayen.The closed dashed black curve encircles the mean trajectory of the Iceland hotspot when viewed on a very large spatial scale.

Figure 11 .
Figure 11.(a) Trajectories of Eifel since 70 Myr (with dot every 10 Myr) for the three referentials described in the text (in white for O'Neill et al. (2005), in black forMüller et al. (1993) and in purple forDoubrovine et al. (2012)) superimposed to the topography (meter) and intraplate volcanism.The present-day position of this hotspot is plotted as a green dot.(b) Scale-Orientation diagrams of the YY Bouguer gravity gradient component over central part of Europa (longitude between 5°and 20°a nd latitude between 48°and 53°N).ZZ isotropic vertical component (c) and YY directional horizontal component (d) of Bouguer gravity gradients computed at a 300 km spatial scale in a local reference frame rotated about 50°with respect to the vertical axis.Thin black lines: contour lines of the topography with a contour interval of 300 m.

Figure 12 .
Figure 12.(a) Trajectories of the so-called "hidden hotspot" and of Bermuda for the three referentials described in the text (in white for O'Neill et al. (2005), in black for Müller et al. (1993) and in purple for Doubrovine et al. (2012)) superimposed to the topography (meter) and intraplate volcanism.The present-day position of Bermuda is plotted as a green dot.Red triangles denote the 75-million-year old kimberlite in Kentucky and the 100-million-year old kimberlite in Texas.The gray box contours the lower lithosphere seismic anomalies (see Figure 3 of Chu et al. (2013)).(b) ZZ isotropic vertical component and (c) YY directional horizontal component of the Bouguer gravity gradients computed at a 350 km spatial scale in a local reference frame rotated by ∼110°around the vertical axis.Thin black lines: contour lines of the topography with a contour interval of 250 m.

Table 1
Present Day Positions of Hotspots and Mantle Plume Source Depth