Gravity, Topography, and Melt Generation Rates From Simple 3‐D Models of Mantle Convection

Convection in fluid layers at high Rayleigh number (Ra ∼106) have a spoke pattern planform. Instabilities in the bottom thermal boundary layer develop into hot rising sheets of fluid, with a component of radial flow toward a central upwelling plume. The sheets form the “spokes” of the pattern, and the plumes the “hubs.” Such a pattern of flow is expected to occur beneath plate interiors on Earth, but it remains a challenge to use observations to place constraints on the convective planform of the mantle. Here we present predictions of key surface observables (gravity, topography, and rates of melt generation) from simple 3‐D numerical models of convection in a fluid layer. These models demonstrate that gravity and topography have only limited sensitivity to the spokes and mostly reflect the hubs (the rising and sinking plumes). By contrast, patterns of melt generation are more sensitive to short‐wavelength features in the flow. There is the potential to have melt generation along the spokes but at a rate which is relatively small compared with that at the hubs. Such melting of spokes can only occur when the lithosphere is sufficiently thin ( ≲80 km) and mantle water contents are sufficiently high ( ≳100 ppm). The distribution of volcanism across the Middle East, Arabia, and Africa north of the equator suggests that it results from such spoke pattern convection.


Introduction
What is the planform of mantle convection?The largest, and most obvious, planform in the convection system is that associated with plate motions, which involves horizontal scales as large as 10,000 km (e.g. the Pacific plate).However, it is also clear that convection takes place at shorter horizontal scales.This scale of convection manifests in hot-spot volcanism, and the swells and troughs in gravity and topography observed at wavelengths of around 1,000 to 2,000 km (McKenzie, 1994;Crosby et al., 2006;Crosby & McKenzie, 2009).A well-cited example of this comes from Africa (Holmes, 1965;McKenzie & Weiss, 1975;Burke, 1996;Jones et al., 2012), which shows a clear pattern of swells and troughs across the continent.The strong correlation between the patterns of gravity and topography, and in particular the characteristic ratio of around 50 mgal km −1 between the two, has been used to strongly argue for convective support of the topography (McKenzie, 1994;Crosby et al., 2006;Jones et al., 2012).
In theory, maps of gravity and topography should provide information on the planform of mantle convection.But extracting that information is challenging.While gravity provides information on density variations within the Earth, the process of inverting gravity data for density is highly non-unique.However, density variations within the Earth on length scales of hundreds of kilometers are not arbitrary, but are controlled by the fluid dynamics of convection.
The aim of the present manuscript is to get a better understanding of the surface expressions of mantle convection from a series of the simplest possible numerical simulations of convection, and to compare these with geophysical and geological observations.The fluid dynamical problem is one that has been extensively studied and discussed in the literature: Rayleigh-Bénard thermal convection in a 3D rectangular box with a fluid of constant viscosity at Rayleigh numbers around 10 6 .This problem was first studied using laboratory experiments in large aspect ratio tanks, using a layer of silicone oil whose depth was a few centimetres (Busse & Whitehead, 1971,9;Richter & Parsons, 1975).These experiments showed that the planform of the convection was a spoke pattern, with hot rising plumes joined to each other by hot sheets near the lower thermal boundary layer, and cold sinking plumes by colds sheets near the top thermal boundary layer.White (1988) carried out similar experiments using a fluid whose viscosity was a strong function of temperature, and showed that at high Rayleigh number the convective planform was also a spoke pattern.
We carried out our calculations with an infinite Prandtl number fluid, and used the Boussinesq approximation throughout.We used large aspect ratio boxes, similar to those previously used in tank experiments, to allow the convective circulation to determine its own planform rather than being dominated by the lateral boundaries.We use our experiments to show which features of the observations are readily explained by the simplest models, and which features are not.We do not attempt to construct a realistic model of the Earth.In this respect our aim is the same as that of the early tank experiments, and, in addition, to extend them to encompass the observables: the gravity field, surface deformation and melt generation.In particular, by allowing the density variations to arise naturally from fluid dynamics, such experiments allow an exploration of short wavelength (< 100 km) temperature variations, which are likely to be most clearly expressed by volcanism.
Two effects that are known to be important in the Earth are not taken into account in the simple constant viscosity Boussinesq model we use.The first is the variation of viscosity with temperature.It is this effect that produces plates, and therefore our modelling does not include the dynamics of plate motions.The other effect is viscous dissipation, which is intimately related to vertical density variations that result from lithostatic pressure (Spiegel & Veronis, 1960;Jarvis & McKenzie, 1980;Schubert et al., 2001).Viscous heat generation has little effect on the circulation even when the relevant term in the equations cannot be neglected.Such heating occurs in boundary layers where temperature gradients are large.As a result the entropy and potential temperature (appendix A.5) is little affected (Jarvis & McKenzie, 1980), though the convection becomes more time dependent.
The approach taken here complements a popular alternative approach to modelling gravity and topography using information from seismic tomography (Hager & Richards, 1989;Flament et al., 2013).In such studies, estimates of density variations within the mantle are inferred from tomography and used to make predictions of gravity and dynamic topography.Though these studies have had some success at predicting the very long wavelength (> 6,000 km) features of Earth's gravity, they depend on knowing the relationship between density ρ and the seismic velocities V P and V S , and also on a rheological model of the mantle.Furthermore, since these calculations are based on seismic tomography, they are limited by its resolution, which is not yet sufficient to map rising and sinking plumes in the upper mantle.Our approach also departs from the common assumption of many Earth Scientists, who believe the convective planform of mantle convection consists solely of plumes and the plates.This assumption arises from the work of Wilson (1963) and Morgan (1971), who showed that the relative motion between major volcanic centres beneath plate interiors was sufficiently slow that they could be used to define a single world-wide reference frame.Their ideas have been enormously influential.But they are based on an intuitive conception about the planform of high Rayleigh number convection, rather than on fluid dynamical experiments.They also predate our understanding of polybaric melt generation.
The manuscript is organized as follows.Section 2 describes the fluid dynamical simulations, and how they are scaled to parameter values appropriate for the Earth's mantle.Section 3 discusses the predicted gravity, topography and their spectral properties.Section 4 discusses melt generation.Section 5 compares the results from the fluid dynamical experiments with the observed gravity and topography, and with the volcanism of Africa and the Middle East, and conclusions follow in section 6.An appendix provides further technical details on the simulations and data processing.

Numerical experiments
We ran 12 numerical experiments of isoviscous thermal convection in a rectangular box.Temperature was fixed at the top and bottom boundaries.To examine the influence of dynamical boundary conditions, runs were made for all combinations of freely slipping or rigid boundary conditions on the top and bottom boundaries.Reflection boundary conditions were applied at the side boundaries.
All convection simulations were performed using v2.01 of the ASPECT mantle convection code (Dannberg & Heister, 2016;Heister et al., 2017;Bangerth et al., 2018).The code was used to solve the dimensionless versions of the Boussinesq governing equations of thermal convection in an 8 × 8 × 1 rectangular box, through a small modification of the "convection-box" example discussed in the ASPECT manual.The governing equations in dimensional form are where v is the velocity, P is the difference in pressure from hydrostatic, θ is potential temperature, ρ 0 is the reference mantle density, α is the thermal expansivity, η is the viscosity, and κ is the thermal diffusivity.We assume constant thermal conductivity, constant heat capacity, constant viscosity, and constant thermal expansivity.In dimensionless form the governing equation are where all lengths have been scaled by the layer depth d, time by the diffusion time d 2 /κ, pressure by ηκ/d 2 , and potential temperature by the potential temperature difference ∆T p across the layer.Just one dimensionless parameter describes this simple system, the Rayleigh number, defined by Runs were performed at three different Rayleigh numbers: Ra = 10 5 , 3 × 10 5 , and 10 6 .A uniform resolution of 32 cells in the vertical and 256 cells in the horizontal was specified for all simulations.Quadratic finite elements were used for temperature and velocity, and linear finite elements for pressure.Simulations were run until the system reached a quasi-steady state, which was monitored by examining the behaviour of the mean temperature and RMS velocity over time.Each simulation ran for a minimum of six times the thermal time constant for the layer (= 6 d 2 /(π 2 κ)).Snapshots of the 12 experiments are shown in Figure 1.The images are made to mimic the shadowgraph visualization technique commonly used in laboratory experiments (Busse & Whitehead, 1971,9;Richter & Parsons, 1975;Whitehead & Parsons, 1977;White, 1988).To make a shadowgraph in the laboratory, light is shone from below through the layer of fluid and projected onto a screen.Refraction causes the light to focus and defocus according to the temperature variations within the fluid through which it passes.The effect is to make upwellings appear dark and downwellings appear bright.Mathematically, a shadowgraph produces a plot of the Laplacian of the vertically averaged temperature (Jenkins, 1988;Travis et al., 1990), and this is how the images in Figure 1 were computed.The planform of convection at these Rayleigh numbers is spoke pattern (Busse & Whitehead, 1974).Hot and cold plumes develop from the top and bottom of the layer: the "hubs" of the planform.Hot and cold sheets, the "spokes" of the planform, radiate from the plumes.These sheets are formed by instabilities in the top and bottom boundary layers, and the shearing associated with the plumes suppresses instabilities with other geometries.The planform is time-dependent.
Several well-known features stand out from the shadowgraph pictures.The first is that as the Rayleigh number increases the thicknesses of the upwellings and downwellings become narrower.For the free-free and rigid-rigid cases the upwellings and downwellings are symmetric, but there is a notable asymmetry in the nature of upwellings and downwellings in the free-rigid and rigid-free cases (Kvernvold, 1979;Weinstein & Christensen, 1991).Also notable is the change in horizontal distance between upwelling and downwelling with boundary condition: this distance is shorter for rigid-rigid simulations than for free-free simulations.
Our main interest here is the surface observables associated with mantle convection.The simulations were performed using dimensionless variables, with the dimensionless Rayleigh number as the only control parameter.To make predictions about observable quantities, these simulations must be scaled appropriately to Earth-like values.For most quantities we can simply use typical mantle estimates, and these values are given Table 1.We consider here just upper mantle convection, and so choose as an appropriate layer thickness d = 600 km to scale all lengths.
Choosing an appropriate scaling of temperature is less straightforward.The numerical simulations essentially provide dimensionless potential temperatures.The potential temperature is the temperature that the mantle material would have if it were moved to the Earth's surface isentropically and without melting (McKenzie (1970), appendix A.5).
The scaling of the dimensionless temperature is required to satisfy two conditions: The first is that the average interior temperature must correspond to a mantle potential temperature of 1315 • C.This choice ensures that the thickness of the oceanic crust is 7 km, generated by isentropic decompression beneath a spreading ridge using the parametrisation of Katz et al. (2003).The second condition arises from the top of the convecting system not being at the Earth's surface.We envisage that there is a rigid mechanical boundary layer (MBL) which separates the Earth's surface from the top of the convecting system.The temperature near the base of the  4); ∆g RMS , root mean square of the gravity anomaly at the top of the convecting region (as plotted in Figure 3).
MBL is where T /T s , where T s is the melting temperature, is highest, and therefore the viscosity is lowest (e.g.Frost & Ashby, 1982).For this reason we used both free-slip and rigid boundary conditions on the top surface of the convecting box.A stress-free boundary is probably the better approximation to the behaviour of the real Earth.Heat transfer through the MBL is purely by conduction, and we assume the thermal profile is linear through this region, and equal to 0 • C at the Earth's surface.At the top of the convecting region we assume that both the horizontallyaveraged temperature and heat flux are continuous with that in the MBL.Finally we determine the thickness of the MBL by prescribing a lithospheric thickness, which we define as the intersection of the linear conductive profile of the MBL with an isentropic profile at the interior mantle potential temperature.In what follows the lithosphere thickness is fixed at 100 km, except when considering melt generation where we have varied this parameter, as melt generation is particularly sensitive to it.This choice of interior potential temperature and lithospheric thickness fixes the heat flux for all simulations to be 50 mW m −2 , which is similar to that through old sea floor (Hasterok, 2013).The resulting associated parameters, which include the upper mantle viscosity, are given in Table 2.The horizontally-averaged potential temperature profiles after scaling are shown in Figure 2.With Rayleigh numbers in the range 10 5 to 10 6 , Table 2 shows mantle viscosity values vary from 1.4 × 10 20 Pa s to 5.3 × 10 21 Pa s, which are around the range expected for the upper mantle (∼ 10 21 Pa s).The higher the Rayleigh number, the lower the inferred viscosity, the smaller the potential temperature difference across  the layer, and the thicker the MBL.Appendix A.4 discusses further the behaviour of the parameters with Rayleigh number.

Gravity and topography
Examples of the gravity and topography at the top of the convecting box when there is no lithosphere present are shown in Figures 3 and 4, calculated from the expressions in Parsons & Daly (1983) (appendix A.1).The gravity field takes account of the contribution from the thermal expansion of the fluid and that from the topography (assuming deformable top and bottom boundaries).Both shadowgraphs and the plots of gravity and topography are filtered versions of the temperature field in the box, but they involve distinctly different kinds of filters.Since the shadowgraph represents the Laplacian of the vertically-averaged temperature, it is equally sensitive to temperature variations at all depths, and acts to emphasise short-wavelength features.By contrast, gravity and topography act to attenuate short-wavelength features in the temperature field, particularly those at depth (Parsons & Daly, 1983).Before comparing the results from the numerical models with the observations a further filter needs to be be applied, to account for the effect of the overlying mechanical boundary layer.This layer has a number of effects.First, gravity anomalies are attenuated by the thickness of the layer.Second, the elastic properties of the overlying plate acts to filter out the topographic expression of short-wavelength features, with the magnitude of this effect dependent on the assumed effective elastic thickness (T e ).Figures 5 and 6 show the expected gravity and topography at the surface, after applying a filter for the mechanical boundary layer, assuming an elastic thickness T e = 30 km.The effect of all this filtering is to remove much of the short wavelength information that is visible in the shadowgraph images.In particular, the filters emphasise the "hubs" of the convective pattern, and suppress the expression of the "spokes".As will be seen in the next section, the effect of the filter on melt generation is even more important.As the thickness of the MBL increases it first suppresses melt generation in the spokes and then in the hubs.
The images in Figures 5 and 6 demonstrate three key points: First, the predicted amplitude of gravity anomalies is comparable with the observed anomalies (compare scale of Figure 5 with Figure 13(c)).Second, the spatial patterns of gravity anomalies and topography are sensitive to the assumed boundary conditions.Third, there is a good correlation between the topography and the gravity in plots of the values of gravity and topography at each spatial location (Figure 7), which closely resemble similar plots that have been made using observed values of gravity and residual topography (e.g.see Figures 6 and 7 of Crosby & McKenzie (2009)).The slope of these plots represents a characteristic average value of the admittance (ratio of gravity to topography) associated with the convection.A characteristic air-loaded admittance between 43-53 mgal km −1 is inferred from Figure 7, corresponding to values of 30-37 mgal km −1 when overlain by water.These values are similar to those previously reported (Parsons & Daly, 1983).x (1000 km) free-free x (1000 km) free-free x (1000 km) free-free x (1000 km) free-free  An alternative way to assess the predictions of gravity and topography is to work in the frequency domain rather than the spatial domain.Figure 8(a) and (b) show the power spectral density of gravity and topography at the top of the convecting box (appendix A.2), along with estimates for the Earth.The figures show that the convection experiments have most power over a range of wavelengths from around 500 to 2000 km, where the estimates of gravity anomalies are comparable in magnitude with those of the Earth.This behaviour reflects the fact that the convection organises into cells where the distance from one upwelling to the next is around twice the layer depth, corresponding to a wavelength around 1200 km.Away from this broad peak the power from the convection experiments decays.At the short wavelengths shown in the plot, this decay in the power spectral density scales roughly as k −3 .The different boundary conditions in the convection experiments lead to subtle differences in the gravity spectra in Figure 8.For example, there is greater power at long-wavelengths for those simulations with free-slip bottom boundary conditions than rigid.There is also greater power at short wavelengths for the simulations with rigid top boundary conditions than free-slip top boundary conditions.
The power spectral density of gravity for the Earth is notably different from that of the convection experiments.Overall, the Earth's power spectral density decays broadly with wavenumber as k −2 , a power law which is referred to as Kaula's rule.In detail, the slope of decay flattens slightly in the 2000 km to 500 km wavelength band, before becoming steeper at wavelengths shorter than 500 km, but not decaying as steeply as in the convection experiments.The difference between the convection experiments and the Earth for wavelengths shorter than 500 km is to be expected: on Earth surface loading causes short wavelength topography that is supported by elastic stresses in the plate.The short wavelength power in the Earth's gravity field arises from surface loading, not mantle convection.At long wavelengths (>2000 km) mantle convection should play an important role in determining the gravity field on Earth, but the convection experiments here have notably less power than the Earth.Thus the upper mantle convection modelled here does not account for the magnitude of the long-wavelength portion of the Earth's gravity field.
Figure 8(b) shows the corresponding spectra for predictions of dynamic topography.The shape of the topography spectra of the convection runs is different from that for the gravity.The topography spectra are broadly flat up to around a wavelength of 1000 km, and then decay steeply at shorter wavelengths.Thus dynamic topography has relatively more long wavelength power than does gravity, and this can be seen in the space domain plots of Figure 3 and Figure 4 -the topography plots look smoother than the gravity plots (as noted by Craig & McKenzie (1987)).
It is more difficult to compare predictions of dynamic topography with observations.The reason for this is that there are significant long-wavelength features in the Earth's topography that are not associated with mantle convection e.g. the difference in elevation between the oceans and the continents, and other topography associated with variations in crustal thickness.A recent attempt has been made to estimate the power spectrum of dynamic topography by Hoggard et al. (2016), by fitting spherical harmonics to point observations of residual depth in the oceans, and a fixed scaling of long-wavelength gravity anomalies in the continents.Their estimate of power indicates a Kaula-rule-like decay in the power spectral density as k −2 , and is plotted in Figure 8(b).At wavelengths longer than 2000 km the convection experiments are not able to explain the Hoggard et al. (2016) estimate of dynamic topography.At shorter wavelengths spectra become more comparable, although only a limited comparison can be made because the Hoggard et al. (2016) estimate is limited to spherical harmonic degree 30 (a wavelength of 1300 km).The Hoggard et al. (2016) estimates are based on air-loading in the continents and water-loading in the oceans.The power spectra of the convection experiments in plotted in Figure 8(b) assume air loading.The effect of water loading would be to increase the amplitude of the power spectral density of the topography by a factor of (1 − ρ w /ρ 0 ) 2 ≈ 2, which would represent only a small shift on the log-scale plot of Figure 8(b).
In addition to comparing the observed and calculated individual spectra of gravity and topography, it is also of interest to look at their relationship to one another (Parsons & Daly, 1983).This relationship is typically characterised in terms of the admittance and coherence of the two signals, with the topography taken as input and the gravity as output (appendix A.3). Figure 8(c) shows the admittance in the spectral domain, calculated from the numerical experiments.The value of the admittance increases with increasing wave number (decreasing wavelength), approximately as log(k) for the range of wave numbers shown.The behaviour of admittance with wave number is similar for the different boundary conditions over the range of interest, with significant departures only noticeable at long wavelengths (>2000 km).Indeed, at long wavelengths the admittance for the free-rigid case becomes negative around a wavelength of 3700 km (as noted by Parsons & Daly (1983), see their Figure 5).Figure 8(c) is consistent with the behaviour in the cross-plots of Figure 7.The gravity anomalies have most power at wavelengths around 1000 km, corresponding to an air-loaded admittance of around 45 mGal km −1 in Figure 8(c), which is broadly the slope obtained from the cross-plots.The corresponding coherence in Figure 8(d) is close to 1, reflecting the good correlation between the two observables that can be seen in the space domain plots.Only at wavelengths longer than 2000 km is a weak coherence between the signals seen, and then only significantly for the free-rigid case.This behaviour mirrors the admittance at long-wavelengths for the free-rigid case, where gravity and topography correlate positively except for wavelengths longer than 4000 km when they correlate negatively.
The frequency domain plots in Figure 8 illustrate the spectral properties of the gravity and topography at the top of the convecting box, as shown in the space domain in Figures 3 and 4. When considering signals at the Earth's surface, the spectral properties will be further modified by the filtering effect of the MBL on top.How significant this effect is depends on the effective elastic thickness and the thickness of the MBL.For an elastic thickness T e = 30 km as used in Figures 5 and  6, the wavelength at which the Fourier coefficients of topography are reduced by a factor of 2 is λ 1/2 flex = 330 km.Wavelengths shorter than this will be significantly attenuated by the flexural filtering; wavelengths longer than this will not.The topography spectra in Figure 8 would only be significantly different at wavelengths shorter than λ 1/2 flex were a flexural filter to be applied.The effect of flexural filtering on the gravity is more complicated, and acts to produce a modest increase in the gravity signal in a wavelength band around that associated with the MBL thickness and that associated with flexure (see appendix A.1 and Figure 15 for further discussion).

Melt generation
The generation of melt, its separation from its source regions and the time τ required for it to move from its source to the surface, have all been extensively studied, both theoretically using two-phase flow equations, and observationally, using a variety of geochemical approaches.The two-phase flow equations show that basaltic melt separates from its source regions when the melt fraction by volume φ 0 exceeds ∼ 0.5% (McKenzie, 1985a).Studies of the composition of abyssal peridotites (Johnson et al., 1990;Warren, 2016) show that the incompatible elements that were present in the source before melting occurred have been removed by the melt.Estimates of the melt fraction present during the melting are between 0.2 and 0.7% (Slater et al., 2001;Liang & Peng, 2010).Estimates of φ 0 and τ can also be obtained from measurements of U-series disequilibria (McKenzie, 1985b;Kokfelt et al., 2003;Stracke et al., 2006;Koornneef et al., 2012;Turner et al., 2016) .Most estimates give φ 0 ≤ 0.5% and τ ≤ 1 ka.Perhaps the strongest constraint on the values of both φ 0 and τ comes from modeling the generation of melt by deglaciation.When most of the ice covering Iceland melted at the end of the last glaciation the melt production rate suddenly increased (Maclennan et al., 2002;Eason et al., 2015).The thickness of the ice that melted was about 2 km, increasing the melt fraction present in the source region by only about 0.2% (Jull & McKenzie, 1996;Eksinchol et al., 2019).This increase was sufficient to generate large shield volcanoes within about 1 ka of the removal of the ice.These models and observations all show that melt generated by decompression melting in the upper mantle rapidly moves to the surface, and that no appreciable volume remains in the source region.We therefore calculated the rate of melt production by simply vertically-integrating the melting rate over the thickness of the layer in which melt was being produced.
Figure 9 shows the calculated rate of melt production for the four Ra = 10 6 simulations, with a lithospheric thickness of 80 km, and assuming various water contents.The hydrous melting parametrisation of Katz et al. (2003) was used (appendix A.5). Melt production in Figure 9 only occurs where the mantle is upwelling, where the gravity anomalies and topography are positive (Figures 5 and 6).However, the spatial extent of the regions of high melt rate are smaller than the regions of positive gravity anomaly.Most of the melting takes place in the regions directly above the narrow (∼ 60 km diameter) upwelling plume conduits, where the decompression rate is greatest.However, melting also takes place, albeit to a lesser degree, in a broader region around the conduits and above some of the rising sheets connecting neighbouring plumes (i.e. on the spokes as well as the hubs of the spoke-pattern).Melt generation along some of the spokes is particularly clear for the free-free simulation, with linear bands of relatively low melt production connecting concentrated centres of relatively high melt production.
The rate of melt production is particularly sensitive to the thickness of the lithosphere and the water content of the mantle.Figure 10 illustrates the effect of varying these two parameters for the free-rigid Ra = 10 6 simulation.At low water contents, x (1000 km) free-free Figure 9: Vertically integrated rate of melt production beneath an 80 km thick lithosphere with 100 ppm water.Note that the colourscale for melt rate is logarithmic.Numbers in white give the total rate of melt production in km 3 ka −1 for each contiguous zone of melting.Regions of total melt production less than 1 km 3 ka −1 are not labelled.and beneath thick lithosphere, melting is restricted to the hubs in the spoke-pattern, if indeed melting happens at all.However, high water contents and thin lithosphere result in melting along the spokes.As Figure 10 illustrates, for lithosphere as thick as 180 km, melting is suppressed unless the water content is sufficiently high (∼ 1000 ppm).For lithosphere as thin as 70 km, melting along both the hubs and the spokes can be seen even for modest (∼ 100 ppm) water contents.The water content of the convecting upper mantle is probably between 100 and 200 ppm (Michael, 1995;Saal et al., 2002) .In contrast the water content is likely to be considerably greater where the base of the lithosphere has been enriched by metasomatism.It is not straightforward to estimate the water content of such regions using that in the nodules brought up by magmas such as kimberlites from depths of 100-200 km, because they are often infiltrated by the host magma.Protons are especially mobile.
More reliable estimates can be obtained from the Ce concentration, because Ce and H have similar bulk partition coefficients between magma and peridotite (Aubaud et al., 2004), and the Ce concentration in nodules is less affected by infiltration than is that of H (Erlank et al. (1987) p 283).The Ce concentration in the commonest class of nodules is ∼ 10 ppm (Erlank et al., 1987), or about 10× that of the convecting upper mantle.Therefore the metasomatically enriched region at the base of thick old lithosphere probably has a water concentration of ∼ 1000 ppm.Where the lithosphere is thin, Figure 10 shows that such high water concentrations will lead to widespread melting along spokes, and at the hotter hubs melting at depths as great as 180 km.

Terrestrial Observations
The numerical experiments described above show that the planform of mantle convection will be most obviously expressed in the surface observables when both the elastic thickness and lithospheric thickness are small.Figure 11 shows that the lithospheric thickness exceeds 120 km over large regions of western and southern Africa.In these regions the volcanism consists of small-volume alkalic eruptions, such as kimberlites, that contain high concentrations of carbonates and hydrous minerals.Where kimberlites are diamond-bearing, Figure 11 shows that the lithospheric thickness generally exceeds 150 km.The limited spatial resolution of surface wave tomography, of ∼ 250 km, probably accounts for the few diamond-bearing locations in Figure 11 that appear to have thinner lithosphere.
The lithospheric thickness within much of the rectangle marked by the thick continuous black line in Figure 11 is less than 80 km, and the horizontal extent of this region is similar to that of the numerical experiments.In the eastern and northern parts of the area there are sufficient surface gravity measurements to allow T e to be estimated from the transfer function between the free air gravity anomalies and the topography (McKenzie & Fairhead, 1997), giving values of 3-4 km (see supplementary material).The coherence method and Bouguer anomalies provides only an upper bound, not an estimate, of the value of T e (McKenzie, 2016).Within the box marked by the dotted lines in Figure 11 there are few surface gravity measurements.Instead the satellite gravity field DIR-R5 can be used to show that the admittance is about 50 mGal km −1 between wavelengths of 200 and 1000 km, and that the elastic  11 marked by the dotted line.The gravity field was calculated from DIR-R5 (Bruinsma et al., 2014) by setting the coefficients from l = 2 to 7 to 0, and applying a taper f, = (l − 7)/5, to those from l = 8 to 11 (wavelengths 3333 to 5000 km).A low pass filter falling to 1/2 at 250 km was applied to remove the short wavelength anomalies associated with elastic flexure.The topography is taken from ETOPO5.The admittance was calculated using the topography as input, gravity as output.The two dotted lines show the flexural admittance for two values of the elastic thickness T e .(c)-(e) show corresponding plots for the Pacific (see supplementary material for maps) and the Indian Oceans (see Figure 14).The gravity anomalies were calculated from the DIR-R5 coefficients with those of degree 2 set to 0. The admittances in (c) and (e) were calculated from the ratio of the spectral coefficients.The dotted lines show the flexural admittance for T e = 20 km.thickness is probably less than 4 km (Figure 12(a) and (b)).Therefore in this region and wavelength band both the gravity and topography are controlled by convection.The same is not the case in southern Africa, where the elastic thickness is about 30 km (McKenzie et al., 2015).
Figure 13 shows maps of gravity, topography, and subaerial volcanism across Africa and Arabia.The volcanism beneath the Red Sea and in the Afar results from upwelling of the mantle between separating plates.However, in Ethiopia, Kenya and Kivu the upwelling from the limited extension is insufficient to account for the extensive volcanism.There is no obvious orientation of the topography and gravity anomalies, probably because Africa is almost stationary with respect to the hotspot frame.The correspondence between the features in Figure 13 and the maps from the numerical experiments is striking.Volcanism is almost entirely restricted to regions where the lithospheric thickness is less than 70 km.The only clear exception is Kivu, where the spatial resolution of the surface wave tomography is probably insufficient to resolve the thickness of the thin lithosphere beneath the Western Rift.The linear volcanic feature extending from Kenya to the Kars Plateau resembles similar linear features in Figure 9, with localised regions of concentrated upwelling being associated with positive gravity and elevation, and with enhanced volcanism.Like the numerical experiments, the volcanism is more localised than are the associated positive gravity and topographic anomalies.What is less clear is which of the four combinations of boundary conditions best fits the observations.Figures 5 and 6 show that the observed horizontal scales of the convective features are probably larger than those of the rigid-rigid, and smaller than those of the freefree, experiments.The scales of the anomalies in the other two experiments are similar, both to each other and to the observed scales.The rigid-free experiment has broadly circular patterns of positive gravity anomalies surrounded by linear negative anomalies, whereas the free-rigid case has the opposite.The viscosity of the lower mantle is greater, and that of the asthenosphere immediately below the lithosphere less, than that of the upper mantle.These viscosity variations suggest that the free-rigid experiment should match the observed patterns better than the rigid-free case.However, the patterns in Figure 13 are not obviously more like the free-rigid features than the rigid-free ones.Furthermore the variation of viscosity with temperature, which has been ignored, may have a strong influence on the geometry, and in particular whether the planform is dominated by rising or sinking fluid in the hubs.At low Rayleigh numbers the answer to this question is controlled by the sign of dη/dT , with the flow in the plumes being in the direction of increasing viscosity (Segel & Stuart, 1962).If the same is true at large Rayleigh numbers the hubs in both these experiments will consist of hot rising material.
The Cameroon Line forms a curve, similar to features in Figure 9.As the lithosphere thickness increases to the SW, where the Line lies beneath Atlantic lithosphere, the volcanism decreases.The association of positive gravity anomalies, elevated topography and volcanism is clearly expressed even in relatively small features like Aïr and Darfur.Many volcanic centres, such as Aïr and Hoggar, and Haruj and Tibesti, are linked to each other by lines of positive gravity and topography, where there is limited volcanism.An especially obvious feature extends from S. Arabia to Anatolia, where the volcanism is beneath Western Arabia, not the Red Sea.Such lines are most clearly visible in Figure 13(c) and (e) where the lithosphere is thin beneath NE Africa, Arabia and Anatolia.All these features are similar to those of the numerical experiments, and all are consistent with a spoke pattern of convection existing beneath the region.In particular, and as expected from Figure 9, the extent of the volcanism is controlled by variations in lithospheric thickness, and is limited in the south and east by the thick lithosphere of the Congo Craton.
The Cameroon Line and the line of active volcanism that extends from Afar to eastern Anatolia have long puzzled geophysicists, because their geometry is difficult to reconcile with a planform of mantle convection consisting of plumes.Sleep (1997Sleep ( , 2008) ) and Ebinger & Sleep (1998) argued that these linear volcanic features were produced by lateral flow from plumes in channels beneath the lithosphere.One problem with this proposal is that the composition of the volcanics along the Cameroon Line shows so little variation (Fitton, 1987;Lee et al., 1994).
A different model was proposed by Milelli et al. (2012) which emphasised the location of the volcanism, which has remained in the same region of Africa as the continent has moved.They argued that this behaviour required the volcanism to result from thermal instabilities in the lower part of the lithosphere, rather than being the surface expression of convective upwellings in the upper mantle below the plates.The numerical experiments described above show that the expected planform of upper mantle convection is that of hubs joined by spokes, both of which can generate melt if the lithosphere is sufficiently thin.The observed linearity of the Cameroon Line and other features in NE Africa therefore requires no special explanation.The experiments also show that Milelli et al. (2012)'s observations may also have a simple explanation, since the melting rate, and not the planform, is controlled by the lithospheric thickness, and the volcanism of the Cameroon Line lies along the northern edge of the thick lithosphere of the Congo Craton.In contrast the line of volcanism from the Afar to Kars lies within a region of relatively uniform lithospheric thickness.It is therefore unlikely that such linear features all form from edge convection like that discussed by King & Anderson (1998).
In general it is not possible to compare the melt generation rates in Figure 9 with those observed because they are so rarely estimated by the geologists who map the volcanics.An exception is Mount Cameroon, which is the most active volcano in Africa.Its eruption rate was estimated by Suh et al. (2003) to be about 700 km 3 /Ma.When the lithospheric thickness is 80 km the larger hubs in the free-rigid experiment in Figure 9 produce about 6 × 10 4 km 3 /Ma and the smaller ones 1 × 10 3 km 3 /Ma.The rates of melt generation in the numerical experiments can therefore easily account for the observed rates.But they are quite inadequate to account for the production rates that occur during major flood volcanism, which commonly exceed 1 × 10 6 km 3 /Ma.Like the long wavelength gravity anomalies discussed below, simple isoviscous upper mantle convective models cannot account for such events.
The box in Figure 11 is too small to be used to study the long wavelength components (wavelengths > 1000 km) of the Earth's dynamic topography and gravity.These are best studied in the Indian and Pacific Oceans.The elastic thickness of old oceanic lithosphere is about 20 km (e.g.McKenzie et al. (2014)).Therefore wavelengths greater than about 800 km are little affected by the thickness and elastic properties of the lithosphere (Figure 12(c) and (e)).In many continental regions the surface topography is dominated by isostatically compensated variations in crustal thickness.Such structures are less common in oceanic regions, where the bathymetry is principally controlled by the age of the lithosphere.The effect of plate cooling can be removed by using a depth-age model, and regions of thick crust removed by hand (e.g.Crosby et al., 2006;Crosby & McKenzie, 2009;Hoggard et al., 2017).
The resulting residual topography is then largely supported by convection.Figure 12 shows a comparison between free air gravity, calculated from DIR-R5 with the coefficients of degree 2 set to zero, and Crosby's values of residual depth in the Indian and Pacific Oceans.Those obtained from Hoggard's values of residual depth are similar and are shown in the supplementary material.Figure 14 shows maps of such anomalies in the Indian Ocean: those for the Pacific are illustrated in the supplementary material.The admittance between wavelengths of 1000 and 2000 km is about 30 mGal/km, in agreement with the values from the numerical experiments in Figure 8.However, at wavelengths greater than about 2000 km the gravity and residual depth anomalies in both oceans cease to be coherent.This incoherency is particularly striking in the Indian Ocean, where the large negative gravity anomaly covering the NE part of the Ocean (Figure 14(a)) has no expression in the residual depth (Figure 14(b)).At wavelengths greater than 2000 km the observed power spectrum also differs from that calculated from the box models (Figure 8(a)).The power in the observed gravity field continues to increase at wavelengths longer than 2000 km, unlike that from the numerical experiments.This behaviour shows that simple isoviscous convective models cannot account for the longest wavelength part of the Earth's gravity field, and suggest that it is not maintained by upper mantle convection.Though the size of the boxes used for the numerical experiments is too small to determine the power at such wavelengths accurately, there is no suggestion in the planforms that the presence of lateral boundaries governs the scale of the convection.

Conclusions
The numerical experiments described above show that the observed gravity and topographic anomalies are reproduced by the simplest isoviscous fluid dynamical model of thermal convection.The wavelength at which convective support dominates elastic support is controlled by the elastic thickness T e , and varies from about 200 km where T e ≤ 4 km in NE Africa to ∼ 500 km in southern Africa, where T e ∼ 30 km.Melt generation occurs where mantle material moves upwards, and is therefore controlled by the lithospheric thickness, and not by the value of T e .The correspondence between the volcanism and the gravitational and topographic anomalies in NE Africa is striking, and shows that they all result from the convective circulation.
There is no similar correspondence between the results from the numerical experiments and the gravity and residual depth anomalies at wavelengths greater than 2000 km.Furthermore the absence of correlation between residual depth and gravity anomalies with wavelengths greater than 2000 km in the Pacific and Indian Oceans is unlike the behaviour observed at shorter wavelengths in these oceans.The sim-  ple numerical models discussed here cannot account for the long wavelength gravity anomalies with spherical harmonic degrees l ≤ 20.The close correspondence between the calculated and observed topography, gravity and volcanism suggests that it should be possible to use the surface observables where the lithosphere is thin and T e is small, together with the isoviscous convective equations, to map the convective circulation in the upper mantle.8, but at the surface after filtering through the MBL and an elastic layer with thickness T e = 30 km.Notice the significant decrease in the power of gravity anomalies at short wavelengths and a moderate increase in the power at intermediate wavelengths.
variations at depth, which will be attenuated by a factor e −ktm if there is a MBL of thickness t m on top.The gravity anomalies in Figure 5 are shown after this filtering process, which is calculated assuming the MBL is laterally uniform.The corresponding spectra are shown in Figure 15a.If ∆g 0 and h 0 represent the gravity and topography estimated at the top of the convecting box, the corresponding gravity ∆g 1 at the surface can be calculated in the Fourier-domain from assuming air-loading.This equation splits the gravity anomaly into two components: The term on the far right represents the gravity due to the surface topography, which is attenuated according to the factor F (k) in ( 8).The other term represents the density variations at depth, which are attenuated as e −ktm due to the finite thickness of the MBL.At long wavelengths gravity anomalies are unchanged by this additional filtering.At wavelengths significantly shorter than both the flexural wavelength λ flex and the wavelength λ m = 2πt m associated with the MBL, gravity anomalies are strongly attenuated, as both the e −ktm and F (k) filters tend to zero.At intermediate wavelengths, particularly in the wavelength band around λ m and λ 1/2 flex , gravity anomalies actually slightly increase in magnitude due to this additional filtering.This behaviour occurs because a MBL acts to separate mass excesses at the surface (associated with positive gravity anomalies) from mass deficits at depth (associated with negative gravity anomalies).The thicker the MBL, the greater the attenuation of the negative anomalies, and the larger the net positive anomaly.Correspondingly, there is a modest increase in the admittance due to the addition of a MBL (Figure 15c; see section 2.3 of Crosby & McKenzie (2009) for further discussion).It should be noted that having a laterally uniform thermal structure in the MBL is a poor approximation, as it implies a discontinuity in heat flux between the MBL and the top of the convecting box (where the heat flux varies laterally).However, we have made such an approximation here because the convection simulations fix the temperature at the top of the convecting layer, rather than at the Earth's surface.A better approach would model the temperature structure of the MBL during the convection simulations.However, such a modification is unlikely to make more than minor changes to the results of the calculations, because the thickness of the MBL is small compared to that of the convecting layer (Table 2), and the temperature on the upper boundary of the MBL is fixed.If the temperature structure of the MBL is included in the calculations it is then no longer accurate to obtain the topography and gravity for variable thicknesses of the MBL simply by filtering in the spectral domain.

A.2 Power spectral density
Power spectral density (PSD) estimates were calculated using the method described by Rexer & Hirt (2015).The initial data (either gravity or topography) is a regularly spaced grid of points representing a region of dimensional extent L x by L y .Let the number of grid points in the x-direction be N , and the number in the y-direction be M .The initial data is given as a matrix of values d rs where r = 0, 1, . . .N −1, and s = 0, 1, . . .M − 1.The convection simulations have reflecting boundary conditions at the sides, so the natural Fourier representation to use is a Discrete Cosine Transform of the first type (DCT-I), which is equivalent to a discrete Fourier Transform of a 2N − 2 by 2M − 2 extended grid of data exploiting the even symmetry.The discrete Fourier coefficients f pq are defined by equation ( 8) of Rexer & Hirt (2015), The even symmetry extends the data such that for N ≤ r ≤ 2N − 2 the value is taken from the original grid at r = 2N − 2 − r, and for M ≤ s ≤ 2M − 2 the value is taken from original grid at s = 2M − 2 − s.In dimensional variables, the corresponding wavenumbers of the transform are k x p = πp/L x and k y q = πq/L y .Owing to the reflection boundary conditions, the discrete Fourier coefficients f pq are real and even in both directions.A 2D grid of power spectral density can then be computed from equation (10) of Rexer & Hirt (2015), Finally, the 2D-PSD were then azimuthally averaged in wavenumber space with a bin-width of 1.3 × 10 −3 km −1 , to produce the 1D-PSD profiles that are shown in Figures 8(a) and Figure 8(b) (with units of either mgal 2 km 2 or km 2 km 2 ).Data for the Earth is typically given in terms of spherical harmonic coefficients, which need to be manipulated before they can be compared directly with the 1D-PSD profiles calculated for the Cartesian geometry of the convection simulations.This process is also described by Rexer & Hirt (2015).The spherical harmonic degree l can be related to the Cartesian wavenumber k by the Jeans relation approximation where R is the radius of the Earth.An estimate of the 1D Cartesian PSD can be obtained from φ PSD (k) = 4πR 2 P l 2l + 1 ( 14) where P l is the power at spherical harmonic degree l (see equation ( 13) of Rexer & Hirt (2015)).

A.3 Admittance and coherence
The admittance in Figure 8(c) was computed as where •, • represents the cross-power of the signals as a function of wavenumber k, calculated in the same way as the power spectra by multiplying the Fourier coefficients and then azimuthally averaging.Since the Fourier coefficients are real, the admittance is also real.The coherence in Figure 8(d) was calculated similarly as

A.4 Rayleigh number scalings
Boundary layer theory suggests that there should be systematic power-law scalings with Rayleigh number for properties of the convecting system, such as the thickness of boundary layers, and the heat flux.Table 3 shows such approximate scaling laws that have been obtained from the 12 numerical runs presented here.These scaling laws should be used with caution: they were determined from the properties of the system at a single snapshot in time, for a limited range of Rayleigh numbers.However, they illustrate the broad trends.As expected, Nusselt number increases with Rayleigh number: the behaviour for the free-free system as Nu ∝ Ra 0.3 is close to to the 1/3 power law expected from boundary layer theory (McKenzie et al., 1974).Since the dimensional scaling used in the main text essentially fixes the heat flux, the potential temperature difference across the layer given in Table 2 scales as the inverse as the Nusselt number scaling, i.e. approximately as Ra −0.3 .Table 3 also shows scalings for dimensionless topography and dimensionless gravity, which all scale roughly as Ra −0.3 although there are some differences in detail.This scaling can be understood in broad terms from the expectation that the dynamic topography should be proportional to the boundary layer thickness (Parsons & Daly, 1983).Since the scaling used in Table 3 includes a ∆T p factor, the scalings of the dimensional RMS gravity and topography in Table 2 go roughly as Ra −0.6 .The behaviour of the power spectral density with Rayleigh number is illustrated in Figure 16 for the free-rigid simulations.The principal effect of changing Rayleigh number is to change the amplitude of the power spectral density.This can be understood from the scalings of RMS gravity with Rayleigh number above, and the Parseval's theorem result which relates the square of the RMS value ∆g RMS to its power spectral density φ ∆g PSD (k).Since RMS gravity in Table 2 scales roughly as Ra −0.6 , the power spectral density would be expected to scale as Ra −1.2 .This effect can be seen in Figure 16: reducing the Rayleigh number by a factor of 10 leads to a little over a order of magnitude shift in the amplitude of the spectra.In addition to the change in amplitude, there are some more subtle changes in the spectra associated with changing the Rayleigh number.The lower Rayleigh number simulations appear to have relatively higher power at long wavelength than short wavelengths.This is as expected from the nature of the boundary layers, which are thicker for the lower Rayleigh number runs.

A.5 The relationship between temperature and potential temperature
To perform melting calculations it is necessary to convert from dimensionless potential temperature back to real temperature.In this section we describe this conversion, and justify the approximate form of energy conservation that has been used in (6).In dimensional variables, conservation of energy can be written where S is the specific entropy of a fluid parcel, k is the thermal conductivity (assumed constant), and Ψ is the viscous dissipation.The potential temperature θ can be defined in differential form as where c p is the specific heat capacity at constant pressure (also assumed constant).The energy equation ( 18) becomes where κ = k/(ρc p ) is the thermal diffusivity, and viscous dissipation has been neglected.The use of (6) as a dimensionless governing equation for potential temperature is justified provided the approximation is accurate, and that viscous dissipation is sufficiently small to be neglected.

A.5.1 No melting
In the absence of melting, the differential of potential temperature can be related to those of temperature and pressure through the standard relationship where α is the thermal expansivity.The principal variation in pressure is hydrostatic.Writing dP = −ρg dz, ( 22) can be written in terms of temperature and depth as where h a is the adiabatic scale height, defined by h a = c p /(αg) ≈ 3, 300 km.Integration of ( 23) yields the relationship between temperature and potential temperature in regions which are not partially molten, where z ref is a reference depth, the depth at which potential temperature is chosen to be equal to real temperature.This depth is chosen to be the Earth's surface in this work.
From (23) it follows that The magnitude of the second and third terms on the right hand side relative to the first term scales approximately as l/h a and (l/h a ) 2 where l is a typical scale over which the temperature varies.If that scale l were the whole of the convecting layer then l/h a = d/h a = 0.18 (a parameter known as the Dissipation number), which is relatively small.In fact, the length scale of the vertical temperature variations will be much smaller than the layer depth, with boundary layer thicknesses on the order of 100 km or less, giving l/h a = 0.03.Thus both the second and third terms on the right hand side of (25) are sufficiently small that the approximation in ( 21) is well justified in the regions that are not partially molten (McKenzie, 1970).The small Dissipation number for upper mantle convection also justifies the neglect of viscous dissipation term in ( 18).An additional approximation has been made in writing the buoyancy term on the right-hand side of the Stokes equation in ( 2) in terms of the potential temperature θ.
Formally, density variations in the fluid are determined by the actual temperature, not potential temperature, and the right-hand side of (2) should be −ρ 0 gαT ẑ.The convective flow is driven by horizontal gradients in the actual temperature, not the potential temperature.From (24) it follows that the horizontal gradients in temperature are related to the horizontal gradients in potential temperature by with a similar expression for the y-derivative.The horizontal gradients of potential temperature and actual temperature differ by an exponential factor whose magnitude is at most the exponential of the Dissipation number.For the upper mantle convection we consider here, this is a relatively small difference, and justifies the approximation made in using the potential temperature in (2).

A.5.2 Melting
The convection simulations provide a 3D grid of potential temperature (entropy) within the box.To turn this into melting rate, the hydrous melting parametrisation of Katz et al. (2003) was used to calculate the expected degree of melting F at each grid point assuming isentropic decompression melting to the given potential temperature and pressure at each grid point.The original parametrisation of degree of melting in Katz et al. (2003) is given with pressure and temperature as the thermodynamic variables.This parametrisation can be recast in terms of pressure and entropy (or potential temperature) by numerically integrating the relevant differential expressions.The differential expression for entropy when melting is where F is the degree of melting, ∆S is the specific entropy difference between the two phases, and The Katz et al. (2003) parametrisation accounts for the different thermal expansivities α s , α f ; and densities ρ s , ρ f of the two phases (solid and liquid respectively), but the specific heat c p is assumed identical for both phases.All parameter values used here are identical to those in Table 2 of Katz et al. (2003), with the exception of the specific entropy difference between the two phases which we set as ∆S = 400 J kg −1 K −1 .The parametrisation provides the degree of melting F as a function of temperature T and pressure P which can be expressed in differential form as Given a parcel of material that is subsolidus at a given potential temperature θ and depth z, (24) gives the relationship between temperature and depth (or pressure) throughout the subsolidus region.Once the material crosses the solidus the relationship between temperature and pressure at constant entropy can be obtained by numerically integrating equations ( 27), (28), and ( 29) with dS = 0. Knowing the temperature and pressure then allows F to be calculated.From a series of constant entropy integrations for different potential temperatures a parametrisation of F as a function of entropy and pressure was generated.Using this entropy parametrisation, the grid of potential temperature and depth values in the box were converted to a grid of F values.This grid of F was then converted to a melting rate Γ using where v is the velocity.The time derivative was calculated using a first-order accurate finite difference approximation.The spatial gradient was calculated using a second-order accurate finite difference approximation.The melting rate Γ was vertically integrated to produce the plots shown in Figures 9 and 10.Only those regions where the melt rate was positive (i.e.melting) were included in the vertical integral.In the calculation of these figures, the advective term (v • ∇F ) was larger by more than an order of a magnitude than the time dependent term (∂F/∂t), and an excellent approximation to the melting rates can be obtained from the advective term alone.
The calculation of melting rates were performed here as a postprocessing operation after running standard single-phase convection simulations.We assume that all the melt that is generated moves to the surface, and that none remains in the source region to freeze as the mantle material cools.It should be noted that this calculation neglects potentially important back-effects that melt can have on the flow, e.g.arising from the buoyancy of the melt, and the thermal effects of the consumption of latent heat.Indeed, in regions where melt is present, the approximation in (21) can cease to be good approximation.However, the melting regions are only a small proportion of the overall domain, and the changes in temperature due to melting are small.For the free-rigid case with a lithospheric thickness of 80 km the average temperature change in the regions undergoing melting is 10 • C, and the maximum change in the whole box is 90 • C. Furthermore such changes in temperature occur within the thermal boundary layer where heat is transported by conduction and where the temperature variations are large whether or not melting occurs.The effect of the temperature changes resulting from melting on the large-scale dynamics is therefore negligible.
There are other back-effects of melt extraction on the convective flow that have been neglected.When melt is extracted from the mantle, the remaining residue has a different density than it had before the melt was extracted.There is thus the potential for this depletion by melting to change buoyancy forces, and hence the flow.However, the density changes in the residue are small.Even for 20% melt extraction the relative density changes on depletion are on the order of -0.5% (Schutt & Lesher, 2006), equivalent to a density change from temperature variations of 125 • C.
Another potential back-effect that has been neglected arises from the effect of volcanic loading on melt production rates.That changes in loading at the Earth's surface can influence melt production rates is well-known from studies that have looked at the volcanic response to changes in ice cover in Iceland (Jull & McKenzie, 1996;Maclennan et al., 2002;Eksinchol et al., 2019).Beneath Iceland most melt generation occurs in the upper 100 km of the mantle, where the upwelling rate, of about 10 mm/a, is driven by the separating plates.The thickness of the ice on Iceland reached about 3 km, equivalent to a thickness of rock of about 1 km, over about 10 5 years, corresponding to an equivalent rock accumulation rate of 10 mm/a.Melt generation therefore ceases during the construction of the icecap, and all the melt that would normally have been generated during 10 5 years is instead released when the ice melts, in about 10 3 years.The behaviour of melt generation within the convecting region beneath the lithosphere is very different.For the free-rigid case with a lithospheric thickness of 80 km shown in Figure 9, the upwelling rate of the solid mantle where the melt generation rate is fastest is 27 mm/a.The accumulation rate of melt at the surface is about 300 m/Ma, or 0.3 mm/a.This rate is therefore about 1/100th of the upwelling rate, and will have no significant effect on the melt generation rate.Moreover, volcanic loads at the surface will be eroded over time, and their influence on the mantle beneath is further attenuated by the finite elastic strength of the overlying lithosphere.S3) from the admittance, taking the topography as input, gravity from Eigen6c as output.

Figure 1 :
Figure 1: Artificial shadowgraphs of the 12 numerical experiments of thermal convection in an 8 × 8 × 1 rectangular box.Each column is at a given Rayleigh number, and each row is at a given choice of boundary conditions for the Stokes flow.In each case the first word corresponds to the top boundary, the second word to the bottom boundary; so free-rigid refers to a free-slip top and rigid bottom boundary condition.

Figure 2 :
Figure 2: Profiles of horizontally-averaged potential temperature.Zero depth represents the top of the convecting region (the base of the mechanical boundary layer).The four panels show the four choices of boundary condition.The dimensional scaling is such that the interior potential temperature is 1315 • C and the lithospheric thickness is 100 km.

Figure 3 :
Figure3: Gravity anomalies at the top of the convecting box for the Ra = 10 6 simulations.The region shown is approximately the same size as the region outlined by the thick black line in Figure11.Each panel shows a different combination of boundary conditions.In each case the first word corresponds to the top boundary, the second word to the bottom boundary; so free-rigid refers to a free-slip top and rigid bottom boundary condition.

Figure 4 :
Figure 4: Dynamic topography at the top of the convecting box for the Ra = 10 6 simulations.

Figure 5 :
Figure5: Gravity anomalies expected at the Earth's surface when the lithospheric thickness is 100 km.Plots are as in Figure3, except the attenuation of gravity anomalies through the mechanical boundary layer and an elastic plate with T e = 30 km has been taken into account.

Figure 6 :
Figure 6: Topography expected at the surface, after flexural filtering through an elastic plate with T e = 30 km.

Figure 7 :Figure 8 :
Figure 7: Cross plots of gravity against topography at the top of the convecting box for the four Ra = 10 6 simulations.In each case a geometric regression line has been calculated, and marked on each plot is the slope Z of that line, and r 2 , the square of the correlation coefficient.If the boxes are overlain by water the values of Z are reduced by a factor of (ρ 0 − ρ w )/ρ 0 0.7

Figure 10 :
Figure10: Vertically integrated rate of melt production for the free-rigid simulation with Ra = 10 6 , showing the effect of varying the water content (as labelled horizontally in ppm), and different lithospheric thicknesses (as labelled vertically in km).

Figure 11 :Figure 12 :
Figure11: Lithospheric thickness calculated from surface wave tomography(Priestley et al., 2018).The thick square box indicates a region that is of the same horizontal extent as the convection simulations.The elastic thickness of the region within the dotted lines is likely to be less than 4 km (Figure12).

Figure 13 :
Figure 13: Maps of the volcanism, lithospheric thickness, gravity, and topography within the region marked by the heavy black line in Figure 11.The gravity and topography are shown without (Figure 13(c) and (e)) and with (Figure 13(d) and (e)) the regions covered by volcanics, most of which are Miocene or younger, taken from Thorpe & Smith (1974), Ball et al. (2019), and earthwise.bgs.ac.uk.

Figure 15 :
Figure15: Spectral characteristics as in Figure8, but at the surface after filtering through the MBL and an elastic layer with thickness T e = 30 km.Notice the significant decrease in the power of gravity anomalies at short wavelengths and a moderate increase in the power at intermediate wavelengths.

Figure 16 :
Figure 16: Power spectral density of gravity anomalies at the top of the convecting box as in Figure 8(a), but showing the variation with Rayleigh number for the free-rigid simulations.

Figure
Figure S1: (a) Gravity field for the Pacific from DIR-R5, with coefficients l = 2 set to 0 and a filter applied, falling to 1/2 at 250 km, to remove the short wavelength components.(b) Residual depths, averaged over 2 • × 2 • boxes (Crosby et al. 2006).The dots show the locations of the resulting averages.Oblique Mercator projection with axis 40 • N, −50 • E.

Figure S3 :
Figure S3: Boxes used to estimate the elastic thickness T e of different parts of the Middle East, superimposed on gravity anomalies from Eigen6c (Förste et al. 2011),with the coefficients from l = 2 to 7 to 0, and applying a taper f, = (l − 7)/5, to those from l = 8 to 11.A low pass filter falling to 1/2 at 50 km was applied to the coefficients to remove the short wavelength anomalies from uncompensated topography.

Figure S4 :
FigureS4: Estimate of the elastic thickness for S and W Arabia (see FigureS3) from the admittance, taking the topography as input, gravity from Eigen6c as output.

Figure S5 :
FigureS5: Estimate of the elastic thickness for E Turkey and NW Iran (see FigureS3) from the admittance, taking the topography as input, gravity from Eigen6c as output.
Estimate of the elastic thickness for Anatolia (see Figure