The Role of Slab Remnants in Modulating Free Subduction Dynamics: A 3‐D Spherical Numerical Study

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 underexplored. 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 toward 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 toward 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 reorganization events.

The existence of remnants requires slab detachment or break-off, with several mechanisms proposed.In plate reorganization events, trenches may be abandoned, with the connected slab eventually detaching from the surface plate (e.g., Matthews et al., 2012;Müller et al., 2016;Whittaker et al., 2007).The arrival of an active spreading center or buoyant continental lithosphere at the trench can also induce subduction termination (e.g., Burkett & Billen, 2009; J. H. Davies & von Blanckenburg, 1995;Duretz et al., 2014;Faccenna et al., 2006;Seton et al., 2015;Wong A Ton & Wortel, 1997;Wortel & Spakman, 2000).Furthermore, the subduction of faults or buoyant anomalies can facilitate slab tearing (e.g., Abratis & Wörner, 2001;Pallares et al., 2007;Thorkelson & Taylor, 1989), which can also occur to accommodate changes in plate geometry, as postulated for STEP faults (e.g., Govers & Wortel, 2005;Obayashi et al., 2009;Thorkelson, 1996).In addition, slab remnants can be generated directly at 660 km depth, through so-called slab orphaning, where the slab tip breaks off at the top of the lower mantle and is abandoned by its parent slab (Grima et al., 2020).Following break off from their active subduction systems, the negative buoyancy and rheological properties of remnants make it likely that they continue to influence mantle convection and the dynamics of nearby active subduction systems.The dynamics of such interactions, however, has not yet been systematically examined.
In this paper, we investigate how actively subducting slabs interact with, and are influenced by slab remnants, above, within, and below the MTZ.We focus on this depth range for two reasons: (a) whilst slabs sink through the upper mantle in a few million years, they can stagnate within the MTZ for tens of millions of years in response to phase buoyancy, rheological complexities, and a likely viscosity increase at these depths (e.g., Christensen & Yuen, 1984;Čížková & Bina, 2019;Garel et al., 2014;Goes et al., 2017;Karato et al., 2001;Ringwood, 1975;Tackley et al., 1993); and (b) slab-transition zone interaction strongly affects the subduction mode and thus controls the surface expressions of subduction systems, for example, through plate motions and the evolution of trench shape (e.g., Cerpa et al., 2022;Garel et al., 2014;Goes et al., 2017;Ribe, 2010;Stegman et al., 2010;Torii & Yoshioka, 2007).
Slab remnants are frequently imaged around the MTZ, most in close proximity to active subduction systems (Figure 1).Remnants from the extensive Farallon and adjacent plates have been imaged beneath North and South America (e.g., Sigloch et al., 2008;van Benthem et al., 2013;van der Lee & Nolet, 1997).Proto-Caribbean slab fragments have been identified at ∼1,000 km depth below northeastern South America by combining seismic tomography and plate reconstructions (e.g., Braszus et al., 2021;van Benthem et al., 2013).Tethyan slab fragments have been found over a large area, extending from beneath the Mediterranean to India, across a range of depths (e.g., Hafkenscheid et al., 2006;Wortel & Spakman, 2000).South-east Asia has a complex subduction history resulting in a multitude of interacting slabs and remnants (e.g., Hall & Spakman, 2015;Li et al., 2008;van der Meer et al., 2018).Proto-South China Sea Plate remnants have been imaged in the lower mantle beneath Northern Borneo (e.g., Hall & Spakman, 2015;Pilia et al., 2023;Wu & Suppe, 2018;Zahirovic et al., 2014); Wu and Suppe (2018) further suggest that northern remnants of the Proto-South China Sea plate are stagnating at the transition zone under the South China Sea, adjacent to the slab currently subducting at the Manila trench.To the east of the region, subducted remnants of the Philippine Sea plate can be inferred from tomography at ∼800 km depth, in agreement with the reconstructed locations of major subduction zones in the region at 30 Ma (e.g., Hall & Spakman, 2015;Widiyantoro et al., 2011).A relic of the Kula slab is interpreted to be present above the transition zone beneath the Bering Sea, just to the north of where the Pacific Plate subducts at the Aleutian trench (e.g., Gorbatov et al., 2000;van der Meer et al., 2018).Flat lying slab remnants in the transition zone, proposed to be from Paleogene Pacific subduction at the Melanesia Arc, have been imaged in close proximity to current Tongan subduction (e.g., Hall & Spakman, 2002;Pysklywec et al., 2003).There are, therefore, many examples of slabs subducting into a region of the mantle with pre-existing remnants.
The most often discussed effect of remnants on subduction systems is tectonic plate reorganization.Pysklywec and Ishii (2005) studied the effect of slab remnants on the dynamics of active subduction using a suite of 2-D numerical models, finding that they can either drive trench retreat or induce slab detachment and a reversal of subduction polarity, depending on remnant location relative to the trench.This mechanism may contribute to tectonic plate reorganization events, where Pysklywec et al. (2003) suggest that the subducted Tonga slab may have caused detachment of the Vitiaz slab and initiation of subduction at the New Hebrides Trench.It has also been suggested that slab-induced downwelling flow may play a role in subduction initiation (Baes & Sobolev, 2017;Crameri et al., 2020), that mantle flow induced by slab remnants influences the dip of adjacent subduction systems (Hu & Gurnis, 2020), and that the flow-field driven by remnants may facilitate the opening of intra-plate basins (e.g., Capitanio et al., 2009;Pysklywec & Mitrovica, 1999;Yang et al., 2018).Remnants also play a role in driving and modulating large-scale mantle flow, which will ultimately influence the dip and orientation of subducting slabs, trench migration rates and upper plate deformation (e.g., Baes & Sobolev, 2017;Chertova et al., 2018;Conrad & Lithgow-Bertelloni, 2002;Ficini et al., 2017;Hager & O'Connell, 1979;Holt & Royden, 2020;Husson, 2012;Stotz et al., 2018).Despite this, a systematic study into how the dynamics of subduction are influenced by slab remnants (and vice versa) has not yet been undertaken.Furthermore, the effects of a 3-D spherical shell domain, which constrains the space available for mantle flow (F.Chen et al., 2022a), on slab-remnant interaction, have not previously been considered.
In this study, we fill this gap by investigating the effect of slab remnants on the dynamics of subducting slabs and the evolution of subduction systems in a 3-D spherical shell domain.We build on seismic tomography observations of remnant slabs (e.g., Li et al., 2008;van der Meer et al., 2018) and previous models (e.g., Pysklywec & Ishii, 2005), and simulate the interaction between remnant slabs and subduction zones through a suite of 3-D spherical shell cases, using the modeling approach developed in F. Chen et al. (2022aChen et al. ( , 2022b)).
As with freely subducting slabs (Capitanio et al., 2007;Ribe, 2010), remnants act as Stokes sinkers, as they have viscosities that are orders of magnitude higher than the surrounding mantle, limiting their deformation (e.g., Jarvis & Lowman, 2007;Quéré et al., 2013).The velocity of Stokes sinkers depends on their shape, density and the viscosity of surrounding material (e.g., Capitanio et al., 2007).Accordingly, we vary remnant density, shape and size, allowing us to examine the impact of remnants with different sinking velocities and associated mantle flow fields.Subducting slabs set up significant mantle flow around them, to distances of approximately slab half width, up to 1,000-1,500 km in along-strike direction and on the order of their slab length in trench-perpendicular directions (distances dependent on mantle viscosity) (F.Chen et al., 2022b;Király et al., 2017;Piromallo et al., 2006;Schellart et al., 2007).Sinking remnant slabs are expected to set up mantle flow regimes of similar strength and scale.The flow regime induced by a sinking remnant will thus exert forces on nearby subduction systems, the impact of which not only depends on remnant properties, but also on remnant location relative to the subducting plate.Accordingly, we also vary remnant location and orientation relative to the trench, placing remnants at different distances from the slab, below the mantle wedge, the slab tip, and the sub-slab regions.In addition, to isolate the role of remnant buoyancy, we examine whether the presence of a purely viscous remnant (i.e., a remnant with no buoyancy anomaly relative to adjacent mantle) can alter the flow regime and dynamics of an active subduction system.
In the following sections, we first describe the setup of our numerical models, including the governing equations, the initial geometric configurations, boundary conditions and material properties, and summarize the different cases examined.We then quantitatively analyze our models to reveal: (a) how remnant density influences the dynamics of adjacent subduction zones; (b) how remnant location and orientation modulate their effects; and (c) how remnants of different dimensions influence trenches of different width.We end by discussing the implications of our results, and explore their applicability for understanding interactions between slab remnants and active subduction zones on Earth.

Governing Equations and Numerical Strategy
We simulate free subduction of a composite visco-plastic plate into ambient underlying mantle inside a 3-D spherical shell domain using a multi-material approach.In cases where we can exploit the symmetry of the system, our simulations are undertaken in a hemispherical shell, following F. Chen et al. (2022b).Assuming incompressibility, the governing equations for this problem are the continuity equation, the conservation of momentum equation at infinite Prandtl number, and an advection equation for tracking different materials, In the above equations, u represents velocity, p dynamic pressure, μ viscosity, ρ density, g gravitational acceleration, k the unit vector in the direction opposite to gravity, and Γ the material volume fraction (Γ = 1 in a region occupied by a given material and Γ = 0 elsewhere).At material interfaces, the average viscosity is calculated through a geometric mean, where μ i is the viscosity of material i, and Γ i is the relative volume fraction of material i in the vicinity of the finiteelement node at which the effective viscosity μ ave is needed.As demonstrated by many previous studies, such multi-material mechanical models capture the key primary characteristics of subduction dynamics, without the need to model thermal evolution (e.g., Bellahsen et al., 2005;Capitanio et al., 2007;F. Chen et al., 2022b;Ribe, 2010).
Simulations are undertaken using Fluidity (e.g., D. R. Davies et al., 2011Davies et al., , 2016;;Kramer et al., 2012;Kramer, Davies, & Wilson, 2021;Le Voci et al., 2014), an adaptive, anisotropic, unstructured-mesh finite element and control volume computational modeling framework, capable of efficiently simulating multi-material wholemantle visco-plastic (Tosi et al., 2015) subduction in spherical shell geometries (F.Chen et al., 2022aChen et al., , 2022b).Fluidity's adaptive mesh capabilities allow our simulations to achieve a local resolution of ∼3 km in regions of dynamical significance (i.e., in regions of high solution curvature within and adjacent to the subducting slab, the remnant and at the interface between upper and lower mantle), with coarser resolution of up to ∼300 km elsewhere.
The importance of sphericity in simulating the dynamics of subduction, particularly for wider subduction systems, has been demonstrated by F. Chen et al. (2022a), who identify two key limitations of Cartesian compared to spherical models: (a) the presence of sidewall boundaries in Cartesian models, which modify the flow regime; and (b) the reduction of space with depth in spherical shells, alongside the radial gravity direction, which cannot be captured in Cartesian domains.This motivates the use of spherical models herein.

Geometry, Boundary Conditions, and Material Properties
The configuration of our reference models without remnants (cases with full widths-W-of W2400 and W4800) follows the setup of F. Chen et al. (2022a) for a subducting plate with plate buoyancy and thickness in the midrange for natural subduction (cases S_W2400 and S_W4800, respectively).Simulations are undertaken in a hemispherical shell domain (thus exploiting the symmetry of the system to halve the computational domain's extent) with outer and inner radii that correspond to Earth's surface and CMB (Figure 2a).A free-surface boundary condition is applied on the outer surface, with free-slip conditions on the symmetric mid-plane and CMB.Gravity points radially toward the center of the sphere.A factor of 50 viscosity increase is included at 660 km depth.Parameters common to all simulations are listed in Table 1.
The subducting plate has length L = 2,200 km, thickness h = 70 km and a density contrast with adjacent mantle of Δρ = 80 kg m 3 .Highly viscous side plates that cover the entire domain adjacent to the subducting plate are required to prevent the undesired narrowing effect on the subducting plate from shallow lateral flow (as in Holt et al., 2017).The initial slab tip geometry is prescribed with a bending radius of 250 km and an angle of 77°( Figure 2b), following F. Chen et al. (2022b).The subducting lithosphere is a composite plate comprising a core isoviscous layer (thickness h c = 30 km) embedded in upper and lower visco-plastic layers with viscosities following a von Mises law, building on OzBench et al. (2008).Upper and lower plates are assigned the minimum viscosity between the Newtonian viscosity μ Newt and an effective von Mises viscosity μ vM , such that purely viscous deformation occurs when the second invariant of the stress tensor τ II = 2μ εII (where εII is the second invariant of strain rate tensor) is below the critical yield stress, τ yield .The effective viscosity of visco-plastic layers is given by: With the parameters specified in Table 1, this yields a slab-mantle viscosity contrast of 100 within the core of the slab, and capped at this value in the visco-plastic upper and lower layers.This is within the inferred range of slab strength estimates, albeit toward the lower end.For example, based on inferences from slab shape and trench retreat, Funiciello et al. (2008) and Ribe (2010) suggest slab-mantle viscosity contrasts of 100-500, although recent geoid based studies tend towards the lower end of this range (Mao & Zhong, 2021).
We simulate a wide-range of remnant cases, assuming a region of pre-existing isoviscous plate-like material in the mantle with a thickness (h rem ) of 70 km and a length (L rem ) of 400 km.For symmetric cases, remnants are designed to either have the same, or half, the subducting plate width, as illustrated in Figure 2c.The remnant can be horizontally (Figure 2d) or radially oriented (Figure 2e).When horizontally oriented, remnants can occupy three positions: (a) sub-slab, where the remnant is offset 500 km laterally from the initial trench location beneath the downgoing plate (at a horizontal distance of ∼700 km from the slab tip); (b) under, where the remnant is placed directly below the initial trench location; and (c) wedge, where the remnant is offset 500 km from the initial trench in the mantle wedge direction (at a horizontal distance of ∼300 km from the slab tip: Figure 2d).In asymmetric cases, remnants are horizontally oriented in either the under or wedge location, with the key difference being that the center of the remnant is not aligned with the center of the subducting plate, as illustrated in Figures 2g and 2h.Unless otherwise specified, simulations are run for 12 Myr, which provides sufficient time for slabs to sink through the upper mantle under the influence of remnant-induced flow, and interact with the upperlower-mantle transition.

Cases Examined and Quantitative Model Diagnostics
We investigate a total of 17 cases, 15 of which are symmetric, with varying plate width (w), remnant width (w rem ), remnant orientation, remnant position, remnant depth (D rem -denoting its deepest point), and density contrast between remnant and ambient mantle (Δρ rem ).In the following sections, the widths specified refer to the full widths of the plate or remnant, but in practice, for symmetric cases, we only simulate half of the width.Our choice of reference plate widths of 2,400 and 4,800 km is motivated by results from F. Chen et al. (2022b), which show that at a width of 2,400 km the subducting slab behaves relatively uniformly along-strike, but at a width of 4,800 km the trench develops a "W"-shaped curvature and slab morphology varies along-strike (see Figures 8 and 10 in Section 3 below).For horizontal remnants, the majority of the cases have an initial remnant depth of 660 km, sitting on the MTZ where slab stagnation is commonly observed (e.g., Fukao et al., 2009;Gorbatov et al., 2000;van der Meer et al., 2018;Wu & Suppe, 2018).An additional case with a remnant depth of 1,000 km is also examined to investigate the effect of deeper remnants (e.g., Braszus et al., 2021;Fukao & Obayashi, 2013;van der Meer et al., 2018;Zahirovic et al., 2014).The radially oriented remnant is placed between 300 and 700 km depth, in close proximity with the slab.While most remnants have the same density contrast with ambient mantle as our subducting plates, we have also simulated a neutrally buoyant remnant to investigate the effect of a purely viscous anomaly.All cases, and their associated parameter values, are listed in Table 2.
As Earth's subduction systems and their interactions with remnant slabs are rarely symmetric, we design two additional asymmetric cases that utilize the full spherical shell domain and illustrate the complexity of slab-remnant interactions under these scenarios.Case Asym_Under_Neutral has a neutral horizontally oriented remnant at 660 km depth, beneath the initial trench location under half of the 4,800 km wide subducting plate, as illustrated in Figure 2g.Case Asym_Wedge_Negative has a negatively buoyant remnant of 2,400 km width in the wedge location, with the center of the remnant aligned to the edge of the 4,800 km wide subducting plate, and thus extending beyond the slab edge (Figure 2h).While a more complete and systematic examination of different combinations of remnant properties and locations is required to fully quantify the impact of remnants on subduction zones, due to the computational costs of our simulations (typically requiring 96 hr on 1,680 cores, on Gadi super-computing facilities at Australia's National Computational Infrastructure (NCI)), we focus on two demonstrative asymmetric cases.Nonetheless, these cases show how the asymmetric interaction between slabs and remnants alters the subduction process, providing a basis for future studies.
We calculate several diagnostic outputs to quantify slab-remnant interactions.When doing so, the boundary of the slab is defined by the 0.5 contour of the plate material volume fraction (material volume fraction = 1 when the material is plate, 0 otherwise).Based on this contour, we extract the slab tip depth, the trench location and the trailing edge position, as well as rates of slab descent, trench retreat and plate advance.We calculate the average slab dip in the upper mantle from the surface to 650 km depth by approximating the geometry of the slab to be linear and comparing the difference in lateral positions over this depth range, with respect to the direction of gravity at the slab center at 325 km depth.Similarly, remnant depth and velocity are extracted using the 0.5 contour of the remnant material volume fraction.The angular distance between slab and remnant is defined as the minimum difference in longitude between the deepest point of the remnant and the portion of the downgoing plate that is below 250 km depth.Measurements are taken at the symmetry plane unless otherwise specified.We also trace the evolution of trench geometry relative to the initial trench shape.

Results
In the following sections we first investigate the role of remnant density by comparing cases W2400_Neutral and W2400_Wedge with their W2400 reference case.Next, we explore the influence of different remnant locations and orientations, before examining the effect of remnant size on narrow (2,400 km) and wide (4,800 km) subduction zones with a focus on the evolution of trench shape.With a general understanding of how remnant properties influence the evolution of subducting systems from symmetric cases, we evaluate two asymmetrical cases with different remnant location and density, highlighting some of the main effects that remnants can have on subduction dynamics.

Role of Remnant Density
In the reference case (W2400) the plate width is close to Earth's mean present-day trench length (e.g., F. Chen et al., 2022b;Heuret et al., 2007;Müller et al., 2016).At this width, the slab remains reasonably uniform in its along-strike morphology as it descends, with its temporal evolution and corresponding cross-sectional flow field at the symmetry plane displayed in Figure 3a.During the initial phase of subduction, the slab tip steepens, and two poloidal flow cells develop in the mantle wedge and sub-slab regions, over a distance similar to the depth of the upper mantle and length of the slab (as in Piromallo et al. (2006)).These poloidal cells flow from regions of high dynamic pressure toward regions of low dynamic pressure, as illustrated in Figure 4a.During the upper mantle sinking phase, the trench retreats steadily and accounts for ∼30% of total subduction, with the remainder accommodated via trailing edge advance (Figures 5b-5d).Following interaction with the viscosity jump at 660 km depth, flow velocities in the wedge poloidal cell diminish as the slab tip deflects and sinks into the more viscous lower mantle (Figure 3a).At the same time, mantle material flows across the transition zone into the lower mantle, into a region of higher dynamic pressure beneath the slab tip (Figure 4a).At this point, the slab sinking rate reduces from its free upper-mantle Stokes velocity of ∼9 cm/yr to a substantially lower velocity of ∼3 cm/yr (Figures 5a and 5e), as the upper mantle portion of the slab steepens and develops buckling folds (Figures 3a  and 5f).
When a neutrally buoyant remnant is placed in the wedge region at 660 km depth (case W2400_Neutral, red dashed line in Figure 5), we find that slab tip depth, trailing edge advance, slab sinking rate and upper mantle dip angle diagnostics are similar to the W2400 case (Figures 5a, 5c, 5e, and 5f).W2400_Neutral displays a slightly increased trench retreat velocity (Figures 5b and 5d), which is also evident from the increased velocity amplitude of the sub-slab poloidal flow cell (Figures 3a and 3b).Upon slab transition-zone interaction, the viscous remnant in front of the slab tip hinders slab tip advance, leading to more retreat as the slab adjusts to maintain its sinking rate (Figure 5b), resulting in a shorter deflected slab tip.The neutrally buoyant remnant remains close to its initial location until it is pulled into the lower mantle and toward the subducting slab via slab induced flow (Figures 5g, 5h, and 5i), initially sinking at the side closest to the subducting plate (Figure 3b).Despite these subtle differences in the flow regime, the dynamic pressure field in the case with a neutrally buoyant remnant is similar to the reference case (Figure 4b).
Case W2400_Wedge (continuous red line in Figure 5) illustrates the effect of a negatively buoyant remnant in the mantle wedge region above the transition zone.The remnant's negative buoyancy drives downwelling flow, leaving in its wake a region of low dynamic pressure of similar size and strength (Figure 4c) as that of the slab once it penetrates into the lower mantle (Figure 4a).The slab is pulled toward this low pressure area, increasing the slab descent velocity and trailing edge advance rate relative to the reference case (Figures 5a, 5c, and 5e), and forcing the slab into an advancing regime (negative trench retreat) with a steep upper mantle dip angle (Figures 5b, 5d, and 5f).Remnant induced downwelling pulls material across the MTZ long before the arrival of slab material at 660 km depth, which differs from the two previous cases analyzed (Figure 4c), whilst also modifying the mantle wedge and sub-slab poloidal flow cells (Figure 3c).The remnant's initial sinking velocity is similar in magnitude as that of the slab once it enters the lower mantle, as expected from the similar size and density contrast of the slab and remnant.Over the duration of our simulation, the remnants' sinking velocity increases from 3 to 3.5 cm/yr.Beyond the initial phase of subduction, the angular distance between the plate and remnant decreases over time, increasingly aligning the flow fields induced by the slab and remnant (Figures 5h  and 5i).
These cases demonstrate some of the impacts of remnants on the evolution of subduction systems.Even in the absence of buoyancy anomalies, highly viscous material can modulate the flow regime, altering the dynamics and morphology of adjacent descending slabs.Negatively buoyant remnants drive downward flow, aiding the descent of adjacent downgoing slabs.At the same time, slabs aid the descent of remnants.Over time, the remnant and slab move toward each other, resulting in an increase in velocity of both as the lateral distance between them reduces (Figures 5e, 5h, and 5i).

Role of Remnant Location and Orientation
We next compare simulations where identical remnants (geometry, volume, density, and viscosity) are placed in different positions with respect to case W2400, to investigate the effect of remnant location and orientation.Cases W2400_Wedge, W2400_Subslab, and W2400_Under simulate horizontal remnants at varying lateral distances from the initial trench position, and at different positions within the asymmetric flow field set up behind and in front of the subducting slab.We compare cases W2400_Under and W2400_Under_1000 to investigate how the depth of the remnant modulates its influence on the subducting slab.Comparisons between cases W2400_Wedge (horizontally oriented remnant) and W2400_Radial (radially oriented remnant) highlight how remnant orientation influences interaction with adjacent subducting systems.Remnants of similar mass and shape should set up similar Stokes-sinking mantle flow fields around them, but radially oriented remnants would be expected to experience somewhat lower drag than horizontally oriented remnants in the gravity driven flow.Case W2400_Subslab, where the negatively buoyant remnant is placed in the sub-slab region above the transition zone, displays the greatest rate of trench retreat across all 2,400 km wide simulations considered, accounting for ∼40% of total subduction (Figures 5b and 5d).The upper mantle portion displays minimal vertical folding, unlike the reference case where vertical folding is clearly visible (Figure 6d).Although the remnant sinks at a velocity close to the remnant of case W2400_Wedge (Figures 5g and 5h), the slab displays the slowest sinking and trailing edge advance velocities of all 2400-km wide cases considered with negatively buoyant remnants (Figures 5a, 5c, and 5e).The remnants in cases W2400_Subslab and W2400_Wedge are equidistant from the initial trench location, but in case W2400_Subslab it is further away from the slab tip (Figure 5i).Nonetheless, the remnant is able to pull the slab laterally, evidenced through enhanced trench retreat rates (toward the remnant) and a reduced angular distance of ∼4°between remnant and slab (Figure 5i).A negatively buoyant remnant is placed directly below the initial trench location in case W2400_Under.In this scenario, we observe the fastest sinking velocities for both the subducting plate and the remnant across all 2,400km cases considered (Figures 5a, 5e, 5g, and 5h) as the mantle flow fields set up by the slab and remnant most effectively enhance each other.Trench retreat is reduced relative to the reference case, but trailing edge advance is the largest of all cases considered, accounting for more than 80% of total subduction (Figures 5b, 5c, and 5d).The slab has a steep angle (greater than 69°throughout the simulation, Figure 5f) and, as a result, slab tip deflection at the transition zone is reduced (Figure 6e).The angular distance between the slab tip and the remnant reduces initially, as the slab tip rotates toward the remnant; once their angular distance is 0°(i.e., the descending slab is radially above the deepest point of remnant), the slab and remnant move in the same radial direction.We note that when the remnant is placed deeper, 1,000 km below the trench (case W2400_Under_1000), its influence on the slab is reduced, although general trends for each diagnostic are similar to Case W2400_Under.
In case W2400_Radial, the remnant is oriented radially in the mantle wedge between 300 and 700 km depth, 500 km away from the initial trench position.The orientation and location of the remnant in closer proximity to the descending slab and a somewhat stronger remnant flow field drive an increase in trench advance and the slab's upper mantle dip angle compared with the W2400_Wedge case.Indeed, of all 2,400-km wide cases considered, this displays the most trench advance and the steepest upper mantle slab dip angle (Figures 5b, 5d, and 5f), albeit with a similar slab sinking velocity and plate advance velocity to its horizontal counterpart, W2400_Wedge (Figures 5a, 5c, and 5e).The angular distance between the slab and remnant remains stable after slab transition zone interaction, whereas in the horizontal W2400 Wedge case, the angular distance continues to reduce (Figure 5i).At this stage of model evolution, the radial remnant and slab are descending in their respective radial directions toward the center of the sphere, hence there is no longer any effective lateral flow pulling them toward each other (Figure 6b).Horizontal remnants in W2400_Wedge and W2400_Subslab cases rotate toward the descending slab so that their angular distances reduce while the slab also moves toward the remnant either through advance or retreat.
The distance between the remnant and the subducting slab plays an important role in determining remnant impact on the subduction process.As shown in Figure 7a, negatively buoyant remnants can increase slab sinking velocities even when displaced horizontally ∼1,000 km from the slab tip, with an increase in sinking velocity of ∼0.5 to ∼1.0 cm/yr relative to the reference case in Case W2400_Subslab, at an angular distance of 8°.Cases W2400_Wedge, W2400_Under, and W2400_Under_1000 demonstrate that the increase in slab sinking velocity relative to the reference case is proportional to the proximity of a negatively buoyant remnant, with relative sinking velocities generally increasing as the slab-remnant distance decreases, as slab and remnant flow fields align more closely.Where distances between the remnant and slab are less than 100-200 km, sinking velocities can be enhanced by up to a factor of 2 (see also Figure 5e).A comparable effect is also measured for the remnant, as shown in Figure 7b, where all remnants exhibit higher sinking velocities as they approach the subducting slab.
These cases demonstrate that, for a given buoyancy anomaly, the proximity of a remnant to a subducting slab (both laterally and radially) controls its ability to accelerate slab descent.Furthermore, more energy is expended on horizontal transportation and/or rotation of the downgoing slab when the lateral distance between remnant and slab increases.If the remnant is located in front of a subducting slab, the slab tends to advance and steepen, but if the remnant is located behind the subducting slab, beneath the downgoing plate, the trench retreats to reduce the distance between slab and remnant.Downgoing flow from the descending slab also rotates horizontal remnants toward the slab, reducing the angular distance between them.These motions lead to increasing alignment of the flow set up by the slab and remnant.As a result, when the remnant and slab are oriented radially in the direction of gravity, further rotations or lateral movements are minimal, with their angular distance unchanged for the remainder of their descent.

Role of Remnant Size and Slab Width
Subduction zones develop curvatures at their edges due to toroidal flow around the slab.The evolution of trench shape depends on plate age and width: narrow plates tend to develop "C"-shaped trenches, where the trench retreats most at the center, while wider plates tend toward developing "W"-shaped curvatures, where trench retreat is lower at the center and the edges of the slab (e.g., F. Chen et al., 2022b;Schellart et al., 2007).Here, we examine the influence of remnants for cases with plate widths of 2,400 and 4,800 km, which develop "C" and "W"-shaped trenches in their reference cases, respectively (W2400 and W4800: Figures 8 and 10).Geochemistry, Geophysics, Geosystems 10.1029/2023GC011180 CHEN ET AL.

Influence on 2,400-km Wide Slabs
Our previous work demonstrates that trenches consistently develop an elongated "C"-shape for 2,400 km wide slabs (F.Chen et al., 2022b).In cases W2400_Wedge and W2400_Subslab, the addition of a full-width remnant does not modify predicted "C"-shaped trench curvature (Figure 8), despite influencing trench motions, as discussed in Section 3.2.
To examine the influence of smaller remnants, we investigate two additional cases with remnant widths that are half of the plate width, positioned beneath the center of the trench.The half-width remnant in Case W2400_Wedge_half has a lower sinking velocity (maximum of ∼3.2 cm/yr) than its full-width remnant counterpart (maximum of ∼3.6 cm/yr-Figure 9h) and, hence, drives less of an increase in slab descent velocity.As a result, relative to its full-width equivalent, the slab sinks slower with less trailing edge advance, at a shallower dip angle (Figures 9a, 9c, and 9f).The half-width remnant is also less able to pull the slab toward it.As a result, despite significantly reduced trench retreat rates compared to case W2400, the trench remains in a retreating regime, rather than entering the advancing regime of case W2400_Wedge (Figures 9b and 9d).While the angular distance between the half-width remnant and the center of the plate reduces (Figure 9i), no significant curvature develops along the trench (Figure 8).
The half-width remnant in Case W2400_Subslab_half drives a similar response to Case W2400_Wedge_half, with its influence on the subducting plate less than its corresponding full-width case, W2400_Subslab.Due to the lower remnant sinking velocity (maximum of ∼3.3 cm/yr vs. ∼3.9cm/yr in W2400_Subslab), the slab in W2400_Subslab_half sinks slower, whilst trench retreat and trailing edge advance also reduce (Figure 9), with a negligible influence on trench curvature.
Thus, smaller remnants have less of an influence on the descending slab than full-width remnants.This is in line with expectations from their Stokes sinking velocities.On 2,400-km wide cases that tend to develop "C"-shaped trenches, such remnants are unable to generate sufficient pull to overcome the bending resistance of the subducting plate, with trench shape remaining similar to corresponding reference cases.

Influence on 4,800-km Wide Slabs
At a plate width of 4,800 km, the reference case (W4800) develops a "W"-shaped trench ("S"-shaped in the halved domain, Figure 10) with the center of the trench stagnating at its initial location.In this section, we examine the effects of full-width remnants via cases W4800_Under, W4800_Wedge, and W4800_Subslab; and half-width remnants via cases W4800_Wedge_half and W4800_Subslab_half.
Similar to case W2400_Under, where the remnant is placed directly beneath the trench, case W4800_Under also displays the fastest descent velocities for both the subducting slab and the remnant over all 4,800 km wide cases considered (Figures 11a,11e,11g,and 11h).As in its 2,400 km wide equivalent, this case displays elevated trailing edge advance and a steeper upper mantle dip angle compared to its reference case (Figures 11c and 11f).
The trench advances from its initial location (Figure 11b), with the "W"-shaped curvature less prominent compared to the reference case (Figures 10, 12a, and 12b).The full-width remnant in case W4800_Wedge drives significant trench advance and steepening of the subducting slab (Figures 11b,11d,and 11f).Trench curvature is reduced relative to the reference case (Figure 10), and slab morphology becomes more uniform along strike (Figure 12c).The half-width remnant (case W4800_Wedge_half) is less able to drive trench advance at the symmetry plane (Figures 11b and 11d).However, as this anomaly is placed below the stagnating region of the trench, it enhances the contrast in trench curvature between the advancing center and the retreating edges of the slab (Figure 10), driving substantial along-strike variations in slab morphology (Figure 12d).
Sub-slab remnants in cases W4800_Subslab and W4800_Subslab_half increase trench retreat while reducing the upper mantle slab dip angle (Figures 11b,11d,and 11f).On the symmetry plane, full-width and half-width remnants have similar effects on slab dynamics, except that the angular distance between the slab and remnant reduces more for the full-width remnant.The full-width remnant reduces "W"-shaped trench curvature, but not to the same extent as the half-width remnant, which almost generates an elongated "C"-shaped trench with a straight center (Figure 10).The localized retreat-driving flow induced by the half-width sub-slab remnant counteracts the subducting plate's natural tendency to stagnate at the center, limiting along-strike variations in trench retreat velocity.
In summary, the trench shape of wider subduction systems are more strongly influenced by remnants than narrower plates.Because of the larger width, the effect of toroidal flow around the slab is limited to regions within 1,000-1,500 km from the slab edges.When additional descent-driving forces act uniformly along strike, the relative impact of toroidal flow at the edge reduces, resulting in less trench curvature.When a half-width remnant drives the slab center to retreat, counteracting the natural tendency of the trench to stagnate, the trench remains relatively straight.Conversely, when the half-width remnant acts to enhance slab stagnation or advance at the center, the trench is able to develop increased curvature.For narrower plates, these different driving forces act in close proximity along the trench, and are typically unable to overcome the plate bending resistance and, hence, these slabs remain in a "C"-shape, regardless of remnant size and location.

Asymmetric Remnants
While the hemispherical cases allow us to systematically explore how remnants with different properties interact with subducting slabs, the full spherical models considered in this section allow us to demonstrate asymmetric slab interactions with neutrally and negatively buoyant remnants via cases Asym_Under_Neutral and Asym_-Wedge_Negative. Case Asym_Under_Neutral is designed to examine the effect of an active slab partially subducting on to a stagnant remnant at the transition zone, where the remnant is not actively sinking and driving mantle flow.Case Asym_Wedge_Negative allows us to illustrate a more dramatic influence of a remnant, where it actively drives a downwelling flow cell in front of the subduction zone segment that displays the most trench retreat in the corresponding reference case.
For Case Asym_Under_Neutral, where half of the subducting slab descends into a neutrally buoyant remnant, slab descent is reduced locally, as evidenced by the along-strike variations in slab tip depth illustrated in Figures 13a, 14a, and 14c: the remnant's high viscosity hinders slab sinking, with the slab deflecting and resting above the remnant.At the edge of the remnant, the slab tip's depth transitions to the same state as the case without a remnant, penetrating deeper.These variations drive small changes in trench morphology, with more bulging observed toward the trench center, and the region of most advance/stagnation shifted toward the remnant side when compared to the reference case (Figure 13c).
In case Asym_Wedge_Negative, the negatively buoyant remnant drives an increase in slab sinking velocities (Figures 13b, 14b, and 14d) and strongly enhances trench advance along the proximal segment of the slab (Figure 13c).This highlights the asymmetric response: not only is trench advance enhanced, but the location of most advance shifts from the center of the subducting plate (as in case W4800) toward the remnant (Figure 13c), illustrating the remnant's ability to alter trench shape and drive substantial along-strike variations in slab morphology (Figures 14e and 14f).To demonstrate how such a complex system evolves, we extended the simulation time of this case to 21 Myr: the system continues to develop more complex trench curvature and associated slab morphology as subduction evolves.

Applicability of Models
Our models are a simplified representation of Earth's subduction systems, but they provide fundamental insight into the first-order effects that slab remnants have on active subduction zones, and vice-versa.They are executed  and 21 Myr (dashed lines).Distance along trench is defined as the distance from the symmetry plane in the direction parallel to the initial trench location.Blue boxes mark the initial trench section that is underlain by the remnant (they share the same latitude), note that the width of the remnant extends beyond the plotted region in (b, c) trench shapes and the amount of retreat for cases W4800, Asym_Under_Neutral, and Asym_Wedge_Negative at 12 and 21 Myr.Case W4800 was simulated in the hemispherical domain, thus for comparative purposes, the trench profile is mirrored at the symmetry plane. in a 3-D spherical shell domain, which is important for simulating Earth's subduction systems, particularly those that exceed 2,400 km in width (F.Chen et al., 2022a).They also better capture the mechanical properties of subducting slabs than previous global models of mantle convection that incorporate slab remnants (e.g., Becker & Faccenna, 2011;Bunge et al., 2002;D. R. Davies et al., 2012D. R. Davies et al., , 2015;;Lowman, 2011;Yanagisawa et al., 2010).
There are, however, simplifications of our model setup that may influence the evolution of slab morphologies and velocities.For example, the MTZ is simulated solely as a viscosity jump that hampers flow into the lower mantle (e.g., Hager & Richards, 1989).We do not account for phase transitions, which will (temporarily) further resist material transfer across this boundary (e.g., Goes et al., 2017;Ito & Yamada, 1982;Ringwood, 1975).Such a simplification may lead to our models predicting less transition-zone slab stagnation.
Another limitation of our model set-up is the lack of an overriding plate, which will influence predicted trench migration rates and associated trench geometries and slab morphologies (e.g., Capitanio, Stegman, et al., 2010;Garel et al., 2014;Heuret et al., 2007;Jarrard, 1986;Lallemand et al., 2005;van Dinther et al., 2010).Our neglect of a composite rheology in surrounding mantle (particularly the influence of dislocation creep) may also impact some of the properties analyzed, for example, affecting slab dips and trench retreat (e.g., Billen & Hirth, 2007;Holt et al., 2017).In addition, given the computational requirements of our simulations, we have only been able to examine a set of cases where slabs have a consistent viscosity contrast relative to background mantle.As the force imposed from a sinking remnant would induce more deformation in a weaker slab (lower bending resistance) it is important to analyze the sensitivity of results to this choice in future work, in terms of both the evolution of slab morphologies and along-strike variations in trench curvature.
Finally, our model is incompressible and lacks a temperature field.The thermal structure of a subducting plate controls its thickness, density and rheology.While we use visco-plastic upper and lower layers of the subducting plate to mimic the deformation predicted by thermo-mechanical subduction models with non-linear temperature and stress dependent rheology, our simulations do not capture adiabatic heating, viscous dissipation, the diffusion of temperature, and variations in thickness and buoyancy from ridge to trench (e.g., Agrusta et al., 2017;Capitanio, 2013;F. Chen et al., 2022b;Garel et al., 2014;Lee & King, 2009;OzBench et al., 2008;Stegman et al., 2010).Furthermore, on Earth, the rheology of slabs and remnants changes with temperature, such that remnants at depth could become less plate-like than those considered herein (e.g., Andrews & Billen, 2009;Stadler et al., 2010;van Hunen et al., 2001).Nonetheless, assuming a thermal diffusivity of 1.0 × 10 6 m 2 s 1 , the time it takes for a 70 km thick slab to cool substantially would be ∼40 Myr.At an upper mantle sinking speed of 5-10 cm/yr (consistent with predictions herein and observational constraints on Cenozoic plate motions, Goes et al., 2011), a slab crosses the upper mantle in 6-12 Myrs and thus retains most of its rheological contrast with background mantle.Buoyancy is modified even less, because while the thermal signature diffuses as the slab sinks, the integrated negative buoyancy only decreases if parts of the slab are sufficiently weakened that they can be eroded or removed by flow.Consequently, over the depth and time range investigated here, the impact of thermal diffusion on slab-remnant interaction will be secondary to the effects of density and strength contrasts that our models capture (e.g., Jarvis & Lowman, 2007;Kundu & Santosh, 2011;Quéré et al., 2013).We are therefore confident that our chosen model design, alongside the systematic comparisons across our model cases, isolates how different remnant properties (e.g., density, size, location) influence slab dynamics.

Potential Influences of Remnants on Subduction Zones
Based on the importance of slab pull for driving surface plate motions (e.g., Becker & O'Connell, 2001;Forsyth & Uyeda, 1975;Lithgow-Bertelloni & Richards, 1998), alongside the fact that regional subduction models predict an important role for incoming plate age and strength in dictating the dynamics of subduction (e.g., Bellahsen et al., 2005;Garel et al., 2014;Stegman et al., 2010;Suchoy et al., 2021), many studies have examined the correlations between different subduction parameters (e.g., Cruciani et al., 2005;Doglioni et al., 2007;Heuret & Lallemand, 2005;Lallemand et al., 2005;Verard, 2019;Vérard et al., 2015) in an attempt to reveal what drives the observed diversity of subduction tectonics.While global correlations between some observed subduction features can be established, such as slab dip and back-arc deformation with upper plate motion and upper plate thickness, be it with notable exceptions (e.g., Heuret & Lallemand, 2005;Lallemand et al., 2005), many other parameters have poor correlations (e.g., Cruciani et al., 2005;Doglioni et al., 2007;Vérard et al., 2015).Our results demonstrate the ability of slab remnants to significantly enhance sinking velocities and affect trench motion.It is therefore likely that the effect of mantle flow driven by subducted remnants contributes to an explanation for these poor correlations.
One of our key findings is that the downwelling flow generated by negatively buoyant slab remnants can enhance the descent of nearby subducting plates.Numerical models have suggested that slab sinking velocity should increase with plate age, which enhances slab pull (e.g., Agrusta et al., 2017;F. Chen et al., 2022b;Garel et al., 2014;Goes et al., 2011).Nonetheless, despite there being some positive velocity-age trend for larger plates (e.g., Goes et al., 2011;Lallemand et al., 2005), overall there is a poor correlation between subducting plate velocities and age (e.g., in the recent global plate reconstructions of Müller et al., 2016Müller et al., , 2019)).This suggests that there are factors, other than age, impacting the velocities of subducting plates, one of which may be the interaction with nearby remnants or ongoing subducting slabs.As shown in Figures 5, 9, and 11 (panel e), a nearby negatively buoyant remnant can increase slab sinking velocities.Such slab-remnant interactions would complicate the trend of measured plate velocities with subducting plate age.Likewise, interaction with remnants may provide a plausible explanation for slabs currently residing at a deeper depth than expected from tectonic reconstructions, such as the Arabian and Kalimantan slabs, as noted by van der Meer et al. (2018).
In addition to increasing slab sinking velocities, the downwelling flow driven by negatively buoyant remnants could potentially help to sustain or initiate subduction.Our results suggest that slabs and remnants tend to move toward each other, which may help to anchor and sustain subduction at a given location for a prolonged period of time.There is a strong correlation between subduction zone initiation events and the presence of nearby subduction (Baes & Sobolev, 2017;Crameri et al., 2020).Pysklywec and Ishii (2005) used 2-D numerical models to demonstrate that slab remnants may trigger slab detachment, and subsequent initiation of subduction of the opposite polarity.Nikolaeva et al. (2010) modeled subduction initiation at passive margins, and suggested the criteria for spontaneous initiation are hard to achieve naturally.Even in other settings, it seems subduction initiation usually requires additional forcing (e.g., Lallemand & Arcay, 2021).Capitanio and Replumaz (2013) modeled slab break-off, and showed that while break-off episodes provide short-lived and localized large stresses in the upper plate interior, the detached slab remnants sustained the subduction dynamics that drive convergent motion.Such dynamics could explain long-lived under thrusting at the India-Asia convergence zone and episodic lithospheric faulting in the Asian continent.Significant downwelling from detached slabs may also facilitate the continuation of subduction.For example, the remnant of the Izanagi plate may have provided the driving force for Pacific Plate subduction at the East Asian margin after the sub-parallel subduction of the Izanagi-Pacific ridge at ∼50 Ma (e.g., Seton et al., 2015;Wu & Wu, 2019).
Our models also show that remnants can influence trench shape.Observed along-strike variations in trench shape have been attributed to different factors: (a) the influence of upper plate heterogeneity (e.g., Arnulf et al., 2022;Capitanio, Stegman, et al., 2010); (b) the impact of subducting a buoyant anomaly, such as ridges and oceanic plateaus, which resist subduction and trench retreat (e.g., Mason et al., 2010;Suchoy et al., 2022); and (c) plate width, where the center of a wider plate tends to retreat less than the edges, or even stagnate, resulting in a "W"shaped trench (e.g., F. Chen et al., 2022b;Schellart et al., 2007).Our models suggest that depending on a remnants' location relative to the subducting plate, it can either enhance trench curvature (as in case W4800_Wedge_half, Figure 10), reduce trench curvature (most notably in case W4800_Subslab_half, Figure 10), induce trench rotation toward the remnant (case Asym_Wedge_Negative, Figure 13c), or influence the location of the protruding stagnation zone along the trench (case Asym_Wedge_Negative, Figure 13c).Whilst both slab remnants and buoyant ridges could lead to localized convex curvature at the trench, there is a clear distinction between them: at the protruding center of the trench induced by a remnant, the slab is advancing forward with a steep dip (Figures 11b and 11f).This differs to the stagnating trench shape induced by the subduction of a buoyant anomaly, which is associated with shallower dip angles (Suchoy et al., 2022).
Our models suggest that the coupled flow from negatively buoyant remnants and the downgoing subducting slab acts to reduce the distance between them, primarily via rotation.This rotation is seen in the increasing slab dip (panel f in Figures 5, 9, and 11) and remnant rotation in Figures 6d and 6e.However, such a rotation is limited by the direction of gravity: when the remnant and slab are descending in a radial orientation aligned with gravity, they cannot rotate any further, as illustrated by case W2400_Radial, and they do not further approach.This contributes to an explanation of why many distinct slab fragments are still imaged in the mantle (e.g., Braszus et al., 2021;Ren et al., 2007;van der Meer et al., 2018;Wu & Suppe, 2018).
Finally, it has long been recognized that super-continents assemble and break up episodically throughout Earth's history, and this cycle is intimately linked to whole-mantle convection (e.g., Mitchell et al., 2021;Nance et al., 1988Nance et al., , 2014;;Rogers & Santosh, 2003;Rolf et al., 2014).In particular, the assembly stage of super-continent cycles is heavily influenced by subduction.Collins (2003) suggested that long-lived slab pull forces controlled Pangaea assembly and dispersal, whilst Santosh et al. (2009) suggested that the process of super continent assembly is driven by super-downwelling through double-sided subduction as seen in the Western Pacific.Our models suggest that the downwelling flow generated by slab remnants can concentrate the locations of slabs, leading to the potential development of super-downwellings.The negative buoyancy of these accumulated remnants may also aid subduction initiation of overlying oceanic lithosphere, which further concentrates negative buoyancy and feedbacks into potential super-downwelling systems.Therefore, based on our results, we speculate that the location and volume of remnant slabs in the mantle may be a crucial factor controlling ongoing and new subduction zones, global plate reorganization events and super-continent cycles.

Potential Examples of Slab-Remnant Interactions
Seismic tomography studies (e.g., Li et al., 2008;van der Meer et al., 2018) show an abundance of remnants in the mantle within close proximity to active subduction zones (see Figure 1).These regions typically have complex tectonic histories, and the results from our models may contribute toward an improved understanding of how the evolution of these subduction systems have been shaped by interactions with nearby remnants.Given that such remnants are scattered around the globe and are capable of affecting subduction zones that range in scale from the large (e.g., Tethyan and Farallon subduction) to the small (e.g., South-East Asian subduction zones), it is important to understand slab-remnant interactions.Below, we highlight a few examples of where slab-remnant interaction may have affected tectonic evolution.
Closure of the Tethys Oceans left substantial subducted remnants below a region extending from the Mediterranean subduction zones to the India-Eurasian collision.Tomographic imaging shows that slab remnants from Tethyan subduction under India are located near the ancient locations where they began subduction, and that Tethyan slab remnants are largely above previous subducted fragments, implying only small amounts of lateral movement, which is an indication of an anchored and ongoing subduction system (e.g., Hafkenscheid et al., 2006;van der Voo et al., 1999) (Figure 1a).This is similar to our W2400_Under model, where the downwelling flow from the pre-existing remnant reduces the distance between the remnant and the subducting slab and effectively pins the location of subduction (Figure 5i).Capitanio, Morra, et al. (2010) examined the force balance required to drive reconstructed plate motions associated with India-Asia convergence, demonstrating that these are three times higher than what can be explained by this subduction system in isolation.They suggest that a component of this convergence is driven by subduction of the dense Indian continental slab, and also invoke forcing related to lower mantle suction created by sinking Paleo-Tethys slabs.Becker and Faccenna (2011) further suggest that this slab induced lower mantle suction combines with active mantle upwelling to create a "conveyor belt," facilitating the ongoing indentation of the Indian and Arabian plates into Eurasia.Furthermore, the flow patterns generated by the relics of west Tethyan subduction may have aided Mediterranean subduction and its rollback (Faccenna et al., 2014).These studies, considered alongside our results, demonstrate a key role for remnant-induced mantle flow in the tectonic forcebalance.
Seismic tomography models illuminate remnants of the Farallon and possibly other plates in and below the MTZ below the North America plate (Goes & van der Lee, 2002;Ren et al., 2007;Sigloch et al., 2008;van der Lee & Nolet, 1997) (Figure 1d).These remnants could have driven downwelling flow that aided scenarios where multiple basins subducted on top of each other, generating large-scale tomographic anomalies interpreted as vertical "slab walls," and facilitating proposed subduction polarity flips (Sigloch & Mihalynuk, 2013).Remnants currently in the lower mantle below the northwestern United States have been attributed to subduction that started at an intra-oceanic arc west of the North American continental margin (e.g., Ren et al., 2007;Sigloch & Mihalynuk, 2013) and, although other factors likely contributed to the westward motion of the North America plate (e.g., Müller et al., 2019), the flow generated by these slab remnants could have aided such motion.
The Antilles and South American subduction zones are tectonically complex (Figure 1c).The Antilles slab appears to be highly fragmented (e.g., Bezada et al., 2010;Braszus et al., 2021;van Benthem et al., 2013) and is affected by motions of the major surrounding plates.It is suggested that some Antillean slab fragments are currently located in the lower mantle beneath northeastern South America (e.g., Braszus et al., 2021;van Benthem et al., 2013).Other subducted remnants, likely from Farallon subduction, have also been identified in the MTZ and lower mantle below southeastern North America (e.g., Bunge & Grand, 2000;Ren et al., 2007;Sigloch et al., 2008).Flow associated with these relic slabs could well have facilitated subduction and the penetration of Farallon and Nazca slabs into the lower mantle beneath Central America and north-central South America (e.g., van der Meer et al., 2018).Furthermore, in conjunction with the effects of plate width (e.g., F. Chen et al., 2022b;Schellart, 2017) and the subduction of buoyant anomalies (e.g., Gutscher et al., 1999;Suchoy et al., 2022), slab fragments below northern South America could have enhanced oroclinal bending of the trench at Bolivia, similar to the enhanced trench curvature predicted for case W4800_Wedge_half herein.
South East Asia has a complex tectonic history and, as a result, a complex underlying mantle structure (e.g., Hall & Spakman, 2015;Pilia et al., 2023;Replumaz et al., 2004;van der Meer et al., 2018) (Figure 1f).Material subducted along the former Banda trench forms a large flat lying remnant just below the MTZ (Spakman & Hall, 2010); and despite the "stalled" appearance of the slab, its stagnation at 660 km discontinuity may well be transient (Goes et al., 2017), in which case its negative buoyancy likely enhances subduction and sinking of the many small slabs below northeastern Indonesia and the Philippines.The remnant in the sub-slab region of the lower mantle beneath Sumatra-Andaman has been associated with Tethyan subduction (e.g., Hafkenscheid et al., 2006;Spakman & Hall, 2010;Widiyantoro & van der Hilst, 1997); and based on our results, a remnant in such a location relative to the subducting plate may facilitate trench retreat and back-arc opening behind Andaman.
Below the Tonga-New Hebrides plate boundary, seismic tomography models suggest distinct 3-D structures that represent slab remnants from past subduction at the Melanesia Arc (e.g., W.-P. Chen & Brudzinski, 2001;Hall & Spakman, 2002) (Figure 1e).As occurs in our models, active subduction in this region is likely affected by the negative buoyancy of these remnants (Pysklywec & Ishii, 2005;Pysklywec et al., 2003).Variable sinking velocities along the subduction zone, suggested based on tectonic reconstructions (Schellart & Spakman, 2012), could be an indicator of a non-uniform, asymmetric influence from remnants, similar to the case Asym_-Wedge_Negative analyzed herein.

Conclusions
We have presented a suite of subduction models in a 3-D hemispherical shell domain, investigating how different slab remnant properties-density, size, location and orientation-influence interaction with a nearby subduction zone.We have also presented two full spherical models with an asymmetric setup, to illustrate the effect of remnants in driving asymmetric along-strike variations in trench shape and slab morphology.
Our models show that downwelling flow generated by negatively buoyant slab remnants can be of similar scale and magnitude as that of the subducting slab and when located within a few 100-1,000 s of km from the slab tip, this flow enhances the sinking velocity of nearby actively subducting slabs by up to a factor 2 (depending on remnant size and location).The joint effects of remnant downwelling and slab pull may explain the observed poor correlations between subducting plate velocities and the age of the subducting lithosphere.
The location of a remnant relative to a nearby subduction system can be an important factor controlling the evolution of trench shape.Sinking remnants and subducting slabs move toward each other via rotation, leading to increasing alignment of the flow set up by the slab and remnant.Remnants located in the sub-slab region, under the incoming plate, tend enhance trench retreat, whereas remnants on the mantle wedge side facilitate trench advance.For wide subduction zones in particular, which develop a convex trench stagnation/advance zone at their center, the mantle flow field generated by a sinking remnant can influence along-strike trench shape variation and potentially drive the convex stagnation point away from the center of the trench toward the remnant.This process may have contributed to trench evolution in regions such as the Bolivian Orocline above the Nazca subduction system, alongside other influences, such as subduction of buoyant anomalies (e.g., Suchoy et al., 2022).
With the abundance of relic subduction fragments that have been identified in mantle tomography, and examples of interactions with subduction zones around the globe, we suggest that slab-remnant interaction is an important process that influences subduction dynamics on Earth, significantly affecting plate velocities and trench shapes.
Remnants likely also help to anchor and sustain subduction systems, and facilitate subduction initiation events that drive large-scale plate reorganizations and super-continent cycles.

Figure 2 .
Figure 2. Simulation setup for (a-e) symmetric cases modeled in a hemispherical shell domain and (f-h) asymmetric cases modeled in a full spherical shell domain: (a) Hemispherical shell domain setup, where the domain is bounded by the symmetry plane of the system, whilst bottom and top (inner and outer) boundaries approximate Earth's core-mantleboundary and surface, respectively; (b) Initial slab tip geometry of our layered visco-plastic plates; (c) Front view of fullwidth and half-width remnant configurations with respect to the width of the subducting plate; (d) Side view of positions of horizontally oriented remnants, in the sub-slab region, under the initial trench location, and in the mantle wedge region; (e) Side view of the configuration of a radially oriented remnant; (f) Spherical shell domain setup, with two side plates covering the domain adjacent to the subducting plate; (g) Front view of remnant location in case Asym_Under_Neutral; (h) Front view of remnant location in case Asym_Wedge_Negative.

Figure 3 .
Figure 3. Snapshots illustrating the spatio-temporal evolution of slab morphology through the viscosity field (top), and poloidal flow cells in the mantle wedge and subslab regions (bottom) at simulation times of 4 8 and 12 Myr for cases: (a) W2400; (b) W2400_Neutral; and (c) W2400_Wedge, respectively.The largest arrow in the bottom panels represents a velocity magnitude of 9.6 cm/yr.

Figure 4 .
Figure 4. Evolution of dynamic pressure field at the symmetry plane at simulation times of 2, 4, 8, and 12 Myr, respectively, for cases: (a) W2400; (b) W2400_Neutral; and (c) W2400_Wedge.Arrows indicate the direction and magnitude of velocity up to 9.6 cm/yr.

Figure 5 .
Figure 5. Comparisons between simulations of a 2,400 km wide plate with 2,400 km wide remnants with different properties, locations and orientations: (a) slab tip depth as a function of time, where the upper-lower mantle boundary is indicated by the black dotted line at 660 km depth; (b) amount of trench retreat; (c) amount of plate advance, measured at the plate's trailing edge; (d) ratio of trench retreat to total descent, which is the sum of trench retreat and trailing edge advance; (e) slab sinking velocity; (f) average slab dip in the upper mantle, with the black dashed line indicating a vertical slab with dip angle of 90°; (g) remnant depth; (h) remnant sinking velocity; and (i) the minimum angular distance between subducting slab tip and the deepest point of the remnant.Triangles indicate the time of first slab tip interaction with the viscosity jump at 660 km depth.All measurements are taken at the symmetry plane.

Figure 7 .
Figure 7. Comparisons of: (a) differences in sinking velocity relative to the reference W2400 case (measured when the slab tip is at the same depth)-values above 0 indicate that the presence of the remnant increases slab sinking velocities; (b) remnant sinking velocities, as a function of angular distance between remnant and slab tip.

Figure 9 .
Figure 9. Comparisons between simulations of a 2,400 km wide plate with remnants of different widths and positions: (a) slab tip depth as a function of time, where the upper-lower mantle boundary is indicated by the black dotted line at 660 km depth; (b) amount of trench retreat; (c) amount of plate advance, measured at the plate's trailing edge; (d) ratio of trench retreat to total descent, which is the sum of trench retreat and trailing edge advance; (e) slab sinking velocity; (f) average slab dip in the upper mantle, with the black dashed line indicating a vertical slab with dip angle of 90°; (g) remnant depth; (h) remnant sinking velocity; and (i) the minimum angular distance between subducting slab tip and the deepest point of the remnant.Triangles indicate the time of first slab tip transition-zone interaction.All measurements are taken at the symmetry plane.

Figure 11 .
Figure 11.Comparisons between simulations of a 4,800 km wide plate with remnants of different widths and positions: (a) slab tip depth as a function of time, where the upper-lower mantle boundary is indicated by the black dotted line at 660 km depth; (b) amount of trench retreat; (c) amount of plate advance, measured at the plate's trailing edge; (d) ratio of trench retreat to total descent, which is the sum of trench retreat and trailing edge advance; (e) slab sinking velocity; (f) average slab dip in the upper mantle, with the black dashed line indicating a vertical slab with dip angle of 90°; (g) remnant depth; (h) remnant sinking velocity; and (i) the minimum angular distance between subducting slab tip and the deepest point of the remnant.Triangles indicate the time of first slab tip transition-zone interaction.All measurements are taken at the symmetry plane.

Figure 13 .
Figure 13.Along-strike variations of asymmetric spherical remnant models: (a, b)-variations in slab tip depth along the subducting slab for cases Asym_Under_Neutral and Asym_Wedge_Negative, respectively, at times of 12 Myr (solid lines)and 21 Myr (dashed lines).Distance along trench is defined as the distance from the symmetry plane in the direction parallel to the initial trench location.Blue boxes mark the initial trench section that is underlain by the remnant (they share the same latitude), note that the width of the remnant extends beyond the plotted region in (b, c) trench shapes and the amount of retreat for cases W4800, Asym_Under_Neutral, and Asym_Wedge_Negative at 12 and 21 Myr.Case W4800 was simulated in the hemispherical domain, thus for comparative purposes, the trench profile is mirrored at the symmetry plane.

Figure 14 .
Figure 14.Slab and remnant morphology at 12 Myr for cases Asym_Under_Neutral (a, c) and Asym_Wedge_Negative (b, d), and at 21 Myr for the case Asym_Wedge_Negative (e, f) from different orientations.The subducting slabs are colored by depth, with different colorbars used for at 12 Myr (a-d) and 21 Myr (e, f) to highlight variations in slab tip depth.The remnant is colored in blue with the lower mantle in green.

Table 1
Parameters Common to All Simulations

Table 2
Simulations Examined and Associated Parameter Values Note. w = slab width; W rem = remnant width; D rem = remnant depth; Δρ rem = density contrast between remnant and ambient mantle.