Control of Subduction Zone Age and Size on Flat Slab Subduction

Flat slab subduction is an enigmatic style of subduction where the slab attains a horizontal orientation for up to several hundred kilometers below the base of the overriding plate. It has been linked to the subduction of buoyant aseismic ridges or plateaus, but the spatial correlation is problematic, as there are subducting aseismic ridges and plateaus that do not produce a flat slab, most notably in the Western Pacific, and there are flat slabs without an aseismic ridge or plateau. In this paper an alternative hypothesis is investigated which poses that flat slab subduction is associated with subduction zones that are both old (active for a long time) and wide (large trench-parallel extent). A global subduction zone compilation is presented showing that flat slabs preferentially occur at old (>∼80–100 Myr) and wide (≥∼6000 km) subduction zones. This is explained by the tendency for wide subduction zones to decrease their dip angle in the uppermost mantle with progressive time, especially in the center. A set of numerical subduction models confirms this behavior, showing that only the central parts of wide slabs progressively reduce their slab dip, such that slab flattening, and ultimately flat slab subduction, can occur. The models further show that a progressive decrease in slab dip angle for wide slabs leads to increased vertical deviatoric tensional stresses at the top surface of the slab (mantle wedge suction), facilitating flat slab subduction, while narrow slabs retain steep dip angles and low vertical deviatoric tensional stresses. The results provide a potential explanation why present-day flat slabs only occur in the Eastern Pacific, as only here subduction zones were old and wide enough to initiate flat slab subduction, and why Laramide flat slab subduction and South China flat slab subduction were possible in the geological past.


INTRODUCTION
The cross-sectional geometry of upper mantle slabs varies considerably on Earth, with variations in slab dip angle, bending curvature and potential deflection of the slab in the mantle transition zone (e.g., Jarrard, 1986;Lallemand et al., 2005). For most active subduction zones on Earth, the uppermost ∼200 km of the slab is defined by one convex-upward slab hinge located close to the trench (Figures 1A,B) (e.g., Kuril-Kamchatka, Izu-Bonin-Mariana, Sunda, Tonga-Kermadec-Hikurangi, New Hebrides). Some subduction segments, however, show two or three slab hinges in the uppermost 200 km. The former generally has a very gentle, convex upward, slab hinge near the trench and a second, more pronounced, convex upward hinge several hundred kilometers downdip, with a very low angle slab segment in between ( Figure 1C) (e.g., Alaska, Nankai, Cascadia). The latter, with three slab hinges, has one convex-upward hinge near the trench, one concave upward hinge that marks the start of a flat slab segment dipping ≤10 • , and one convex-upward hinge that marks the end of the flat slab segment ( Figure 1D). It is this subduction geometry, with three slab hinges that is most enigmatic and that is the subject of this study. In this contribution, only the subduction geometry with three slab hinges will be referred to as flat slab subduction. There are only a few current cases of such flat slab subduction, including central Chile, southern Peru, central Peru and Mexico, as has been discussed recently (Manea et al., 2017) and as can be seen in global slab models (e.g., Hayes et al., 2012). Several examples have also been proposed for the geological past, such as western North America, Central Andes and South China (Henderson et al., 1984;Li and Li, 2007;Ramos and Folguera, 2009).
The origin of flat slab subduction remains enigmatic. Previous work proposed that flat slabs occur through subduction of buoyant ridges or plateaus such as the Juan Fernandez Ridge and Nazca Ridge at the South American subduction zone (e.g., Pilger, 1981;Gutscher et al., 2000;van Hunen et al., 2002). But the spatial correlation between ridge/plateau subduction and flat slab subduction has many exceptions, with regions of aseismic ridge subduction lacking flat slab subduction (e.g., Tonga with Louisville Ridge, New Hebrides with d'Entrecasteaux Ridge, Mariana with Marcus-Necker Ridge, Kamchatka with Emperor Ridge) and some flat slab subduction segments lacking an aseismic ridge/plateau (e.g., Mexico) Clayton, 2011, 2013;Manea et al., 2017). Others have proposed that flat slab subduction might result from forced trench retreat (e.g., van Hunen et al., 2004;Schepers et al., 2017), strong suction forces in the mantle wedge (e.g., Tovish et al., 1978), or slab-plume interaction (e.g., Betts et al., 2009). The latter mechanism has not generally been applied to Cenozoic examples of flat slab subduction, and recent geodynamic models of slabplume interaction for present-day Earth-like settings indicate that plumes generally do not affect slab geometry as their upward buoyancy flux can be more than two orders of magnitude smaller than the downward slab buoyancy flux (Mériaux et al., 2016). Wedge suction forces are generally considered to play a role in enhancing flat slab subduction (Tovish et al., 1978;van Hunen et al., 2004), but the question is, what might enhance wedge suction? A number of recent modeling works indicate that a relatively thick or far-field cratonic overriding plate will enhance wedge suction and thus facilitate flat slab subduction (Manea et al., 2012;O'Driscoll et al., 2012;Taramón et al., 2015). Forced trench retreat can also reproduce flat slab subduction in geodynamic models (van Hunen et al., 2004), but this forcing is generally externally imposed (mostly as a velocity boundary condition on the overriding plate) and so the question remains what the source of this forcing is and if sufficient forcing can arise in a buoyancy-driven geodynamic environment in nature (and model).
A number of relatively recent works argue that flat slab subduction requires a combination of physical ingredients. Proposed combinations include forced trench retreat and enhanced suction due to a thick cratonic overriding plate (Manea et al., 2012), buoyant plateau subduction with forced trench The southern Kuril slab segment as an example of normal subduction with one convex-upward slab hinge and a normal slab dip angle [geometry derived from Slab1.0 model (Hayes et al., 2012)]. (B) The southern New Hebrides slab segment as an example of normal subduction with one convex-upward slab hinge and a steep slab dip angle [geometry from Slab1.0 model (Hayes et al., 2012)]. (C) The Alaska slab segment as an example of unusual subduction with two convex-upward slab hinges and a low-dip-angle slab segment in between [geometry derived from Ohta et al. (2006)]. (D) The Central Peru slab segment as an example of unusual subduction with two convex-upward slab hinges (trench hinge and convex flat slab hinge) and one concave upward slab hinge (concave flat slab hinge) in between [geometry derived from Phillips and Clayton (2014)]. A (sub)horizontal slab segment is located between the concave and convex flat slab hinges. It is only this slab geometry with three hinges that is referred to as flat slab subduction in this work and that is the focus of this study. Note that all images show the geometry of the top of the slab. and plate motion (Liu and Currie, 2016), buoyant ridge/plateau subduction, enhanced suction, young oceanic plate subduction, and forced trench and plate motion , or ridge/plateau subduction, suction and forced trench retreat (van Hunen et al., 2004;Antonijevic et al., 2015).
Here I build on these earlier findings and investigate the hypothesis that flat slab subduction (with three slab hinges as in Figure 1D) preferentially occurs at old subduction zones (those that have been active for a long time) that are very wide (large trench-parallel extent), because long-term subduction facilitates lower mantle slab penetration, which -through deepmantle return flow -forces trenchward overriding plate motion, overriding plate thickening, upper mantle slab dip reduction and trench retreat (Schellart, 2017), with trench retreat being slow due to the wide-slab setting (Schellart et al., 2007), while wide-slab subduction enhances wedge suction, especially in the center (Dvorkin et al., 1993). A compilation from geological data of subduction zones in nature is presented, as well as results from buoyancy-driven geodynamic subduction models, that both provide support for the hypothesis that subduction zone width and age play an important role in flat slab subduction.

Data Sources for Subduction Zones in Nature
A major compilation of subduction zone age and size of the active subduction zones on Earth has been constructed, either with or without flat slab subduction, as well as a number of subduction settings from the geological past for which flat slab subduction has been reported (Figure 2). The data sources that have been used are listed in Tables 1, 2, and consist of previous works on the geological and tectonic evolution of subduction zones. Table 1 lists those subduction zones that have a presentday flat slab or a flat slab in the geological past as defined in Figure 1D in the paper (i.e., a subduction zone with a trench hinge and two flat slab hinges). Table 2 lists those subduction zones that have no reported flat slab and subduction segments of relatively large subduction zones for which that particular segment itself currently has no flat slab. For both tables, slab widths are listed. For Table 1 these are the slab widths during (or close to) subduction zone formation and the slab widths at the time when flat slab subduction commenced. For Table 2 these are the slab widths during (or close to) subduction zone formation and the present-day slab widths. All the present-day slab widths have been derived from Schellart et al. (2007). For slab width values during the geological past, a variety of published plate tectonic reconstructions have been used. Both tables also list the age of the subduction zone (the time of subduction initiation), while Table 1 also lists the time of flat slab inception and the subduction zone age at the time of flat slab inception. Details and justifications for data selection and uncertainties associated with the data are listed in the Appendix.

Numerical Modeling Method
The hypothesis that flat slab subduction occurs for old and wide subduction zones is tested with geodynamic models of long-lived, buoyancy-driven, progressive subduction in which slab width is varied. The models use the Underworld numerical modeling code (Stegman et al., 2006;Moresi et al., 2007) and the model approach is the same as discussed in Schellart (2017). The models are run non-dimensionally and are later scaled to natural values. The reader is referred to the earlier works for FIGURE 2 | (A) Diagram showing the dependence of flat slab subduction in nature on slab width (trench-parallel extent) and subduction zone age. Diamonds represent subduction zones/segments without a flat slab. Filled circles and white circles with an outline represent subduction zones with a present-day flat slab and flat slab in the geological past, respectively. For the black circles and white circles with black outline subduction zone age and slab width represent the age and width at the time of flat slab inception. For the black diamonds subduction zone age and slab width represents their age and width at present. For the gray circles, white circles with a gray outline and gray diamonds slab width represents the width at subduction initiation. See Appendix and Tables 1, 2 for data choices, justification and uncertainties. Black bars and gray bars are not error bars but indicate range of age estimates and slab width estimates, where data point location represents the average or the best estimate (see Appendix and Tables 1, 2). Thick horizontal gray lines connect subduction initiation data point with present-day subduction zone data point (diamonds) or flat slab subduction initiation data point (circles) of the same subduction zone (segment). Subduction zone data points: 1-South America with 1a-Southern Peru segment (Nazca ridge), 1b-Central Peru segment, 1c-Central Chile segment (Juan Fernandez ridge), 1d-Colombia segment, 1e-Northern Chile segment, 1f-Southern Chile segment; 2-Mexico-Central America with 2a-Mexico segment, 2b-Central America segment; 3-South America (Central Andes segment); 4-South China; 5-Farallon (Laramide segment); 6-Calabria; 7-Gibraltar; 8-Hellenic; 9-Scotia; 10-Manila; 11-Sangihe; 12-Puysegur; 13-Halmahera; 14-North Sulawesi; 15-Cascadia; 16-Nankai-Ryukyu; 17-Lesser Antilles-Puerto Rico; 18-Tonga-Kermadec-Hikurangi; 19-Aleutians-Alaska; 20-Melanesia (New Britain-San Cristobal-New Hebrides); 21-Northwest Pacific (Kamchatka-Kuril-Japan-Izu-Bonin-Mariana) with 21a-Kamchatka-Kuril-Japan segment, 21b-Izu-Bonin-Mariana segment; 22-Sunda (Burma-Andaman-Sumatra-Java-Banda) with 22a-Burma-Andaman segment, 22b-Burma-Andaman-Sumatra-Java segment. The light gray zone indicates the approximate location of the boundary separating the domains where flat slab subduction is not possible (lower left) and possible (upper right). Note that for the individual subduction zone segments (1a-f and 2a-b) the entire width of the slab that these segments form part of has been plotted.  100-120 * Slab width is based on the average of the estimated minimum and maximum slab width at the time of inception of flat slab subduction. It is derived from tectonic reconstructions presented in Gordon and Jurdy (1986);Collins (2003), Schellart et al. (2007Schellart et al. ( , 2010, Domeier and Torsvik (2014) and Seton et al. (2012). ∧ Trench-parallel distance from flat slab segment to closest lateral slab edge at the time of flat slab inception. † Subduction zone age (time of subduction zone initiation) based on the best estimate (be) or average (av) of the estimated minimum and maximum age. ‡ Age of the subduction zone at the time of inception of flat slab subduction. § Southern Peru flat slab spatially correlated with the subducting Nazca Ridge (Pilger, 1981;Gutscher et al., 2000). @ Central Peru flat slab possibly related to subduction of the Inca plateau (Gutscher et al., 2000). ¶ Central Chile flat slab spatially correlated with the subducting Juan Fernandez Ridge (Pilger, 1981;Gutscher et al., 2000). # Mexican flat slab, not spatially associated with any ridge or plateau (Skinner and Clayton, 2011). Note that a flat slab subduction initiation age of 30-25 Ma is chosen, based on the timing of migration of arc magmatism from the Sierra Madre del Sur to the Trans-Mexican Volcanic Belt (Ferrari et al., 1999;Morán-Zenteno et al., 1999;Kim et al., 2010), and not a younger initiation age of 16-9 Ma proposed more recently (Manea et al., 2017). A consequence of adopting the older age is that the slab width then includes the South American segment and Mexico-Central America segment (and possibly the Baja California segment), because the Nazca and Cocos plates still formed one single plate until ∼23 Ma (Lonsdale, 2005). & Central Andes (southern Peru-northernmost Chile) flat slab in the latest Eocene-Oligocene (Ramos and Folguera, 2009). + South China flat slab in the Middle Triassic-Early Jurassic (Li and Li, 2007). % Farallon flat slab in western North America, thought to be responsible for the Laramide orogeny and eastward migration of the magmatic arc (Dickinson and Snyder, 1978;Henderson et al., 1984). The numbers in between the brackets refer to the following references: 1 Ferrari et al. (1999), 16 Ramos and Folguera (2009), 17 Domeier and Torsvik (2014), 18 Collins (2003), 19 Li and Li (2007), 20 Burchfiel andDavis (1975), 21 DeCelles (2004), 22 Henderson et al. (1984), 23 Liu et al. (2010), 24 Copeland et al. (2017).

Mexico-Central America
Central America 2b 3100 5800 (10) /11,000 ± 1000 (18) 220 (be) 174-201 (37) 220 (38) The numbers in between the brackets refer to the following references: 1 van Hinsbergen et al. (2014), 2 Séranne (1999), 3 Rosenbaum et al. (2002), 4 Lonergan and White (1997), 5 Wortel and Spakman (2000), 6 Dilek and Sandvol (2009) details on the numerical method and approach. The model setup involves buoyancy-driven subduction of a higher-density, viscously stratified, subducting plate below an overriding plate into a lower density, lower viscosity, stratified mantle that extends down to the core-mantle boundary (Figure 3). Each model has free-slip conditions along all its boundaries and starts with a short, 206-km-long slab perturbation dipping at 29 • that triggers buoyancy-driven subduction. The model domain is 10,000 km long and 2900 km deep. Four models are presented that are all exactly the same except for their subduction zone width (W) and box width (BW): a narrow slab model (W600 with W = 600 km, BW = 4000 km), an intermediate-width slab model (W2000 with W = 2000 km, BW = 6000 km), a wide slab model (W6000 with W = 6000 km, BW = 12,000 km), and a very wide slab model (W-infinite, a model with a 2D spatial set-up, so W = infinite, BW = infinite). The models do not contain lateral side plates, assuming that the transform faults along the sides of the plates are very weak. Furthermore, the trailing edges of the plates are free, mimicking spreading ridges, offering minimal resistance to lateral motion, following earlier works on free subduction (e.g., Kincaid and Olson, 1987;Chen et al., 2016).
The subducting plate has a homogeneous density that is 60 kg/m 3 higher than that of the sub-lithospheric mantle, and is viscously stratified into three layers, including a viscoplastic top layer with a von Mises rheology to allow for decoupling of the subducting plate from the top surface, a high-viscosity Newtonian central layer and a low-viscosity Newtonian bottom layer (Figure 3). The effective viscosity averaged over the thickness of the slab falls within the range 270-640η UM−Max (with η UM−Max being the maximum sub-lithospheric upper mantle viscosity), which is very comparable to viscosity estimates from earlier works, such as Ribe (2010) with an estimated slab/upper mantle viscosity ratio of 140-510, and Schellart (2008) with an estimated ratio of 100-700. The overriding plate has several lateral domains with each domain having a constant viscosity. It further has a vertical density stratification with a 30-km-thick top layer with a density that is 480 kg/m 3 less than that of the sub-lithospheric upper mantle, mimicking the density contrast between continental crust and sub-lithospheric upper mantle, while the density of the lithospheric mantle is 30 kg/m 3 higher, mimicking continental lithospheric mantle with a density that is moderately higher than that of the sublithospheric mantle.
The models include a sub-lithospheric upper mantle domain with a non-linear stress-dependent viscosity η UM down to a depth of 660 km, with a stress exponent n = 3.5 (Mackwell et al., 1990), a minimum viscosity η UM−Min and maximum viscosity η UM−Max , such that η UM−Min = 0.1 η UM−Max . The variation in sub-lithospheric upper mantle viscosity is therefore limited to one order of magnitude to facilitate reasonable convergence rates in the numerical calculations and to ensure that the scaled time of individual time steps does not drop significantly below the average scaled time. The dimensionalized η UM−Min = 5 × 10 19 Pa s and η UM−Max = 5 × 10 20 Pa s. Such values are within the estimated range of values for the sub-lithospheric upper mantle (10 19 -10 21 Pa s) in nature (Artyushkov, 1983;Ranalli, 1995). Note, however, that uncertainty in mantle viscosity values in nature directly affect dimensionalized velocity values in the models.
The 2240 km thick lower mantle domain has a Newtonian viscosity η LM = 100 η UM−Max . As such, a minimum viscosity jump of a factor of 10 2 was applied between the upper and lower mantle. This viscosity step implementation represents all the effects of the 660 km discontinuity (viscosity changes, density changes, phase transitions) and captures the discontinuity's geodynamic essence through reducing the slab velocity in the lower mantle, as implied by earlier studies on mineral physics and phase transitions in the deep mantle (Torii and Yoshioka, 2007;Ganguly et al., 2009). The choice for this simple implementation is further supported by results from recent subduction modeling, which show that, unless a complete treatment of compositional layers and phase transitions is implemented, large-scale deformation of slabs is better approximated by a model with no phase transitions rather than including an incomplete approximation that over-predicts slab folding (Arredondo and Billen, 2016).
The isothermal conditions in the model, and thereby the absence of slab warming, will not affect the slab morphology, slab viscosity and slab-mantle density contrast in the upper mantle significantly due to the relatively rapid rate of subduction (up to 7.5 cm yr −1 ) and the slow rate of conductive slab warming. Lower mantle slab warming is likely more significant, producing a weaker lower mantle slab, which would thereby likely result in stronger slab folding and tighter slab folds. However, the density contrast of the lower mantle folded slab pile is not significantly affected, because it sinks as a whole and includes the entrained mantle material enclosed within the folds. Considering that warming of the folded slab segment coincides with the cooling of the entrained ambient mantle in the folds (Schellart, 2017), the thermal buoyancy contrast between folded slab pile and ambient mantle does not diminish and disappear on a timescale representing the duration of the numerical models. Therefore, the driving mechanism of the lower mantle slab is not significantly affected.
The model resolution in cross-section is 512 (horizontal) × 192 (vertical) elements with 20 particles per element. A mesh refinement has been implemented such that a 3,000 km (length) by 290 km (depth) domain around the subduction zone has a maximum resolution with cells that are 9.8 km long and 7.6 km deep. This spatial refinement is required for properly resolving the subduction zone interface. The model resolution in the trench-parallel direction is 48 (W600), 80 (W2000), and 160 (W6000) elements. The models with a 3D spatial set-up were run on the Australian national supercomputer (Raijin) using 240-496 cores and individual models took between 3 months (W600) and more than 2 years (W6000) to complete, the latter model requiring more than 90 restarts to complete a total of more than 5000 time steps (with one time step scaling to, on average, ∼0.04 Myr, so about 40,000 years). Considering the large amount of time and computational resources required to complete these models, it has not been possible to do an extensive parametric investigation other than investigating the influence of slab width. A detailed parametric study is beyond the scope of this work but is part of future studies using models with a 2D spatial set-up, in which the effect of various parameters and boundary conditions on flat slab subduction are investigated, including the subducting plate thickness, subducting plate viscosity, overriding plate viscosity, viscous stratification of the sub-lithospheric mantle, temperature-dependent viscosity and temperature-dependent density, as well as free-slip, no-slip and periodic boundary conditions.
The numerical models use Cartesian geometry. It could be conceived that the usage of (more realistic) spherical geometry would facilitate the formation of flat slabs, especially for wide slabs, because the reduction in spherical surface area with increasing depth requires a sinking slab to accommodate horizontal shortening strains and compressive stresses during progressive sinking, while the total required horizontal shortening increases linearly with slab width. The slab resistance to such strains and the associated in-plane stress arching would provide an additional lifting force to the slab that would promote slab flattening. However, the slab strain magnitude associated with the subduction of a spherical cap in spherical geometry (in contrast to a flat plate in Cartesian geometry which lacks such strain) is independent of the width of the slab, as it merely depends on the vertical sinking distance (with respect to the surface) and the radius of the sphere. Furthermore, it can be shown that strain rates in nature associated with such strains are relatively low. For example, if one considers a subducting plate that sinks from the surface down to 100 km depth at a vertical sinking rate of 5 cm/yr, this would amount to a trench-parallel shortening strain and strain rate of ∼0.016 and ∼2.5 × 10 −16 s −1 , respectively, due to the decreased radial distance to the Earth's center and the associated reduction in trench-parallel line length. If we compare such a shortening strain rate with strain rates within subducting plates near the trench in nature resulting from plate convergence and bending, which reach values up to ∼10 −14 s −1 (e.g., Kreemer et al., 2003;Stadler et al., 2010), then it is clear that the trench-parallel strain rates are subordinate, being up to 40 times smaller than the convergence/bending related strain rates. Considering these relatively low trench-parallel strain rates in spherical geometry, implying low trench-parallel within-slab compressive stresses, it is here argued that the simplification of using Cartesian geometry does not have a significant effect on flat slab subduction formation and evolution. Figure 2 shows a compilation of the active subduction zones on Earth, either with or without flat slab subduction, and a number of subduction settings from the geological past for which flat slab subduction has been reported. Age since subduction initiation has been plotted against slab width. For those subduction zones without a flat slab (e.g., Scotia, Tonga-Kermadec-Hikurangi, Sunda, Nankai-Ryukyu, Aleutians-Alaska, Hellenic), their current age and width (their trenchparallel extent) have been plotted. Note again that, following the definition outlined in the introduction, Nankai, Alaska and Cascadia are not counted as flat slab settings because their uppermost slab geometry is not defined by three slab hinges (as in Figure 1D with one concave-upward hinge) but only by two convex-upward hinges (as in Figure 1C). For those subduction zones with a flat slab (as in Figure 1D), their age and width have been plotted as it was at the onset of flat slab subduction. For all zones, their width at or close to the time of subduction initiation has also been plotted. The data generally plot in two clusters:

Flat Slab Subduction Relation With Subduction Zone Age and Width in Nature
(1) One cluster that plots toward the bottom left where most subduction zones are relatively narrow and young (e.g., Scotia, Manila, North Sulawesi, Nankai-Ryukyu, Puysegur), some are older but narrow (e.g., Hellenic, Cascadia, Lesser Antilles-Puerto Rico) and some are wide but young (e.g., Melanesia, Sunda, Northwest Pacific); and (2) Another cluster that is located toward the top right where subduction zones are both old and wide (e.g., South America). Those subduction zones that plot in the first domain all lack flat slab subduction as defined in Figure 1D (but note that low angle subduction as in Figure 1C can occur). Those subduction zones that plot in the second domain in the upper right region do have one or more segments that are characterized by flat slab subduction (e.g., segments 1a-c) but also segments that are not (e.g., segments 1d-f).

Geodynamic Model Evolution
The four models all show subduction through trenchward subducting plate motion and through trench retreat. The first 35-40 Myr of the model runs are characterized by upper mantle subduction, which are followed by whole mantle subduction as the deepest part of the slab sinks into the lower mantle. Each of the models shows a similar, relatively steep, slab geometry during upper mantle subduction (Figure 4). During whole mantle subduction the different models develop variability in slab geometry.
The minimum slab dip angle in the depth range 100-200 km is tracked to identify flat slab subduction. This minimum dip angle is comparable for the four experiments in the early subduction stage and of the order 60-70 • (Figure 5A). Once the slab tip approaches the 660 km discontinuity at ∼35 Ma, the minimum dip angles for the different experiments start to diverge. All models show an overall, long-term, decrease in dip angle, but the magnitude of decrease depends on the slab width, with the largest decrease occurring for the widest slab. The narrow and intermediate-width slabs only show a moderate decrease in slab dip (with slab dip minima of 48 and 38 • , respectively), while the wide and very wide slabs show a large decrease (with minima of 21 and 2 • , respectively). For the wide slab model, slab flattening starts in a late stage at ∼171 Myr (Figures 4C, 5A). Here, the start of slab flattening is defined as the onset of formation of the concave-upward slab hinge between ∼50-150 km depth in between the two convex-upward hinges. This slab flattening only occurs in the central part of the wide subduction zone (central ∼700 km) and does not occur in segments closer toward the lateral slab edges. Indeed, the local slab dip angle near the lateral slab edges is consistently steeper than in the center with observed differences in dip angle up to 27 • toward the end of the model run ( Figure 5B).
For the very wide slab model slab flattening starts at ∼134 Myr and flat slab subduction starts at ∼150 Myr, lasting ∼25 Myr (Figures 4D, 5), although briefly interrupted by a ∼4 Myr period where the local dip angle exceeds 10 • . Flat slab subduction is followed by approximately normal subduction with dip angles of ∼20-30 • until the end of the model run at ∼185 Myr. Note again that flat slab subduction is defined as the moment when the flat slab segment downdip of the concave hinge has a dip angle ≤10 • . Also note that time in the models depends on the scaling of viscosity and velocity. Average subduction rates in the models (3-4 cm/yr) are somewhat slower than the average present-day subduction rate on Earth (5.5 cm/yr, Schellart et al., 2007). If we scale the average subduction velocities in the models to those in nature (implying a η UM−Max = 2.7-3.6 × 10 20 Pa s in nature rather than η UM−Max = 5 × 10 20 Pa s, with all values falling Colors indicate the non-dimensional effective viscosity field. Note that the somewhat fuzzy-blocky appearance of the slab in the lower part of the upper mantle in the bottom panels of (B,C) is the result of the implemented spatial adaptive mesh, which gives a lower resolution in the lower part of the upper mantle compared to the upper part, and the many times that repopulation of cells with particles was required to avoid the occurrence of empty cells. The high number of repopulation exercises, in particular for W = 2000 and W = 6000, caused the boundaries between different particle fields to become progressively fuzzier. within the natural range), then flat slab subduction in the very wide slab model would be reached earlier, at ∼80-110 Myr rather than ∼150 Myr. But note that stresses are not affected by any rescaling of velocities.

Stress in the Mantle Wedge and Sub-Slab Domain
The vertical deviatoric normal stress (σ YY ) in the mantle wedge and in the sub-slab domain has been plotted in Figure 6 for models W600, W6000 and W-infinite. For all models it can be observed that σ YY is mostly tensile in the mantle wedge near the top of the slab and compressive in the sub-slab domain near the base of the slab. Furthermore, σ YY values are relatively small for the narrow slab model (Figures 6A,B), in particular in an advanced stage of subduction. Here, σ YY is in the range −1.0 to 2.7 MPa at the slab top from the tip of the wedge down to 200 km depth and in the range −1.2 to 0.4 MPa at the base of the slab within the 100-200 km depth range. At an early stage of subduction σ YY values are somewhat higher, with a maximum of ∼4.0 MPa (tension) at the slab top and a minimum of about −1.9 MPa (compression) at the base of the slab.
For the very wide slab model, σ YY values at an early subduction stage ( Figure 6G) are comparable to those of the narrow slab model at the slab top with a maximum of ∼4.0 MPa (tension), and are close to neutral (minimum of −0.2 MPa) at the base of the slab above 200 km depth. At an advanced subduction stage, however, σ YY is generally much higher compared to the narrow slab model, with a maximum of ∼7.5 MPa (tension) at the slab top and a minimum of about −2.0 MPa (compression) at the base ( Figure 6H). For the wide slab model W6000, the stresses in the central section (Figures 6E,F) are comparable to those of W-infinite (Figures 6G,H), but the stresses in the section close to the lateral slab edge in a late stage of subduction ( Figure 6D) are smaller.

Mechanism of Flat Slab Subduction in the Geodynamic Models
The major dip angle decrease for the wide and very wide slab models results from the large-scale whole mantle circulation they produce once the slabs enter the lower mantle. This flow drags the overriding plate trenchward, pushing the subduction hinge backward and forcing the overriding plate to thicken (Schellart, 2017) (forced hinge rollback or pushback, Figures 4C,D), while the slab itself resists rollback due to its great width (Schellart et al., 2007) and its partial anchoring in the lower mantle, thus attaining a progressively lower dip angle. The slab flattening that occurs in both models is facilitated by a relatively high effective suction force in the tip of the mantle wedge in a late stage of subduction. The very wide slab model with the 2D set-up, and thus an implied infinitely wide slab, develops the highest suction force in the tip that increases with time (Figures 6G,H), because lateral mantle inflow is absent, thereby facilitating flat slab formation. The wide slab model, despite its 6000 km width, allows for (minor) lateral mantle inflow around its lateral slab edges, reducing Note that the short-wavelength structures in panel (D) (the dark green and red blobs at ∼300-600 km depth in the mantle wedge and just below the base of the overriding plate) are formed by slab material from the lateral slab edge that has been advected into the mantle wedge by the rollback-induced toroidal upper mantle return flow. This toroidal flow, which is directed from the sub-slab and around the lateral slab edge into the mantle wedge (e.g., Schellart et al., 2007), has curved and sheared the slab edge into the mantle wedge and has also sheared off small fragments from the slab edge. Additionally, a consequence of the curvature is that the local trench orientation strikes at an oblique angle with respect to cross-section (D). The slab edge material has a higher viscosity than the ambient upper mantle material and therefore shows higher deviatoric stresses than the ambient upper mantle material. the wedge suction, and thus slab flattening develops more slowly (Figure 5A).
In case the slab is narrow or intermediate in width, slab rollback is directly driven by the negative buoyancy of the slab (slab-driven hinge rollback or pullback, Figures 4A,B) rather than being forced by the trenchward motion of the overriding plate, as is evident by the absence of overriding plate thickening but the presence of overriding plate thinning in these models. Rollback is facilitated by the presence of nearby lateral slab edges, which allow for efficient toroidal return flow from the sub-slab region toward the mantle wedge region (Dvorkin et al., 1993;Schellart et al., 2007), sustaining a relatively steep slab dip angle and lowering the suction in the mantle wedge. Indeed, σ YY just above the top surface of the slab is significantly lower for a narrow slab (Figure 6B), which has two nearby lateral slab edges, than for the center of a wide slab (Figure 6F), which has two distant lateral slab edges, and for an infinitely wide slab (Figure 6H), which has no lateral slab edges. Furthermore, for trench segments of wide subduction zones located close to lateral slab edges, which only have one nearby lateral slab edge, σ YY is intermediate to that of narrow slabs and wide slabs (Figure 6). What is more, the vertical slab lifting force (F L ) that results from σ YY is about an order of magnitude larger for a very wide slab than for a narrow slab. For W600, this lifting force, calculated from the tip of the mantle wedge down to 200 km depth, is ∼1.5 × 10 11 N per meter trench length. Here, ∼66% results from the vertical suction at the top of the slab, while ∼34% results from vertical compressive stresses at the base of the slab. For the very wide slab model the vertical lifting force is ∼1.4 × 10 12 N per meter trench length, with ∼76% resulting from the vertical suction at the slab top and ∼24% from the vertical compressive stress at the base of the slab. F L is also more significant for wide slabs when compared with the negative buoyancy force (F Bu ) of the same slab segment for which the lifting force has been calculated. For the very wide slab model F L /F Bu = ∼8.7%, while for the narrow slab model F L /F Bu = ∼1.7%.
The proposed important role of subduction zone size and age, and the proposed mechanism of slab-induced whole-mantle return flow in driving trenchward overriding plate motion, slab dip angle decrease and eventually slab flattening at wide subduction zones, can be reconciled with earlier modeling studies. Indeed, a number of earlier modeling works have advocated the important role of forced trench retreat (e.g., van Hunen et al., 2004;Manea et al., 2012), large wedge suction (e.g., Tovish et al., 1978;van Hunen et al., 2004;Manea et al., 2012;Rodríguez-González et al., 2012) and a thick overriding plate (e.g., Rodríguez-González et al., 2012) or far-field thick craton in the overriding plate (e.g., Manea et al., 2012;O'Driscoll et al., 2012;Taramón et al., 2015) in promoting flat slab subduction. The current work shows that, in a buoyancy-driven environment, forced trench retreat, large wedge suction and a thick overriding plate will develop in a wide subduction zone setting after a long period of subduction (of the order 100 Myr) once the slab tip has reached a significant depth in the lower mantle. The current work also shows that such features will not develop in a narrow or intermediate-width subduction zone setting, or close to lateral slab edges.
The mechanism for flat slab generation as proposed above requires a slab in nature to behave as a strong entity that is generally continuous from the surface to the lower mantle without major holes or tears in the upper mantle part of the slab. Indeed, it is generally thought that slab tearing can facilitate the collapse and removal of a flat slab (e.g., Schepers et al., 2017). For a strong and continuous slab, this requires natural slabs to have a greater viscosity than the ambient upper mantle, which is generally thought to be the case with slabs that are, on average, about 140-510 times (Ribe, 2010) or 100-700 times (Schellart, 2008) more viscous than the ambient mantle. The current models adopt comparable average viscosity ratios and slab deformation is accommodated through stretching, thickening, bending and buckling, not by slab tearing, thereby facilitating flat slab subduction. Another requirement is for slabs to be continuous from the surface down to the lower mantle, which is the case for a number of slabs in nature that include flat slab segments, as implied by mantle tomography studies. In particular, several tomography studies indicate that part of the Nazca slab subducting below South America continues down to ∼2400 km (e.g., van der Meer et al., 2018), and that part of the Cocos slab subducting along the Mexico-Central America subduction zone continues down to ∼2600-2800 km (e.g., Boschman et al., 2018).

Present-Day Flat Slabs in the Eastern Pacific
Flat slab subduction occurs in the Eastern Pacific in Peru, Central Chile, and Mexico and the proposed flat slab mechanism, involving a wide subduction zone (≥∼6000 km) and prolonged (>∼80-110 Myr) subduction, applies to these examples (Figure 2). Indeed, they formed at the Farallon/South American subduction zone at a time when it had been active for a long period (since the Middle Triassic-Early Jurassic) (Burchfiel and Davis, 1975;Coira et al., 1982;Boschman et al., 2018) and when it was extremely wide (7000-13,000 km). The Altiplano flat slab that formed in the central Andes at ∼35 Ma and terminated at ∼25 Ma (Ramos and Folguera, 2009) can also be explained in this manner. Although the Mexico flat slab is currently part of an intermediate-width subduction zone (Mexico-Central America), at the time when this flat slab formed at 30-25 Ma, as implied by the inboard migration of arc magmatism (Ferrari et al., 1999;Morán-Zenteno et al., 1999;Kim et al., 2010), the subduction zone was very wide (∼11,000 km) ( Table 1).
There are also several subduction segments of the South American subduction zone and Mexico-Central America subduction zone that do not contain a flat slab, including the Colombia, Northern Chile, Southern Chile and Central America segments. This observation implies that an old and wide subduction zone is a requirement for flat slab subduction to occur, but it is not a guarantee that a flat slab will always be present. The observation is consistent with the geodynamic modeling results, which show that for the very wide slab model flat slab subduction is a transient phenomenon, for this model lasting some 25 Myr, and is preceded and succeeded by normal subduction. The transient nature of flat slab subduction is at least partly explained by the periodic folding of the slab at the 660 km discontinuity and the initiation of flat slab subduction is likely triggered by a new slab folding phase. Such folding is also responsible for the periodic dip angle changes (Schellart, 2017) that characterize each model (Figure 5). In nature, flat slab subduction at an old and wide subduction zone could also be triggered by a slab folding phase at the 660 km discontinuity, or it could be triggered by another process such as the initiation of subduction of an aseismic ridge or plateau.

Latest Cretaceous-Early Cenozoic Flat Slab Subduction in Western North America
Flat slab subduction has been proposed for Western North America to explain the eastward migration of the magmatic arc and the eastward migration of the deformation front during cordilleran mountain building in the Late Cretaceous and Early Cenozoic (Henderson et al., 1984;DeCelles, 2004). Formation of the flat slab has been ascribed to the trenchward motion of the overriding North American plate and the subduction of aseismic ridges or oceanic plateaus (Henderson et al., 1984;Liu and Currie, 2016). The current work implies that flat slab subduction would not have been possible in the first place were it not for the extreme width (∼12,000 km) and longevity (∼110 Myr) of the subduction zone (Figure 2), and the large separation of the flat slab segment from lateral slab edges (Table 1).

Mid Triassic-Early Jurassic Flat Slab Subduction in South China
A model of flat slab subduction has been proposed to explain the broad intracontinental Mesozoic orogen in South China and the postorogenic magmatism in the region (Li and Li, 2007), with flat slab subduction starting in the mid Triassic (∼230 Ma) (Li and Li, 2007) along a subduction zone in the paleo Western Pacific. Reconstructions indicate that the subduction zone has been active since the mid Carboniferous (330 Ma) (Domeier and Torsvik, 2014) or Middle Devonian-Early Carboniferous (390-356 Ma) (Collins, 2003;Domeier and Torsvik, 2014). As such, the subduction zone had been active for 100-160 Myr when the flat slab developed along a subduction zone that was likely very wide (Figure 2 and Table 1), which is consistent with the proposed mechanism.

Absence of Present-Day Flat Slabs in the Western Pacific and Elsewhere
The general absence of current flat slabs (i.e., with a slab geometry as in Figure 1D) in the Western Pacific, Indian Ocean, Atlantic Ocean, and Mediterranean can be explained by the lack of at least one of two crucial ingredients, namely subduction zone longevity and large slab width (Figure 2). There are three large subduction zones (W ≥ 4000 km) in these regions, Sunda, Northwest Pacific and Melanesia (Schellart et al., 2007), but Sunda has been active only since the Middle Eocene (Hall, 2012), Northwest Pacific came into existence only after the Izu-Bonin-Mariana segment formed in the Early Eocene (Arculus et al., 2015), while Melanesia only formed in the Miocene (Hall, 2002; Table 2). Apart from these three subduction zones, the remaining ones in the Western Pacific, Atlantic, Indian Ocean and Mediterranean all have small or intermediate widths (Schellart et al., 2007). The lack of subduction zone longevity and/or large width did not allow for a significant decrease in slab dip angle (Figure 5) nor did it enhance suction forces (Figures 6A-C,E), thereby preventing subduction zones in these regions to form a flat slab. The lack of subduction zone longevity and/or large width also prevented those Western Pacific subduction zones subducting an aseismic ridge or plateau (e.g., d'Entrecasteaux Ridge at the New Hebrides subduction segment, Marcus-Necker Ridge at the Mariana-Bonin segment) to form a flat slab. This implies that buoyant ridge/plateau subduction is not sufficient on its own to cause flat slab subduction, in agreement with earlier work on Nazca Ridge subduction (Antonijevic et al., 2015). It is also consistent with recent geodynamic experiments of buoyancy-driven subduction with an aseismic ridge, indicating that subduction of a large ridge (e.g., Carnegie Ridge or Nazca Ridge) reduces the slab dip angle in the shallow part of the upper mantle by a mere ∼10 • and does not produce flat slab subduction on its own (Flórez-Rodríguez et al., 2019).

CONCLUSION
The observational constraints (Figure 2) and models (Figures 4-6) imply that flat slab subduction dominantly initiates at wide subduction zones (≥ ∼6000 km) and only after a prolonged period of subduction (≥ ∼80-110 Myr). This finding is consistent with earlier works advocating the role of forced trench retreat, large wedge suction and a thick or cratonic overriding plate in promoting flat slab subduction (e.g., van Hunen et al., 2004;Manea et al., 2012Manea et al., , 2017O'Driscoll et al., 2012;Rodríguez-González et al., 2012;Taramón et al., 2015;Schepers et al., 2017), because only in the center of wide subduction zones do forced trench retreat and overriding plate thickening develop in a late subduction stage, and is the effective wedge suction high. The conceptual model works for the proposed Laramide flat slab in western North America (Henderson et al., 1984;DeCelles, 2004), the Early Mesozoic South China flat slab (Li and Li, 2007), and the current flat slabs in the Eastern Pacific, as they developed at old and wide subduction zones (Figure 2). The finding can also explain the absence of present-day flat slabs in the Western Pacific, Indian Ocean, Atlantic Ocean and Mediterranean, because subduction zones here are either relatively young (< ∼80-110 Ma), less than ∼6000 km wide, or both. Subduction of a buoyant aseismic ridge or plateau in these regions will not induce flat slab subduction, as slabs in these regions will not have had sufficient time to reduce their dip angle at shallow depth and/or will not have a significant mantle wedge suction force due to the limited slab width. It thereby solves a longstanding debate on the lack of spatial correlation between ridge subduction and flat slab subduction in regions outside the Eastern Pacific.

DATA AVAILABILITY STATEMENT
All data and methods necessary to understand, evaluate, replicate, and build upon the reported research are presented in this article.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.

FUNDING
This work has been funded by a Vici Fellowship (016.VICI.170.110) from the Dutch National Science Foundation (NWO), and has been supported by computational resources from the NCI National Facility in Australia through the National Computational Merit Allocation Scheme (project qk0).

APPENDIX Justification for Data Selection and Uncertainties Associated With the Data
There are uncertainties associated with the data presented in Tables 1, 2 that are used as a basis for Figure 2. For the presentday slab width values listed in Table 2, the uncertainties are relatively small, in general less than 5-10%, which will therefore have no discernable impact on the data distribution as shown in Figure 2A. Uncertainties with estimated slab width values for the geological past will be larger and are generally hard to quantify, as they mostly result from differences in plate tectonic configurations in the published literature. For the relatively narrow subduction zones listed in Table 2 (numbers 6-14) the uncertainties can be up to several hundred kilometers, but this will have no discernable impact on the general distribution of data points in Figure 2A. For the slab width data of the other subduction zones, which have a moderate to very large width (numbers 1-5 and 15-22) there are larger uncertainties. These subduction zones generally also have a larger associated uncertainty in subduction zone age, or a subduction zone age that is more controversial. Below, each of these subduction zones will be briefly discussed and a justification for the choice of subduction zone age and width is provided.
South American Subduction Zone (Points 1a-f and 3 in Figure 2 and Tables 1, 2) The age of the South American subduction zone has long been interpreted as very old, around 200 Ma, and thus dating from the Late Triassic/Early Jurassic (e.g., Coira et al., 1982;Vásquez et al., 2011;Maloney et al., 2013). This old age has very recently been contested by Chen et al. (2019), who argue for an 80 Ma subduction zone in the north, becoming progressively younger southward. The young age for the South American subduction zone proposed by Chen et al. (2019) is in contradiction with geological data showing a continuous arc magmatic record from the Jurassic to the present (∼200-0 Ma), such as reported for northern Chile (Scheuber et al., 1994). It is also in contradiction with seismic tomographic data presented in van der Meer et al. (2018), which show a continuous slab from the trench down to 2400 ± 200 km depth in the Peru segment of the South American subduction zone, which they interpreted as subduction being active from 200-173 Ma (Jurassic) to the present. Recent work on sedimentary systems in the retroarc and arc region in the southern Central Andes also imply continuous subduction activity since the Late Triassic/Early Jurassic (e.g., Mackaman-Lofland et al., 2019). Based on these different lines of evidence, I adopt an old age (∼200 Ma) for the South American subduction zone.
Tectonic reconstructions show that at ∼200 Ma the entire western margin of the South American continent was an eastdipping subduction zone and generally also show it to be continuous with the Mexico-Central America segment (e.g., Seton et al., 2012). This gives a subduction zone width of the order 11,000 ± 1000 km.
For data points 1a, 1b and 1c, the flat slab formed at ∼4, ∼8, and ∼20 Ma and so formed when the Nazca plate had already separated from the Cocos plate, which occurred at ∼23 Ma (Lonsdale, 2005). Thus the width of the subduction zone was roughly the north-south length of the western margin of South America, which is of the order 7000 km. The flat slab for data point 3 formed at ∼35 Ma, so before Nazca-Cocos separation, and so the slab width at that time was roughly the sum of the South American subduction zone width and the width of the Baja California-Mexico-Central America segment.
Mexico-Central America Subduction Zone (Points 2a and 2b in Figure 2 and Tables 1, 2) The Mexican flat slab formed at ∼30-25 Ma (Ferrari et al., 1999;Morán-Zenteno et al., 1999;Kim et al., 2010), and the Cocos spreading ridge separating the Cocos plate from the Nazca plate formed only at ∼23 Ma (Lonsdale, 2005). Thus, the Nazca + Cocos single plate subducted eastward as one very large slab, most likely along one continuous subduction zone at 30-25 Ma that would be ∼10,500-12,000 km wide (some 7000 km for the South American subduction segment and some 3000-5000 km for the Mexico-Central America segment). This is consistent with tectonic reconstructions (e.g., Gordon and Jurdy, 1986;Seton et al., 2012;Müller et al., 2016).
Recent work indicates that at least the Mexico-Central America subduction segment has been active for a very long time. Indeed, according to Boschman et al. (2018) there is a continuous record of subduction since 220 Ma along the Mexico-Central America subduction zone, as indicated by a continuous slab geometry from the surface down to 2500-2800 km depth in seismic tomography models. According to van der Meer et al. (2018) this old subduction zone had a large lateral extent, stretching from the Cocos ridge in the south to the Mendocino triple junction in the north (their page 351), which would imply a subduction zone extent of ∼5800 km. This width is used as a lower limit of the slab width of the subduction zone during its formation. As an upper limit the reconstruction of Seton et al. (2012) is used, which shows a continuous subduction zone from Mexico in the north to Patagonia in the south at 200 Ma, which would thus imply an ∼11,000 km wide subduction zone at a time close to its inception.
Another reconstruction, that of Müller et al. (2016), shows a different tectonic setting in the Mexico-Central America-Caribbean region with a subduction flip and ∼west-dipping subduction at ∼130-90 Ma, which would not be consistent with the current work and that of Seton et al. (2012). The recent work of Boschman et al. (2018), however, can be used to refute the proposed tectonic scenario from Müller et al. (2016), as it shows a long and continuous subducted slab, from the surface down to ∼2500-2800 km depth, which implies long-lived continuous east-dipping subduction along the Central America-Mexico subduction zone. Figure 2 and Table 1) The age of inception and the width of the South China subduction zone have been derived from the reconstruction papers of Collins (2003) and Domeier and Torsvik (2014). It is clear that there are large uncertainties associated with the estimated age and width, because we are concerned with subduction that was active in the Paleozoic and Triassic, and flat slab subduction that started in the Triassic (Li and Li, 2007). Despite these large uncertainties, the results are consistent with the conceptual model proposed in the paper (Figure 2A).

South China Subduction Zone (Point 4 in
Farallon Subduction Zone (Points 5 and 15 in Figure 2 and Tables 1, 2) The Farallon subduction zone is generally thought to have been a very large, east-dipping subduction zone along the west coast of the Americas that subducted the Farallon oceanic plate during part of the Mesozoic and the Early Cenozoic (Burchfiel and Davis, 1975;Gordon and Jurdy, 1986;DeCelles, 2004). Earlier work indicates that the subduction zone is very old and formed in the Triassic or Jurassic (160-210 Ma) (Burchfiel and Davis, 1975;DeCelles, 2004), and relatively recent reconstructions imply a comparably old age (e.g., Seton et al., 2012). Such reconstructions also imply a very large subduction zone extent during subduction zone infancy exceeding 10,000 km (Seton et al., 2012). Only at ∼45-30 Ma, significantly after the inception of Farallon flat slab formation at ∼85-65 Ma (Henderson et al., 1984;Liu et al., 2010;Copeland et al., 2017), do reconstructions show that the Farallon plate, slab and subduction zone segment due to the formation of slab windows (Gordon and Jurdy, 1986;Schellart et al., 2010), resulting in the formation of a relatively small Farallon/Juan de Fuca plate subducting along the Cascadia subduction zone and a much larger Nazca-Cocos plate.
Nankai-Ryukyu Subduction Zone (Point 16 in Figure 2 and Table 2) The Nankai-Ryukyu subduction zone likely formed in the earliest Miocene (∼18 Ma) through subduction of the Philippine Sea plate below Eurasia (Lallemand et al., 2001). When the subduction zone formed, it was initially narrower (∼1600 km) and from ∼8 Ma it grew southwestward to include the southern Ryukyu subduction segment (Lallemand et al., 2001).
Lesser Antilles-Puerto Rico Subduction Zone (Point 17 in Figure 2 and Table 2) The age and the subduction zone width at the time of subduction zone formation of the Lesser Antilles subduction zone are somewhat controversial. The minimum age of 45 Ma is based on the oldest ages of Lesser Antilles arc magmatism (Burke, 1988). There is a ∼50 Ma age based on a reconstruction of Boschman et al. (2014), a 71-94 Ma age based on reported ages of arc magmatism from Hispaniola (Hastie et al., 2013), and a ∼120 Ma age as argued and shown in a reconstruction in Pindell et al. (2006). van der Meer et al. (2018) estimate the age using tomography and geological data and propose subduction initiation at 45-55 Ma for the lesser Antilles. Estimates of subduction zone width at subduction initiation vary significantly, as shown in tectonic reconstructions, including ∼1400 km (Boschman et al., 2014) and ∼3300 km (Pindell et al., 2006). Considering the controversy, the average is taken for both the subduction zone age and width.
Tonga-Kermadec-Hikurangi Subduction Zone (Point 18 in Figure 2 and Table 2) The subduction zone age and the width at the time of the Tonga-Kermadec-Hikurangi subduction zone formation are somewhat controversial as well. The minimum age is generally based on the oldest ages of fore-arc magmatism in the Tonga arc (49-51 Ma) (Meffre et al., 2012) and reconstructions that adhere to this relatively young formation age imply that the subduction zone width during initiation was relatively small (1800 ± 800 km) (Meffre et al., 2012). More recently, an even younger minimum age has been suggested for the formation of the subduction zone (∼30 Ma) based on the termination of the New Caledonia subduction zone around that time (van de Lagemaat et al., 2018). This very young age is not likely, though, as backarc spreading in the North Loyalty Basin and South Fiji Basin was already active before 30 Ma. Reconstructions from Schellart et al. (2006) show a subduction zone that was already active at 90 Ma, as implied by convergence between the Australian plate and the Pacific plate, with a large width (∼5000 km) that incorporates the Tonga-Kermadec-Hikurangi segment and the Vitiaz-Solomon segment. Due to the uncertainty and controversy regarding the slab width at subduction initiation and time of subduction initiation, an average between that of Meffre et al. (2012) and Schellart et al. (2006) has been taken. Figure 2 and Table 2) It is generally agreed upon that the Aleutians-Alaska subduction zone formed sometime in the Early Cenozoic as based on ages of the oldest arc magmatic rocks in the Aleutian Islands, with an older age of ∼55 Ma suggested earlier (Scholl et al., 1986) and a younger age (∼46 Ma) suggested more recently (Jicha et al., 2006). The younger age is here taken as the best estimate of the true age. The Aleutians segment is thought to have formed in an intra-oceanic setting and was limited in the west by the Shirshov Ridge, after which the subduction zone grew westward to form a cusp with the Kuril-Kamchatka arc (Schellart et al., 2003). Hence, the subduction zone had a smaller width (by ∼500 km) at subduction initiation than its present width.

Aleutians-Alaska Subduction Zone (Point 19 in
Melanesia Subduction Zone (Point 20 in Figure 2 and Table 2) It is generally agreed upon that the Melanesia subduction zone is very young and formed in the Late Cenozoic, around 10-15 Ma (Hall, 2002;Schellart et al., 2006;Richards et al., 2011). In the short time since 15 Ma it has grown to a relatively wide subduction zone with a width of ∼4400 km (Schellart et al., 2007).
Northwest Pacific Subduction Zone (Points 21, 21a, and 21b in Figure 2 and Table 2) The width of the Northwest Pacific subduction zone, or Kamchatka-Kuril-Japan-Izu-Bonin-Mariana subduction zone, is very large, ∼6550 km, but it formed rather recently due to the conjoining of the Izu-Bonin-Mariana segment and the Kamchatka-Kuril-Japan segment. The former segment formed around 52 Ma, as implied by dated forearc magmatic rocks (Ishizuka et al., 2011), while the latter segment formed at ∼60 Ma, after accretion of the Okhotsk terrane to Eurasia (Schellart et al., 2003). A maximum age of subduction zone formation is thus 52 Ma (assuming that the two segments conjoined immediately, which is not very likely due to the difference in strike of the two segments at that time), while a minimum age of subduction zone formation is based on a reconstruction of Sdrolias et al. (2004), who show that the two segments connected at ∼25 Ma.
Sunda Subduction Zone (Points 22, 22a, and 22b in Figure 2 and Table 2) The width of the slab subducting at the Sunda subduction zone is very large, ∼7850 km, as it includes the Burma-Andaman segment, the Sumatra-Java-Bali-Nusa segment and the Timor-Banda segment (Schellart et al., 2007). These individual segments, however, formed at different times. The Andaman subduction segment likely started around 94 Ma based on suprasubduction-zone ophiolites of this age, as dated by Sarma et al. (2010) and interpreted as such by Advokaat et al. (2018), but this subduction segment is only ∼1400 km wide (∼3150 km if one includes the Burma segment). The central segment of the Sunda subduction zone (Sumatra-Java-Bali-Nusa), which is ∼3300 km wide, formed at ∼50-45 Ma (Hall, 2012), while the easternmost segment (Banda segment), which is ∼1400 km wide, formed only at ∼15 Ma (Spakman and Hall, 2010). So from ∼94-50 Ma the slab was ∼3150 km wide, from ∼50-15 Ma it was ∼6450 km wide, and from ∼15-0 Ma it was ∼7850 km wide. As such, it acquired its large width only very recently.