Seismic tomography of Earth’s mantle images abundant slab remnants, often located in close proximity to active subduction systems. The impact of such remnants on the dynamics of subduction remains under explored. Here, we use simulations of multi-material free subduction in a 3-D spherical shell geometry to examine the interaction between visco-plastic slabs and remnants that are positioned above, within and below the mantle transition zone. Depending on their size, negatively buoyant remnants can set up mantle flow of similar strength and length scales as that due to active subduction. As such, we find that remnants located within a few hundred km from a slab tip can locally enhance sinking by up to a factor 2. Remnant location influences trench motion: the trench advances towards a remnant positioned in the mantle wedge region, whereas remnants in the sub-slab region enhance trench retreat. These motions aid in rotating the subducting slab and remnant towards each other, reducing the distance between them, and further enhancing the positive interaction of their mantle flow fields. In this process, the trench develops along-strike variations in shape that are dependent on the remnant’s location. Slab-remnant interactions may explain the poor correlation between subducting plate velocities and subducting plate age found in recent plate tectonic reconstructions. Our results imply that slab-remnant interactions affect the evolution of subducting slabs and trench geometry. Remnant-induced downwelling may also anchor and sustain subduction systems, facilitate subduction initiation, and contribute to plate reorganisation events.

Ocean Island Basalts (OIBs) are generated by mantle plumes, with their geochemistry controlled by a combination of source composition, temperature, and thickness of overlying lithosphere. For example, OIBs erupting onto thicker, older oceanic lithosphere are expected to exhibit signatures indicative of higher average melting pressures. Here, we quantitatively investigate this relationship using a global dataset of Neogene and younger OIB compositions. Local lithospheric thicknesses are estimated using theoretical plate-cooling models and Bayes factors are applied to identify trends. Our findings provide compelling evidence for a correlation between OIB geochemistry and lithospheric thickness, with some variables SiO2,  Al2O3, FeO, Lu, Yb and λ2) showing linear trends that can be attributed to increasing average melting pressure, whereas others (λ0 and λ1, CaO) require a bi-linear fit with a change in gradient at ~55 km. Observed variations in highly incompatible elements are consistent with melt fractions that decrease with increasing lithospheric thickness, as expected. Nevertheless, at thicknesses beyond ~55 km, the implied melt fraction does not decrease as rapidly as suggested by theoretical expectations. This observation is robust across different lithospheric thickness estimates, including those derived from seismic constraints. We interpret this result as weak plumes failing to effectively thin overlying lithosphere and/or producing insufficient melt to erupt at the surface, in combination with a 'memory effect' of incomplete homogenisation of melts during their ascent. This view is supported by independent estimates of plume buoyancy flux, indicating that OIB magmatism on older lithosphere may be biased towards hotter plumes.

Rapid plate motion, alongside pronounced variations in age and thickness of the Australian continental lithosphere, make it an excellent location to assess the relationship between seismic anisotropy and lithosphere-asthenosphere dynamics. In this study, SKS and PKS shear-wave splitting is conducted for 176 stations covering the transition from the South Australian Craton to eastern Phanerozoic Australia. Comparisons are made with models of lithospheric thickness as well as numerical simulations of mantle flow. Splitting results show uniform ENE-WSW aligned fast directions over the Gawler Craton and broader South Australian Craton, similar to the orientation of crustal structures generated during an episode of NW-SE directed compression and volcanism ~1.6 billion years ago. We propose that heat from volcanism weakened the lithosphere, aiding widespread lithospheric deformation, which has since been preserved in the form of frozen-in anisotropy. Conversely, over eastern Phanerozoic Australia, fast directions show strong alignment with the NNE absolute plate motion. Overall, our results suggest that when the lithosphere is thin (<125 km), lithospheric contributions are minimal and contributions from asthenospheric anisotropy dominate, reflecting shear of the underlying mantle by Australia’s rapid plate motion above. Further insights from geodynamical simulations of the regional mantle flow-field, which incorporate Australian and adjacent upper mantle structure, predict that asthenospheric material would be drawn in from the south and east towards the fast-moving continental keel. Such a mechanism, alongside interactions between the flow field and lithospheric structure, provides a plausible explanation for smaller-scale anomalous splitting patterns beneath eastern Australia that do not align with plate motion.

The effects of sphericity are regularly neglected in numerical and laboratory studies that examine the factors controlling subduction dynamics. Most existing studies have been executed in a Cartesian domain, with the small number of simulations undertaken in a spherical shell incorporating plates with an oversimplified rheology, limiting their applicability. Here, we simulate free-subduction of composite visco-plastic plates in 3-D Cartesian and spherical shell domains, to examine the role of sphericity in dictating the dynamics of subduction, and highlight the limitations of Cartesian models. We identify two irreconcilable differences between Cartesian and spherical models, which limit the suitability of Cartesian-based studies: (i) the presence of sidewall boundaries in Cartesian models, which modify the flow regime; and (ii) the reduction of space with depth in spherical shells, alongside the radial gravity direction, which cannot be captured in Cartesian domains. Although Cartesian models generally predict comparable subduction regimes and slab morphologies to their spherical counterparts, there are significant quantitative discrepancies. We find that simulations in Cartesian domains that exceed Earth’s dimensions overestimate trench retreat. Conversely, due to boundary effects, simulations in smaller Cartesian domains overestimate the variation of trench curvature driven by plate width. Importantly, spherical models consistently predict higher sinking velocities and a reduction in slab width with depth, particularly for wider subduction systems, enhancing along-strike slab buckling and trench curvature. Results imply that sphericity must be considered when simulating Earth’s subduction systems, and that it is essential for accurately predicting the dynamics of subduction zones of width ~2400 km or more.

A significant component of Earth’s surface topography is maintained by stresses induced by underlying mantle flow. This ‘dynamic’ topography cannot be directly observed, but it can be approximated — particularly at longer wavelengths — from measurements of residual topography, which are obtained by removing isostatic effects from the observed topography. However, as these measurements are made at discrete, unevenly-distributed locations on Earth’s surface, inferences about global properties can be challenging. In this paper, we present and apply a new approach to transforming point-wise measurements into a continuous global representation. The approach, based upon the statistical theory of Gaussian Processes, is markedly more stable than existing approaches — especially for small datasets. We are therefore able to infer the spatial pattern, wavelength and amplitude of residual topography using only the highest-quality oceanic spot measurements within the database of Hoggard et al. (2017). Our results indicate that the associated spherical harmonic power spectrum peaks at l=2, with power likely in the range 0.46–0.76 km^2. This decreases by over an order of magnitude to around 0.02 km^2 at l=30. Around 85% of the total power is concentrated in degrees 1–3. Our results therefore confirm previous findings: Earth’s residual topography expression is principally driven by deep mantle flow, but shallow processes are also crucial in explaining the general form of the power spectrum. Finally, our approach allows us to determine the locations where collection of new data would most impact our knowledge of the spectrum.

Many of the factors expected to control the dynamics and evolution of Earth’s subduction zones are under-explored in an Earth-like spherical geometry. Here, we simulate multi-material free-subduction of a complex rheology slab in a 3-D spherical shell domain, to investigate the effect of plate age (simulated by covarying plate thickness and density) and width on the evolution of subduction systems. We find that the first-order predictions of our spherical cases are generally consistent with existing Cartesian studies: (i) as subducting plate age increases, slabs retreat more and subduct at a shallower dip angle, due to increased bending resistance and sinking rates; and (ii) wider slabs can develop along-strike variations in trench curvature due to toroidal flow at slab edges, trending towards a ‘W’-shaped trench with increasing slab width. We find, however, that these along-strike variations are restricted to older, stronger, retreating slabs:. Younger slabs that drive minimal trench motion remain relatively straight along the length of the subduction zone. We summarise our results into a regime diagram, which highlights how slab age modulates the effect of slab width, and present examples of the evolutionary history of subduction zones that are consistent with our model predictions.

Geodynamical simulations underpin our understanding of upper-mantle processes, but their predictions require validation against observational data. Widely used geophysical datasets provide limited constraints on dynamical processes into the geological past, whereas under-exploited geochemical observations from volcanic lavas at Earth's surface constitute a valuable record of mantle processes back in time. Here, we describe a new peridotite-melting parameterization, BDD21, that can predict the incompatible-element concentrations of melts within geodynamical simulations, thereby providing a means to validate these simulations against geochemical datasets. Here, BDD21's functionality is illustrated using the Fluidity computational modelling framework, although it is designed so that it can be integrated with other geodynamical software. To validate our melting parameterization and coupled geochemical-geodynamical approach, we develop 2-D single-phase flow simulations of melting associated with passive upwelling beneath mid-oceanic ridges and edge-driven convection adjacent to lithospheric steps. We find that melt volumes and compositions calculated for mid-oceanic ridges at a range of mantle temperatures and plate-spreading rates closely match those observed at present-day ridges. Our lithospheric-step simulations predict spatial and temporal melting trends that are consistent with those recorded at intra-plate volcanic provinces in similar geologic settings. Taken together, these results suggest that our coupled geochemical-geodynamical approach can accurately predict a suite of present-day geochemical observations. Since our results are sensitive to small changes in upper-mantle thermal and compositional structure, this novel approach provides a means to improve our understanding of the mantle's thermo-chemical structure and flow regime into the geological past.

Several of Earth’s intra-plate volcanic provinces are difficult to reconcile with the mantle plume hypothesis. Instead, they exhibit characteristics that are better explained by shallower processes involving the interplay between uppermost mantle flow and the base of Earth’s heterogeneous lithosphere. The mechanisms most commonly invoked are edge-driven convection (EDC) and shear-driven upwelling (SDU), both of which act to focus upwelling flow, and the associated decompression melting, adjacent to steps in lithospheric thickness. In this study, we undertake a systematic numerical investigation, in both 2-D and 3-D, to quantify the sensitivity of EDC, SDU and their associated melting to several key controlling parameters. Our simulations demonstrate that the spatial and temporal characteristics of EDC are sensitive to the geometry and material properties of the lithospheric step, in addition to the depth-dependence of upper mantle viscosity. These simulations also indicate that asthenospheric shear can either enhance or reduce upwelling velocities and predicted melt volumes, depending upon the magnitude and orientation of flow relative to the lithospheric step. When combined, such sensitivities explain why step changes in lithospheric thickness, which are common along cratonic edges and passive margins, only produce volcanism at isolated points in space and time. Our predicted trends of melt production suggest that, in the absence of potential interactions with mantle plumes, EDC and SDU are viable mechanisms only for Earth’s shorter-lived, lower-volume intra-plate volcanic provinces.

Several of Earth’s intra-plate volcanic provinces are hard to reconcile with the mantle plume hypothesis. Instead, they exhibit characteristics that are more compatible with shallower processes that involve the interplay between uppermost mantle flow and the base of Earth’s heterogeneous lithosphere. The mechanisms most commonly invoked are edge-driven convection (EDC) and shear-driven upwelling (SDU), both of which act to focus upwelling flow and the associated decompression melting adjacent to steps in lithospheric thickness. In this study, we undertake a systematic numerical investigation, in both 2-D and 3-D, to quantify the sensitivity of EDC, SDU, and the associated melting to key controlling parameters. Our simulations demonstrate that the spatio-temporal characteristics of EDC are sensitive to the geometry and material properties of the lithospheric step, in addition to the magnitude and depth-dependence of upper mantle viscosity. These simulations also indicate that asthenospheric shear can either enhance or reduce upwelling velocities and the associated melting, depending upon the magnitude and orientation of flow relative to the lithospheric step. When combined, such sensitivities explain why step changes in lithospheric thickness, which are common along cratonic edges and passive margins, only produce volcanism at isolated points in space and time. Our predicted trends of melt production suggest that, in the absence of potential interactions with mantle plumes, EDC and SDU are viable mechanisms only for Earth’s shorter-lived, lower-volume intra-plate volcanic provinces.

Several of Earth's intra-plate volcanic provinces occur within or adjacent to continental lithosphere, with many believed to mark the surface expression of upwelling mantle plumes. Nonetheless, studies of plume-derived magmatism have generally focussed on ocean-island volcanism, where the overlying rigid lithosphere is of uniform thickness. Here, we investigate the interaction between mantle plumes and heterogeneous continental lithosphere using a series of geodynamical models. Our results demonstrate that the spatio-temporal magmatic expression of plumes in these continental settings is complex and strongly depends on the location of plume impingement, differing substantially from that expected beneath oceanic lithosphere. Where plumes ascend beneath thick continental cratons, the overlying lid locally limits decompression melting. However, gradients in lithospheric thickness channel plume material towards regions of thinner lithosphere, activating magmatism away from the plume conduit, sometimes simultaneously at locations more than a thousand kilometres apart. This magmatism regularly concentrates at lithospheric steps, where it may be difficult to distinguish from that arising through edge-driven convection, especially if differentiating geochemical signatures are absent, as implied by some of our results. If plumes impinge in regions of thinner lithosphere, the resulting asthenospheric flow regime can force material downwards at lithospheric steps, shutting off pre-existing edge-related magmatism. In addition, under certain conditions, the interaction between plume material and lithospheric structure can induce internal destabilisation of the plume pancake, driving complex time-dependent magmatic patterns at the surface. Our study highlights the challenges associated with linking continental magmatism to underlying mantle dynamics and motivates an inter-disciplinary approach in future studies.

The role of Earth’s spherical geometry in modulating the evolution of subduction zones is poorly understood. Here, we simulate multi-material free-subduction in a 3-D spherical shell domain, to investigate the effect of plate thickness, density (combined approximating age) and width on the evolution of subduction systems. To isolate the role of sphericity, we compare results with equivalent Cartesian models. The first-order predictions of our spherical cases are generally consistent with existing Cartesian studies: (i) slabs retreat more, at a shallower dip, as plate age increases, due to increased bending resistance and sinking rates; and (ii) wider slabs can develop along-strike variations in trench curvature, trending towards a ‘W’-shape, due to toroidal flow at slab edges. We find, however, that these along-strike variations are restricted to older, stronger, retreating slabs. When compared to slabs in Cartesian models, in a spherical domain: (i) slabs descend faster, due to the convergence of downwelling material with depth; (ii) these faster sinking rates reduce the time available for bending at the trench, resulting in effectively stronger slabs; (iii) the curvature of slabs increases their effective strength; and (iv) the curvature of the transition zone tends to enhance slab stagnation. These differences between spherical and Cartesian cases become more prominent as slab width increases. Taken together, our results suggest that Cartesian models are suitable for simulating narrow subduction zones, but spherical models should be utilised when investigating subduction zones wider than ~ 2000 km: at such length-scales, the consequences of Earth’s curvature cannot be ignored.