Hotspots and mantle plumes revisited: Towards reconciling the mantle heat transfer discrepancy

Abstract Mantle convection is the principal mechanism by which heat is transferred from the deep Earth to the surface. Cold subducting slabs sink into the mantle and steadily warm, whilst upwelling plumes carry heat to the base of lithospheric plates where it can subsequently escape by conduction. Accurate estimation of the total heat carried by these plumes is important for understanding geodynamic processes and Earth's thermal budget. Existing estimates, based upon swell geometries and velocities of overriding plates, yield a global heat flux of ∼2 TW and indicate that plumes play only a minor role in heat transfer. Here, we revisit the Icelandic and Hawaiian plumes to show that their individual flux estimates are likely to be incorrect due to the assumption that buoyancy is mainly produced within the lithosphere and therefore translates at plate velocities. We develop an alternative methodology that depends upon swell volume, is independent of plate velocities, and allows both for decay of buoyancy through time and for differential motion between asthenospheric buoyancy and the overlying plate. Reanalysis of the Icelandic and Hawaiian swells yields buoyancy fluxes of 4.0 ± 0.5 Mg s−1 and 2.9 ± 0.6 Mg s−1, respectively. Both swells are used to calibrate a buoyancy decay timescale of ∼45 Myr for the new volumetric approach, which enables buoyancy fluxes to be estimated for a global inventory of 53 swells. Estimates from magmatic hotspots yield a cumulative lower bound on global plume flux of 2 TW, which increases to 6 TW if amagmatic swells are also included and if all buoyancy is assumed to be thermal in origin. Our results suggest that upwelling plumes play a significant role in the transfer of heat into the uppermost mantle.

The mantle acts as a giant heat engine, with various sources, sinks, and transport mechanisms (Figure 1). Using measurements of the conductive geothermal gradient in continental and oceanic lithosphere, the amount of heat currently escaping through the Earth's surface is estimated to be 43±3 TW (Jaupart et al., 2007). The magnitudes of the principal heat sources at the present day are less well constrained, but include heat supplied through the core-mantle boundary, decay of radiogenic nuclides (typically estimated by assuming the Earth has a bulk chondritic composition), and residual mantle heat from formation of the Earth and gravitational segregation of the core (Lay et al., 2008).
Transport of this heat to the Earth's surface occurs by vigorous convection and is the principal mechanism that drives mantle dynamics and plate tectonics (Turcotte & Schubert, 2002). Midoceanic ridge spreading causes passive upwelling of mantle material to the surface to form oceanic lithosphere, whilst small-scale convection at the base of the plates facilitates heat exchange with the underlying asthenosphere (Ballmer, 2017). Cold, dense lithospheric slabs sink into the mantle to different depths where they steadily warm, whilst active upwellings are buoyant and advect heat into the shallow mantle. The relative magnitudes and total contained heat flux of these transport mechanisms have significant ramifications for geodynamical problems. In the case of upwelling plumes, for example, their heat flux has implications for the generation, supply timescales, and geochemistry of melts at volcanic hotspots, the efficiency of mantle convection, core-mantle boundary heat flux, growth of the inner core, and the history of geodynamo activity (e.g. Vidal & Bonneville, 2004;Bourdon et al., 2006;Yoshida & Ogawa, 2005;Lay et al., 2008;Nimmo, 2015;Olson, 2016). For these reasons, considerable research effort has been expended on attempting to constrain the magnitude of heat supplied to the surface by active upwellings.
The morphology of plume swells results from interaction between upwelling convection currents and overlying rigid lithospheric plates. Therefore, these topographic and bathymetric anomalies are often used to constrain plume heat flux and to compare the relative strengths of plumes. Davies (1988) suggested a calculation based on the rate at which new surface topography (i.e. buoyancy) is generated, assuming swell topography is supported by simple thermal isostasy. Pioneering studies developed standardised methodologies for estimating buoyancy flux, with a view to producing a global inventory of individual plumes (Sleep, 1990). Associated cumulative estimates of global plume heat flux are 2.0 ± 0.3 TW, which represents a small fraction of heat escaping from the upper mantle and implies that upwelling plumes play only a minor role in heat transfer (Hill et al., 1992;Turcotte & Schubert, 2002). Over the intervening 30 years, however, many geochemical and geophysical studies have shed new light on the dynamics of plume swell formation that call into question some of the early methodological assumptions. Existing estimates of buoyancy flux are therefore worth revisiting, starting with one of the most well-studied examples beneath Iceland.

The Icelandic plume
A notable feature of global oceanic residual topography is that the North Atlantic region is anomalously elevated by up to 2 km over an area with diameter ∼ 2500 km (Figure 2a; . Iceland lies at the center of this swell, and it is generally agreed that the pattern of gravity anomalies, sub-plate shear-wave speeds, excess crustal thickness, and active volcanism are consistent with the presence of a major upwelling mantle plume (Vogt, 1971;Davies, 1988;Sleep, 1990;Ito, 2001;Rickers et al., 2013).
The Icelandic plume has had a significant influence on the stratigraphic, crustal, and lithospheric architecture of fringing continental margins throughout the Cenozoic Era (Howell et al., 2014). Arrival of the plume head is associated with Late Cretaceous opening of the North Atlantic Ocean and emplacement of considerable volumes of volcanic rocks from Baffin Island to Lundy (White & McKenzie, 1989). Rare Earth element geochemical signatures suggest initial excess temperatures of 100-200 • C, peaking at ∼ 56 Ma (White & McKenzie, 1989). Continued upwelling and enhanced mid-oceanic ridge melting throughout Cenozoic times has generated thickened oceanic crust along the Greenland-Iceland and Iceland-Scotland rises. Variations in the extent of smooth (i.e. plume-influenced) oceanic basement versus highly fractured (i.e. cooler) crust along and on either side of the Reykjanes Ridge imply that the plume-head radius has evolved with time (Figure 2b and 2d; Poore et al., 2009). The estimated radius is > 1000 km at break-up time, rapidly shrinking to < 400 km at ∼ 42 Ma before steadily growing again to ∼ 1000 km at the present day. Complementary cyrstallisation temperature estimates from olivine-spinel aluminium exchange thermometry on basalt samples within the Northern Rift Zone suggest that excess plume temperatures of ∆T ≈ 160 • C occur at the present-day conduit (Matthews et al., 2016). : Schematic cartoon of mantle heat budget and potential role of actively upwelling plumes. Blue boxes = heat flux lost through Earth's surface; green boxes = mechanisms for transporting heat through mantle; red boxes = sources of heat; black dashed line = 660 km discontinuity. Total heat absorbed from re-warming subducted lithospheric slabs is relatively well constrained but locations where it occurs within mantle depends on slab penetration depth, on sinking velocity, and on time taken to conductively warm slab interiors. Other flux estimates taken from Jaupart et al. (2007); Lay et al. (2008); Nimmo (2015); and . Cartoon not to scale.
Transient temperature anomalies within the plume head are manifest by diachronous V-shaped ridges visible in seafloor bathymetry and in free-air gravity anomalies southwest of Iceland ( Figure 2b). These nested chevrons of thicker crust converge along the Reykjanes Ridge and are thought to track the outward flow of hot mantle ripples that spread out from the conduit within an asthenospheric channel (Vogt, 1971). Hotter temperatures within the core of each ripple generate increased melt fraction beneath the oceanic spreading center (Ito, 2001). This asthenospheric flow is inferred to be at least semi-radial based upon the existence of episodic phases of clastic deposition in surrounding basins and of transient uplift events on fringing continental margins (White & Lovell, 1997). An Eocene ephemeral landscape identified in the Faroe-Shetland Trough records a particularly striking example of this perturbation, undergoing > 0.5 km of transient vertical motion on timescales shorter than 3 Myr (Hartley et al., 2011). Neogene vertical motions have also modified the pattern of oceanic circulation, with periodic uplift of the Denmark Straits and Faroe-Shetland Trough shutting off overflow of Northern Component Water from the Arctic Ocean and reducing sedimentation rates in contourite drifts (Poore et al., 2009;Parnell-Turner et al., 2015). Full-waveform tomographic studies suggest that the plume head has a highly irregular planform with fingers of slow shear wave velocity material protruding beneath Scotland, southern Norway and Greenland (Rickers et al., 2013). Fluid dynamical experiments indicate that the fingering may be a large-scale manifestation of the Saffman-Taylor instability -behaviour that strongly depends upon a combination of asthenospheric thickness, the viscosity contrast between plume-head material and ambient mantle, and the plume buoyancy flux (Schoonman et al., 2017).

Traditional on-axis estimates of Icelandic buoyancy flux
The standard buoyancy flux approach for plumes located beneath mid-oceanic spreading centres was developed and applied to the Icelandic plume by Sleep (1990). Volume flux, Q V , is gauged by assuming that material carried upward by the plume is in steady state with that being removed by seafloor spreading. Removal occurs via diverging lithospheric plates overlying a low-viscosity asthenospheric channel, within which Couette flow dominates such that lateral velocity drops linearly from the plate half-spreading rate at the top to zero at the base of the channel (Figure 3a). Thus where z l is lithospheric thickness, za is thickness of the asthenospheric channel, zy is the along-axis ridge length supplied with excess heat, and vs is the full-spreading rate (see Table 1). A key assumption is that Couette flow dominates within the channel, which means that the average asthenospheric velocity is one   corrected for age-depth subsidence using RHCW18 plate model . Circles = measurements including isostatic correction for crustal thickness; upwards/downwards pointing triangles = minimum/maximum constraints; grid = residual bathymetry calculated from ship-track measurements; red/black/blue contours = positive/zero/negative values of GGM03C free-air gravity anomalies spaced every 10 mGal and band-pass filtered 9000 > λ > 730 km (Tapley et al., 2007). (b) V-shaped ridges straddling Reykjanes mid-oceanic spreading center, visible within free-air gravity anomalies that are highpass filtered using cosine taper over range of 100-400 km (Sandwell et al., 2014). (c) Swell planform from degree-30 spherical harmonic fit to oceanic residual bathymetric measurements with degree 0-3 components removed; dashed line = approximate extent of plume swell based upon along-Reykjanes Ridge extent of smooth-to-rough basement transition. (d) Interpretation of crustal structure based on short-wavelength gravity anomalies (pale background shading); black tramline = mid-oceanic spreading center; long dashed lines = V-shaped ridges; short dashed lines = prominent fracture zones; green line = boundary between smooth and rough basement.
half that of the overriding plates. Lithospheric and asthenospheric channel thicknesses of z l = za = 100 km and an along-axis ridge length of zy = 800 km were used, based upon bathymetric gradients and the extent of axial geochemical anomalies. The presentday full spreading rate was estimated as vs = 16.5 mm yr −1 , yielding Q V = 62.7 m 3 s −1 . Volume flux can subsequently be converted into buoyancy flux, Q B . The simplest approach is to estimate the mantle's thermal expansion using inferences of the mean excess temperature of plumehead material, ∆T , such that where ρm = 3300 kg m −3 is the density of mantle at the surface and α = 3 × 10 −5 K −1 is the thermal expansivity. ∆T can be estimated from independent constraints such as crustal thickness produced by melting or major/trace element geochemistry. Sleep (1990) adopted ∆T = 225 • C for Iceland, based upon melt thickness estimates, yielding Q B = 1.40 Mg s −1 . More recent estimates from Matthews et al. (2016) revise the excess temperature down to ∆T ≈ 160 • C. In addition, a modern, regional full-waveform tomographic model suggests that slow shear-wave speeds associated with the plume head occur within a ∼ 125 ± 25 km thick layer directly beneath the plates (Rickers et al., 2013). Updating these values of ∆T and za yields Q B = 1.08 Mg s −1 .  (Sleep, 1990); vs = full spreading rate; z l = lithospheric thickness; za = asthenospheric channel thickness; zy = along-axis length of ridge supplied with plume material; blue/black arrows = plate motion; red arrow labelled Q V = plume volume flux. (b) Plate velocity-based approach for plume beneath single plate (Sleep, 1990); vp = plate velocity; A = cross-sectional area of elevated region, perpendicular to direction of plate motion; red arrow labelled Q B = plume buoyancy flux. (c) V-shaped ridge approach for Icelandic plume; grey = crust; red = plume; dashed black line = smooth-to-rough basement transition; Rmax = finger propagation radius; r• = interstitial plume-head radius; ∆t = duration of V-shaped ridge propagation (from magnetic reversal history); za = asthenospheric channel thickness; Q V = plume volume flux. (d) Volumetric approach; V = swell volume; Q B = plume buoyancy flux. Cartoon not to scale.

Revised estimates of Icelandic buoyancy flux
One key assumption of the spreading-rate calculation is that plume-head material horizontally translates at an average of half the velocity of the overriding plate. However, evidence in the North Atlantic Ocean from prominent V-shaped ridges, ephemeral landscapes and off-axis uplift of oceanic gateways suggests that this assumption is probably incorrect. The angle and curvature of each "V" with respect to the Reykjanes Ridge is controlled by interaction between spreading rate and the along-axis velocity of each hot ripple as it propagates away from the plume conduit. Spreading rates are independently known from magnetic reversal histories. Therefore, along-axis asthenospheric velocities can be shown to vary between 87 and 282 mm yr −1 , which is an order of magnitude faster than the full plate-spreading rate and occurs at a highly oblique angle (Poore et al., 2009). Additional evidence for differential motion between the plate and asthenosphere comes from the rapid reburial rates of ephemeral landscapes, which are too fast for conductive decay of a thermal anomaly. Instead, horizontal advection within an asthenospheric channel at similarly fast velocities of ∼ 350 mm yr −1 is required (Hartley et al., 2011). This mode of flow, whereby material in the sub-plate channel travels faster than the overlying lithosphere, is known as Poiseuille flow and is driven by lateral pressure gradients. Numerical experiments of thermal convection beneath moving plates also exhibit this phenomenon. Yoshida & Ogawa (2005) found in their two-dimensional models that heads of hot upwelling plumes spread out laterally an order of magnitude faster than overlying plate velocities. Measuring the heat flow using plate velocitybased methods accounted for less than 10% of the real plume heat flux for all of their experiments. Thus, inferring asthenospheric velocities to always be slower than plate spreading results in substantial underestimation of buoyancy flux for the Icelandic plume. We therefore require an alternative approach that does not assume Couette flow within the plume head.
One option is to exploit the architecture of V-shaped ridges in combination with the distance of the transition from smooth to rough crustal fabric along the Reykjanes spreading center (Figure 3c). The pattern of volcanism, gravity anomalies, and slow sub-plate shear wave speeds has recently been used to infer that the plume head is highly irregular, with several fingers protruding away from Iceland within a sub-plate asthenospheric channel, including one aligned with the spreading ridge (Schoonman et al., 2017). The down-axis smooth-to-rough transition marks the maximum extent, Rmax, of this finger. The relationship between each "V" and the magnetic reversal history allows us to extract the time taken, ∆t, for pulses of hot material to spread down the length of this finger to the outer margin of the plume head. Adopting a simple geometric representation of the outer extent of the fingering plume head in Figure 3c, it can be shown that the volume flux can be approximated by where r• is the minimum radius of the plume head in the interstices between fingers (H. Galbraith-Olive, personal communication, 2019). The present-day smooth-to-rough transition occurs at a distance of 1200 ± 100 km down the ridge axis, and the "V" that reaches this point left the plume conduit at ∼ 15 ± 1 Ma (Parnell-Turner et al., 2014). Adopting za = 125 ± 25 km from the Rickers et al. (2013) tomography and assuming that the minimum radius is r• = 200 ± 50 km, Equation (3) yields a geometrically derived Qv = 16.1 ± 4.2 km 3 yr −1 for Iceland, where uncertainties in za, Rmax, r•, and ∆T are assumed to be uncorrelated and have been propagated according to the variance formula. Notably, this volume flux is a factor of two to three smaller than values that have been obtained by assuming that hot ripples within the plume head spread within a circular disk (e.g. Poore et al., 2009;Parnell-Turner et al., 2014). On the other hand, it is a factor of two to three larger than numerical simulations of the plume, which are tuned to fit first-order hotspot characteristics, such as melt generation rates inferred from crustal thickness measurements (e.g. Ito, 2001;Howell et al., 2014). If the finger aligned with the Reykjanes Ridge is longer than other fingers, or if the ripple duration actually records the travel time for solitary waves that traverse faster than the kinematic flow rate of hot material, then this new geometric volume flux represents an upper bound. However, there are also significant uncertainties in models of plume dynamics arising from poorly constrained knowledge of the appropriate mantle rheology, particularly when variable source composition, generation of large quantities of melt, its extraction to the surface, and any associated effects on the residuum must be included. For example, numerical models have thus far been unable to reproduce the along-axis velocities required by V-shaped ridges (Ito, 2001). Reconciling these issues is likely to be a fruitful endeavour for further research, but given the wide scope of model uncertainties associated with the rheological effects of significant melt removal, we favour the geometrically derived volume flux determined here. Adopting ∆T ≈ 160 • C from Matthews et al. (2016) and assuming that the average excess temperature of plume material rising up the conduit may be approximately half of this value, Equation (2) yields Q B = 4.0 ± 1.0 Mg s −1 for Iceland. An alternative method, that is also independent of assumptions regarding plate-spreading rates, is to exploit the total volume of the plume swell. The volume and the duration over which it has grown provide information on the flux of material supplied by the plume. Swell volume is a balance between enlargement, caused by the input of new buoyancy at the conduit, and reduction due to buoyancy loss, such as cooling of hot material within the plume head (Crosby & McKenzie, 2009). This heat loss occurs by transport of melt products to the surface and conduction both into the surrounding ambient mantle and through the lithosphere, which is aided by small-scale convection and thermal erosion at the base of the plate (Ballmer, 2017). These processes all occur on characteristic timescales.
Assuming that the swell is in approximate steady-state, the buoyancy flux of upwelling material will be balanced by the rate of buoyancy loss. In this case, the buoyancy flux can be estimated using a volumetric approach given by where ρa = 3200 kg m −3 is the density of asthenospheric mantle, ρ i is the density of displaced surface fluid (i.e. air or water), V is the volume of the swell, and τ is a characteristic timescale of buoyancy loss (Crosby & McKenzie, 2009). This approach implies that if a plume were to suddenly cease upwelling, topographic effects at the surface would decay away over this time period. Therefore, during this time, all buoyancy must be resupplied by the plume to maintain steady-state. An important corollary of using this volumetric method is that individual plumes younger than τ or increasing in magnitude through time will yield underestimates of present-day buoyancy flux, whilst those that are waning or dying out will provide overestimates. The steady-state approximation is likely to be least applicable during the earliest phases of a plume's lifetime when buoyant plume-head material first impacts the base of the plate, often generating large igneous provinces (Hill et al., 1992;White, 1993). The volumetric approach can be applied to estimate buoyancy flux of the Icelandic plume. Here, the swell geometry is obtained using the oceanic residual topography database of , which has been corrected for the effects of sedimentary loading, crustal thickness variations, and age-depth cooling using the RHCW18 plate model . The longest wavelengths of residual topography are primarily supported by flow in the deeper mantle rather than asthenospheric isostasy . We therefore fit spherical harmonic basis functions to this database up to and including degree 30 (wavelengths approximately ≥ 1300 km) and subtract the degree 0-3 components from the field (Figure 2c). Defining the extent of the plume swell is subjective, but the contour that coincides with the present-day extent of V-shaped ridges down the Reykjanes Ridge provides a useful guide and yields a water-loaded swell volume of V = 2.6 × 10 6 km 3 . We note that the limited spatial resolution of the residual topography, particularly following spherical harmonic filtering, does not currently fully resolve the individual plume-head fingers visible in the seismic tomography and gravity datasets (Schoonman et al., 2017). Future studies with denser data coverage may start to recover these features and allow for better determination of the swell extent. The possible additional hotspot  . located at the junction between the Kolbeinsey and Mohns ridges beneath the volcanically active island of Jan Mayen is included within this swell.
The geometrically derived buoyancy flux of Q B = 4.0 Mg s −1 obtained from Equation (3) and V-shaped ridge propagation rates can be used to calibrate the buoyancy decay timescale, yielding τ = 45 Myr ( Figure 4a). For this timescale, systematic ±100 m errors within the swell topography grid give a range of Q B =3.6-4.4 Mg s −1 , whilst adopting τ =45 +10 −5 Myr for the buoyancy decay time gives Q B =3.3-4.5 Mg s −1 . In order to recover the plate spreading rate-based value of 1.1 Mg s −1 from Section 1.1, the decay timescale would need to be ∼ 165 Myr.

The Hawaiian plume
The Hawaiian plume located in the central Pacific Ocean is a widely studied mantle upwelling. Since at least ∼ 76 Ma, this hotspot has produced basaltic volcanism recorded by an ageprogressive chain of seamounts stretching southeastward across the Pacific plate (Clouard & Bonneville, 2005;Jicha et al., 2018;Figure 5a). Modern basalts at Kilauea are erupted on top of Cretaceous seafloor and the present-day planform is visible in both residual depth measurements, which peak at ∼ 1 km, and longwavelength free-air gravity anomalies (Figure 5b). Sleep (1990) described a method for constraining buoyancy flux for plumes located away from plate margins by measuring the rate at which new swell topography is generated and assuming that it is supported by thermal isostasy. The effects of advective stresses are neglected, which is probably valid if upward flow is concentrated within a plume conduit that is much narrower than the width of plume-head material within the asthenosphere. If excess buoyancy is generated within the lithosphere, it moves downstream at the velocity of the overriding plate. The buoyancy flux can then  Clouard & Bonneville (2005). Dark grey polygons = land above sea level. (b) Oceanic residual depth anomalies from  corrected for age-depth subsidence using RHCW18 plate model . Circles = measurements including isostatic correction for crustal thickness; upwards/downwards pointing triangles = minimum/maximum constraints; grid = ship-track residual topography; red/black/blue contours = positive/zero/negative values of GGM03C free-air gravity anomalies spaced every 10 mGal and band-pass filtered 9000 > λ > 730 km (Tapley et al., 2007). (c) Swell planform from degree-30 spherical harmonic fit to oceanic residual depth anomalies with degree 0-3 components removed; dashed line = approximate swell extent predominantly following zero contour; circles = seamount ages in Ma.

Traditional off-axis estimates of Hawaiian buoyancy flux
be calculated by combining plate velocity with the cross-sectional area of the swell above the conduit, measured perpendicular to the direction of plate motion (Figure 3b). Thus where A is the cross-sectional area perpendicular to plate motion and vp is the plate velocity. Critically, buoyancy flux obtained by using this off-axis approach does not require an estimate of the plume excess temperature, ∆T . For plumes beneath fast-moving plates, Sleep (1990) predicts that the majority of excess topography downstream of the conduit is generated by thermal isostasy within the lithosphere, and that only a small fraction of the swell is underlain by hot asthenosphere. Hawaii is considered to be a typical example of this situation. The cross-sectional area of the Hawaiian swell was estimated as A ≈ 1430 km 2 and local velocity of the Pacific plate as vp = 83 mm yr −1 . Assuming ρm = 3300 kg m −3 and ρ i = 1000 kg m −3 , a buoyancy flux of Q B = 8.7 Mg s −1 was calculated. This plume flux is the largest within Sleep's database and is six times the value for Iceland. Notably, despite a significant increase in the quantity and quality of residual topography observations around Hawaii over the last three decades, the location of the swell nose and its general shape have remained broadly unchanged . Updating the plate velocity to ∼ 87 mm yr −1 from Jicha et al. (2018) increases the buoyancy flux estimate to 9.1 Mg s −1 . A significant assumption behind this method is that buoyancy is generated within the lithosphere and so travels at plate velocities. This value is therefore potentially erroneous if the swell is supported by asthenospheric material that is not perfectly coupled to the plate.

Revised estimates of Hawaiian buoyancy flux
In the absence of a mid-oceanic spreading ridge that independently tracks the velocity of asthenospheric material at Hawaii, we must turn instead to numerical simulations of the plume. The swell morphology has been the focus of several studies of plume-head dynamics (e.g. Watson & McKenzie, 1991;Ballmer et al., 2011). For example, Ribe & Christensen (1994) exploited a three-dimensional Cartesian box with temperature-and pressure-dependent viscosity. A layer of cold, high-viscosity lithosphere is imposed at the top with a fixed lateral velocity. A thermal anomaly on the base of the box represents the plume conduit, whose size and temperature excess (i.e. buoyancy flux) can be varied to optimise the fit between the observed and predicted swell. Their results suggest that there has been only minor thermal erosion of lithosphere, and that the swell is mainly supported by buoyancy variations within the asthenosphere. These factors are corroborated by analysis of the Hawaiian geoid-to-topography ratio and downstream rates of seamount subsidence (Cadio et al., 2012;Huppert et al., 2020). Updating the model to include effects of partial melting and associated depletion of plume-head material suggests a present-day thermal buoyancy flux of Q B =2.2-3.5 Mg s −1 . For comparison, the model of Ballmer et al. (2011) also successfully matches the firstorder characteristics of the swell geometry and melt production observed at Hawaii. They fixed Q B = 4.0 Mg s −1 and tested the effect of a significantly more temperature-dependent rheology for the mantle material. This rheology leads to the development of vigorous small-scale convection, which initially causes lithospheric erosion of up to ∼ 15 km, resulting in an increase in swell topography (i.e. dynamic rejuvenation) and affecting the geoid-totopography ratio downstream of the conduit in a manner that is compatible with the observations (Cadio et al., 2012). Small-scale convection subsequently leads to more rapid decay of the swell as heat escapes more efficiently into the surrounding ambient mantle and by conduction through this mildly thinned lithosphere.
Notably, these numerical model values are a factor of two to four smaller than buoyancy flux derived from the traditional platevelocity approach. This difference arises because a significant component of the buoyancy provided by hot plume-head material is generated within the asthenosphere that, as in the Icelandic case, can have strong differential motions with respect to the overlying Pacific lithosphere. For the model of , Couette-style flow dominates and the average downstream velocity of plume-head material is only ∼ 65% of the plate velocity. Conversely, Ballmer et al. (2011) find that the plume head spreads faster than the plate in a Poiseuille-type flow, with differential velocities initially 200 mm yr −1 greater in the vicinity of the conduit, reducing to only ∼ 25 mm yr −1 (∼ 25%) faster than the plate at a distance of 900 km downstream from the conduit. This contrast emphasises the uncertain, yet crucial, role of rheology in controlling plume-head dynamics.
We can use the results of these numerical models to attempt a calibration of the volumetric buoyancy flux approach for Hawaii. The total water-loaded volume of excess material contained within the swell in Figure 5c is 1.79 × 10 6 km 3 , which is 10% larger than the value of 1.60 × 10 6 km 3 obtained by Vidal & Bonneville (2004) using the PS77 plate model. Application of Equation (4) yields a range of buoyancy fluxes that depends upon the value of the decay timescale. Calibration using Q B = 2.9 Mg s −1 from the model of  suggests that τ = 43 Myr (grey bars in Figure 4b). Varying τ =43 +10 −5 Myr gives a buoyancy flux range of 2.4-3.3 Mg s −1 , whilst systematic swell topography errors of ±100 m with τ = 43 Myr gives a range of 2.2-3.7 Mg s −1 . Within the confines of the steady-state assumption, this result suggests that the processes of buoyancy loss from the Hawaiian swell appear to occur over timescales of ∼ 45 Myr.
Using seamount ages as a proxy for conduit location allows reconstruction of plate position through time and suggests that the majority of excess topography has indeed decayed over 45-55 Myr since the lithosphere overlay the plume conduit ( Figure 5c). This value agrees with the swell decay observations of Vidal & Bonneville (2004) and Crosby & McKenzie (2009), but is shorter than the ∼ 65 Myr obtained from linearly extrapolating the decay of uplift in the numerical model of , which is likely to be an upper bound reflecting the absence of smallscale convection in their plume head (Supplementary Material). In the case where plume material travels slower than the plate (Couette flow), this seamount-derived age will be an underestimate of the buoyancy decay timescale, whereas when it travels faster (Poiseuille flow), the asthenosphere supporting downstream topography will be comparatively 'younger' than the volcanism of the overlying plate, making 45-55 Myr an upper bound. This latter case, with a significant contribution from Poiseuille flow, is supported both by the most recent rheological modelling of the Hawaiian plume and by a Pacific-wide analysis of the forces driving plate motion over the last 15 Ma (Ballmer et al., 2011;Stotz et al., 2018).
A volumetrically derived Q B = 2.9 Mg s −1 for Hawaii is one third of the buoyancy flux of Q B = 9.1 Mg s −1 estimated using the plate velocity-based approach. For the volumetric approach to agree with this larger value, the buoyancy decay timescale would have to be ∼ 11 Myr, indicating exceptionally rapid loss of excess buoyancy and heat. An alternative explanation is that the volumetric method underestimates the true buoyancy flux because the Hawaiian plume is increasing in strength through time. The steady-state assumption would therefore be invalid. One possible test of this hypothesis comes from calculations of the melting rate through time estimated from the volume of igneous material contained within seamounts. Intrusive melt products contained within crustal roots must be included, requiring estimation of the densities and elastic parameters associated with loading of the lithospheric plate. Analysis of the Hawaiian-Emperor chain indicates that longterm melt productivity has probably increased by a factor of two to three over the last 70 Ma (White, 1993;Vidal & Bonneville, 2004;Crosby & McKenzie, 2009). This observation provides compelling evidence for an increase in plume flux, although caveats exist when linking melt productivity to buoyancy flux due to the requirement for independent constraints on temporal variations of excess temperature, melting depth range and source fusibility.
A second, oft-cited piece of observational evidence in support of a general increase in buoyancy flux at Hawaii since 30 Ma comes from direct analysis of the swell geometry (e.g. Davies, 1992;Vidal & Bonneville, 2004;Crosby & McKenzie, 2009). These studies reconstruct the buoyancy flux at earlier times using the plate velocity-based approach, taking cross-sectional areas at multiple positions along the swell axis and using the plate velocity to infer the age at which each location originally overlay the conduit. These calculations therefore assume that excess buoyancy travels at the velocity of the lithospheric plate and that it does not decay through time. Given that numerical studies of Hawaiian plumehead dynamics appear to invalidate these two assumptions, this line of evidence for increasing plume flux through time is question- able. In particular, if excess buoyancy does diminish with time, downstream parts of the swell will have had longer to decay and their cross-sectional areas will reduce. The analysis will therefore produce an apparent increase in Hawaiian buoyancy flux through time that may be erroneous. Nevertheless, if the buoyancy flux at Hawaii (or Iceland) has indeed steadily increased, as suggested by Hawaiian melt production arguments, then τ ≈ 45 Myr likely represents a lower bound on the buoyancy decay timescale.

Interpretation of the buoyancy decay timescale, τ
Analysis of the Icelandic and Hawaiian swells suggests that an effective buoyancy decay time of ∼ 45 Myr provides results that are consistent with the available independent constraints. However, physical interpretation of this value is complicated because a combination of different processes may be occurring that operate on a range of decay timescales.
The largest buoyancy source is likely to be thermal in origin, and its decay depends upon the rate at which excess heat is lost by conduction from hot plume-head material into the surrounding ambient mantle and through the overlying lithospheric plate. The thickness of a layer, z, that can lose heat by conduction over a given duration is given by where τc is the characteristic conductive timescale and κ = 1 mm 2 s −1 is the thermal diffusivity of silicate mantle rocks (Turcotte & Schubert, 2002). The buoyancy decay timescale of τ ≈ 45 Myr obtained for Iceland and Hawaii yields a characteristic conductive layer thickness of 120 km, whilst a range of τ = 30-65 Myr givens a thickness range of 100-140 km. This estimate raises the question -what physical feature might such an effective thickness relate to? There are two potential options that are most apparent. First, oceanic lithosphere reaches a maximum thickness of ∼ 120 km at older ages, similar to that in non-cratonic continental regions (e.g. . Therefore, plumes upwelling in the center of plates, such as Hawaii, will need to lose heat by conduction through a lithospheric lid of approximately this thickness, whilst plumes beneath a mid-oceanic spreading ridges, such as Iceland, will undergo corner flow, drawing material into a layer of this thickness that then cools conductively to become oceanic lithosphere. A second explanation is that the plume head itself may spread out beneath the plate within a ∼ 120 km thick asthenospheric layer, resulting in conductive cooling on ∼ 45 Myr timescales. It is important to also consider the role of small-scale convection in increasing the efficiency of heat loss. Numerical models of asthenospheric convection using rheologies with a realistic temperature dependence illustrate that this phenomenon is likely to be ubiquitous in low-viscosity regions close to thermal boundary layers, including at the top of the mantle (Ballmer, 2017). The model of the Hawaiian plume by Ballmer et al. (2011) suggests that small-scale convection increases the rate of heat loss in two ways. First, deformation inter-fingers plume-head material with colder lithosphere and ambient mantle, reducing the effective conductive lengthscale and shortening the required cooling time according to Equation (6). Secondly, it erodes the base of the plate, which thins the lithosphere and reduces the timescale for conductive heat loss to the surface. Evidence in support of this thermal erosion and general rejuvenation at the base of the lithosphere in some specific cases has been obtained from the analysis of observed geoid anomalies, including examples at Bermuda, Cape Verde and Hawaii (Sleep, 1990;Cadio et al., 2012).
In addition to thermal buoyancy supporting swell topography, another potentially significant buoyancy source results from compositional depletion of residual plume material during melting near the conduit. For example, melt generation rates at Hawaii have been estimated to be ∼ 5 m 3 s −1 based on seamount volumes (White, 1993;Vidal & Bonneville, 2004;Crosby & McKenzie, 2009). An axisymmetric numerical model of the plume suggests that upwelling material may undergo an average melt removal of 6.6% (Watson & McKenzie, 1991). Three-dimensional numerical experiments suggest that plume-head depletion could contribute up to a quarter of total buoyancy at Hawaii . Depletion is therefore potentially a significant buoyancy source, particularly in locations with thin lithosphere such as Iceland where plumes can upwell to shallow depths and undergo larger amounts of fractional melting (White, 1993). However, it should be noted that some studies have also suggested that a negative feedback exists that limits large quantities of decompression melting, whereby removal of melt and volatiles results in a viscosity increase in the residuum that prevents further upwelling to shallow depths (e.g. Ito, 2001). The decay timescale associated with loss of depletion buoyancy is likely to be significantly longer than ∼ 45 Myr, being dominantly controlled by lateral asthenospheric flow and mechanical mixing rates. This effect would lead to our volumetrically derived buoyancy fluxes for Iceland, and perhaps for Hawaii, representing upper bounds. Detailed analysis of the melt productivity of individual plumes is required to better constrain the relative contributions of thermal and depletion buoyancy. Later studies using plate velocity-based methodology have exploited revised bathymetric, topographic, sedimentary thickness and crustal age grids, updated plate reconstruction models, and improved methods for estimating swell cross-sectional areas and volumes (e.g. Vidal & Bonneville, 2004;Adam et al., 2005). The greatest uncertainties are usually encountered for the smallest swells, whilst buoyancy fluxes for larger plumes such as Hawaii are generally unchanged. King & Adam (2014) directly calculated buoyancy flux for 54 hotspots using the off-axis plate velocitybased approach. They compared different methods for estimating swell geometries, including the MiFil technique, which was designed to remove topographic contributions arising from volcanic crustal thickness variations and flexural features (Adam et al., 2005). Their average global buoyancy flux is 47 ± 7 Mg s −1 , which is similar to both Q B = 41 Mg s −1 from Davies (1988) and Q B = 55 Mg s −1 from Sleep (1990).
The issue surrounding differential velocities between the asthenosphere and an overlying plate can be made most apparent by considering the contribution of each plate in isolation. Here, individual plumes are assigned to their overlying plate, with any that fall close to a plate boundary being divided equally between each plate (N.B. Azores and Bowie have been split into three). Summing the total contribution and dividing by plate area yields the areal concentration of plume buoyancy flux (Figure 6c).
Upwellings are initiated by buoyant instabilities generated within the lower thermal boundary layer of a convecting system (Turcotte & Schubert, 2002). Their locations and relative strengths might therefore be expected to be relatively independent of overlying conditions, particularly for plates larger than typical convective planforms, which are less biased by spatial location. However, Figure 6c shows that these normalised fluxes are not independent of average plate velocities. Larger flux concentrations tend to occur on faster moving plates. The Pacific plate, in particular, dominates the global budget since it has both the largest plate area and the fastest average velocity. Figure 7 demonstrates that Pacific buoyancy flux values systematically vary with plate velocity/spreading rate, with large fluxes for Hawaii, the South Pacific Superswell, Samoa and Caroline where plate velocities exceed 70 mm yr −1 . The on-axis Easter hotspot is also large and occurs where full plate spreading peaks at ∼ 150 mm yr −1 on the East Pacific Rise between the Pacific and Nazca plates. In contrast, Louisville, Baja, Bowie and Cobb have smaller buoyancy fluxes, partly because plate velocities decrease to both the north and south.
The general relationship between plate velocity and plume flux may represent an actual feedback between the volume of material supplied by plumes and the forces driving plate motion, but given the co-dependence of these variables in this methodology, drawing this conclusion is potentially problematic and casts doubt on the use of plate-velocity terms in calculating plume buoyancy flux (King & Adam, 2014). The flux beneath slow-moving plates is almost certainly under-represented. This issue can be taken to the extreme case of a plume upwelling beneath the center of a stationary plate -irrespective of the resultant swell size, calculated buoyancy flux will be zero due to the fact that vp = 0 mm yr −1 . Calculations that depend upon plate velocity break down if plumehead material has markedly different lateral velocities compared with the overriding plate. It is therefore desirable to attempt application of the volumetric method, which does not make use of plate velocities.

Revised swell volume-based estimates
As with the comprehensive analyses of Sleep (1990) and of King & Adam (2014), it is important to calculate volumetrically-derived buoyancy flux associated with a global suite of swells. To constrain swell volumes, we use the recent residual topography measurements of  that have been corrected for age-depth cooling using the RHCW18 subsidence relationship . Spherical harmonic basis functions are fitted to these measurements using the approach described by . The longest wavelength degree 0-3 components have been subtracted from the spherical harmonic representation, since these components are probably dominated by flow within the lower mantle. The maximum spherical harmonic degree included is l = 30, which corresponds to minimum wavelengths of ∼ 1300 km. This cut-off should help to filter out features associated with small-scale convection and lithospheric flexure. The resulting swell geometry contains a combination of buoyancy supplied either by plumes or by mid-scale upper mantle convection. The original grid is built from accurate oceanic residual depth measurements in the offshore, and more controversial continental constraints based on the linear scaling of free-air gravity anomalies. Buoyancy fluxes associated with oceanic swells are therefore likely to be more reliable.
Swell outlines have been digitised and total excess volumes calculated. Picking swell outlines is a subjective undertaking. Away from Iceland, we do not have the smooth-to-rough basement transition to demarcate swell extent. In general, we have therefore relied upon identification of the zero contour line. Where swells overlap or continue into probably unrelated features, approximate extents have been assigned (see Supplementary Information). As pointed out by Sleep (1990), some swells overlie more than one putative plume, and in these cases the composite swell has been equally split. In the Pacific Ocean, these swells include Bowie and Cobb, and Juan Fernandez and San Felix. Composite swells in the Atlantic Ocean include the Azores and Great Meteor, Ascension and St. Helena, Discovery and Tristan, Meteor and Bouvet, and Marion and Crozet on the Southwest Indian Ridge. As things stand, St. Paul, Eastern Australia, Lord Howe and the Tasmantid hotspots do not appear to be associated with swells in the global topography model, although better resolved regional studies may recover such features in the future. The total volume of excess topography for each swell has been computed and converted into buoyancy flux using Equation (4), assuming τ = 45 Myr in all cases (Table 2; Figure 8a). The resulting relationship between buoyancy flux concentration and plate velocity for these magmatic swells is less obvious than for the Sleep (1990) estimates, and most plates have values of 1-6 kg m −2 yr −1 (Figure 8c). Larger values of ∼ 10 kg m −2 yr −1 are obtained for the Arabian and Cocos plates, which is probably due to bias resulting from the small spatial extent of these plates with respect to the underlying convective planform. Table 2: Plume buoyancy flux estimates, Q B , derived using plate velocity-based and volumetric approaches. 48 magmatic swells are considered, including all 37 hotspots from Sleep (1990) and 47 of 54 considered by King & Adam (2014). 5 amagmatic African swells are also included. NS = no associated swell visible in global swell topography map.
b Plate velocity-based values using MiFil volume from King & Adam (2014).
c Volumetric-based values of this study using τ = 45 Myr.

Global plume flux
Previous studies have estimated a global buoyancy flux of 40-55 Mg s −1 by summing up contributions from individual magmatic swells (Davies, 1988;Sleep, 1990;King & Adam, 2014). By adopting this approach, we obtain a similar value of Q B = 46.2 Mg s −1 for the volumetric values, where τ = 40-55 Myr gives a range 38-52 Mg s −1 . However, it is important to note that the process of picking swell outlines is subjective and that identification based on volcanic activity alone results in amagmatic swells being overlooked. For example, in addition to the continental magmatic hotspots of Afar, Hoggar, and Tibesti, several prominent anorogenic swells occur throughout Africa that have also been linked to the underlying planform of mantle convection (e.g. Burke & Gunnell, 2008;. These features include the Atlas, Fouta Djallon, Angolan, Namibian, and South African domes, which each have a buoyancy flux of ∼ 1 Mg s −1 using the volumetric approach, but have no volcanic activity since thicker lithosphere inhibits decompression melting at shallow depths (Table 2). In the Pacific Ocean, where the global swell topography map is based on accurate residual depth measurements, there is evidence for additional swells with wavelengths of ∼ 1000 km that are not associated with known hotspots. The total volume of all amagmatic features is double that of the magmatic swells, and estimating the buoyancy flux associated with total excess topography gives Q B ≈ 150 ± 25 Mg s −1 for τ = 45 +10 −5 Myr (Figure 9). Whether the value of τ obtained by calibration at Iceland and Hawaii is globally applicable will require further analysis. However, the assumption of a steady-state planform, which is inherent to the volumetric approach, is probably more realistic from a global perspective. Individual plumes may be forming or dying out over the chosen time window, leading to under-and over-predictions of buoyancy flux, respectively. On longer timescales, however, the mantle convects in a quasi-steady state such that despite individual plumes evolving, the general morphology is more likely to be consistent.
The excess swell topography grid is built from a combination of residual depth measurements from the oceanic realm and scaled long wavelength free-air gravity anomalies onshore. It is therefore worth emphasising that the offshore buoyancy flux contribution, which covers approximately two-thirds of the Earth's surface, represents 59% of the total amount. These values suggest that there is limited onshore versus offshore bias arising from the observational constraints and that swell locations are probably independent of the overlying type of lithosphere.

Mantle heat budget and the role of upwelling plumes
Investigating the implications of these buoyancy flux estimates requires an understanding of the mantle heat budget (Figure 1). The total conductive heatflow escaping through continental regions is thought to be ∼ 14±1 TW (average of 65 mW m −2 , based on thousands of borehole measurements). Approximately 6.5 ± 0.5 TW of this value is generated by decay of radionuclides within the crust and a further 0.5 ± 0.5 TW is generated within the lithospheric mantle (Jaupart et al., 2007). Thus, the remaining 7 ± 1 TW must be basal heat supply from the convecting mantle, which is thought to be dominated by contributions from tectonically active regions and from passive margins (Jaupart et al., 2007). In contrast, oceanic crust contains only minor quantities of radioactive elements and so internal heat production is insignificant. In the oceanic realm, the principal mechanism of heat loss is passive upwelling of hot asthenospheric mantle beneath mid-oceanic spreading centres and its subsequent conductive cooling to produce lithospheric plates. Additional basal heat is supplied by the upper mantle through a mixture of conduction and small-scale convection, yielding a steady state geotherm and constant bathymetry in older oceanic seafloor (Turcotte & Schubert, 2002;Ballmer, 2017). Directly measured oceanic heatflow tends to have large scatter due to the combined effects of hydrothermal circulation and measurement bias (Jaupart et al., 2007). Observations are particularly unreliable for young oceanic crust, which unfortunately has the most significant heatflow contribution. To sidestep this issue, models of the thermal structure of oceanic lithosphere are fitted to reliable heatflow and subsidence data from older regions and subsequently extrapolated into young lithosphere. These extrapolations are then weighted for the age-distribution of oceanic lithosphere in order to calculate total heatflow through ocean floor. A recent analysis yields a total surface flux of 29 ± 2 TW, ∼ 5 TW of which is basal resupply beneath older lithosphere . Combining these oceanic estimates with those from continents yields a total heat flux through the Earth's surface of ∼ 43 ± 3 TW (an average of 84 mW m −2 ). Subtracting the contribution from radionuclide decay in the continental lithosphere leaves 36 ± 4 TW being supplied from below.
The two principal mechanisms by which this heat is delivered to the top of the upper mantle are through passive mantle upwelling (in response to sinking slabs and mid-oceanic ridge spreading) or through actively buoyant hot plumes. The heat flux of an upwelling plume, Q H , can be related to its buoyancy flux by assuming that excess buoyancy is generated by thermal expansion. Thus where Cp is the specific heat per unit mass of mantle. Using Cp = 1250 J kg −1 K −1 and α = 3 × 10 −5 K −1 , a useful rule of thumb is that heat flux in TW is equal to 1 24 of the buoyancy flux measured in Mg s −1 .
Previous global flux estimates derived using plate velocitybased methodologies obtain a total cumulative heat flow of 2.0 ± 0.3 TW carried to the surface by plumes, which is in good agreement with the new volumetrically derived value of 1.9±0.3 TW for the magmatic swells alone, and supports the notion that passive upwelling is the dominant mode of heat transfer within the upper mantle (Davies, 1988;Sleep, 1990;King & Adam, 2014). However, we argue that this value represents a lower bound for three reasons. First, using all excess topography (including amagmatic swells) suggests that, for the endmember case where all of these features represent active upwellings, the total heat flux carried into the upper mantle could be as high as 6.3±1.1 TW. Secondly, the removal of the degree 0-3 components from the global swell topography grid may remove too much excess topography from Africa and the Pacific Ocean, which particularly penalises fluxes for magmatic swells. Although these long-wavelengths undoubtedly include contributions from lower mantle flow, this process may also remove long-wavelength support from shallow buoyancy sources, such as the proposed plume-fed sub-oceanic asthenosphere that is thought to be particularly well developed beneath the Pacific Ocean (Morgan et al., 1995). Thirdly, it has been suggested that there may be significant chemical heterogeneity within plumes sourced from the lower mantle. This intrinsically dense basal material reduces plume buoyancy, which must therefore be compensated by a thermal buoyancy contribution that is up to ∼ 50% larger in order to facilitate rising into the upper mantle (Ballmer et al., 2013;Dannberg & Sobolev, 2015).
An additional question concerns the depth at which these buoyant features originate. The existence of plumes that traverse the entire mantle versus those that are generated from intermediate depths, such as the mantle transition zone, is much debated (e.g. Davaille, 1999). Assuming whole-mantle transfer, previous studies have used plume buoyancy flux estimates to directly infer coremantle boundary (CMB) heat flux (e.g. Davies, 1988;Sleep, 1990). If the CMB is in thermal steady state with no significant build-up or loss of heat through time, basal flux into the mantle from the core must be balanced by removal. Two mechanisms exist to facilitate this removal -buoyant plumes that form in the thermal boundary layer and advect heat upwards away from the CMB, and conductive reheating of cold, subducted oceanic slabs that have sunk through the mantle until they lie upon the CMB (Silver et al., 1988). The slab heat flux entering the mantle at subduction zones can be gauged using the thermal structure of oceanic lithosphere. In a non-expanding Earth, the mean areal flux of subducting material is equal to the ridge production rate of ∼ 3 km 2 yr −1 , and the linear decrease of areal extent of ocean floor as a function of age indicates that there is an equal probability of subducting any seafloor out to a maximum of ∼ 180 Ma (Rowley, 2002). Adopting the RHCW18 plate model, plate thickness can be approximated using the base of the mechanical boundary layer, which coincides with the 1175 • C isothermal surface . By equally subducting all ages of this plate at a rate of 3 km 2 yr −1 , a total heat flux of 23 ± 3 TW is required to rewarm this lithosphere up to a potential temperature of 1333 • C.
The depths at which these refrigerative effects are concentrated depends upon the sinking rate, thickness and depth of penetration of the slab. Anderson (2002) suggests that 3.8 ± 0.6 TW of heat may be absorbed at the CMB by slabs that have traversed the whole mantle, with the remaining heat removed by a mixture of strong and weak plumes. Numerical simulations of mantle convection have produced contrasting results on the importance of CMB heat removal by slabs versus upwellings. Some studies suggest that warming of slab material dominates, with plumes typically carrying ∼ 40% of CMB flux and possibly as little as 10% (Labrosse, 2002;Mittelstaedt & Tackley, 2006). Other studies have found that 80-90% of CMB heat flux is removed by upwelling plumes, although this flux could potentially halve during rise to the upper mantle as a result of diffusive and sub-adiabatic cooling (Bunge, 2005;Zhong, 2006;Leng & Zhong, 2008).
Independent estimates of CMB heat flux have been obtained from a range of seismic and mineral physics observations, as well as geodynamo simulations (see Lay et al., 2008;Nimmo, 2015;Olson, 2016 and references therein). Despite the range of suggested values, consensus has begun to emerge for a value of ∼ 12 ± 5 TW. While the heat sink provided by warming of cold subducted slabs remains relatively unconstrained, it is possible that a large fraction of this heat is removed by upwelling plumes, and is not inconsistent with an upper mantle plume flux that is significantly larger than the purely magmatic swell-derived value of 2 TW. Indeed, tanta-lising evidence from seismology suggests that large lower mantle plumes do rise from the CMB, although these structures are at the limits of tomographic resolution French & Romanowicz, 2015). Estimates of the total heat carried by these features are 10-30 TW .

Conclusions
Two different methods for determining buoyancy flux from plume swells have been compared. Both approaches neglect the effect of normal stresses arising from vertical flow within the conduit and assume that excess topography is supported by thermal isostasy. The traditional approach exploits plate velocities and spreading rates. Its key assumption is that excess buoyancy is primarily within the lithosphere and that the asthenosphere is dominated by Couette flow, such that buoyancy travels away from the conduit at velocities determined by the overriding plate. However, there is growing evidence for asthenospheric support of swells and for Poiseuille flow of sub-plate material, particularly in the North Atlantic region where lateral velocities in the asthenosphere can be an order of magnitude faster than overriding lithospheric plates. Thus it is desirable to implement an alternative volumetric approach that is free of plate velocity terms. The key assumption of this method is that the plume is in quasi-steady state, where influx of new buoyancy is balanced by buoyancy loss, such as conduction of heat into the surrounding ambient mantle and through the overlying lithospheric plate.
The volumetric approach has been applied in two case studies. The Icelandic plume has a buoyancy flux of 3.9 ± 0.5 Mg s −1 and Hawaii has 2.9 ± 0.6 Mg s −1 . Both estimates assume a buoyancy decay timescale of 45 Myr, which may represent conductive cooling of a ∼ 120 km thick layer of hot plume-head material. Extending this analysis to 48 major magmatic swells yields a total upper mantle heat flux of ∼ 2 TW, which is almost certainly a lower bound for plume flux. This value could be as high as 6.3 ± 1.1 TW if amagmatic swells are also included. Episodic arrival of large plume heads that generate large igneous provinces will supply pulses of heat in addition to this steady-state background flux. These greater values are more consistent with recent CMB heat flux estimates of 12 ± 5 TW, suggesting significant basal heat may be transported by upwelling plumes throughout the mantle. These results accord with seismological evidence for large plumelike features within the lower mantle.
Future measurements of excess swell morphology will undoubtedly require plume buoyancy fluxes to be revised. Continental constraints on residual topography are particularly difficult to isolate at present, such that oceanic swell fluxes are likely to be more reliable. In future studies, the accuracies of plate velocity-based and volumetric approaches for estimating plume flux could be tested using numerical models of mantle convection. Further work is also required to explore physical implications of the buoyancy decay timescale, including whether our assumed value of 45 Myr is appropriate for all plume swells.
Continued analyses will help to illuminate complex interactions between plume buoyancy flux, rheology, excess temperature, melt fraction, melting depths and isotopic signatures. Improved understanding of these relationships is necessary to constrain the distribution and lengthscales of chemical heterogeneities within Earth's mantle. Finally, we note that our volumetric approach for deriving plume buoyancy flux can be applied to other planets with convecting interiors, where velocity-based analysis may not be possible due to the absence of plate tectonics.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Appendix: Summary of Sleep (1990) plate velocitybased plume flux catalogue Sleep (1990) produced a global inventory of 37 buoyancy flux estimates using plate velocity-based methodologies. In addition to Iceland, the on-axis approach (Section 1.1) was applied at the Afar, Azores, Bouvet, Cobb, Easter, Galapagos and Tristan plumes. The affected along-axis ridge lengths were separately estimated in each location, based upon the extent of elevated bathymetry and presence of sub-aerial plateaux. Spreading rates were taken from plate reconstructions or local measurements, and the lithospheric and asthenospheric thicknesses were assumed to consistently be z l = za = 100 km. Each volume flux was subsequently converted into buoyancy flux using Equation (2) and the same excess temperature of ∆T = 225 • C obtained from Iceland.
The off-axis approach (Section 2.1) was applied to the East Australian, Hoggar, Reunion and Yellowstone plumes. A slightly modified version was used for plumes on very slow-moving plates, as sub-lithospheric temperature anomalies were anticipated to build up, resulting in a component of the swell being supported by asthenospheric isostasy. Careful analysis of heat flow and geoid anomalies allows quantification of this asthenospheric contribution. In the case of Cape Verde, Sleep (1990) estimates that approximately one third to half of the excess topography is generated within the asthenosphere. Thus only the lithospheric component of excess swell topography is used in the calculation of buoyancy flux of Q B = 1.6 Mg s −1 . Bermuda is the other location where this asthenospheric contribution to the swell was removed.
The buoyancy flux of the Crozet plume was quantified using a third approach based upon stagnation distances. The basic principle is that radial flow of hot asthenosphere away from the conduit interacts with the drag of the overriding plate. A stagnation point occurs upstream of the conduit where the velocities are equal and opposite. The closest approach of this point to the conduit can therefore be used to quantify the volume flux, assuming that the average asthenospheric velocity is half of the plate velocity when unperturbed by a plume. Sleep (1990) fits a stagnation distance of 140 km to the Crozet geoid anomalies, giving a volume flux of 22 m 3 s −1 . Application of Equation (2) with an excess temperature of ∆T = 225 • C yielded a buoyancy flux of Q B = 0.5 Mg s −1 .
To compile a global catalogue, a range of additional buoyancy flux values were inferred based upon comparisons to aforementioned swells. Canary was considered to be approximately twothirds the size of Cape Verde, whilst the Great Meteor, Fernando, St Helena, Trindade, Discovery and Meteor were thought to be about one third the size. Tasmantid and Lord Howe were assigned the same value as the East Australia hotspot. Several swells overlap to generate the South Pacific Superswell, which has a total Q B = 13.3 Mg s −1 using the off-axis approach. This flux was equally divided amongst the Macdonald, Marquesas, Pitcairn and Tahiti hotspots. The Caroline, Samoa, San Felix, and Juan Fernandez were considered to be half the strength of these superswell hotspots. The Louisville swell was considered similar to those in the Tasman Sea, whilst Bowie and Baja were inferred based on the Cobb buoyancy flux. Kerguelen was assigned the same value as Crozet. Finally, there are a handful of notable omissions including the Comoros, St Paul, Eifel, Tibesti and Cameroon hotspots. The plate velocity-based buoyancy flux estimates for all 37 swells considered in the analysis of Sleep (1990) are listed in Table 2  2. Buoyancy decay timescale from the  Hawaiian numerical model.

Swell Outlines
Swell topography is obtained using the oceanic residual topography database of Hoggard, Winterbourne, Czarnota, & White (2017), which has been corrected for the effects of sedimentary loading, crustal thickness variations and age-depth cooling using the RHCW18 plate model . The longest wavelengths of residual topography are primarily supported by flow in the deeper mantle rather than asthenospheric isostasy .
We therefore fit spherical harmonic basis functions to this database up to and including degree 30 (wavelengths down to ∼ 1300 km) and subtract the degree 0-3 components from the field. Swell outlines for individual volumetric analysis were digitised by eye and are shown in Figure S1. Figure S1: Swell polygons and topography. Air-loaded onshore, water-loaded offshore. 1

swell decay timescale
The decay of swell topography for the  model of the Hawaiian plume is shown in their Figure 2a. Here, we have digitised their results and converted the x axis into time using their plate velocity of 86 mm yr −1 ( Figure S2). The upper green line is for the full model including both thermal and depletion buoyancy sources, whilst the lower blue line is for a purely thermal run. Unfortunately, the model only extends to ∼ 20 Myr downstream of the conduit. Assuming that the decay is linear, extrapolating the purely thermal run suggests that the swell topography will fully subside by the time the lithosphere has been downstream of the conduit for 54 Myr. Note that the model with depletion includes a contribution to swell uplift of ∼ 400 m arising from depletion buoyancy (difference between the two curves) that will almost certainly be permanent, and therefore not compatible with this simplified linear swell decay. Nevertheless, this value is likely to be an upper bound on the thermal buoyancy decay time, due to the absence of small-scale convection in the Ribe & Christensen (1999) model. Figure S2: Hawaiian swell decay in the model of . Blue line = swell uplift for a model including purely thermal buoyancy; green line = same for a model including depletion buoyancy; red dashed lines = linear fits to swell decay.