The impact of along-strike variations in overriding plate thickness on subduction dynamics and mantle flow

Subduction zones are inherently three-dimensional systems and exhibit pronounced trench-parallel variability in key observables, including the deformation style of the overriding plate, trench migration rates, slab geometry, and mantle flow patterns. Geodynamic models typically invoke external mantle flow and/or along-strike variations in the properties of the subducting slab to explain this variability, often neglecting the influence of the overriding plate, despite growing evidence of its strong control on subduction dynamics. In this study, we use self-consistent three-dimensional numerical models to explore how along-strike heterogeneities in the overriding plate structure can generate significant variations in subduction dynamics and mantle flow. Our results demonstrate that trench-parallel variations in overriding plate thickness produce large along-strike differences in trench retreat velocities, leading to strongly arcuate trench geometries. We further conduct a suite of models incorporating a mechanically weak zone in the subducting plate, representing the subduction of a transform fault oriented perpendicular to the trench. These experiments show that along-strike variations in overriding plate thickness can promote vertical slab tearing and segmentation of the subduction system into distinct slab segments. Slab tearing facilitates focused mantle upwelling through the tear, potentially triggering tear-related magmatism during slab rollback. Natural examples of subduction zones characterized by vertical slab tears include the Melanesian subduction system, the South Shetland margin and the Tyrrhenian-Apennines collision system.

Subduction is the core process of Earth’s plate tectonics, driving the surface plate movement and mantle flow through geological time scales (e.g., Forsyth and Uyeda, 1975; Davies and Richards, 1992; Becker and O’Connell, 2001). Numerical and analogue models of symmetric subduction zones, in which the subducting and the overriding plates have no heterogeneities along the trench-parallel direction, have greatly contributed to the understanding of subduction zone dynamics (e.g., Funiciello et al., 2004; Piromallo et al., 2006; Schellart et al., 2007; Yamato et al., 2009; Stegman et al., 2010; Billen and Jadamec, 2012; Chen et al., 2015; Holt and Becker, 2016; Schliffke et al., 2022). Assuming laterally homogeneous plates has made possible to isolate the effects of other parameters on subduction dynamics, such as the slab width (e.g., Stegman et al., 2006; Schellart et al., 2007; Guillaume et al., 2010; Strak and Schellart, 2016; Chen et al., 2022) and thickness (e.g., Bellahsen et al., 2005; Funiciello et al., 2008; Stegman et al., 2010; Capitanio and Morra, 2012; Gea et al., 2024) or the mantle rheology (e.g., Holt and Becker, 2016; Pusok et al., 2018). However, subduction zones on Earth present a large number of variabilities along the trench-parallel direction, such as the trench shape, the mantle flow, the slab depth and dip angle and the deformation style and structure of the overriding plate. For example, the Melanesia subduction zone (Schellart et al., 2006) shows relatively little overriding plate deformation at its centre (in front of the San Cristobal Trench), but the edge segments (New Britain and New Hebrides) are characterised by backarc extension (Schellart et al., 2007; Müller et al., 2019). Similarly, the Middle America subduction zone shows regions of both extension and compression along-strike in its overriding North American plate (Müller et al., 2019).

To explain these variations along-strike, more recent numerical and analogue models have relied on lateral heterogeneities in the properties of the subducting plate. Most of these studies have focused on lateral variations on slab buoyancy, caused by oceanic plateaus and aseismic ridges (e.g., Martinod et al., 2013; Strak and Schellart, 2018; Flórez-Rodríguez et al., 2019; Suchoy et al., 2022; Zhao and Leng, 2023) or contiguous continental plates (van Hunen and Allen, 2011; Duretz et al., 2014; Magni et al., 2014; Menant et al., 2016; Pusok et al., 2018; Hertgen et al., 2020), often ignoring the role of lateral heterogeneities in the overriding plate. Yet, in natural settings, subducting plates frequently subduct beneath strongly variable overriding lithosphere. For example, the nature (oceanic or continental) and thickness of the overriding lithosphere can vary across different segments of the same subduction zone. Notable examples of along-strike variations of continental OP mechanical thickness (as defined by the depth of the 1100 C isotherm) include the change from 66.4 km in Alaska to 51.6 in the Aleutians; from 98.1 km in Mexico to 58.5 km in North Costa Rica, and from 96.3 km in Nankai to 83.2 km in Kyushu (Kerswell et al., 2021).

The numerical modelling by Gea et al. (2024) showed that the influence of variable OP thickness is modulated by the imposed lateral boundary conditions. Fixed trailing edges maximize this influence and are representative of narrow slabs attached to larger slowly moving plates, as Mediterranean (Calabrian, Hellenic), South Scotia and Makran subduction zones. Free trailing edges better represent wide subduction zones bordering fast moving plates. Prominent examples are South America, Northwest Pacific and Sunda subduction zones.

Recognising and studying lateral heterogeneities in the overriding plate is essential for advancing our comprehension of subduction systems. The effects of overriding plate thickness on subduction dynamics and mantle flow have been investigated in detail in previous studies (e.g., Meyer and Schellart, 2013; Sharples et al., 2014; Hertgen et al., 2020; Gea et al., 2024). However, just a few studies have attempted to study the impact of lateral variabilities in the structure of the overriding plate on trench-parallel flow (Rodríguez-González et al., 2014) and subducting plate velocity (Rodríguez-González et al., 2016), but those models used a stationary trench that prevented them from studying features such as the trench shape evolution and trench kinematics or the evolution of the induced mantle flow. Rodríguez-González et al. (2014) also showed that OP heterogeneities may result in highly contorted slab geometries, which may ultimately develop slab tears in the presence of intraslab weak zones. The subduction of fracture zones and transform faults represented as weak zones in the SP has been investigated in previous 3D geodynamic modelling (e.g., Burkett and Billen, 2010; Manea et al., 2014; Liu and Pysklywec, 2023), but they did not take into consideration OP heterogeneities.

In this study, we have performed self-consistent 3D numerical models with a mobile trench to investigate how along-strike variations in the OP thickness can induce along-strike variations in subduction dynamics and mantle flow. Additionally, in one set of models we have included a narrow weak zone in the middle of the subducting plate simulating a fracture zone, in order to investigate the combined effects of OP heterogeneities and the subduction of fracture zones on the occurrence of slab segmentation. Finally, we compare our model results with observations from natural subduction zones.

Basic equations and approximations

The simulations have been performed with version 2.4.0 of the finite-element code ASPECT (Advanced Solver for Problems in Earth’s ConvecTion) Kronbichler et al. (2012); Heister et al. (2017); Gassmöller et al. (2020); Bangerth et al. (2022). We have used the Boussinesq approximation to solve the coupled conservation equations for mass, momentum and energy (Equations 1–3), which assumes that the density is constant in all equations except in the buoyancy term of the momentum equation (right hand term in equation (Equation 2)): where ε ̇ = 1 2 ∇ u + ∇ u T is the strain rate tensor, ρ is the density, μ is the viscosity, P is the pressure, u is the velocity, g is the gravitational acceleration, c p is the specific heat, T is the temperature, k is the thermal conductivity, ρ o is the density at a reference temperature T o, and H is the radiogenic heating, which is neglected in our models.

Besides these equations, ASPECT solves the advection of compositional fields c i, which are Eulerian finite element fields used to track materials and their properties throughout the simulations. The compositional fields are advected on the mesh with the velocity field (Equation 4):

Model setup

We have adopted an initial setup similar to Gea et al. (2024), where self-driven subduction is modelled in a three-dimensional space by including a 220 km-long proto-slab with a dip angle of 40 ° (Figure 1). We do not include external forces or velocities, and progressive subduction and mantle flow are entirely driven by the negative buoyancy of the slab. The subducting plate (SP) has a width of 1,400 km and is laterally homogeneous, representing a relatively mature oceanic lithosphere with a thickness of 95 km. It consists of a 10 km thick weak oceanic crust and an 85 km thick lithospheric mantle. The subducting oceanic crust has a low viscosity of 3 ⋅ 1 0 20 Pa s and acts as a decoupling layer between the subducting and the overriding plates. In the reference models (RM), we use an uniform thickness of 70 km for the overriding plate (OP), corresponding to 30 km thick continental crust and a 40 km thick lithospheric mantle. The OP thickness of 70 km is consistent with the average thickness of 60 km (depth to the 1,200 °C isotherm) for backarcs obtained by Currie and Hyndman (2006), on the basis of surface heat flow, seismic velocity, and xenolith thermobarometry indicators. Following Kerswell et al. (2021), we assume that backarcs are representative of initial steady-state continental lithosphere (unless affected by strong crustal thinning or magmatic activity). In the asymmetric models (AM), we include two portions of OP (700 km wide each) along the trench-parallel direction with thicknesses of 70 and 95 km (Figure 1a). This pairing is representative of the variation from an average thickness of a backarc lithosphere to a unperturbed continental lithosphere thickness. We have performed additional tests with different pairings of OP thickness. Given that the lithosphere is thermally (not compositionally) defined, the abrupt initial change in lithospheric thickness partially relaxes thermally and evolves to a more gradual thickness change. The model also includes two lateral plates, one on each side of the SP, separated by a 20 km wide zone of 1 0 21 Pa s. We use a moderate viscosity of 1 0 21 Pa s to prevent large amounts of trench shortening along the trench-parallel direction caused by lateral subduction along the sides of the SP (e.g., Yamato et al., 2009). For simplicity and computational convenience, the lateral plates are made of 95 km thick lithosphere, without any compositional stratification. In a suite of models, AM-WZ, we have introduced in the SP a weak zone 20 km wide perpendicular to the trench (red line in Figure 1). Table 1 summarizes the main characteristics of the models shown in this study.

We use a Cartesian geometry in which the box size is 3,000 km long, 3,000 km wide and 660 km deep. With a box width of 3,000 km and a SP width of 1,400 km, we achieve a lateral distance between the slab edges and the lateral sidewalls of 800 km, which meets the criteria proposed by Schellart et al. (2007) that the lateral distance should be at least half of the slab width (700 km in this case) and larger than the upper mantle depth (660 km). This distance prevents the box edges from affecting the subduction dynamics and mantle flow. We use an adaptive mesh refinement that allows for an increased resolution in areas of strong gradients of temperature, viscosity, and/or composition, resulting in a high resolution of 3.9 km (length) by 3.9 km (width) by 2.6 km (depth) within the subducting and overriding plates, while maintaining a coarse mesh of 125 × 125 × 82.5 km far from the subduction zone. A high resolution is required to resolve the subduction zone interface and to ensure realistic subduction evolution.

We present models in which the far-field boundary conditions (free or fixed) of the subducting and overriding plates represent two end-members. In models RM-fixed and AM-fixed, the subducting and overriding plates are fixed at their trailing edges (Figure 1b), representing a scenario where the subduction zone is part of much larger and relatively stationary plates. In this case, both the subducting and the overriding plates are 1,500 km long. In models RM-free and AM-free, the subducting and overriding plates are free to move at their trailing edges. To achieve that, we have implemented an oceanic ridge on the edge of the SP, while for the OP we have reduced its length to 1,000 km and have added a weak corner with a constant viscosity of 1 0 20 Pa s at its trailing edge (Figure 1c). The weak corner is 95 km thick and extends to the lateral side of the model box, and it helps to ensure the free motion of the OP, mimicking ridge-like thermal structures. Previous subduction modelling (Billen and Arredondo, 2018) has shown that adding a ridge instead of a weak zone on the overriding has only a small effect on slab dynamics.

The temperature is set to 293 K at the surface, and it increases linearly to 1,643 K at the lithosphere-asthenosphere boundary. The bottom boundary is fixed at the mantle temperature (1643 K). For the free models, the thermal structure in the SP is described by the plate model (Fowler, 2004) to assign an age-dependent temperature, with a bottom temperature equal to that of the mantle. Mechanical boundary conditions are free-slip for all boundaries, except for the bottom of the box, which is modelled as a rigid no-slip boundary.

Rheology

We use a visco-plastic rheology in which the viscous deformation is given by a composite viscosity that takes into account both diffusion creep and dislocation creep (Equations 5, 6). In this way, we make a huge computational effort to improve the rheology of the 3D models of Gea et al. (2024) by including dislocation creep. The viscosities associated with both creep mechanisms can be expressed in a general way using the same relationship (Karato and Wu, 1993; Hirth and Kohlstedt, 1996): where the subindex i denotes the deformation mechanism (diffusion or dislocation), A is a prefactor of the equation, d is the grain size, m is the grain size exponent, n is the stress exponent, ε ̇ i i is the square root of the deviatoric strain tensor second invariant, E and V are an activation energy and volume respectively and R is the gas constant. Each of these parameters varies depending on the mineral (except for the grain size) and the deformation mechanism. For the diffusion case (n = 1 and m ≠ 0) viscosity depends on the grain size but there is no dependence on the strain rate tensor. We use a constant over time grain size of 1 0 − 2 m (Billen and Hirth, 2007). In the dislocation case (n > 1 and m = 0), there is a power-law dependence between viscosity and strain rate and there is no dependency on grain size.

Both types of viscous deformations can act simultaneously and the contributions of diffusion μ diff and dislocation μ disl to the effective viscosity are harmonically averaged into a composite viscosity μ comp:

Plastic deformation will occur if the applied stress exceeds the yield stress σ y. We use the expression for σ y given by the 3D formulation of the Drucker-Prager yield criterion (Equation 7): where C is the cohesion and ϕ is the internal friction angle. In this visco-plastic approach, the final viscosity (Equation 8) is dominated by the minimum between the plastic and the creep mechanisms:

We use rheological parameters from wet olivine (Hirth and Kohlstedt, 2003) for the upper mantle (asthenosphere and lithosphere). For the crust of the OP, we use a diffusion rheology from wet anorthite feldspar (Bürgmann and Dresen, 2008). All model parameters are listed in Table 2. Finally, the viscosity is limited between minimum and maximum values of 1 0 19 and 1.57 ⋅ 1 0 23 Pa s to preserve numerical stability. Additional tests indicated that using a lower minimum viscosity cutoff of 1 0 18 Pa s makes the simulations much more time-consuming while leading to the same qualitative results.

Subduction evolution and trench kinematics

As previously mentioned, the far-field boundary conditions of the models represent two end-members. The fixed case (Figure 1b), in which the trailing edges of both the SP and OP are fixed, maximises slab rollback, which is entirely accommodated by OP extension. In the free case (Figure 1c), where the SP and OP are free at their trailing edges, slab rollback is reduced and the subduction process is accommodated by trenchward motion of both the subducting and overriding plates.

Figures 2, 3 show the subduction evolution of reference models with a homogeneous OP (RM-fixed and RM-free) and models with an asymmetric OP (AM-fixed and AM-free), respectively. Supplementary Figures S1-S3 show trench-perpendicular cross sections at the same timesteps. The subduction dynamics of RM models show the typical evolution of a subduction zone, as described in previous works (e.g., Funiciello et al., 2003; Schellart, 2004; Magni et al., 2012; Guillaume et al., 2021; Gea et al., 2024). First, the slab freely descends into the upper mantle due to the density contrast between the SP and the sub-lithospheric mantle (Figures 2a,b). The free sinking phase is defined here as the period before the tip of the slab enters the viscous deeper part of the model, which is assigned a depth of 400 km for comparison purposes. During this phase, the trench retreats towards the SP and its velocity increases as the amount of slab entering the upper mantle increases (Figure 4a). When the slab tip approaches the bottom of the box (z > 400 k m), there is a short period of interaction in which the dip angle of the slab decreases and the slab lies flat over the 660-km discontinuity. After the interaction between the slab and the bottom of the box, the slab continues to flatten on the 660-km discontinuity, only driven by its negative buoyancy (Figures 2e,f and Supplementary Figure S1).

In the free case, RM-free, the SP is allowed to move towards the trench. Close to the left boundary, new thin thermal lithosphere is created over time due to the fixed cold temperature top-boundary condition. The trenchward motion of the SP causes the slab to sink faster into the upper mantle and to enter earlier the viscous deeper part of the box (compare Figures 2c,d). The tip of the slab reaches a 400 km depth after 5.5 Myr in RM-free model and 7.5 Myr in the RM-fixed model. Initially, both models exhibit similar trench retreat velocities, with model RM-free reaching a maximum V T of ∼ 2 cm/yr after 6–7 Myr of evolution (Figure 4a), just after the slab reaches the deep viscous zone (z > 400 k m). Afterwards, trench retreat velocities are much reduced compared to the model with fixed trailing edges. This reduction of trench retreat velocity occurs because the slab continues to sink vertically, developing a folded geometry (Figures 2d,f and Supplementary Figure S1). This result agrees with previous studies showing that free-trailing edges tend to develop quasi-stationary trenches with slabs forming folded piles atop the lower mantle (Stegman et al., 2010). Model RM-fixed does not exhibit the decrease or stabilization of trench retreat velocity during the phase of sinking into the viscous deeper portion of the box (t > 7,5 Myr) obtained in previous studies (e.g., Faccenna et al., 2001; Schellart, 2004; Gea et al., 2024; Sharma et al., 2026). Instead, trench retreat velocity continues to increase, albeit at a reduced rate. The reason for this discrepancy is that the decrease in sinking rate is counteracted by localised weakening and extension of the OP (Supplementary Figure S1), which ultimately leads to OP opening at 11 Myr. This opening is accompanied by a sharp increase in retreat velocity (Figure 4a). The trench geometry remains practically linear at its centre during the free sinking phase in both models (Figures 4c,d). After the OP opening, the trench in model RM-fixed continues to retreat and it seems to adopt a slight “w”-shaped geometry after 14 Myr (Figure 4c). In contrast, model RM-free shows an essentially linear trench that evolves towards a stationary location (Figure 4d), with the SP moving trenchwards and resulting in a folding geometry (Figure 2f and Supplementary Figure S1). This model does not exhibit OP breakup because the OP moves laterally trenchwards with little internal deformation (Supplementary Figure S1).

In the AM models, the OP presents two portions of varying thickness, 70 km and 95 km respectively, along the trench-parallel direction (Figure 1a). This asymmetrical OP structure results in more complex subduction dynamics (Figure 3). Below the thicker part of the OP, the mantle wedge is more viscous because there is less space available for slab-induced high strain-rate mantle flow (which would reduce the viscosity). In turn, this reduced corner flow reduces thermal erosion of the OP lithospheric mantle above the mantle wedge and, as a consequence, the amount of effective coupling between the SP and the OP is higher (e.g., Rodríguez-González et al., 2012; Hertgen et al., 2020). This stronger coupling reduces the slab sinking velocity. During the free sinking phase, the difference in sinking rate along-strike causes the slab to sink deeper into the upper mantle below the thin part of the OP (Figure 3), which in turn further increases the difference in slab pull. The difference of trench retreat velocity along the thick and thin parts of the OP gradually increases during the free sinking of the slab (Figure 4b). This progressive differentiation of retreat velocities along the trench-parallel direction results in a clockwise trench rotation (Figures Figures 4e,f), producing an arcuated trench, which is much more pronounced when the plates are fixed at their trailing edges (Figure 4e). For this AM-fixed model, the trench retreat velocity above the thin part of the OP quickly increases after 10 Myr of model evolution due to the opening of the upper plate, which is not observed above the thicker part (Figure 4b). The reason for this difference is that a thicker (colder and stronger) OP opposes more resistance to extension and therefore trench retreat is reduced compared to a thinner and weaker OP (e.g., Hertgen et al., 2020). The tip of the slab in the AM-fixed model reaches a depth of 400 km after 8 Myr and 10 My of evolution below the thin and thick portions of the OP, respectively. The AM-free model shows a faster slab sinking rate, reaching a depth of 400 km after 5.5 Myr and 6.5 Myr beneath thin and thick OP, respectively. Both AM models show that beneath thick OP, the deep flat segment of the slab is shorter and that the slab subducts more steeply due to the reduced trench retreat velocity (Figures 3e,f and Supplementary Figures S2,S3).

Slab-induced mantle flow

In order to visualize the slab-induced mantle flow, we show in Figure 5 the horizontal (arrows) and vertical (color palette, with positive values for upward direction) components of the velocity field at a depth of 250 km. The poloidal component of the mantle flow shows the classic pattern of a retreating slab, with two flow cells developing through the subduction zone (e.g., Funiciello et al., 2006; Piromallo et al., 2006): one clock-wise in the sub-slab region and one counter-clockwise in the mantle wedge. We analyse toroidal flow using horizontal vorticity derived from the horizontal components V X and V Y.

The RM-fixed model shows the typical toroidal flow pattern around a subduction zone, with a returning flow from the sub-slab region to the mantle wedge (Figure 5a) (e.g., Funiciello et al., 2006; Strak and Schellart, 2016; Goldberg and Holt, 2024). Two symmetrical toroidal flow cells develop around the edges of the slab. Both in the sub-slab region and in the mantle wedge, the mantle flow transitions from slab-perpendicular in the centre to slab-parallel at the edges of the slab (Figure 5a), resulting from along-strike pressure gradients (Goldberg and Holt, 2024). The toroidal flow starts to develop during the free sinking of the slab into the upper mantle, but it becomes more pronounced in advanced stages of subduction. Model RM-free also shows a similar flow pattern with two symmetrical flow cells surrounding the slab (Figure 5b). However, both the maximum horizontal velocities and the extent of the mantle flow cells are reduced due to the reduced slab rollback.

In contrast, the AM-fixed model exhibits a much more complex toroidal flow pattern. The thinner OP allows for efficient toroidal mantle flow, which causes the returning mantle flow surrounding both edges of the slab to converge within the mantle wedge domain beneath the thin portion of the OP (Figure 5c). This pattern results in a dominantly trench-parallel flow below the thick portion of the OP and a mainly trench-perpendicular flow beneath the thinner upper plate, exhibiting V X up to ∼ 6 cm/yr during the free sinking phase. The extent of the toroidal flow cells is greater at the edge of the subduction zone corresponding to the thinner OP. During the last stages of slab sinking, the trench-parallel extent (defined here by horizontal velocities greater than 1.5 cm/yr) of the horizontal mantle flow is ∼ 450 km at the slab edge located beneath thin OP, while the toroidal flow just extends ∼ 200 km from the slab at its edge beneath thick OP. The trench-perpendicular extent of the toroidal flow also varies, changing from ∼ 700 km below the thin OP to ∼ 500 km below the thicker part of the upper plate. Allowing the plates to freely move at their trailing edges (AM-free model) results in a similar flow pattern, but much less pronounced since the difference in trench retreat velocities along-strike is smaller compared to model AM-fixed (Figure 5d).

The effect of localising deformation in the subducting plate

To investigate the possibility of the occurrence of slab tearing related to OP heterogeneities, we have included in the asymmetric model AM-fixed a weak zone 20 km wide in the middle of the SP (red line in Figure 1a) model AM-WZ). This weak zone has a viscosity of 1 0 20 Pa s and acts as a fracture zone that localises deformation and facilitates the process of tearing. We tested different values for the thickness of the weak zone and obtained that 20 km gave a good compromise between properly reproducing the viscosity decrease and, at the same time, allowing for a computationally convenient grid resolution (Gea et al., 2024). We have kept the subducting and overriding plates fixed at their trailing edges to maximise subduction rollback and trench retreat. For simplicity, we do not introduce any variation of the subducting lithosphere across the modelled fracture zone.

The evolution of this model AM-WZ is illustrated in Figure 6 and in vertical cross sections in Supplementary Figure S4. The dynamics is initially similar to that of model AM-fixed (Figures 3a,c,e), with the slab sinking deeper and retreating faster at its edge below the thin part of the OP (Figure 6a). The increase in the difference in trench retreat velocities along-strike coupled with the presence of a weak zone causes a slab tear to form after 10 Myr of model evolution (Figure 6a), which continues to grow as subduction evolves (Figures 6b-d).

Figure 7 shows the horizontal and vertical components of the mantle flow during the process of slab tearing at a depth of 250 km. At 10 Myr, when the slab tear is not developed, the mantle flow follows the same pattern as that described for model AM-fixed (Figure 5c). Once the slab tearing develops, mantle flow coming from the sub-slab region rapidly infiltrates through the tear (Figure 7b). This results in a predominant clockwise toroidal flow through the tear that converges in the mantle wedge beneath the OP thin portion (Figure 7c). A minor part of mantle material coming through the tear attempts to form a counterclockwise toroidal flow cell towards the wedge beneath thick OP, but this contribution is counterbalanced by the toroidal flow around the slab edge beneath thick OP, converging at the centre of the subduction zone. The process of slab tearing at the centre of the subducting slab also leads to strong upwellings of mantle material through the tear (Figures 7c,d). In contrast, a control test introducing a weak zone in the reference model with a homogeneous OP exhibits the same evolution as the RM-fixed model, the slab sinks with an along-strike constant dip and there is no significant intraslab stress to produce the tearing. Similarly, introducing the weak zone in the asymmetric OP model with free trailing edges does not lead to slab tearing because the contorted geometry of the slab is much less pronounced than in the fixed edges case, and therefore not enough intraslab stress arises to develop slab tearing. We performed additional simulations by placing the weak zone in the prolongation of thin and thick OP (models AM-WZ-Thin and AM-WZ-Thick in Table 1). The results are illustrated in the Supplementary Figure S5. When the weak zone is subducting either beneath thin or thick OP it starts developing tearing about 3–4 Myr later than in the model with the weak zone placed right at the location of OP thickness contrast. In this case the weak zone coincides with the area of maximum slab rollback velocity gradient (Figure 4f); this produces a markedly larger aperture in a much shorter interval (t = 10–12 Myr) than in models with displaced weak zones (Supplementary Figure S5).

Effects of along-strike variations in overriding plate thickness

Along-strike heterogeneities in slab dynamics and trench kinematics are commonly observed in subduction zones. Previous numerical and analogue models have extensively studied the role of laterally varying properties in the SP, such as slab buoyancy due to subduction of aseismic ridges (e.g., Martinod et al., 2005; Martinod et al., 2013; Flórez-Rodríguez et al., 2019; Suchoy et al., 2022; Zhao and Leng, 2023) or lateral changes in slab age (Morra et al., 2006). Simultaneous subduction of adjacent continental and oceanic lithospheres has also been shown to result in along-trench variations between the oceanic subduction side and the continental collision side (van Hunen and Allen, 2011; Duretz et al., 2014; Magni et al., 2014; Menant et al., 2016; Pusok et al., 2018). Slab tearing propagating from the collision side leads to differential trench retreat velocities along-strike and associated trench deformation (e.g., Duretz et al., 2014; Sternai et al., 2014).

The numerical models of this work focus on the lateral changes of the upper plate structure, rather than on the SP, thus complementing the above-mentioned studies by demonstrating that along-strike variations in OP thickness lead to strong along-strike changes in slab morphology and trench retreat velocities. The length of the subduction interface and thus the coupling between plates increases below the thicker OP, significantly reducing the slab-sinking rate and resulting in shallower slab depths below the thick part of the OP. In fact, the progressive larger slab volume submerged in the upper mantle below the thin part of the OP enhances the slab pull difference along-strike during the free sinking of the slab. In addition, the thin/weak OP portion provides a lower resistance to extension and therefore promotes trench retreat (e.g., Meyer and Schellart, 2013; Hertgen et al., 2020; Gea et al., 2024). The combination of both effects, the enhanced slab pull and thin OP extension, generates an along-strike gradient of trench retreat velocity which leads to trench rotation. These results are in agreement with the 3D subduction models of Zuhair et al. (2022). Using numerical modelling, they studied the role of along-strike slab pull gradients by initially varying the slab depth along the trench, while lateral variations of the slab depth are naturally achieved in our models due to the variable thickness of the OP. Their results showed that the difference in slab depth along the trench-parallel direction drives a rotation in the trench, with the retreat of the long slab and the advance of the trench above the short slab, comparable to the asymmetric model presented in this work. The main difference is that we do not observe significant trench advance above the short slab segment of the models (below the thick OP), but rather a nearly stationary trench. This can be explained by the fact that in our models the slab depth slowly increases below the thick OP, providing sufficient slab pull to facilitate subduction and prevent trench advance. Previous models from Rodríguez-González et al. (2016) have also shown that along-strike variations of the OP thermal structure can lead to along-strike variations in SP velocity, with the SP being faster beneath the warmer (thinner) portion of the OP. Both models with variable OP thickness (AM-fixed and AM-free) exhibit the same qualitative results. Still, the arcuate trench geometry and the along-strike heterogeneous mantle flow are much more pronounced for the fixed trailing edges case. Another significant difference is that models with free trailing edges never develop slab tearing. Supplementary Figures S6,S7 show the results of additional tests (with fixed trailing edges) considering different pairings of OP thickness. As expected, stronger OP thickness contrasts result in more contorted slab geometries (Supplementary Figure S6) and more arcuate trench shapes (Supplementary Figure S7). These tests reveal that differential trench retreat does not induce slab tearing in the absence of a pre-defined oceanic weak zone. The trench geometry is modulated by the OP thickness absolute values (compare 40/70 km vs. 80/110 km OP thickness pairings in Supplementary Figure S7). In contrast, additional tests using mixed boundary conditions (SP-fixed, OP-free and SP-free, OP-fixed) reveal a negligible effect of OP thickness variations on subduction dynamics, due to the small OP internal deformation. For example, with a fixed OP but a free SP, the subduction process is accommodated by a trenchward motion of the SP with almost null trench retreat, as observed in the Cocos subduction zone.

The asymmetry of mantle flow around the slab edges is a feature that can be compared with seismic anisotropy data. Assuming that fast SKS-wave directions primarily map subduction-induced mantle flow, a possible example of this asymmetric toroidal flow is the Calabrian subduction zone, where fast directions only turn around the SW edge of the Ionian slab in the western Sicily side (e.g., Civello and Margheriti, 2004; Baccheschi et al., 2024). Similarly, for the Gibraltar Arc subduction system fast directions only turn around the SE end of the Betic-Rif chain (Diaz et al., 2010). However, caution must be taken when attributing asymmetric toroidal flow to OP variable plate thickness, as other reasons (e.g., plate rotations, asynchronous continental collision) cannot be ruled out.

Present numerical models have also shown that in the presence of an asymmetric OP, the slab-induced mantle flow deviates from the idealised flow pattern with symmetrical toroidal flow cells developing around each slab edge. Instead, our results show that differential trench retreat velocities along-strike induce a much more complex mantle flow. In the sub-slab region, the mantle flow transitions from slab-perpendicular beneath the slab with highest rollback (in the thinner OP side) to slab-parallel in the edges of the slab, resulting in a large region of trench-parallel flow below the slab. Similarly, the returning mantle flow coming from the sub-slab region tends to converge beneath the mantle wedge domain beneath the thin OP (Figure 5c), resulting in dominantly trench-parallel flow beneath the thick OP and trench-perpendicular beneath the thin OP. Rodríguez-González et al. (2014) previously investigated the mantle flow in models with along-strike variations in the OP thermal structure. In their models, large regions of trench-parallel flow were observed both in the sub-slab domain and in the mantle wedge. However, they had a fixed trench that prevented slab rollback and therefore underestimated toroidal flow (e.g., Faccenda and Capitanio, 2012; Meyer and Schellart, 2013; Schellart and Moresi, 2013; Chen et al., 2016).

Comparison with natural subduction zones

General comparison

The numerical results presented in this study contribute to a better understanding on the origins of along-strike heterogeneities observed in natural subduction zones, such as the slab morphology, trench retreat velocities, slab-induced mantle flow or slab segmentation. We emphasise that these models include several simplifications, corresponding to an idealised subduction process, and do not include some factors that influence subduction zones on Earth such as interactions with external slabs (Király et al., 2016), external mantle flow coming from other processes (Chertova et al., 2018) or phase transitions (Li et al., 2025). More specifically, we acknowledge that the slab width of our numerical models (1,400 km) represents relatively narrow subduction zones. This choice has been motivated by the fact that most of the Earth’s subduction zones are relatively narrow, with a median width of 1,050 km (Schellart, 2023). However, some subducting plates might extend for several thousands kilometres along-strike, including some of the examples that we discuss in this section. The precise comparison with natural examples is challenging because of different reasons. First, the present-day OP thickness in the selected location may differ significantly from the initial value. Second, lithosphere thicknesses obtained from different methods (flexure, seismic data, heat flow, mantle xenoliths) can provide significantly different estimates. Third, external controlling process, for instance, along-strike variation of the convergence velocity or heterogeneities in the SP, can play a major role on the observed along-strike variability. Therefore, natural comparisons are restricted to qualitative trends and must be approached with caution.

The Tonga-Kermadec-Hikurangi subduction zone is a very complex system where the Pacific plate subducts beneath the Indo-Australian plate. Tectonic reconstructions and tomographic studies have shown a considerable variability in the style of subduction along the trench (Schellart et al., 2006). The Tonga trench is characterised by high trench retreat velocities and rapid trenchward SP motion, with high sinking velocities and the slab both flattening at the mantle transition zone and penetrating into the lower mantle. The Kermadec-Hikurangi segment, however, presents slower trench retreat velocities and trenchward plate velocities, with slow sinking velocities and a steep slab almost reaching the mantle transition zone (Schellart and Spakman, 2012). The fast trench retreat of the northern edge of the Tonga-Kermadec-Hikurangi subduction zone is thought to be caused by two main factors. First, the Samoa plume-slab interaction has resulted in the migration of plume materials around the northern edge of the subduction zone (Chang et al., 2016; Yu et al., 2022). Secondly, the presence of the Hikurangi plateau at the southern segment of the subduction zone resists slab sinking and rapid trench retreat, while further to the north the Tonga trench continues to retreat eastwards (Strak and Schellart, 2018). Overriding-plate lithosphere thickens southward, from roughly 70–90 km in the Tonga-Kermadec oceanic arc to roughly 100–130 km beneath the Hikurangi margin (Fullea et al., 2021). This scenario resembles our results for a model with fixed trailing edges and a heterogeneous OP (AM-fixed), where fast trench retreat is accommodated through OP extension mainly in the thinner (i.e., weaker) portion of the OP lithosphere (Figures 4b,e). We therefore suggest that the thinned Indo-Australian OP at the Tonga-Kermadec segment during the opening of the South Fiji basin in the Oligocene (Schellart et al., 2006) may have also led to fast slab rollback in the northern edge of the subduction zone, combined with the other two mentioned mechanisms to explain the difference in trench kinematics along-strike.

The Ryukyu-Nankai trench marks the subduction of the Philippine Sea Plate beneath the Eurasian Plate. Trench retreat velocities present significant variations along-strike, with velocities ranging from ∼ 5 cm/yr in the Ryukyu segment in the southwest to just ∼ 1 cm/yr in the Nankai trench (Schellart et al., 2007). Similarly, seismic tomography studies in the region indicate that the Ryukyu slab extends deeper into the upper mantle with a flattened shape at ∼ 500 km depth (Hayes et al., 2018), while the slab at the northeast beneath the Nankai trench is steeper and only reaches a depth of 200 km (Lallemand et al., 2001; Hayes et al., 2018). Although there is a difference in slab age along-strike, with the Ryukyu segment characterised by an older (40–50 Ma) and colder slab than in the Nankai trench (≤ 20 Ma) (e.g., Ishizuka et al., 2011), global observational studies on active subduction zones have shown that the correlation between the SP age and the trench velocity is negligible (Heuret and Lallemand, 2005; Schellart, 2010) and the correlation between the age of the SP and its velocity is low (Schellart, 2010). On the basis of results for the AM-fixed model, we suggest that the asymmetric structure of the OP, with the Ryukyu slab subducting beneath the Okinawa backarc (50–70 km thick OP; Fullea et al. (2021) and the Nankai slab sinking beneath the thicker continental Japan lithosphere (90–120 km thick OP; Fullea et al. (2021), contributes to the observed pattern of retreat velocities along the trench.

The model including a weak zone in the SP has also shown that differential trench retreat along-strike, irrespective of its cause, can lead to vertical slab tears across pre-existing weak zones (Figure 6). For example, tomographic studies have imaged two vertical slab tears in the Lesser Antilles (Harris et al., 2018), which the authors proposed to have been caused by differential rollback facilitated by the presence of microplates at the northern boundary of the Caribbean plate. Additionally, Hall and Spakman (2015) propose that a vertical tear in the slab extends through the entire upper mantle to the north of west Sumatra, where a significant bend in the slab’s morphology is observed. We speculate that the OP thickness variations from 65 km below the Andaman Sea to 90–100 km beneath Sumatra to the SE (Yu et al., 2017) can contribute to this folded slab, and that slab tearing was likely facilitated by the subduction of different fracture zones in the region. Similarly, Rodríguez-González et al. (2012) proposed that the large dip variability of the Cocos slab is likely related to the change of the thermal state of the overriding plates, with flat subduction occurring under cold lithosphere in southwestern Mexico and steep subduction under the warmer lithosphere of the Chortis block in the northwestern Caribbean plate. The subduction of the Tehuantepec fracture zone roughly corresponds to the change in Cocos slab dip, and likely enabled the occurrence of a vertical slab tear imaged by Calò (2021). Finally, subduction zone segmentation by means of vertical tears and different trench retreat rates likely occurred in the central and western Mediterranean. The subcrustal structure beneath the Tyrrhenian-Apennines system is characterised by an intricate array of segmented slabs separated by vertical tear faults (e.g., 3-D mapping by Rappisi et al., 2022). Tear faults in this region follow the trend of fracture zones and transform faults inherited from the break-up of Pangea (Pierantoni et al., 2020). Based on our modelling results, we hypothesise that significant heterogeneities in the OP lithosphere (back-arc basins, continental blocks) likely produced significant variations in trench retreat rates. These OP heterogeneities combined with the subduction of weak fracture zones could result in oroclinal bending and in the segmentation of the Apenninic slab along vertical tear faults.

The model with vertical tearing is equally useful to illustrate how the development of a tear can affect the subduction-induced mantle flow. Similar to the flat subduction models from Liu and Pysklywec (2023), we find that the opening of the tear leads to a rapid flow of mantle material infiltrating from the sub-slab region through the tear (Figure 7). The mantle flow through the vertical tear spreads to all directions, but the toroidal flow cell around the slab with faster retreat velocity is more prominent (Figure 7). The upwelling component of the mantle material is also amplified through the slab tear, which could lead to tear-related magmatism during slab rollback and other associated processes. Natural subduction zones that include vertical tears where the above-mentioned mantle flow pattern could be observed include the South Shetland (e.g., Parera-Portell et al., 2023) and the Tyrrhenian-Apennines collision system (e.g., Rosenbaum et al., 2008).

We will now provide a more detailed discussion of the Melanesia subduction zone, a complex subduction zone in the Southwest Pacific that presents a setting comparable to our asymmetric model. We will discuss the recent geodynamic evolution of the subduction zone in light of the numerical results of this study.

Melanesia subduction zone

The Southwest Pacific region is a complex area characterised by subduction and backarc extension processes resulting from the interaction of the Australian and Pacific plates since the Late Cretaceous (e.g., Petterson et al., 1999; Hall, 2002; Schellart et al., 2006; Seton et al., 2016). The present-day Melanesia subduction zone (New Britain-San Cristobal-New Hebrides) (Figure 8), where the Australian Plate subducts northward underneath the Pacific plate, formed during the Early to Middle Miocene (e.g., Hall, 2002; Hall and Spakman, 2002; Schellart et al., 2006). This subduction zone presents a great variability of the slab depth along the trench, with the slab almost reaching the mantle transition zone in the New Britain Trench, but a depth of ∼ 300 km in the New Hebrides Trench and ∼ 100 km in its central part, where the slab subducts beneath the Ontong Java Plateau (Figures 8a,c). The trench retreat velocity also varies strongly along the trench, with the edge segments of the subduction zone retreating much faster than its centre, showing retreat velocities of ∼ 4.5 cm/yr in the northwestern segment of the trench and ∼ 9 cm/yr in the southeastern segment (Figure 8a) (Schellart et al., 2007). Indeed, the rapid rollback along the New Hebrides trench caused the clockwise rotation of the trench and the opening of the North Fiji Basin from ∼ 10 Myr to present, while the southeastward rollback in the New Britain trench triggered the counter-clockwise rotation of the trench together with the opening of the Manus backarc Basin (e.g., Schellart et al., 2006; Wallace et al., 2009). Predictions from AM-fixed model show that the occurrence of back-arc extension in the thinner OP portion allows for rapid slab roll-back and high trench retreat rates (Figures 4b,e and Supplementary Figure S2), leading to trench rotations qualitatively comparable to that of New Hebrides and New Britain trenches.

These observed differences in slab length and trench motion can be explained by our model results showing that the asymmetric thickness of the OP leads to variations in slab depth along-strike, as well as trench kinematic variations and trench rotations. Both the edge segments of the subduction zone (New Britain and New Hebrides) started to subduct below very young oceanic crust (Figure 8b), which triggered fast slab rollback, large slab depths and trench rotation (Figures 8a,c), while the slab at the San Cristobal trench subducted beneath the Ontong Java Plateau, which caused the trench to be near stationary and the slab to reach shallow depths (Figures 8b,c). Thus, we propose that much of the trench kinematics observed in the Melanesia subduction zone is driven by the heterogeneous structure of the overriding Pacific plate, with a seismic LAB depth of 130 km beneath the Ontong Java Plateau, which is about 40–50 km deeper than beneath the surrounding Pacific Ocean seafloor (Isse et al., 2021). It is important to note that we do not discard the effect that other factors might have had on the evolution of this subduction zone. For example, the collision of the West Torres Plateau with the New Hebrides Trench that occurred at ∼ 0.7 Ma (Greene and Collot, 1994) and the arrival of the D’Entrecasteaux Ridge at the trench at ∼ 1.5–3 Ma (Greene and Collot, 1994) may have contributed to the rapid clockwise rotation of the New Hebrides Trench. Slab-slab interactions with the Trobriand subduction zone in the northwest and with the Tonga subduction zone in the southeast may also have influenced both the trench kinematics and the induced mantle flow, since the offset distance between these subduction zones is small (Király et al., 2016; 2021; Zhang et al., 2024). Finally, we acknowledge that the width of the Melanesia subduction zone (∼ 4,400 km) is another important factor affecting the evolution of the subduction zone, since it has been shown that very wide slabs tend to develop a convex geometry towards the OP (e.g., South America), with trench retreat velocities increasing towards the slab edges (Schellart et al., 2007; Stegman et al., 2010; Strak and Schellart, 2016).

Model limitations

The model box has a depth of 660 km and the bottom of the box is modelled as a rigid no-slip boundary simulating an impenetrable upper-lower mantle discontinuity. In nature, some subducting slabs are able to penetrate into the lower mantle, such as in the subduction zones of Hellenic (van der Hilst et al., 1991), Peru (Fukao et al., 2001) or Mariana (van der Hilst et al., 1991). Therefore, the use of an impenetrable barrier at the bottom of the box amplifies the effects of the slab-660 km discontinuity interaction. Future studies including a viscous lower mantle are needed to quantify this effect. Similarly, modelling wide subduction zones requires much larger model domains. Additional limitations stem from simplifications such as not considering phase transitions or grain size evolution. Future work that included a systematic sensitivity analysis of the OP thickness variations combined with different parameters (e.g., lateral OP thickness gradient, SP thickness, subduction zone width, trailing edge conditions) would be helpful to establish a more direct link between the model results and observations from natural subduction systems.

We have performed 3D subduction models to investigate the effects of along-strike variations in the overriding plate thickness on subduction dynamics and mantle flow. The numerical results show that these variations can lead to large variations of trench retreat velocities along-strike, with the trench above the thin overriding plate retreating faster. The difference of velocities along-strike produces arcuated trenches, which tend to rotate with the rotation axis located above the thicker part of the overriding plate. These results are obtained regardless of the far-field boundary conditions of the plates, although they are more evident for fixed plates, which is the scenario most representative for narrow subduction zones that are embedded within much larger and relatively immobile plates. Our models reveal complex mantle flow patterns due to the lateral gradients of trench velocities. The toroidal mantle flow is no longer symmetrical, but becomes more prominent around the slab edge on the side of the thin overriding plate and both toroidal flow cells returning from the sub-slab region tend to converge in the mantle wedge beneath the thin overriding plate, resulting in large regions of trench-parallel flow. Our models with a weak zone within the subducting plate show that along-strike variations in the overriding plate thickness can lead to vertical slab tearing across pre-existing weak zones, leading to slab segmentation. Slab tearing focuses the mantle upwelling through the tear, potentially triggering tear-related magmatism. Finally, we propose that our modelling results can explain to a large extent some aspects of the recent geodynamic evolution of numerous subduction zones.