Overriding Plate Thickness as a Controlling Factor for Trench Retreat Rates in Narrow Subduction Zones

Slab width is a significant factor in controlling subduction zone dynamics, particularly the retreat velocities, which tend to decrease with wider slabs. However, observations of natural narrow subduction zones reveal no correlation between slab width and trench velocities. This suggests that other factors may exert a greater influence. In this study, we employ 3D numerical subduction models to systematically assess the impact of slab width, strength of slab coupling to the lateral plate (LP), and overriding plate (OP) thickness on trench kinematics and geometry. Our models focus on narrow slabs (400–1,200 km), and the results demonstrate that, in the case of narrow subduction zones, the slab width has little effect on trench migration rates and the viscous coupling at the lateral slab edge is only important for very narrow subduction zones (≤800 km). Conversely, the OP thickness emerges as a crucial factor, with increasing plate thickness leading to a strong decrease in trench velocities. These findings provide an explanation for the observed trench velocities in natural narrow subduction zones, where an inverse relationship with OP thickness is evident. Furthermore, our study reveals that not only slab width, but also the OP thickness and the slab coupling to the LP, significantly influence trench geometry. Strong lateral coupling promotes the formation of concave trench geometries, while thick overriding plates favor the development of “w”‐shaped geometries. Overall, a comprehensive understanding of subduction processes necessitates considering the interplay between slab width, OP thickness, and slab coupling to the LP.

Among the factors controlling subduction dynamics, the slab width (W) has been shown to be of major importance (Bellahsen et al., 2005;F. Chen et al., 2022;Di Giuseppe et al., 2008;Guillaume et al., 2010Guillaume et al., , 2021;;Royden & Husson, 2006;Schellart et al., 2007;Stegman et al., 2006Stegman et al., , 2010;;Strak & Schellart, 2016).Earlier geodynamic studies have distinguished three types of trench curvature depending on slab width: concave toward the OP side for narrow slabs (W ≤ 1,500 km), "w"-shaped (also referred to as sublinear) for intermediate width slabs (W ∼ 2,000-3,000 km) and convex toward the OP for wide slabs (W ≥ 4,000 km) (Schellart et al., 2007;Stegman et al., 2010;Strak & Schellart, 2016), although the exact values differentiating these regimes depend somewhat on the specific conditions.In addition, the slab dip angle seems to be controlled by W, increasing (producing steeper slabs) with wider subducting plates (Schellart, 2004;Strak & Schellart, 2016).W also controls the subductioninduced mantle flow (Piromallo et al., 2006;Stegman et al., 2006), causing faster and more localized mantle upwellings near the lateral slab edges for wider slabs (Strak & Schellart, 2016).Regarding the slab kinematics, previous studies have shown that the trench retreat velocity (V T ) decreases as the slab becomes wider (F.Chen et al., 2022;Schellart, 2004;Schellart et al., 2007;Stegman et al., 2006), with wider subduction zones showing negligible trench retreat in their centers (Schellart et al., 2007).All these studies provide useful insights that help to understand the effect of W on subduction processes.However, most of these model setups that specifically focused on slab width did not incorporate an OP, which is known to affect subduction dynamics significantly (e.g., Butterworth et al., 2012;Capitanio et al., 2010;Z. Chen et al., 2015;Duarte et al., 2013;Hertgen et al., 2020;Magni, 2019;Magni et al., 2014;Rodríguez-González et al., 2014;Yamato et al., 2009).Additionally, as we will discuss in this work, the inverse dependence of V T on W predicted by some models is not observed in natural narrow subduction zones, suggesting that other factors may play a more relevant role on trench retreat velocities when narrow subduction zones are involved.Incorporating an OP in subduction models can help to better understand the effect of W on V T and the dominant factors controlling trench retreat rates in narrow subduction zones.
In this study, we have conducted self-consistent 3D numerical models of narrow subduction zones including the subducting and surrounding plates (lateral and overriding) to systematically evaluate the effect of slab width, OP thickness (T OP ) and coupling of the slab with the lateral plate (LP) on trench motion.Based on our geodynamic models and a comparison with observations in nature, this work provides new insights on the factors dominating trench retreat velocities in narrow subduction zones.

Numerical Method
The simulations have been performed with version 2.4.0 of the finite-element code ASPECT (Advanced Solver for Problems in Earth's ConvecTion) (Bangerth et al., 2022;Gassmöller et al., 2018;Heister et al., 2017;Kronbichler et al., 2012).We have used the Boussinesq approximation to solve the coupled conservation equations for mass, momentum and energy, which assumes that the density is constant in all equations except in the buoyancy term of the momentum equation (right hand term in 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 used to track materials and their properties throughout the simulations:

Model Setup
The initial 3-D model setup is shown in Figure 1.The initial geometry has been built using the Geodynamic World Builder version 0.4.0 (Fraters et al., 2019(Fraters et al., , 2021) ) and it measures 2,000 × 1,300 × 660 km in the x, y, and z directions.The width of the model box (1,300 km) is such that the separation between the lateral slab edge and the lateral sidewall is always at least half-width of the slab and larger than the depth of the upper mantle (660 km) in order to minimize the influence of the lateral wall on subduction dynamics (Schellart et al., 2007(Schellart et al., , 2010)).We have tested the effect of the distance between the lateral slab edge and the sidewall by running our model with the widest slab (W/2 = 600 km) with a box width of 2,600 km and have found that the trench retreat velocities just show a maximum increase of 0.2 cm/yr (∼6%) throughout the model evolution, thus confirming that our results are not significantly affected by the sidewall boundary condition.We use an adaptive mesh refinement (AMR) that allows us to locally increase the mesh resolution in areas of strong gradients of viscosity, temperature and/or composition, resulting in a maximum resolution of 7.8 × 5 × 2.6 km and a minimum resolution of 62.5 × 20.3 × 20.6 km.The highest mesh resolution is achieved within the SP and OP to accurately resolve the subduction zone interface, while the broadest meshing is only reached in the deep part of the upper mantle.
The initial geometry consists of an incipient subduction system with a 220 km-long slab and a dip angle of 40°.In this way, no external forces are imposed and the subduction is self-driven by the negative buoyancy of the slab.The subducting plate (SP) has an initial length of 1,000 km and is comprised of weak crust 10 km thick and lithospheric mantle 85 km thick.The weak crust has a low viscosity of 10 20 Pa s to allow for the decoupling of the slab from the OP and to facilitate subduction.This crust mimics a weak layer of sedimentary material at the top of the SP.The OP is also 1,000 km long and has an initial thickness of 70 km in the reference model, consisting of a 30 km thick continental crust and a 40 km thick lithospheric mantle.To avoid unnecessary model complexity, the LP is completely made of lithosphere 95 km thick, without any compositional stratification.We prescribe a weak zone (i.e., transform fault) 20 km wide of 10 20 Pa s in the reference model as a mechanical coupling of the LP with the other plates.The trench is initially located at x = 1,000 km, and we use particles placed along its length to track its movement over time.
The SP and OP are fixed at their trailing edges, simulating a subduction zone that is attached to much larger and relatively immobile plates.This scenario is different from plates with free trailing edges, corresponding with minimum resistance to lateral motion at spreading ridges.Most narrow subduction systems (as we will further discuss) are part of much larger and relatively stationary subducting and OPs.Therefore, the SP-OP-fixed scenario is well representative of most narrow subduction zones (e.g., Z. Chen et al., 2015;Yamato et al., 2009).Additionally, we have carried out some selected experiments with a free OP to examine its effect on subduction dynamics (see the Supporting Information S1).
The initial temperature profile increases linearly from 293 K at the surface to 1,643 K at the lithosphere-asthenosphere boundary.For the top and bottom boundaries we use Dirichlet temperature conditions.Mechanical boundary conditions are free-slip for all boundaries, except for the bottom of the box where a no-slip condition is imposed to model the upper-lower mantle discontinuity.This discontinuity is characterized by a factor of 10-100 viscosity increase (Davies & Richards, 1992;Kaufmann & Lambeck, 2000), which is simulated in our models as an impenetrable barrier to flow.This assumption is motivated by computational constraints and by the fact that narrow slabs that subduct through trench retreat, such as those modeled in this work, tend to stagnate at the 660 km discontinuity (e.g., Funiciello et al., 2003;Goes et al., 2017;Schellart et al., 2007).

Rheology
We use a diffusion creep rheology in which the viscosity is given by: where A is a prefactor of the equation, d is the grain size, m is the grain size exponent, E is the activation energy, V is the activation volume and R is the gas constant.
For the upper mantle (asthenosphere and lithospheric mantle) we use rheological parameters from wet olivine (Hirth & Kohlstedt, 2003).With these parameters, we obtain a viscosity of 1.57 ⋅ 10 20 Pa s at a depth of 150 km.
For the OP crust we adopt rheological parameters from wet anorthite feldspar (Bürgmann & Dresen, 2008).The viscosity is capped by a preset minimum and a maximum viscosity of 10 19 and 1.57 ⋅ 10 23 Pa s (1,000 times the upper mantle viscosity at a depth of 150 km), respectively, to preserve numerical stability.For the oceanic crust, we adopt a low constant viscosity of 10 20 Pa s.A low viscosity of 10 20 Pa s is also used in the transform fault of the reference model as a weak mechanical coupling at the lateral slab edge.All model parameters are listed in Table 1.

Effect of Slab Width and Coupling With the Lateral Plate
In order to focus on narrow subduction zones, we systematically changed the slab width with respect to the reference model (Experiment 2 in Table S1 in the Supporting Information S1) in the range of 400-1,200 km (experiments 1-5 in Table S1 in Supporting Information S1).For each experiment, we varied the mechanical coupling at the lateral slab edge by changing the viscosity of the transform fault (μ TF ).In addition to the value   S1 in Supporting Information S1) and μ TF = 10 22 Pa s (experiments 11-15 in Table S1 in Supporting Information S1).
The subduction dynamics is similar in all the experiments.Initially, the slab freely sinks into the upper mantle due to the density contrast between the SP and the sub-lithospheric mantle and the slab dip angle increases (Figures 2a  and 2c).During the sinking phase, the trench retreat velocities quickly increase (Figure 3d), accompanied by slab rollback and inducing toroidal flow around the slab edges (Figure 2).All models show the maximum V T when the slab tip reaches ∼400 km depth.Thereafter, trench velocities slightly decrease due to the interaction of the slab tip with the deep viscous upper mantle, and approach a roughly constant value (Figure 3d).The high trench retreat rates during subduction lead to the slab lying down horizontally above the deep viscous part of the mantle (Figures 2b and 2d).
Figure 3 shows the influence of SP width on trench shape and kinematics.Models develop two types of trench geometries in the center of the subduction zone depending on W. The trench in models with W ≤ 1,000 km shows a concave geometry toward the OP, with the trench in the center of the subduction zone retreating faster than its edges (Figures 3a and 3b).This characteristic geometry is attained earlier in models with smaller W. For example, the concave geometry for W = 400 km is almost achieved in 2 Myr, while for W = 600 km it is not attained until 10 Myr (Figures 3a and 3b).For the model with W = 1,200 km, the trench geometry is almost rectilinear in its center with concave edges up to ∼8 Myr, and thereafter adopts a "w"-shape, with retreat velocities being higher in between the lateral slab edge and the center of the subduction zone (Figure 3c).
Trench retreat velocities do not change significantly with the slab widths tested in this work.During the first time steps (0-4 Myr), the trench steadily retreats with similar velocities in all cases.Afterward, trench velocities reach a maximum and finally approach a roughly steady value (Figure 3d).The highest steady velocity is found for W = 600 km (∼4.3 cm/yr) and decreases for wider slabs up to a steady velocity of ∼2.9 cm/yr for W = 1,200 km (Figure 3d).The highest maximum trench velocity reached during each simulation is found for W = 400 km (∼4.9 cm/yr) and slightly decreases for wider slabs (Figure 3d).The trench in panel (d) seems not to be connected to the LP due to some amount of subducting plate shortening in the trenchparallel direction, which is not captured in the temperature cutouts.
Concerning the mechanical coupling of the slab with the LP, we find that the model with W = 1,200 km exhibits different trench geometries depending on the viscosity of the transform fault: we obtain a "w"-shaped trench for a weak coupling (μ TF = 10 20 Pa s) (Figure 4a) while a concave trench is developed for a strong coupling (μ TF = 10 22 Pa s) (Figure 4c).For the intermediate case (μ TF = 10 21 Pa s), we obtain a behavior in between those of the two previous models.The trench remains approximately straight with concave edges and only starts to develop a "w"-shape after ∼14 Myr (Figure 4b).A weak lateral coupling leads to some amount of trench shortening in the trench-parallel direction (Figures 3a-3c and 4a), in contrast with models with a strong lateral coupling (Figure 4c; see discussion by Yamato et al. (2009) about trench shortening).
The increase in the viscous resistance at the lateral slab edge causes the trench retreat velocity at the center of the subduction zone to decrease for W ≤ 800 km, but does not change significantly for W ≥ 800 km (Figures 4d and  4e).Decreases of ∼1.7 and ∼1 cm/yr in the maximum trench retreat velocity are observed for W = 400 km and W = 600 km respectively when the viscosity increases from μ TF = 10 20 Pa s to μ TF = 10 22 Pa s, while only maximum differences of ∼0.2 cm/yr are observed for W = 1,200 km (Figure 4e).Additionally, the highest maximum of V T at the center of the subduction zone is found for W = 400 km when the viscosity of the transform fault is μ TF = 10 20 Pa s, but for W = 600 km when μ TF = 10 21 Pa s and for W = 800 when μ TF = 10 22 Pa s (Figure 4e).

Effect of Overriding Plate Thickness
We have varied the OP thickness with respect to the reference model in the range of 40-100 km to study its effect on trench motion (Experiments 2 and 16-21 in Table S1 in Supporting Information S1).During the free sinking phase, thin OPs promote fast slab rollback with small dip angles (Figure 5a) while a thick OP leads to the slab sinking almost vertically without rollback (Figure 5c).Consequently, the toroidal mantle flow around the slab is less intense for thicker OPs.The sustained high retreat velocities during all the simulations with a thin OP lead to an elongated flat-lying slab with a very shallow dip angle (Figure 5b).In contrast, for thick OPs, the slab pull force increases as the slab sinks vertically, progressively increasing the retreat velocity and resulting in a stepped slab with its deepest part becoming flattened (Figure 5d).
As for the trench geometry, the trench develops a concave shape for all OP thicknesses (Figures 6a-6d) except for T OP = 100 km, where a "w"-shape develops after about 10 Myr of evolution (Figure 6b).Significant trench lateral shortening is also observed as T OP increases (Figure 6b).Geochemistry, Geophysics, Geosystems 10.1029/2023GC011345 The trench retreat velocity significantly decreases with thicker OPs (Figure 6e).The maximum V T decreases from 6.6 cm/yr for T OP = 40 km to 4.3 cm/yr for T OP = 80 km (Figure 6e).

Comparison With Previous Studies
Our geodynamic models provide new insights into the effect of W, OP thickness and viscous resistance at the lateral slab edge on trench kinematics, and directly complements previous work on this topic.The numerical models of Schellart et al. (2007) and Stegman et al. (2010) suggest values of W ≤ 1,500 km, W ∼ 2,000-3,000 km and W ≥ 4,000 km leading to concave, "w"-shaped and convex trenches, respectively, while the analog models of Strak and Schellart (2016) indicate values of W = 2,000-2,500 km and W ≥ 3,000 km for "w"-shaped and convex geometries respectively.Király et al. (2017) tested plate widths of 1,000, 1,500, and 2,000 km and their models predict concave geometries for the narrow case (W = 1,000 km) and "w"-shaped trenches for the wider cases (W = 1,500 km and W = 2,000 km).In contrast, the analog modeling of Guillaume et al. (2021) predicts concave and "w"-shaped geometries for much wider slabs (W = 2,000 and W = 4,000 km respectively).Similarly, the recent 3-D spherical shell numerical models of F. Chen et al. (2022) place the transition from concave trenches to "w"-shaped geometries at W = 2,400 km for their reference case.Our results do not show convex geometries, which is expected as we focus only on narrow subduction zones and these geometries are found for very wide slabs (Schellart et al., 2007;Strak & Schellart, 2016).We do, however, identify two of the three types of trench geometries identified by Schellart et al. (2007), Stegman et al. (2010), and Strak and Schellart (2016), but our models show a "w"-shaped geometry for a narrower slab (W = 1,200 km, Figure 3c) than in previous studies (W ≥ 1,500 km).These discrepancies may be partly explained by the far-field boundary conditions applied to the SP.In the models of Schellart et al. ( 2007 2022), the SP was free to move at its trailing edge (which is more appropriate for wide slabs), while in our models it is fixed (a suitable condition for narrow slabs).In the case of a fixed SP, the subduction process is primarily governed by slab rollback, while a free trailing edge allows the SP to advance toward the trench.The models of Király et al. (2017) and Guillaume et al. (2021) both used fixed subducting plates at their trailing edges, but only the findings of Király et al. (2017) ("w"-shaped trenches for W = 1,500 km) are close to the results of this work ("w"-shaped trenches for W = 1,200 km).Discrepancies may come from the use of different viscosity ratios between the slab and the upper mantle, but another likely factor that could explain these variations is the influence of the OP, which was not included in these previous studies: our results show that increasing T OP in models with the same slab width promote "w"-shaped trenches (Figure 6b).
Regarding trench migration velocities, previous modeling studies have found that V T decreases with increasing W (e.g., F. Chen et al., 2022;Schellart et al., 2007;Stegman et al., 2006).However, this behavior is not maintained for narrow subduction zones.Our results show a direct dependence of V T on W for W ≤ 600 km when the mechanical coupling of the slab with the LP is moderate (μ TF = 10 21 Pa s) and for W ≤ 800 km when the coupling is stronger (μ TF = 10 22 Pa s) (Figure 4e).This positive correlation between V T and W has also been suggested in previous works for W ≤ 500-600 km (Schellart, 2004;Stegman et al., 2006) and for W ≤ ∼1,000 km (Strak & Schellart, 2016).Our models using an intermediate coupling (μ TF = 10 21 Pa s) are in agreement with the results of Stegman et al. (2006) finding a maximum retreat velocity for W = 600 km.The discrepancy on the W showing maximum V T between the present work and previous studies suggests a dependence on different model parameters and setup.For example, our models show that the W for which V T peaks increases as the viscous resistance at the lateral slab edge becomes stronger (Figure 4e).More research is needed to clarify the dynamics and factors affecting this phenomenon.Here we propose that a peak V T is observed in Figure 4e due to an energy balance (Magni et al., 2014) between the available gravitational potential energy of the slab E pot and the energy required for frictional dissipation in the mantle E diss,m and at the transform fault E diss,TF .For relatively wide slabs (W ≥ 800 km), E diss,TF is relatively unimportant, E pot linearly increases with slab width, while E diss,m increases faster than E pot .Therefore, trench retreat is slower for wider slabs as the frictional dissipation in the mantle becomes dominant.For narrow slabs (W ≤ 800 km), E pot and E diss,m are both small, and E diss,TF becomes dominant.Since E pot decreases for narrowing slabs, while E diss,TF is independent of slab width, narrow slabs retreat slower than wider slabs for the range of W ≤ 800 km.
The effect of the strength of transform faults on trench kinematics remains poorly studied and the coupling of the slab with the LP is usually assumed to be either weak (e.g., Chertova et al., 2018;Király et al., 2021) or strong (e.g., Yamato et al., 2009).Govers and Wortel (2005) made the first study investigating the dynamics at the subduction zone edges based on observations in nature and numerical models.Later on, Hale et al. (2010) evaluated with 3D numerical subduction models how the strength at the lateral slab edge affects the subduction dynamics.Their models showed that increasing the viscous resistance at the transform fault results in lower trench retreat velocities and more prominent trench curvatures, which was also found by Nijholt and Govers (2015).The recent modeling of Schliffke et al. (2022) further highlights the importance of the transform fault strength in subduction dynamics, showing that back-arc ridge jumps depend on the ratio of transform fault versus OP strengths.The results of Hale et al. (2010) are in agreement with our work showing that a strong coupling with the LP increases viscous resistance at the lateral slab edge, enhancing concave trench geometries (Figures 4a-4c) and reduces the trench retreat velocity for W ≤ 800 km.However, the fact that the three curves in Figure 4e tend to converge with W approaching 1,200 km suggests that effect of the strength of lateral coupling on V T is only significant for very narrow subduction zones.
Finally, the significant decrease of V T with increasing T OP reported in Figure 6e demonstrates the strong impact of the OP thickness on trench retreat velocities.This result is in agreement with previous 3D and 2D models (Butterworth et al., 2012;Gea et al., 2023;Hertgen et al., 2020;Sharples et al., 2014;Yamato et al., 2009) finding that thin/weak OP favors faster trench velocities and rollback compared to a thick/strong OP.This result is enhanced by the fixed boundary condition of the OP at its trailing edge.With a fixed-OP condition, slab rollback is entirely accommodated by OP extension, and thus stronger (thicker) OPs will provide more resistance to extension and, hence, to trench retreat.This is very different from a scenario with a free OP, where slab rollback will be accommodated by both OP extension and OP motion toward the trench.For example, the SP-OP-free models of Meyer and Schellart (2013) found fast trench retreat velocities for both thin and thick OP models and showed that increasing the thickness of the OP just retards trench retreat, which is gradually further accommodated by trenchward OP motion.Other models, though, have found that increasing the thickness of the OP reduces rollback rates, regardless of the condition of the OP at its trailing edge, but with a fixed OP enhancing the dependence of trench retreat velocities on the mechanical properties of the OP (Capitanio et al., 2010;Holt et al., 2015).Our supplementary models with a free OP support these results, showing smaller but still significant decreases in retreat velocities with increasing T OP (see Figure S4 in Supporting Information S1).
The short periods of trench advance found in this work for T OP ≥ 70 km (Figure 6e) are also consistent with the models of Sharples et al. (2014) showing that trench advance only occurs in models with thick and/or strong OPs.The impact of the OP thickness on slab dip angle has also been studied in previous models, some finding that the slab dip angle does not change significantly (e.g., Meyer & Schellart, 2013), while others reporting an increase in slab dip angle at depth with increasing T OP (e.g., Holt et al., 2015;Sharples et al., 2014).Our model results show that T OP does have an effect on slab dip angle through its control on trench retreat velocities: thin OP enhance high trench retreat velocities resulting in small slab dip angles (Figure 5a), while a reduction of retreat velocities due to thicker OP results in an increase of slab dip angle (Figure 5c).The analog models of Meyer and Schellart (2013) also showed that the trench curvature decreases as T OP increases.Although this is not observed in our models (Figure 6), we do find that the thickness of the OP affects the trench geometry, with an OP 100 km thick resulting in a "w"-shaped trench.

Comparison With Nature
Natural subduction zones on Earth are influenced by a large number of factors that are not included in our models, such as slabs of non-uniform age, complex slab morphologies or plate velocities.Similarly, interaction of subduction zones on Earth with mantle flow coming from nearby subduction zones elsewhere or from other processes such as mantle plumes can also affect subduction evolution.Thus, comparing our results with observations on Earth should be approached with caution, being aware of the inherent limitations.Nevertheless, our geodynamic model predictions are useful to explain some observations on natural narrow subduction zones.
The concave trench geometries predicted by our geodynamic models are widely observed on Earth for narrow subduction zones such as Scotia, South Shetland, Calabria or Gibraltar (Figure 7a).Our results have  2019)).(c) V T against W for all subduction zones in Earth with W ≤ 1,000 km (data from Schellart et al. (2007)).(d) V T against T OP for all subduction zones in Earth with W ≤ 1,000 km (data can be found in Table S2 in Supporting Information S1).Note: the T OP plotted values are subjected to significant uncertainty, which in turn depends on the estimation method.The reader is referred to the main text for a more detailed discussion on this point.Gray dashed lines are the least-squares linear best-fit lines and R 2 is the coefficient of determination.

Geochemistry, Geophysics, Geosystems
10.1029/2023GC011345 demonstrated that including an OP promotes the formation of "w"-shaped geometries.For example, the model with W = 1,200 km and μ TF = 10 20 Pa s develops a "w"-shaped trench (Figure 3c) for much smaller W than any of the previous studies not incorporating an OP (e.g., F. Chen et al., 2022;Schellart et al., 2007;Stegman et al., 2010;Strak & Schellart, 2016).Increasing the thickness of the OP has also been shown to affect the geometry of the trench, with an OP 100 km thick leading to a "w"-shaped geometry (Figure 6b).This effect of the OP could explain why some natural subduction zones of narrow to intermediate widths develop "w"-shapes.For example, the Manila trench (W = 1,000 km) has a geometry between concave and "w"-shaped (Figure 7b).The main difference with our results is that the Manila trench is convex toward the OP at its northernmost part, while for our models the trench geometry is symmetric and convex in its center.We attribute these differences to the use of an idealized model setup with uniform subducting and OPs, which results in symmetric trenches.
Additionally, the Manila subduction zone is affected by the Ryukyu subduction system and the convex geometry in its northernmost extremity may be explained by slab-slab interactions with the Ryukyu slab (Király et al., 2016).
The decrease in V T with increasing W observed in the present work for W ≥ 800 km (Figures 3d and 3e) and in previous works (e.g., F. Chen et al., 2022;Royden & Husson, 2006;Schellart et al., 2007;Stegman et al., 2006) is in general agreement with observations in narrow and wide subduction zones.This result explains why narrow subduction zones (e.g., Calabria, South Shetland, Halmahera) exhibit relatively high V T , while wide subduction zones (e.g., Melanesia, South America) are essentially stationary (Schellart et al., 2007).However, when this comparison is restricted to narrow subduction zones, no correlation between W and V T is observed (R 2 = 0.015, Figure 7c).For the same slab width (∼900 km), the Makran subduction zone has a V T of 0.2 cm/yr while the Trobriand system has a V T of 7.6 cm/yr (Schellart et al., 2007).Our results showing that W has little effect on V T for narrow subduction zones, with variations in V T less than ∼2 cm/yr (Figure 3d), are in agreement with the lack of correlation between W and V T observed in nature and provide an explanation for these observations.
Conversely, our models reveal that the effect of OP thickness is much stronger, showing an inverse dependence of V T on T OP (Figure 6e).This negative correlation is also observed in nature for narrow subduction zones (R 2 = 0.643, Figure 7d) and brings an explanation of why two subduction zones with the same W (e.g., Makran and Trobriand) show such different V T .Nevertheless, it is essential to acknowledge the substantial uncertainty associated with estimating OP thickness in natural settings.The accuracy of OP thickness determination in each subduction zone is inherently related to the specific geophysical method used for its assessment.Receiver functions emerge as the best method for lithosphere thickness estimation, but only two data points, namely Gibraltar (Molina-Aguilera et al., 2019) andSouth Shetland (Parera-Portell et al., 2021), have been based on receiver function analysis.For the rest of the subduction zones, the data have been based on other, much less accurate geophysical methods, including tectonic reconstruction studies for Calabria (Rosenbaum & Lister, 2004), numerical models that best fit the observations for Halmahera (Zhang et al., 2017) and North Sulawesi (Dong et al., 2022), seismic crustal studies for Puysegur (Shuck et al., 2021) and Makran (Motaghi et al., 2020), mantle tomography studies for Sangihe (Fan & Zhao, 2018), heat flow studies for Trobriand (Martinez et al., 2001) and the global thermochemical model of the lithosphere and underlying upper mantle WINTERC-G for Scotia and Manila (Fullea et al., 2021).Thus, although the negative correlation between V T and T OP appears evident for narrow subduction zones at first order, the significant uncertainty in the data must be considered.Future studies are needed to better estimate the OP thickness and to further confirm this V T -T OP correlation.
We note that our far-field boundary conditions (SP-fixed and OP-fixed) limit our results to these specific conditions, which represent a scenario where subduction zones are attached to larger and relatively immobile plates.However, our SP-OP-fixed condition is primarily consistent with the natural examples presented in this study, since most narrow slabs are embedded within larger subducting plates and have a relatively stationary OP.The Mediterranean subduction zones (Gibraltar, Calabria, and Hellenic) are good examples of this scenario, as both Eurasia and Africa are moving much slower than the subduction zones inside.The Makran subduction zone is also comparable to our SP-OP-fixed condition as the slab is embedded within the much larger SP (Arabia) and OP (Eurasia).Among the subduction systems presented in this work, those with an OP moving trenchward at high rates (which ideally would be best comparable to a SP-fixed OP-free scenario) include the subduction zones of Puysegur beneath the overriding Pacific plate, Trobriand beneath the Australian plate, Manila below the Philippine plate and Scotia under the South Sandwich plate.Even in these cases, it is not straightforward that an SPfixed OP-free condition would be the most applicable to these scenarios because the high velocities of the OPs may be partly attributed to other nearby subduction zones (e.g., Tonga for Puysegur or the Melanesia subduction zone for Trobriand).In any case, our supplementary models reveal that even including a free OP, the effect of T OP on trench retreat velocities remains more relevant than W for narrow subduction zones (see Figure S4 in Supporting Information S1).
Finally, it is worth noting that although the model setup used in this study includes a continental OP, the qualitative results are also valid for an oceanic OP, as it is the case for some of the above-mentioned natural examples.
The thermal thickness and strength of the oceanic lithosphere are directly related to its age, and therefore a thick (old, cold, and strong) oceanic lithosphere will provide more resistance to internal deformation and trench motion than a thin oceanic OP.

Model Limitations
Our numerical models contain some limitations that should be taken into account.We employ a simplified rheological structure.The advantage of including the three plates and allowing for their relative motions and interaction has the drawback of requiring an increased mesh resolution in the weak zones between the plates.While this is satisfactorily assessed thanks to the AMR feature of ASPECT code, simulations are still expensive in terms of computation time.In order to keep acceptable computation times which allow for a deeper exploration in the parameters space, we used a linear temperature-and pressure-dependent viscosity.The use of simplified rheologies has been shown to satisfactorily reproduce realistic subduction dynamics and slab morphologies (Capitanio et al., 2007;Di Giuseppe et al., 2008;Ribe, 2010;Strak & Schellart, 2018).More complicated (computationally much more expensive) rheologies such as non-linear viscoplastic rheology promote strain localization and back-arc spreading in the OP, resulting in higher rates of trench retreat (Pusok et al., 2018).Still, a thicker and therefore stronger OP would delay or even hinder strain localization compared to a thinner plate and therefore our qualitative results of slower trench retreats for thicker OPs would not change.
We do not consider the effect of nearby subduction zones on our individual subduction models, although the interaction between nearby subduction zones naturally exists in nature.Indeed, it has been shown that slab-slab interactions may influence the trench kinematics, the induced mantle flow, and the deformation of the subduction zones involved (Király et al., 2021).Interactions between narrow slabs in nature include the coupled evolution of the Ryukyu and Manila subduction zones, the divergent double subduction of Halmahera and Sangihe and the recent simultaneous subduction of the Mediterranean slabs.However, we argue that the influence of external mantle flow coming from nearby subduction zones is of secondary importance compared to the negative buoyancy of the slab (Forsyth & Uyeda, 1975).
We fix the trailing edge of the SP and OP (apart from a few free OP models in Supporting Information S1), thereby neglecting LP velocities.In such scenario, the SP subducts almost entirely by slab rollback and is accommodated by extension of the upper plate.This is a simplification of subduction zones in nature, where the plate velocities are not negligible.For example, previous subduction models (and data from subduction zones in nature) have shown that the SP velocity increases with increasing slab width (Schellart et al., 2010).The OP also contributes to the subduction rate, with a global average trench-normal velocity of 0.81 cm/yr (Schellart, 2023).However, most narrow subduction zones, which are the focus of this work, are embedded within much larger and relatively immobile plates, with the subducting and OPs usually having low trench-normal velocities.Therefore, the behavior of narrow subduction zones can be best approximated by the SP-OP-fixed condition that we use in this work (Z.Chen et al., 2015;Yamato et al., 2009).
Our models are restricted to subduction in the upper mantle since the 660-km discontinuity is simulated as an impenetrable barrier, even though it is known that some slabs are able to penetrate into the lower mantle with high dip angles (e.g., Hellenic, Peru) (e.g., Bijwaard et al., 1998;Fukao et al., 2001;van der Hilst et al., 1991).This simplification implies that the interaction of the slab with the upper-lower mantle discontinuity is amplified, for example, by further decelerating the slab sinking.However, previous 3D subduction modeling studies that include a high-viscosity lower mantle have shown that those slabs that subduct primarily by trench retreat tend to stagnate and flatten above the discontinuity (e.g., Funiciello et al., 2003;Griffiths et al., 1995;Schellart et al., 2007;Sharples et al., 2014;Stegman et al., 2006Stegman et al., , 2010)).Therefore, the no-slip boundary condition at the bottom of the models to simulate the 660-km discontinuity seems to be a reasonable approximation for our models in which the subducting slabs are dominated by trench retreat.

Conclusions
The self-consistent 3-D numerical subduction models presented in this work shed light on some factors (slab width, OP thickness and viscous coupling at the lateral slab edge) controlling trench kinematics in narrow subduction zones, and help to explain some observations of natural subduction processes.As opposed to what happens in wide subduction zones, our models show that the slab width has little effect on trench retreat velocities for narrow subduction zones, which explains the lack of correlation between these parameters observed in nature.
In contrast, from our models and observations in nature we conclude that the thickness of the OP is the main controlling factor on trench retreat velocities for narrow subduction zones, with velocities strongly decreasing as the thickness increases.Stronger coupling to the LP has been found to increase the trench curvature, favoring concave geometries, but its effect on decreasing trench retreat velocities is only significant for very narrow subduction zones.Finally, our work reveals that the inclusion of an OP plays an important role in modulating trench geometries, facilitating the formation of "w"-shaped geometries, which are predicted for smaller slab widths than in previous studies.
μ TF = 10 20 Pa s used in the reference model, we tried μ TF = 10 21 Pa s (experiments 6-10 in Table

Figure 2 .
Figure 2. 3D perspective view of the subduction zone for (a, b) W = 400 km and (c, d) W = 1,200 km.The figures are cutouts of the temperature between 293 and 1,500 K. White lines with triangles on the surface mark the trench position.Black arrows show the horizontal velocity field around the slab at 250 km depth.Note that the lateral plate (LP) is shown as transparent.The trench in panel (d) seems not to be connected to the LP due to some amount of subducting plate shortening in the trenchparallel direction, which is not captured in the temperature cutouts.

Figure 3 .
Figure 3. (a-c) Evolution of the subduction trench for three simulations with (a) W = 400 km, (b) W = 600 km, and (c) W = 1,200 km.Note that the entire trench is plotted but we only model half of the subduction zone.(d) V T over time (measured in the center of the subduction zone) for simulations with different W.

Figure 4 .
Figure 4. (a-c) Evolution of the subduction trench for three simulations with W = 1,200 km and (a) μ TF = 10 20 Pa s, (b) μ TF = 10 21 Pa s, and (c) μ TF = 10 22 Pa s.Note that the entire trench is plotted but we only model half of the subduction zone.(d) V T over time (measured in the center of the subduction zone) for simulations with different μ TF .Solid lines indicate models with W = 600 km and dashed lines indicate models with W = 1,200 km.(e) Maximum V T at the center of the subduction zone for different W and different μ TF .
We observe local peaks of trench velocities for the models with T OP = 40 km and T OP = 50 km at 16-18 Myr and 18-20 Myr respectively, which are likely caused by transient changes in the mantle flow.Short periods of trench advance at the beginning of the simulation occur for T OP ≥ 70 km.The time during which these models undergo trench advance also increases with the T OP : 2 Myr for T OP = 100 km and less than 1 Myr for T OP = 80 km (Figure 6e).

Figure 5 .
Figure 5. 3D perspective view of the subduction zone for (a, b) T OP = 50 km and (c, d) T OP = 100 km.The figures are cutouts of the temperature between 293 and 1,500 K. White lines with triangles on the surface mark the trench position.Black arrows show the horizontal velocity field around the slab at 250 km depth.Note that the lateral plate is shown as transparent.See comment about trench representation in the caption of Figure 2.

Figure 6 .
Figure 6.(a-d) Evolution of the subduction trench for four simulations with W = 600 km and (a) T OP = 50 km, (b) T OP = 100 km, (c) T OP = 70 km, and (d) T OP = 90 km.Note that the entire trench is plotted but we only model half of the subduction zone.(e) V T over time (measured in the center of the subduction zone) for simulations with different T OP .

Figure 7 .
Figure 7. (a, b) Subduction fronts at present-day for the (a) Gibraltar subduction system (from Gutscher et al. (2012)) and the (b) Manila subduction zone (from Qiu et al. (2019)).(c) V T against W for all subduction zones in Earth with W ≤ 1,000 km (data fromSchellart et al. (2007)).(d) V T against T OP for all subduction zones in Earth with W ≤ 1,000 km (data can be found in TableS2in Supporting Information S1).Note: the T OP plotted values are subjected to significant uncertainty, which in turn depends on the estimation method.The reader is referred to the main text for a more detailed discussion on this point.Gray dashed lines are the least-squares linear best-fit lines and R 2 is the coefficient of determination.

Figure 1. Three
-dimensional model setup and boundary conditions.The setup includes a subducting plate (SP), an overriding plate (OP), and a lateral plate (LP).The SP is made of oceanic crust 10 km thick and lithosphere 85 km thick.The OP consists of a 30 km thick crust and 40 km thick lithosphere in the reference model.The LP is comprised of 95 km thick lithosphere and it is separated from the OP and SP by a 20 km wide weak zone (i.e., transform fault).W indicates slab width and the red line with triangles marks the initial trench position at x = 1,000 km.The trench is plotted for the whole subduction zone in Figures3, 4, and 6.The temperature is fixed to 293 K at the top boundary and 1,643 K at the bottom boundary.All boundaries are free-slip except the bottom of the box, where a no-slip condition is applied to simulate the upper-lower mantle discontinuity.Note that only half of the subduction zone is modeled due to the symmetry of the problem.

Table 1
Model Parameters