How Slab Age and Width Combine to Dictate the Dynamics and Evolution of Subduction Systems: A 3‐D Spherical Study

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

The age of a subducting slab determines its thermal structure, which controls slab thickness, density and rheology. In turn, these control slab buoyancy and strength, which combine to determine the rate of trench retreat (e.g., Bellahsen et al., 2005;Capitanio et al., 2007;Di Giuseppe et al., 2008;Garel et al., 2014;Goes et al., 2017;Ribe, 2010;Schellart, 2008;Stegman, Farrington, et al., 2010). It is well-established that trench-motion history correlates with slab morphology (e.g., Faccenna et al., 2001;Goes et al., 2017;van der Hilst & Seno, 1993). Goes et al. (2008) suggest that older, colder, oceanic lithosphere is stronger due to the temperature dependence of viscosity, and that this drives significant trench retreat, with slabs more likely to lie flat at the mantle transition zone; conversely, younger lithosphere is weaker and subducts with less trench retreat, tending to buckle at the mantle transition zone. This direct link between slab age and the style of slab-transition zone interaction is supported by laboratory and numerical simulations (e.g., Bellahsen et al., 2005;Capitanio et al., 2007;Funiciello et al., 2008;Garel et al., 2014;Goes et al., 2017;Schellart, 2004).
Slab width also plays an important role in determining the evolution of subduction zones, affecting the shape and curvature of the trench, by influencing the rate of trench retreat (e.g., Bellahsen et al., 2005;Schellart et al., 2007;Stegman, Schellart, & Freeman, 2010;Stegman et al., 2006;Strak & Schellart, 2016). Schellart et al. (2007) advocate an inversely proportional relationship between trench migration rates and the width of subduction zones. In their models, wider slabs develop upper mantle stagnation zones, where the center of the trench exhibits negligible trench retreat relative to its edges. Such a relationship is observed in large subduction systems such as the South American subduction zone, and could be responsible for the varying styles of slab morphology observed along wide trenches.
Plate interactions and lateral variations in plate properties also influence the geometry and evolution of subducting slabs. The structure and motion of overriding plates have an effect on slab dip, trench migration, and slab interaction with the mantle transition zone (e.g., Capitanio et al., 2010;Garel et al., 2014;Heuret et al., 2007;Jarrard, 1986;Lallemand et al., 2005;van Dinther et al., 2010). Interactions with nearby past or ongoing subduction zones may also affect the sinking velocity and slab penetration of the mantle transition zone (e.g., Becker & Faccenna, 2011;Čížková & Bina, 2019;Fukao & Obayashi, 2013;Fukao et al., 2009). The subduction of locally thickened oceanic lithosphere, such as oceanic plateaus, aseismic ridges or seamount chains, has been proposed to influence the shape of the trench and change the geometry of subducting slabs (e.g., Capitanio et al., 2011;Cross & Pilger, 1982;Gutscher, Malavieille, et al., 1999;Martinod et al., 2005;Suchoy et al., 2022). The higher compositional buoyancy of oceanic plateaus and ridges could resist slab sinking into the mantle and, thus, potentially lead to flat slab subduction (e.g., Mason et al., 2010;van Hunen et al., 2002). These factors add to the complexity of observed slab morphology, and may explain why no simple correlations can be made with plate properties.
Existing numerical and laboratory studies in an enclosed Cartesian domain provide valuable insight into the sensitivity of slab morphology to a number of these controlling parameters. However, such models neglect the role of Earth's sphericity, the geometric consequences of which results in curved subduction zones, the build-up of lateral strain, buckling along the trench, and an increase in the geometric stiffness of slabs (e.g., Bayly, 1982;Frank, 1968;Fukao et al., 1987;Laravie, 1975;Mahadevan et al., 2010;Schettino & Tassi, 2012;Tanimoto, 1998;Yamaoka, 1988).
In this paper, our aim is to investigate the effect of subducting plate age and width on slab morphology using 3-D spherical shell numerical models of free-subduction. Only a limited number of spherical subduction models have investigated the dynamics of subduction, and these considered only isoviscous plates (e.g., Butterworth et al., 2012Butterworth et al., , 2014Chamolly & Ribe, 2021;Morra et al., 2009Morra et al., , 2012. To more realistically approximate subduction on Earth, especially for wider plates, we will present the first multi-material subduction models with a composite plate rheology in a spherical shell domain. We consider plate ages (estimated with three combinations of plate thickness and density) ranging from young (∼10 Myr, estimated using a half space cooling model) to old (∼140 Myr), noting that the range of subducting plate ages on Earth at the present day is 0-160 Myr (e.g., Müller et al., 2016). Motivated by a compilation of global trench lengths at the present day (Heuret et al., 2011) and Cenozoic Era reconstructions (Müller et al., 2016), we examine trench widths of 1,200, 2,400, 3,600, and 4,800 km. As illustrated in Figure 1, at the present day, most trenches are less than 5,000 km wide, with mean and median widths of 1,940 and 1,130 km for the data set of Heuret et al. (2011) and 2,000 and 1,230 km for the data set of Müller et al. (2016). At 60 Ma (where there are less data on narrow trenches), mean and median widths are 3,520 and 3,430 km (Müller et al., 2016). It is noteworthy that very few trenches exceed 6,000 km in width: at present, the South America trench is 7,060 km wide, and the Andaman-Sumatra-Java-Timor trench exceeds 6,040 km width; at 60 Myr, the South America and Aleutian trenches were ∼7,100 and ∼6,070 km wide, respectively (Müller et al., 2016).
In this paper, we present, for a systematic set of simulations, a quantitative comparison to demonstrate how slab thickness and density (approximating the buoyancy of slabs of different age) and slab width affects the evolution of subducting slabs. We then discuss the likely mechanisms that can explain our results and their implications for an improved understanding of subduction on Earth.  Heuret et al. (2011) and combined here for commonly recognized continuous trenches; (b and c) trench length histograms based on global tectonic plate reconstructions, at the present day and at 60 Ma, respectively (Müller et al., 2016). The reconstructed trenches are segmented based on changes in lower or upper plate characteristics (green bars with orange outline). Based on lower plate properties, we merged segments that likely subduct as coherent slabs (purple bars). (d) Map of subduction zones from Heuret et al. (2011). (e and f) Map of trench segments based on plate reconstruction at present and at 60 Ma by Müller et al. (2016). In (e and f), black crosses represent edges of subduction zones. Trench segments within each subduction zone were merged when calculating trench length.

Governing Equations and Numerical Strategy
We simulate multi-material free-subduction of a composite visco-plastic plate into an ambient mantle, in a 3-D hemispherical shell domain, which extends from the surface to the core-mantle-boundary (CMB) at a depth of 2,890 km. Assuming incompressibility, the governing equations for this problem are the continuity equation, (1) the conservation of momentum equation for infinite Prandtl number, and an advection equation for composition, where u is velocity, p the pressure, μ the viscosity, Δρ the density difference between plate and mantle, g gravity acceleration, ̂ unit vector in the direction opposite gravity, and Γ the material volume fraction (Γ = 1 in a region occupied by a given material and Γ = 0 elsewhere).
Simulations are carried out using Fluidity (e.g., Davies et al., 2011Davies et al., , 2016Kramer et al., 2012), an anisotropic, adaptive, unstructured mesh computational modeling framework supporting finite element and control volume discretizations. Fluidity has recently been validated for simulations in a spherical shell domain against an extensive set of analytical solutions introduced by . It has also been validated for visco-plastic simulations like those examined herein (e.g., Le Voci et al., 2014;Tosi et al., 2015). In the context of this study, Fluidity has several favorable features. The framework: (a) uses an unstructured mesh, which enables the straightforward representation of complex geometries and materials; (b) dynamically optimizes this mesh, across parallel processors, providing increased resolution in areas of dynamic importance, thus allowing for accurate simulations across a range of length-scales, within a single model; (c) enhances mesh optimization using anisotropic elements; (d) can employ a free-surface boundary condition, which is important for correctly capturing slab decoupling from the surface (Kramer et al., 2012); (e) utilizes the highly scalable parallel linear system solvers available in PETSc (Balay et al., 1997(Balay et al., , 2021a(Balay et al., , 2021b, which can efficiently handle sharp, orders of magnitude variations in viscosity; and (f) has a novel interface-preservation scheme, which conserves material volume fractions and allows for the incorporation of distinct materials (Wilson, 2009). In this study, Fluidity's adaptive mesh capabilities are utilized to provide a local resolution of 3 km in regions of dynamic significance (i.e., at the interface between materials and in regions of strong velocity and viscosity contrasts), with a coarser resolution of up to 300 km elsewhere. It is this adaptive mesh functionality that makes our spherical shell simulations computationally tractable.

Geometry, Boundary Conditions and Material Properties
The configuration of our models is inspired by Stegman, Farrington, et al. (2010) and Garel et al. (2014). Simulations exploit the symmetry of the system to halve the computational domain's extent, modeling half of the plate width in a hemispherical shell. When non-dimensionalized, the hemispherical shell has an outer radius of 2.22 and an inner radius 1.22, thus the computational domain has a thickness of 1. Models have a free-surface boundary condition on the outer surface, and a free-slip condition on the symmetry plane and CMB, as illustrated in Figure 2a. Gravity acts radially, toward the center of the sphere. The subducting plate length (L) is 2,200 km. The initial slab tip geometry is prescribed with a bending radius of 250 km and an angle of 77°, resulting in a 336 km-long slab tip (Figure 2b). The subducting lithosphere comprises a composite plate of constant initial thickness with a core isoviscous layer embedded in upper and lower visco-plastic layers with viscosities that follow a von Mises law. Upper and lower visco-plastic layers are used to approximate the strain-rate weakening that occurs above and below the slab core in thermo-mechanical simulations of subduction (e.g., Garel et al., 2014), following OzBench et al. (2008. Upper and lower layers are assigned the minimum viscosity between the Newtonian viscosity μ Newt and an effective von Mises viscosity μ vM , such that purely viscous deformation occurs as long as the second invariant of the stress tensor II = 2 II (where II is the second invariant of strain rate tensor) does not reach the critical yield stress, τ yield . The effective viscosity of the visco-plastic layers is given by: At material interfaces, the average viscosity is calculated through a geometric mean, where μ i is the viscosity of material i, and Γ i is the relative volume fraction of material i in the vicinity of the finite-element node at which the effective viscosity μ ave is needed.
The subducting plate is surrounded by mantle material, with no overriding or trailing plate. When the plate advances, the mantle material fills in behind the trailing edge. A dome-shaped side plate covers the entire domain adjacent to the subducting plate. It has the same thickness as the subducting plate, and is placed 22 km away from the plate's edge, keeping a constant distance from the symmetry plane around the base of the dome. The side plate is 1,000 times more viscous than adjacent upper mantle material, and prevents lateral flow from narrowing the width of downgoing plate (as in Holt et al., 2017). The lower mantle below 660 km depth is 50 times more viscous than the upper mantle. The viscosity jump is a simplified parameterization of the transition zone's hindering effect of the mantle flow from upper mantle into lower mantle, that in reality arises from not only the viscosity increase with depth but also the endothermic phase transition that is excluded from the model, as in Stegman, Farrington, et al. (2010). Model parameters common to all simulations are listed in Table 1.

Cases Examined and Quantitative Model Diagnostics
We investigate 12 cases across a wide parameter-space. We systematically varied plate width (w) and plate "age," represented by co-varying plate thickness (h), core thickness (h c ), and density contrast between the plate and adjacent mantle (Δρ) to capture the joint increase of negative plate buoyancy and bending strength with age. This allows us to examine how these two factors influence the evolution of subduction and slab morphologies. Our choices are motivated by subduction regime diagrams, as a function of plate age or plate width, from other studies (e.g., Garel et al., 2014;Goes et al., 2017;Stegman, Schellart, & Freeman, 2010). The selected combinations of plate thickness and density contrast produce diverse subduction behaviors, ranging from a vertical-folding type young plate to a retreating and flattening old plate. We note that in the following sections, specified plate widths refer to the full width of the plate (i.e., twice the width simulated, owing to the symmetry plane). Case names, alongside their key parameter values, are listed in Table 2.  To quantify the sensitivity of results to these parameters, we have calculated several diagnostic outputs. When doing so, the boundary of the slab is defined as the 0.5 contour of the mantle material volume fraction (material volume fraction = 1 when the material is mantle, 0 otherwise). Based on this contour, we extract the slab tip depth, the trench location (measured at 15 km depth), and the trailing edge position, as well as rates of slab descent, trench retreat and plate advance. We extract the fastest vertical sinking velocity at the symmetry plane; as well as the maximum sinking velocity throughout the entire domain. We calculate the average slab dip in the upper mantle from the surface to 650 km depth, with respect to the direction of gravity at the slab center at 325 km depth, which is radially toward the center of the sphere from the point of measurement. Measurements are taken at the symmetry plane unless otherwise specified. We also trace the evolution of trench geometry relative to the initial trench shape.

Reference Case
Case W1200_ref is selected as our reference, given its mid-range plate buoyancy and thickness (McKenzie et al., 2005), and width that sits toward the lower end of trench lengths on Earth. The temporal evolution of this case is illustrated in Figure 3a. As subduction starts, the slab tip steepens. During the upper mantle sinking phase (Figure 4a), the trench steadily retreats from its initial position (Figure 4b) with ∼50% of subduction accommodated via this trench retreat (over 60% in the early stages), despite the trailing edge of the plate advancing steadily (Figures 4c and 4d). As the trench retreats, it develops a concave "C" shape, as illustrated in Figure 5b. Following interaction with the viscosity jump at 660 km depth, the slab tip is deflected, the slab sinking rate reduces substantially (Figures 4a and 4e), and the upper mantle section of the slab steepens ( Figure 4f). The slab then slowly sinks into the lower mantle.
Coupling between the sinking plate and the adjacent mantle drives toroidal and poloidal mantle flow (e.g., Funiciello et al., 2006;Schellart, 2004;Stegman et al., 2006). Figures 6a-6c illustrates horizontal flow at 300 km depth at different stages of subduction: trench/slab retreat drives toroidal flow around the edge of the plate, which promotes increased trench concavity (Figure 5b). Figures 6d-6f shows vertical cross-sections through the symmetry plane: two poloidal cells can be identified as the slab sinks in the upper mantle, one above the downgoing plate in the mantle wedge, and the other in the sub-slab region. During the upper-mantle phase of subduction, the mantle wedge cell is more prominent, while flow velocities in this cell diminish as the slab tip deflects and sinks into the more viscous lower mantle.

Influence of Subducting Plate Age
The two cases, W1200_young and W1200_old, were designed to demonstrate how plate age (through its joint effect on buoyancy and strength) modifies subduction dynamics. Our parameter values approximate the properties of younger (decreased Δρ and H) and older (increased Δρ and H) plates, respectively.
The younger slab (W1200_young) stretches and sinks almost vertically through the upper mantle as it subducts, folding upon interaction with the upper-lower mantle interface (Figures 3b and 4e). Trench location remains almost fixed ( Figure 4b) and trench shape does not evolve much over time (Figure 5a). Excluding the initial phase of subduction, trench retreat is minimal: within 5 Myr of the start of subduction, ∼80% of subduction is accommodated by plate advance (Figure 4d).
The older case (W1200_old) exhibits the fastest sinking, trench retreat and plate advance velocities among all cases examined at this width ( Figure 4). The slab tip sinks in the upper mantle at a shallower angle than the younger cases ( Figure 4f). It is deflected at the mantle transition zone, and the sinking rate decreases as the slab moves into the lower mantle (Figures 4a and 4e). Similar to the reference case, after reaching 660 km depth, the upper mantle portion of the slab gradually steepens (Figure 3c). Trench retreat is substantial and accommodates much of the subduction (the trench retreat to total descent ratio remains above ∼55% throughout the simulation- Figures 4b and 4d), with the trench developing a concave curvature over time (Figure 5c).
The cases examined at 1,200 km width clearly display a range of behaviors, with a strong sensitivity to the age, that is, buoyancy and strength, of the subducting slab. The younger plate exhibits the weakest behavior, manifested by a steeper upper mantle subduction angle and minimal trench retreat, with subduction principally accommodated via plate advance. This case falls into the vertical folding regime (e.g., Garel et al., 2014;Goes et al., 2017;Schellart, 2008;Stegman, Farrington, et al., 2010). The older plate is the strongest: it sinks faster,   10.1029/2022GC010597 9 of 21 has a shallower upper mantle subduction angle, and drives significant trench retreat, with more than half of the subduction accommodated via this retreat (Figure 4). This case falls into the weak retreat regime (e.g., Garel et al., 2014;Goes et al., 2017;Schellart, 2008;Stegman, Farrington, et al., 2010). As expected, the reference case has an intermediate strength, with sinking and trench-retreat rates, in addition to the slab dip angle, all between those of the older and younger cases ( Figure 4). As trench retreat accounts for slightly more than ∼50% of the total subduction in the reference case, it also falls into the weak retreat regime.

Effect of Subducting Plate Width
We next examine cases with the same buoyancy and strength combinations as in the previous section, but at larger widths of 2,400, 3,600, and 4,800 km.

Trench Retreat
For cases that share a common plate age, a larger plate width reduces trench retreat. As the slab tries to maintain its sinking rate, this results in stronger bending at the trench. The dynamical behavior can lead to a shift in subduction regime, especially at the center of the plate, where increased slab width causes slabs to steepen at the trench, with the trench sometimes advancing. The behavior at the center of the plate thereby shifts toward a "bending mode," where slab bending at the trench takes up a significant part of the potential energy of the slab, as opposed to a "sinking mode," where bending at the trench uses only 10%-20% of the potential energy, and slab sinking is, in part, achieved through trench retreat (e.g., Capitanio et al., 2007;Ribe, 2010).
For younger cases, at all widths (W2400_young, W3600_young and W4800_young), slabs stretch and sink steeply in the upper mantle, at a dip of >75° (Figure 4f), eventually buckling upon interaction with the transition zone at 660 km depth, like their narrower counterpart (Figure 7). However, as plate width increases, the rate of trench advance also increases. Upon interaction with the upper-lower mantle boundary, the 1,200 km-wide cases display minimal trench motion, whereas the trench has advanced ∼30, ∼50, and ∼50 km for the 2,400, 3,600, and 4,800 km wide cases, respectively (Figure 4b). This folding, with some advance, is a characteristic of a "fold-and-retreat" bending mode (e.g., Goes et al., 2017;Schellart, 2008;Stegman, Farrington, et al., 2010), and the center of wide young slabs display behavior between a vertical folding and fold-and-retreat mode.
For wider cases at the reference age (W2400_ref, W3600_ref, and W4800_ref), slabs retreat prior to interacting with the base of the upper mantle. At the symmetry plane, they steepen and buckle following interaction, thus exhibiting stronger bending at the trench in comparison to the 1,200 km-wide case, which displayed a deflect-andsink behavior ( Figure 3a). As plate width increases, the upper mantle portion of the slab steepens (Figures 4f  and 7). Buckled slabs with a width of 4,800 km have maximum dips that exceed those of the 2,400 km wide case by ∼4°. At the symmetry plane, the trench retreat:total subduction ratio decreases with plate width (∼0.5, ∼0.25, ∼0.1, and ∼0 for 1,200, 2,400, 3,600, and 4,800 km wide cases respectively- Figure 4d), indicating less of a role for trench retreat in accommodating subduction. This is most clearly demonstrated for the W4800_ref simulation, which transitions from retreat at a width of 1,200 km, to stagnation at a width of 4,800 km ( Figure 4d).
For older cases at widths of 2,400 km (W2400_old), 3,600 km (W3600_old), and 4,800 km (W4800_old), slabs sink with shallower angles than corresponding reference age cases in the upper mantle (Figures 4f and 7), deflecting at transition zone depths and, subsequently, sinking through into the lower mantle. As plate width increases from 1,200 to 4,800 km, the maximum upper mantle dip angle increases by ∼14°. The trench retreat:total subduction ratio also decreases as slab width increases, from ∼0.6 to ∼0.3 to ∼0.25 to <0.2 for 1,200, 2,400, 3,600, and 4,800 km wide slabs, respectively. While the older cases across all widths fall into the weak retreat regime, the wider cases (W3600_old and W4800_old) exhibit a significant reduction in trench retreat post interaction with the transition zone ( Figure 4b).
Taken together, our results demonstrate that as plate width increases, slabs display less of a tendency to retreat across all three ages examined. As a consequence, they bend more strongly at the trench.

Trench Curvatures and Along-Strike Variations in Morphology
Different trench shapes are observed across the simulations examined, which can be categorized into 3 types: (a) "I"-type, where the trench is reasonably straight (e.g., Figure 5a); (b) "C"-type, where trench retreat is strongest in the center of the slab relative to its edges (e.g., Figure 5b); and (c) "W"-type, where trench retreat is low in the center of the slab and at the edges, and higher in between ("S" curvature in half-width, as shown in Figure 11b).
We find that "I"-type trenches develop for younger cases across all plate widths: trenches remain reasonably straight, aside from a slight curvature adjacent to the slab edge (Figures 5a, 9a, 10a, and 11a). "C"-type trenches develop for narrow plates that are retreating, for example, in cases W1200_ref and W1200_old (Figures 5b  and 5c). For stronger plates that have moderate width, such as case W2400_old (Figure 9b), the trench develops a gentle curvature close to the edge, but the bulk of the trench remains approximately straight throughout the simulation, in an elongated "C" shape. As slab width increases, "W"-type trenches develop in slabs that would have "C"-type trenches at a narrower width. This can be seen by comparing cases W2400_ref, W3600_ref, and W4800_ref. Case W2400_ref develops a concave curvature at the edges, with the center of the trench retreating slightly less than the edge (Figure 9b), placing it at the transition between "C"-and "W"-type trenches. Conversely, case W3600_ref and case W4800_ref display a "W"-type curvature (Figures 10b and 11b). Similarly, for older slabs, wider trenches develop into a "W" shape ("S" in half-width in Figure 11c). In case W4800_old, the curvature increases following slab transition-zone interaction and the difference in trench retreat between the center and the region of most retreat also increases ( Figure 8b).
Taken together, our results demonstrate that the evolution of trench shape is dependent on both slab age and slab width. Younger and weaker slabs that are in the vertical folding regime develop mostly straight "I"-type trenches, regardless of slab width. For older cases that can drive trench retreat, trench curvatures evolve from a "C"-shape in narrower plates to a "W"-shape in wider plates, with slabs of greater strength transitioning to a "W" shape at a greater width.
Slab morphologies evolve with trench shape. For weaker cases with an "I"-type trench, subducting slab morphology is relatively uniform along strike (Figure 12a). For stronger wide retreating cases that develop a "W"-type trench, along-strike variations in trench retreat translate into morphological variations at depth (Figures 12b, 12c, and 8a): at the symmetry plane, the slab is steep and buckles at the transition zone, with dips up to 9° larger than in the slab at the location of most retreat, where the slab deflects at transition-zone depths (Figures 8d and 8f).

Role of Subducting Plate Age and Width
Our results demonstrate that the evolution of subduction systems is strongly sensitive to slab age (due to its effect on buoyancy and strength), which is consistent with several previous studies (e.g., Capitanio et al., 2007;Garel et al., 2014;Schellart, 2008;Stegman, Farrington, et al., 2010). Slabs with greater negative buoyancy and slab pull promote increased upper mantle sinking velocities (e.g., Garel et al., 2014;Goes et al., 2017;Stegman, Farrington, et al., 2010). Slab thickness determines slab strength (alongside buoyancy), with thicker slabs possessing a higher bending resistance (e.g., Bellahsen et al., 2005;Capitanio & Morra, 2012;Conrad & Hager, 1999;Ribe, 2010) and, accordingly, taking longer to bend at the trench. The regime that a subduction system falls into depends on a delicate balance between the amount of time taken to bend (larger for thicker slabs) and the sinking time (larger for younger slabs) (Ribe, 2010). Taken together, in younger slabs (low slab pull-longer sinking times; low slab strength-shorter bending times), bending dominates over trench retreat, with slabs typically subducting steeply and buckling upon interaction with the upper-lower mantle boundary, owing to the high angle of incidence. Conversely, in older slabs (high slab pull-shorter sinking times; high slab strength-longer bending times), there is insufficient time for substantial bending at the trench, with subduction accommodated principally through trench retreat. As a result, slabs typically exhibit a shallower upper mantle dip angle, which facilitates slab deflection and stagnation at the base of the upper mantle (e.g., Agrusta et al., 2017;Čížková & Bina, 2013;Garel et al., 2014;Torii & Yoshioka, 2007).
The evolution of "C"-and "W"-shaped trenches for our retreating cases are similar to results from Schellart et al. (2007), with curvature at slab edges induced by toroidal flow around and into the slab. Interplay between the size and strength of the toroidal cell, the width of the slab, and the slab's resistance to along-strike bending, dictate how trench shape responds. The size of the toroidal cell determines which segments of the trench experience significant force from adjacent mantle flow and, hence, the location of potential concave curvature development. The strength of the toroidal cell is determined by the amount of trench retreat. The width of the plate relative to the size of the toroidal cell determines the distance between the toroidal cells at both edges. When these factors are coupled with the plate's resistance to bending, the evolution of trench shape can be understood. "C"-shaped trenches are observed for narrow plates, where toroidal cell sizes are close to half of the slab width ( Figure 6). "W"-shaped trenches are observed for wider plates, where the size of the toroidal cell is substantially smaller than the width of the plate: the center of such plates are thus not markedly influenced by toroidal flow (Figure 13a). Plates with less resistance to bending along strike can develop enhanced curvature at the trench (W3600_ref, Figure 10b); conversely, plates with a higher strength can prevent significant curvature development at the trench, remaining in a relatively uniform elongated "C"-shape rather than evolving into a "W"-shape (e.g., W3600_old, Figure 10c).
While the influence of plate width on subduction dynamics has been studied previously (e.g., Di Schellart et al., 2007;Stegman et al., 2006), our results demonstrate that the important role of width is strongly modulated by the age of the plate and its effect on plate strength. The change from a "C"-shaped trench to a "W"-shaped trench with increasing width only occurs for cases that are initially in a retreating regime (i.e., older plates). For younger plates that are in the vertical folding regime, increasing plate width has little impact on along-strike variability, because these plates do not subduct through trench retreat and, hence, cannot generate toroidal cells of the intensity required to induce trench deformation (Figure 13b). Accordingly, the younger cases develop "I"-type trench shapes across all widths examined in this study.
Variations in the amount of trench retreat also translate into along-strike morphological variations at depth. These variations can also be understood by the along-strike influence of toroidal flow, facilitating trench retreat. Vertical folding morphologies correspond to (parts of the) slabs with limited trench retreat or minor advance, while retreating segments lead to shallower dips and deflection at the base of the upper mantle. While all wide slabs display the typical morphology of the vertical folding regime at the center, the young models have tight buckles, whereas older slabs have open folds with larger bending radii. This difference in bending radii illustrates that older slabs have higher bending resistance and strength than younger slabs, despite falling into the same subduction regime. Overall, as plate width increases, the center of the slab shifts from sinking to bending, due to the lack of toroidal flow and its role in driving trench retreat. The competing roles of plate age and plate width in dictating the subduction style are summarized via a regime diagram in Figure 14.

Implications for Subduction on Earth
Our model predictions across different plate ages and widths, can be used to analyze how these aspects combine to control the spatio-temporal evolution of trenches on Earth. While Earth's subduction zones are substantially more complex than those considered in our simulations, due to a multitude of factors including subducting plates of non-uniform age, the subduction of buoyant anomalies, and the influence of overriding plates, the "I," "C," and "W," trench shapes predicted by our models are consistent with present-day trench shapes (e.g., Bird, 2003;Heuret et al., 2011;Schellart et al., 2007) and those in reconstructions of plate motion histories through the Cenozoic Era (Müller et al., 2016(Müller et al., , 2019.

"I"-Type Trenches
Our results demonstrate that "I"-shape trenches typically develop when a young plate subducts with negligible trench motion, regardless of trench width. The tectonic reconstructions of Müller et al. (2019) provide some examples of young plate subduction during the Cenozoic, such as the subduction of the Izanagi plate at the Japan trench at 50-60 Ma and the subduction of the Farallon plate below North America prior to ∼30 Ma ( Figure 15).
As illustrated in Figure 16, the Japan subduction zone was relatively straight (apart from its northern end) with minimal trench motion between 50 and 60 Ma, characteristics typical of "I"-type trenches. The trench measured ∼5,000 km in width, and a young plate (∼10 Myr) was subducted along the whole trench ( Figure 15c). As time advanced past 50 Ma, the main part of the trench evolved from an "I"-type toward a "C"-type example, with increasing trench retreat, trench curvature and trench segmentation (Figure 16a). This is coincident with an increase in subducting plate age, as shown in Figures 15a and 15b, after the subduction of the Izanagi-Pacific ridge at ∼50 Ma (Kimura et al., 2019).
Similarly, the reconstructed North America trench displayed an "I"-type shape prior to 30 Ma, when the very young (∼10 Myr at 30 Ma) Farallon plate was subducting beneath North America, as shown in Figure 15b. Prior to 30 Ma, the trench shape was relatively straight, with very little trench retreat, particularly from 30-50 Ma ( Figure 16). Following breakup of the Farallon Plate in the mid-Cenozoic, into the Juan de Fuca, Cocos and Nazca Plates (e.g., Atwater, 1970;Lonsdale, 2005), the continuity of the 5,900 km-wide trench was lost, and the strike-slip San Andreas Fault developed on the west coast of North America.

"C"-Type Trenches
Our results suggest that "C"-shape trenches should be associated with moderate to old subducting plate ages and moderate slab widths. "C"-shape trenches are the most common trench shape observed on Earth. The Aleutian subduction zone and the Sunda-Java subduction zone are two examples of "C"-shape trenches that developed during the Cenozoic (Figure 16). Although both trenches have also been affected by buoyant structures on the incoming plate (such as the Yakutat terrane below Alaska and Australian continental crust impinging on the Banda part of the Sunda-Java arc), we propose that the combined width and age of the downgoing plate exerted an important control on the evolution of trench shape.
The Aleutian trench extends ∼4,000 km from the south coast of Alaska to Kamchatka (Scholl et al., 1975). Prior to 40 Ma, the Kula plate subducted at the eastern side of the Aleutian trench, and the Pacific plate subducted at the western side. Subduction of young oceanic lithosphere (∼10-40 Myr at 60 Ma) along the trench (Figure 15c) resulted in a gentle curvature, with a shape between "C"-and "I"-types ( Figure 16, Müller et al., 2019). At 40 Ma, Figure 13. Lateral flow patterns at 300 km depth for: (a) case W4800_ref; and (b) case W4800_young. The length and direction of the arrows illustrates the magnitude and direction of tangential velocities (i.e., after the radial component has been removed). In both panels, the largest arrow represents a tangential velocity magnitude of 2.5 cm/yr. For the reference age case in (a), a toroidal cell can be identified at the edge of the slab, which does not extend to the center of the slab. For the younger case in (b), although there is some toroidal flow around the edge of the slab, its magnitude and area of influence is insignificant when compared to the reference case. Figure 14. Schematic diagrams illustrating how slab age and width control subduction styles and trench shape. Regimes: VF-vertical folding; WR-weak retreat (see text for further detail). The three trench shapes, "I," "C," and "W," are separated by gray dashed lines into three approximate domains. Slab behaviors that are in VF at the symmetry plane are represented by the cyan region, and those in WR are represented by the brown region. The rightmost panel illustrates slab morphology at the center of the slab (i.e., the symmetry plane, abbreviated to sym.) and at the location of most trench retreat, which is at the center of the concave curvature (curv.). The 12 cases examined herein are marked by stars on the regime diagram. The young, reference and old plates have estimated ages of 10, 70, and 140 Myr respectively. For young plates, the subduction regime is VF regardless of the width of the plate, and trench shapes are mostly straight, indicated by "I"-type. As plate age increases, it is easier to drive trench retreat and slabs fall into the WR regime; but as width increases, the center of the slab shifts toward the VF regime. Beyond a certain age, the narrower and/or older plates tend to develop "C"-type trenches; wider and/or younger plates tend to develop "W"-type trenches.
when the Pacific plate had started subducting at the Aleutian trench; the age of the subducting plate increased and the Aleutian trench retreated and developed a "C"-shape curvature, with enhanced curvature in the west ( Figure 16). This is consistent with our modeling predictions, and could be related to the non-uniform subducting plate age at the Aleutian trench: the subducting Pacific plate is younger to the east (currently ∼10 Myr) and older to the west (∼120 Myr at present, Figures 15a and 15b), with the older part of the plate driving more retreat, generating the asymmetric "C"-shaped trench.
The Sunda-Java trench also has a complex subduction history. Prior to 43 Ma, the active Wharton ridge was subducting beneath Sumatra, but since then the ridge has become inactive (Whittaker et al., 2007). As a result, the majority of material subducted prior to 43 Ma was young (∼10 Myr old- Figure 15c), leading to minimal trench motion at the Sunda-Java subduction zone, consistent with "I"-type subduction ( Figure 16b). As the Wharton ridge ceased spreading and the subducting plate age increased (Figures 15a and 15b), the trench began to retreat and developed a "C" shape across its ∼5,000 km width, again demonstrating that an "I"-type to "C"-type transition can occur when subducting plate age increases, consistent with our modeling predictions.

"W"-Type Trenches
"W"-type trenches develop with moderate and older subducting plate age and very wide trenches. The South American trench is the textbook example of a "W"-shape (Schellart et al., 2007): it exhibits concave curvature on both edges, with the center of the trench almost stagnant throughout the Cenozoic (Figure 15c). Subduction in the South Pacific also exhibits "W"-type characteristics in the early Cenozoic.
The South American trench is over 6,000 km long, and subducted moderately old material (∼50-80 Myr) throughout the Cenozoic at the center of the trench ( Figure 15). Trench evolution shows increasing oroclinal  bending through the Cenozoic (Schepers et al., 2017) and, hence, more retreat toward the north and south than in the central part at the Bolivian bend ( Figure 16, Müller et al., 2019). The present-day trench shape is typical of our wide plate model predictions, where the Bolivian Orocline is located close to the center of the trench, while sections of the trench close the edges have a concave geometry. The subduction of pre-existing buoyant features on the Nazca Plate and interaction with varying thickness upper plate are further complexities that likely added to its evolution toward the current geometry (e.g., Capitanio et al., 2011;Espurt et al., 2008;Gutscher, Olivet, et al., 1999). However, although it is likely that multiple factors contribute to the shape of the trench at the South America Subduction Zone, our results agree with the previous suggestion by Schellart et al. (2007) that the first-order "W"-shape is dictated by its large width, modulated by its moderate subducting-plate age.
The South Pacific region has a relatively complex tectonic history. In the early Cenozoic, old Pacific plate was subducting under the South Pacific trench, which had length exceeding 6,000 km ( Figure 15c). The trench at 60 Ma has characteristics of a "W"-shape. There is a region of low trench retreat near 20°S, where the oldest (∼100 Myr) part of the plate was subducting. The trench is concave north and south of this point. As the trench evolves, the northern and southern ends of it start to evolve separately, and the shorter trench takes on a "C" type shape (Figures 15c and 16). Post 20 Ma, in response to arrival of the buoyant Ontong-Java plateau (Mann & Taira, 2004;Neal et al., 1997;Stotz et al., 2017) the trench segments further into the New Hebrides, New Britain and Tonga-Kermadec-Hikurangi trenches (e.g., Pelletier et al., 1998). These different stages of evolution and the shape of the resulting trenches were likely preconditioned by the original "W" shape.
Overall, the examples of "I"-, "C"-, and "W"-shape trenches on Earth are in line with our modeling results. "I"-type trenches are associated with very young downgoing plates of ∼10 Myr old; and as plate age increases, some transition into "C"-shape trenches. "W"-shape trenches are observed in subduction zones exceeding 6,000 km width, where older material (greater than 50 Myr old) is being subducted, thus driving trench retreat. There is no doubt that the trench shape at each subduction zone is further modulated by additional complexities, including variable downgoing plate age along strike, subduction of buoyant active or bathymetric ridges, variations in thickness and buoyancy of the upper plate, and large tectonic reconfiguration events such as the India-Eurasian collision. Nonetheless, our results demonstrate the key role that both subducting plate age and width play in controlling the evolution of trench geometry, providing a framework to better understand the evolution of subduction zones.

Conclusions
We have presented new 3-D spherical free-subduction models with a composite visco-plastic plate and viscously layered mantle. We examined the sensitivity of subduction dynamics and trench evolution to plate age (simulated with covarying plate densities and thicknesses) and plate width.
Our results complement previous studies on the effect of age and width on the evolution of subduction zones. As plate age increases, plate strength increases and, as a result, the subduction style transitions from vertically sinking and folding to retreating with a shallower upper mantle dip angle. Our models show a tendency to produce "C" shaped trenches for narrower plates and "W" shaped trenches for wider plates, consistent with the models of Schellart et al. (2007). However, we also find that the effect of width is modulated by the age of the subducting plate. For young plates that are in the vertical folding regime, the trench does not develop a "W" shape, regardless of width. "C" or "W"-shaped trenches are only generated in cases where the slab is able to drive trench retreat. For those cases that do retreat, younger and weaker plates develop enhanced along-strike variability in trench shape.
We now have the means to simulate subduction dynamics in a spherical shell geometry, which is appropriate for Earth's mantle. This opens up new possibilities and will be used in the future to investigate additional factors that affect subduction dynamics and their expression at Earth's surface.

Data Availability Statement
The Fluidity computational modeling framework, including source code and documentation, is available from https://fluidityproject.github.io/; the version used for the simulations presented herein has been archived at . The input files required to reproduce the simulations presented herein have also been made available at Chen (2022). Figures have been prepared using Matplotlib, Cartopy, and Paraview.