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.

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.

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.

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.