Where does subduction initiate and cease? A global, Lecture notes of Topography

Download Where does subduction initiate and cease? A global and more Lecture notes Topography in PDF only on Docsity! EARTH AND PLANETARY SCIENCE LETTERS, VOL. ???, XXXX, DOI:XX.XXXX/, Where does subduction initiate and cease? A global1 scale perspective.2 Martina M. Ulvrova1, Nicolas Coltice2, Simon Williams3, and Paul J. Tackley1 Corresponding author: M. Ulvrova, Department of Earth Sciences, Institute of Geophysics, ETH Zürich, Sonneggstrasse 5, Zürich, 8092, Switzerland. (martina.ulvrova@erdw.ethz.ch) 1Institute of Geophysics, Department of Earth Sciences, ETH Zürich, Zürich, Switzerland. 2Laboratoire de Géologie, École Normale Supérieure, CNRS-UMR 8538, PSL Research University, Paris, France 3EarthByte Group, School of Geosciences, University of Sydney, Sydney, New South Wales, Australia D R A F T July 18, 2019, 4:11pm D R A F T X - 2 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Abstract. The thermo-mechanical evolution of the Earth’s mantle is largely3 controlled by the dynamics of subduction zones, which connect the surface4 tectonic plates with the interior. However, little is known about the system-5 atics of where subduction starts and stops within the framework of global6 plate motions and evolving continental configurations. Here, we investigate7 where new subduction zones preferentially form, and where they endure and8 cease using statistical analysis of large-scale simulations of mantle convec-9 tion that feature self-consistent plate–like lithospheric behaviour and con-10 tinental drift in the spherical annulus geometry. We juxtapose the results of11 numerical modelling with subduction histories retrieved from plate tectonic12 reconstruction models and from seismic tomography. Numerical models show13 that subduction initiation is largely controlled by the strength of the litho-14 sphere and by the length of continental margins (for 2D models, the num-15 ber of continental margins). Strong lithosphere favors subduction inception16 in the vicinity of the continents while for weak lithosphere the distribution17 of subduction initiation follows a random process distribution. Reconstruc-18 tions suggest that subduction initiation and cessation on Earth is also not19 randomly distributed within the oceans, and more subduction zones cease20 in the vicinity of continental margins compared to subduction initiation. Our21 model results also suggest that intra-oceanic subduction initiation is more22 prevalent during times of supercontinent assembly (e.g. Pangea) compared23 to more recent continental dispersal, consistent with recent interpretations24 of relict slabs in seismic tomography.25 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 5 the continents than more recent examples. This raises the question of how important the70 presence of continents is to the life cycle of subduction systems, and whether this influence71 varies between periods of supercontinent assembly and continental dispersal.72 In this paper, we investigate the pattern of subduction initiation and cessation in time73 and space using numerical simulations of mantle convection. Numerical models are de-74 signed in a spherical annulus geometry, and we vary the number of continents, the strength75 of the lithosphere and its structure. We compare aspects of the modelling results to recon-76 structed subduction histories based on both plate kinematics and analysis of slab remnants77 imaged by seismic tomography [van der Meer et al., 2010, 2012; Müller et al., 2016]. We78 focus on global scale models to infer how the continental configuration and the plate layout79 change the distribution of subduction initiation and cessation, and the lifespan of subduc-80 tion zones. The models indicate that subduction initiation is non-randomly distributed81 in the ocean, and cessation happens mostly in the vicinity of continents. Our calculations82 point to di↵erent distributions of subduction inception between phases of supercontinents83 and phases in which continents are dispersed..84 2. Method In order to investigate statistically the spatial relationships between continents and the85 initiation, evolution and cessation of subduction zones, we numerically calculate the solu-86 tion of mantle convection in a spherical annulus [Hernlund and Tackley, 2008] using the87 StagYY code [Tackley, 2008]. The choice of geometry is motivated by the necessity of88 having long temporal series of several billions of years. Employment of the spherical annu-89 lus ensures similar scaling properties compared to the full 3D spherical shell. The model90 features self-consistently generated plate-like surface tectonics and drifting continents.91 D R A F T July 18, 2019, 4:11pm D R A F T X - 6 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION 2.1. Physical and numerical model We determine temperature, velocity, pressure and composition within the mantle by92 solving the equations of conservation of mass, momentum and energy and the advec-93 tion of material composition considering an incompressible mantle under the Boussinesq94 approximation. Below, the equations are given in their dimensionless form.95 r · v = 0 , (1)96 97 r · ⇣ ⌘ ⇣ rv + (rv)T ⌘⌘ rp = Ra (T +B C) er , (2)98 99 @tT + v ·rT = r2 T +H , (3)100 101 @tC + v ·rC = 0 , (4)102 with v the velocity, p the static pressure, ⌘ the viscosity, T the temperature, C the103 composition, H the internal heating rate, Ra the Rayleigh number, B the buoyancy ratio104 and er the radial unit vector. @t is the partial time derivative.105 Viscosity ⌘ follows the Arrhenius law and is strongly temperature and pressure depen-106 dent107 ⌘(T, p) = ⌘A exp ✓ Ea + p Va RT ◆ , (5)108 where Ea = 166 kJmol1 is the activation energy (kept constant for all simulations),109 Va = 6.34 · 107 m3 mol1 the activation volume (constant for all simulations) and R =110 8.314 Jmol1 K1 the gas constant. We give all parameters in Table 1. ⌘A is set such that111 ⌘ matches the reference viscosity ⌘0 at zero pressure and at temperature 1600 K, which112 is the expected temperature at the base of the lithosphere. We apply a viscosity cut o↵113 at 104 times ⌘0 . The viscosity varies over 6 orders of magnitude over the temperature114 variation T , the superadiabatic temperature drop over the mantle. Independently of115 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 7 temperature, viscosity increases exponentially by an order of magnitude with depth. The116 lowest values of the viscosity are in the asthenospheric mantle, while at the core mantle117 boundary viscosity is about 20 times lower than ⌘0.118 To localize deformation in narrow zones and obtain realistic plate boundaries at the119 surface, we use a pseudoplastic rheology [Moresi and Solomatov, 1998; Tackley, 2000a, b].120 After reaching a certain threshold value, the yield stress Y, the rocks undergo plastic121 yielding. Y is depth dependent and follows122 Y = 0 + d 0 Y , (6)123 where 0 is the surface yield stress, d is the depth and 0 Y is the yield stress depth124 derivative. If the stress reaches Y, we calculate the e↵ective viscosity ⌘e↵ on the grid125 ⌘e↵ = min [⌘(T, p), ⌘Y] , (7)126 with ⌘Y as127 ⌘Y = Y 2✏̇II . (8)128 ✏̇II is the second invariant of the strain rate tensor. We vary the surface strength 0 of129 the lithosphere between 7 MPa and 56 MPa while keeping the gradient constant for all130 simulations at 810Pam1.131 Using this kind of rheology results in the self-consistent formation of strong plate in-132 teriors moving with a uniform velocity delimited by narrow plate boundaries character-133 ized by reduced viscosity and an abrupt velocity change [Moresi and Solomatov, 1998;134 Tackley, 2000b]. Importantly, such rheology successfully reproduces seafloor age dis-135 tributions [Coltice et al., 2012] and is suciently realistic to investigate global surface136 tectonics [Coltice et al., 2017; Ulvrova et al., 2019].137 D R A F T July 18, 2019, 4:11pm D R A F T X - 10 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Earth’s Rayleigh number, but lies at the edge of the computational feasibility for the182 given rheology and the presence of sticky air.183 We use a resolution of 128 ⇥ 1024 cells in the radial and horizontal directions, respec-184 tively, except for the case with the free surface when we use a grid with 256⇥ 2048 cells.185 Vertical grid refinement close to the top and bottom limits is employed resulting in 10 km186 and 15 km (5 km and 8 km for the case with the sticky air) thick cells at the surface and187 at the core-mantle boundary (cf. more details in the supplementary material). To track188 composition, we use 4 ⇥ 107 tracers. This means that on average there are around 300189 tracers per cell except in simulations with a free surface, which have higher resolution. In190 this case, the number of tracers per cell is around 75.191 First, we run a model until it reaches a statistically steady state, at which the heat192 budget is balanced and characteristic properties of the system such as mean temperature,193 mean velocity and average surface heat flux fluctuate around some constant values. In this194 initial stage, we do not advect tracers and hence the composition field remains the same195 throughout the initialization period, which means that continents do not move from their196 initial positions. Once statistically steady state is achieved, we use a random snapshot197 from the equilibrated evolution to start the calculation with continents that move freely.198 The choice of the particular snapshot from the statistically steady state does not have any199 influence on the statistics performed on the system. The model is run until a sucient200 number of subduction initiations N is collected. The shortest analyzed period is 1 Gy with201 27 subduction initiation events. The longest analyzed period is 7 Gy with 288 detected202 subduction initiations. The length of the simulation does not have any influence on the203 results as soon as N is large enough and statistically representative of the system. The204 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 11 parameters of the calculations are listed in Table 1. We vary the strength of the oceanic205 lithosphere and number of continents, and we test how the presence of the weak crustal206 layer and presence of a free surface alter the results. In total, we run 13 models having207 di↵erent parameterizations.208 2.2. Model analysis To detect subduction inception and follow its evolution, we analyze the divergence of209 the surface velocity. As soon as a negative peak of the surface velocity divergence appears210 (the threshold of detection being minus one tenth of the surface rms velocity), a new211 subduction zone is formed. The motion of this peak is then tracked through time and we212 record the distance to the closest continental margin until the peak disappears. For each213 case we specifically register this distance at subduction inception and at cessation. These214 values are further analyzed by calculating a cumulative distribution of the distance to215 the closest continental margin specifically when subduction starts and ends. To construct216 these distributions we bin the range of all possible distances into 500 km wide intervals.217 We then count the cumulative number of detected distances falling into a specific bin.218 All subduction zones that initiate at a continental margin fall into the first bin. The219 long duration of each model allows us to collect up to several hundreds of subduction220 initiations, typically several tens of them.221 To investigate systematic patterns within the distributions retrieved from the models,222 we use a Monte Carlo method and calculate the synthetic distribution for subduction223 zones initiating at random positions within oceans. At each model time we generate a set224 of randomly located points within the oceans, which together we take to have a spatial225 distribution equivalent to scenarios where subduction intitiation or cessation locations226 D R A F T July 18, 2019, 4:11pm D R A F T X - 12 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION is also randomly distributed. For models with one continent, the cumulative frequency227 of random point locations as a function of distance to the continent is a straight line228 since the distance between the two continental margins is constant through time and the229 probability that a subduction zone initiates at a certain location is uniform. For more230 than one continent, the distribution is typically curved. This allows us to show how the231 modelled subduction distributions deviate from the random one. However, it is necessary232 that a suciently large number of random subduction initiations is generated at each233 time level, in total at least several thousand. To estimate the variance of the random234 distribution, we calculate distributions of N subduction initiation events that happen235 at random times and random positions. N is the number of initiations detected in the236 particular model.237 3. Results The organization of the system dictates its dynamic evolution: sinking slabs drive238 plate motion while inducing the mantle flow that in turn is at the origin of the plate239 motion (Figure 1). Plate-like behaviour is developed self-consistently using the yielding240 rheology (cf. Section 2.2) with a surface velocity that is constant in plate interiors while241 changing abruptly over plate boundaries. Since the system is heated mainly by internal242 heat sources (and we keep the internal heating rate constant for all simulations), sinking243 slabs and surface plates that compose the upper boundary layer dominate over plumes244 created at the core-mantle boundary. In particular, around 20% of the total surface heat245 flux is due to heating the mantle from the core. . The core derived fraction of the246 total surface heat flux is very similar for all models, with di↵erences smaller than 1%.247 When looking at the temperature structure of the system (Figure 1b,e,h), one can note248 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 15 Subduction zones that are formed next to a continental margin stay glued to the con-294 tinent their whole existence (Figure 5). In contrast, subduction zones that initiate as295 intra-oceanic reside in the oceans where they either merge with another subduction zone296 or migrate in the oceans before reaching a continental margin where they endure for some297 time before ceasing (Figure 5). In some cases, a subduction zone retreats toward the con-298 tinent and once it hits the margin, subduction continues with a reversed polarity (Figure299 2). Rarely, a subduction zone is formed and ceases as intra-oceanic without colliding with300 another convergence zone (Figure 5). In the weak lithosphere case, termination is random301 just like initiation (Figure 3a). Where the continents have an influence, subduction zones302 are more likely to terminate adjacent to continents than to initiate there, which is pre-303 sumably a consequence of intra-oceanic subduction zones being able to migrate freely, but304 once they migrate to a continental margin subduction zones stay there until they cease305 (Figure 3b and c).306 3.2. Influence of number of continental margins Repeated continental assembly and dispersal are observed on Earth in the cycles that307 last for several 100 Myr [e.g. Rogers and Santosh, 2004]. The length of the continen-308 tal margins thus varies according to the continental configuration, being minimal when309 continents are aggregated and maximal when dispersed. To investigate the influence of310 the number of continental margins on the position of subduction initiation within our311 annulus models, we perform a set of calculations with one, two and three continental rafts312 (where each raft has two margins), while keeping the total cover of continents constant at313 30% of the annulus. The number of continental margins in 2D models corresponds to the314 length of the continental margins in 3D. In these models, we keep the yield stress fixed315 D R A F T July 18, 2019, 4:11pm D R A F T X - 16 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION (0 = 35MPa) and include a weak crustal layer at the surface. The number of subduction316 zones fluctuates around four or five (cf. histograms on Figure 6). With increasing num-317 ber of continental margins, the number of intra-oceanic subduction initiations decreases,318 regardless of whether a weak crustal layer is present (cf. Figure 4b). In particular, subduc-319 tion initiation at continental margins increases from 30% for the case with one continent320 to 67% for three continents (cf. Figure 4b and 6). Two fundamental results are consistent321 across all model cases: firstly, all cases significantly di↵er from initiation at random po-322 sition for the chosen lithospheric strength; secondly, we systematically observe that more323 subduction zones cease at a continental margin compared to the initiation position (cf.324 Figure 6). Both these relationships are weakest for the case with a single continent (cf.325 Figure 6a). The threshold distance below which the e↵ect of continents is negligible and326 the distribution follows a random distribution of subduction initiation increases with in-327 creasing number of continental margins and is around 4000 km, 6000 km and 7000 km for328 cases with two, four and six continental margins (Figure 6).329 3.3. Asymmetric subduction zones At the surface of the Earth, during a collision of two oceanic plates, one of them subducts330 into the mantle while the overriding plate stays at the surface. Although numerical mod-331 els still have limited ability to produce such behavior, strong asymmetry of sinking slabs332 is observed in our models as is described in the beginning of this section. The simplest333 model, with a free-slip surface and no weak crustal layer features the least realistic sub-334 duction dipping angles and longer periods of vertical descent. In this case, about 45% of335 subduction inceptions are adjacent to the continental margins (Figure 7a). For more com-336 plex models with more realistic slab dip angles due to lubrication and partial decoupling337 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 17 at the interface between the two colliding plates (models with the weak crustal layer), it is338 more likely that subduction initiates at the continental margin; close to 60% of detected339 subduction zones start in the vicinity of the continent (Figure 7b). This is because the340 lateral density gradient between continental and oceanic lithosphere produces additional341 compressive stresses that would drive continents to spread below the free-slip boundary,342 if the viscosity was low enough. Hence compressive stresses focus at the continent ocean343 boundary [Nikolaeva et al., 2010; Rolf and Tackley, 2011] and favour subduction initiation.344 In combination with a weak crustal layer having lower yield stress and hence localizing345 deformation more easily, more subduction zones initiate at the continental margins. For346 the models with a free, deformable surface that features more realistic slab dips and stress347 states, the number of initiation events at the continental margins drops to around 30%348 (Figure 7c). In this case, the free surface allows the continent to hamper continent spread-349 ing by generating a topography that can accommodate a fraction of the stresses at the350 continent ocean boundary and the lateral density di↵erence between the continents and351 the mantle through isostasy.352 3.4. Reconstructions of subduction, initiation, and cessation The timing and location of past subduction can be reconstructed from geological and353 geophysical constraints, though such reconstructions are subject to large uncertainties354 over the timescales of supercontinent cycles (several hundred Myr). In particular, the355 lengths and locations of plate boundaries within the oceanic realm far from continents356 are poorly known, increasingly so further back in time, since the crust that comprised357 these regions is scarcely preserved at the present day. We compare two reconstructions358 of subduction history (Figure 8) that are derived by di↵erent methods, and use these as359 D R A F T July 18, 2019, 4:11pm D R A F T X - 20 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION that bound the adjacent segments di↵ers between successive plate boundary snapshots 1406 Myr apart (occurring for example in cases of reconstructed subduction polarity reversal,407 which we consider as termination of existing subduction and initiation of a new segment).408 Locations of subduction initiation and cessation defined in this way are highlighted in409 animation S2, and the distances to the nearest continents at these points form the basis410 of the cumulative distribution functions in Figure 9c.411 For subduction history interpreted from seismic tomography, we compute the distance412 to the nearest continent of the line geometries defined by these previous studies at discrete413 reconstruction times being 7 to 17 Myr apart. We make an important distinction between414 slabs tabulated in van der Meer et al. [2010] (Figure 9e) and the longer list of slabs415 considered in van der Meer et al. [2012] (Figure 9d). For the former, the top and bottom416 slab ages are available and we take them as the timings of subduction termination and417 initiation respectively. The histories of additional slabs mapped in van der Meer et al.418 [2012] are not defined with the same level of detail. Consequently we do not include419 pre-Cenozoic subduction interpreted from slab remnants beneath the Pacific in two of our420 distribution plots (Figure 9e,f).421 With these limitations in mind, we ask the question as to whether the patterns of sub-422 duction initiation and cessation contained within current reconstructions follow a similar423 distribution to that observed in our numerical models. The distribution of subduction424 in relation to continents from the alternative subduction histories is illustrated using cu-425 mulative distribution functions (Figure 9), providing some analogy to the annulus model426 results. For both subduction histories, the overall distribution of ongoing subduction is427 typically closer to the continents than expected for randomly distributed points within428 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 21 the oceans (Figure 9b,d,e). This pattern is particularly pronounced for the Müller et al.429 [2016] kinematic reconstruction (Figure 9b), which lacks many intra-oceanic systems in-430 terpreted by van der Meer et al. [2010, 2012] from seismic tomography. When we isolate431 subduction initiation and cessation (Figure 9c,f), a further trend emerges that is apparent432 in both kinematic and tomography-based subduction histories - subduction cessation is433 typically closer to the continents than subduction initiation. This trend is only absent in434 the poorly resolved pre-100 Ma section of the Müller et al. [2016] reconstruction and the435 results do not show an obvious distinction between the periods before and after 100 Ma,436 broadly corresponding to the periods of initial and later stages of dispersal of Pangea.437 However, since the distributions for tomography-based initiation and termination do not438 include additional slabs interpreted to have existed within the middle of the Panthalassic439 ocean (Figure 8), the proportion of cases far from continents while Pangea was assembled440 are likely to be underestimated in our plots.441 4. Discussion Previous studies into the e↵ect of continents on mantle dynamics have shown that con-442 tinents increase the wavelength of the convective flow [e.g. Guillou and Jaupart, 1995] and443 influence heat loss out of the system as they act as thermal insulators [e.g. Lenardic and444 Moresi, 1999; Rolf et al., 2012]. Importantly, numerical simulations and laboratory exper-445 iments suggest that continents change the lithospheric stress distribution and facilitate446 subduction initiation [e.g. Nikolaeva et al., 2010; Rolf and Tackley, 2011]. However, sys-447 tematic study of the locations of subduction initiation and their ensuing evolution taking448 into account global tectonic settings has received very little attention.449 D R A F T July 18, 2019, 4:11pm D R A F T X - 22 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Comparison between the distribution of subduction in numerical simulations and those450 inferred from reconstructions o↵ers insight into the most plausible model parameters.451 Models in which the lithospheric strength is low, such that subduction is e↵ectively ran-452 domly distributed within the oceans, are inconsistent with the inferred patterns of sub-453 duction on Earth during the last >200 Myr (Figure 9). Instead, the distributions from454 reconstructions are more compatible with scenarios in which the lithospheric strength is455 relatively high, such that the continents generate a strain shadow promoting subduction456 initiation closer to the continents. This region adjacent to margins is loaded by the den-457 sity and topography contrast between continent and ocean lithosphere. Therefore, the458 lithosphere more readily yields upon experiencing additional stresses produced by convec-459 tion, whereas the same convective stresses alone are less likely to reach the yield criterion460 further from continents.461 In models where sites of initiation of subduction are biased towards regions close to462 continents, the distribution of subduction throughout its duration and cessation is also463 naturally biased towards these regions. Subduction that initiates at a continental margin464 remains there until it ceases. Subduction systems initiating in the oceans may migrate465 towards or away from the nearby continental margins; those that reach the continent mar-466 gin become continental arcs, and remain so until subduction ceases. Hence, the control467 of continents on patterns of subduction initiation influences the distribution of subduc-468 tion cessation, such that subduction termination along continental margins occurs more469 frequently than subduction initiation, even in the absence of continent-continent collision.470 Comparison between models with di↵erent numbers of continents (and therefore mar-471 gins) o↵ers an insight into the di↵erent distributions of subduction that we may expect472 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 25 plates have memory of previous yielding and can be damaged or undergo healing [e.g519 Bercovici and Ricard, 2016].520 5. Conclusions We present an assessment of where subduction initiates and ceases in global convec-521 tion models with a plate-like surface and continental drift. We compare the results of522 numerical simulations with distributions of subduction initiation and cessation retrieved523 from plate tectonics reconstructions and seismic tomography models. We show that the524 location of subduction initiation and cessation is not randomly distributed within the525 oceans on Earth. Subduction zones that are formed at continental margins tend to stay526 there, while subduction zones formed within the oceans migrate and merge with other527 intra-oceanic subduction zones, or reach continental margins where subduction usually528 continues, changing polarity before eventually ceasing. Hence, we systematically find529 that more subduction zones cease in the vicinity of continental margins compared to sub-530 duction initiation. Numerical models indicate that the critical parameters that influence531 the position of subduction initiation are the lithospheric strength and the number of con-532 tinental margins. Stronger lithosphere (which implies larger plates and fewer subduction533 zones [Tackley, 2000b]) increases the probability of subduction initiation in the vicinity of534 continental margins. With our numerical simulations we also predict that intra-oceanic535 subduction initiation is more likely during the times of supercontinent assembly, while536 continental dispersal favors incipient subduction close to continents. These results favour537 interpretations of intra-oceanic subduction systems within the Panthalassa Ocean during538 the time of Pangea based on seismic tomography [van der Meer et al., 2012; Van Der Meer539 et al., 2014], which are missing from earlier plate tectonic reconstructions.540 D R A F T July 18, 2019, 4:11pm D R A F T X - 26 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Acknowledgments. The support for this research has been provided by the Eu-541 ropean Union’s Horizon 2020 research and innovation program under the ERC grant542 agreement no617588 and the Marie Sk lodowska Curie grant agreement no753755.543 Simulations were performed on the AUGURY super-computer at P2CHPD Lyon.544 The StagPy library was used in this study to process StagYY output data545 (https://github.com/StagPython/StagPy). The pygplates library was used to anal-546 yse plate tectonic reconstructions (https://www.gplates.org/docs/pygplates/).547 References Auzende, J.M., Lafoy, Y., Marsset, B., 1988. Recent geodynamic evolution of the north548 Fiji basin (southwest Pacific). Geology 16, 925–929.549 Baes, M., Sobolev, S.V., Quinteros, J., 2018. Subduction initiation in mid-ocean induced550 by mantle suction flow. Geophys. J. Int. 215, 1515–1522.551 Bercovici, D., Ricard, Y., 2016. Grain-damage hysteresis and plate tectonic states. Phys.552 Earth Planet. Inter. 253, 31 – 47.553 Braz, C., Seton, M., Flament, N., Müller, R.D., 2018. Geodynamic reconstruction of an554 accreted cretaceous back-arc basin in the northern andes. Journal of Geodynamics 121,555 115–132.556 Christensen, U.R., 1985. Heat transport by variable viscosity convection ii: pressure557 influence, non-newtonian rheology and decaying heat sources. Physics of the Earth and558 Planetary Interiors 37, 183 – 205.559 Coltice, N., Gérault, M., Ulvrová, M., 2017. A mantle convection perspective on global560 tectonics. Earth-Sci. Rev. 165, 120 – 150.561 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 27 Coltice, N., Rolf, T., Tackley, P.J., Labrosse, S., 2012. Dynamic causes of the relation562 between area and age of the ocean floor. Science 336, 335–338.563 Crameri, F., Schmeling, H., Golabek, G.J., Duretz, T., Orendt, R., Buiter, S.J.H., May,564 D.A., Kaus, B.J.P., Gerya, T.V., Tackley, P.J., 2012. A comparison of numerical sur-565 face topography calculations in geodynamic modelling: an evaluation of the sticky air566 method. Geophys. J. Int. 189, 38–54.567 Crameri, F., Tackley, P.J., 2014. Spontaneous development of arcuate single-sided sub-568 duction in global 3-D mantle convection models with a free surface. J. Geophys. Res.:569 Solid Earth 119, 5921–5942.570 Dilek, Y., Furnes, H., 2009. Structure and geochemistry of tethyan ophiolites and their571 petrogenesis in subduction rollback systems. Lithos 113, 1–20.572 Doin, M.P., Fleitout, L., Christensen, U., 1997. Mantle convection and stability of depleted573 and undepleted continental lithosphere. J. Geophys. Res.: Solid Earth 102, 2771–2787.574 Domeier, M., Shephard, G.E., Jakob, J., Gaina, C., Doubrovine, P.V., Torsvik, T.H.,575 2017. Intraoceanic subduction spanned the Pacific in the Late Cretaceous–Paleocene.576 Science Advances 3, eaao2303.577 Flament, N., Williams, S., Müller, R., Gurnis, M., Bower, D.J., 2017. Origin and evolution578 of the deep thermochemical structure beneath eurasia. Nature communications 8, 14164.579 Gerya, T., Stern, R., Baes, M., Sobolev, S., Whattam, S., 2015. Plate tectonics on the580 earth triggered by plume-induced subduction initiation. Nature 527, 221–225.581 Gerya, T.V., Connolly, J.A., Yuen, D.A., 2008. Why is terrestrial subduction one-sided?582 Geology 36, 43–46.583 D R A F T July 18, 2019, 4:11pm D R A F T X - 30 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Mortimer, N., Campbell, H.J., Tulloch, A.J., King, P.R., Stagpoole, V.M., Wood, R.A.,628 Rattenbury, M.S., Sutherland, R., Adams, C.J., Collot, J., et al., 2017. Zealandia:629 Earths hidden continent. GSA Today 27, 27–35.630 Müller, R., Landgrebe, T., 2012. The link between great earthquakes and the subduction631 of oceanic fracture zones. Solid Earth 3, 447–465.632 Müller, R.D., Seton, M., Zahirovic, S., Williams, S.E., Matthews, K.J., Wright, N.M.,633 Shephard, G.E., Maloney, K.T., Barnett-Moore, N., Hosseinpour, M., Bower, D.J.,634 Cannon, J., 2016. Ocean Basin Evolution and Global-Scale Plate Reorganization Events635 Since Pangea Breakup. Ann. Rev. Earth Planet. Sci. 44, 107–138.636 Nikolaeva, K., Gerya, T.V., Marques, F.O., 2010. Subduction initiation at passive mar-637 gins: Numerical modeling. J. Geophys. Res.: Solid Earth 115, B03406.638 Reagan, M.K., Ishizuka, O., Stern, R.J., Kelley, K.A., Ohara, Y., Blichert-Toft, J.,639 Bloomer, S.H., Cash, J., Fryer, P., Hanan, B.B., et al., 2010. Fore-arc basalts and640 subduction initiation in the izu-bonin-mariana system. Geochem., Geophys., Geosyst.641 11.642 Rogers, J., Santosh, M., 2004. Continents and Supercontinents. Oxford University Press.643 Rolf, T., Coltice, N., Tackley, P., 2012. Linking continental drift, plate tectonics and the644 thermal state of the Earth’s mantle. Earth Planet. Sci. Lett. 351, 134–146.645 Rolf, T., Tackley, P., 2011. Focussing of stress by continents in 3D spherical mantle646 convection with self-consistent plate tectonics. Geophys. Res. Lett. 38.647 Sigloch, K., Mihalynuk, M.G., 2013. Intra-oceanic subduction shaped the assembly of648 cordilleran north america. Nature 496, 50.649 D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 31 Stern, R.J., 2004. Subduction initiation: spontaneous and induced. Earth Planet. Sci.650 Lett. 226, 275–292.651 Tackley, P.J., 2000a. The Quest for Self-Consistent Generation of Plate Tectonics in652 Mantle Convection Models. American Geophysical Union, Washington, D.C.. volume653 121 of Geophys. Monograph. pp. 47–72.654 Tackley, P.J., 2000b. Self-consistent generation of tectonic plates in time-dependent,655 three-dimensional mantle convection simulations 1. pseudoplastic yielding. Geochem.,656 Geophys., Geosyst. 1.657 Tackley, P.J., 2008. Modelling compressible mantle convection with large viscosity con-658 trasts in a three-dimensional spherical shell using the yin-yang grid. Phys. Earth Planet.659 Inter. 171, 7–18.660 Tackley, P.J., King, S.D., 2003. Testing the tracer ratio method for modeling active661 compositional fields in mantle convection simulations. Geochem., Geophys., Geosyst.662 4, 8302.663 Ulvrova, M.M., Brune, S., Williams, S., 2019. Breakup without borders: How continents664 speed up and slow down during rifting. Geophysical Research Letters 46, 1338–1347.665 Van Der Meer, D.G., Zeebe, R.E., van Hinsbergen, D.J., Sluijs, A., Spakman, W., Torsvik,666 T.H., 2014. Plate tectonic controls on atmospheric co2 levels since the triassic. Proc.667 Nat. Acad. Sci. 111, 4380–4385.668 Williams, S., Flament, N., Müller, R.D., Butterworth, N., 2015. Absolute plate motions669 since 130 ma constrained by subduction zone kinematics. Earth and Planetary Science670 Letters 418, 66–77.671 D R A F T July 18, 2019, 4:11pm D R A F T X - 32 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Nondimensional Dimensional Variable Symbol Value Value Gravitational acceleration g - 9.81m s2 Mantle thickness D 1 2890 km Thermal expansivity ↵0 - 3⇥ 105 K1 Thermal di↵usivity  - 106 m2 s1 Thermal conductivity k - 4Wm1 K1 Gas constant R - 8.314 Jmol1 K1 Reference density ⇢0 1 3300 kgm3 Internal heating rate H 20 5.44⇥ 1012 Wkg1 Reference viscosity ⌘0 1 6⇥ 1022 Pa s Activation energy Ea 8 170 kJmol1 Activation volume Va 3 6.34 · 107 m3 mol1 Surface temperature Ts 0.12 300 K Superadiabatic temperature drop T 1 2500K Rayleigh number Ra 106 - Surface yield stress in oceans 0 103 to 8⇥ 103 7 MPa to 56 MPa Yield stress depth derivative in oceans 0 Y 3.3⇥ 105 810Pam1 Table 1. Dimensional and non-dimensional parameters of the convection model. D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 35 Figure 3. Cumulative distribution of subduction initiation (blue) and cessation (black) as a function of distance from the nearest continent for three di↵erent yield stresses 0 (increasing from left to right). Dark grey line represents the distribution of subduction initiation at random position for a large population of cases. Gray area designates random distributions generated for N subduction initiations with N being the number of initiations detected for a particular model. Number of subduction zones detected is (from left to right) 89, 93 and 78. Histograms in the bottom right corners show the distribution of the number of subduction zones in the respective models. The models have a free slip top boundary but no weak crustal layer. D R A F T July 18, 2019, 4:11pm D R A F T X - 36 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Figure 4. a) Proportion of all subduction zones that initiate (blue) and cease (black) in the vicinity of the continent as a function of the lithospheric strength. One continental raft is present (i.e., two margins) throughout the simulations. b) Proportion of all subduction zones that initiate (blue) and cease (black) in the vicinity of the continent as a function of the number of the continental margins. Solid line is for models with compositionally uniform oceanic lithosphere while dashed line is for runs with weak crustal layer. The yield stress is fixed at 0 = 35MPa. D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 37 Figure 5. Position of the continents (blue) and subduction zones (coloured lines; one colour corresponds to individual subduction zone) through time together with the surface heat flux (gray scale). The model has a weak crustal layer and intermediate yield stress 0 = 35MPa. D R A F T July 18, 2019, 4:11pm D R A F T X - 40 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION a) b) c) d) e) f) g) h) Figure 8. Position of the continents and subduction zones since the Triassic according to two alternative reconstructions (see text), subdivided into 4 distinct time windows from Pangea times to recent. The detailed time-evolution of these reconstructions is illustrated in animations S2-S3. a) Subduction zones and continent positions for the M2016 model between 0 and 50 Ma, plotted at 10 Myr increments; locations of subduction zones are shown in colours corresponding to the color legend, while the continents are shown in gray with darker gray standing for younger positions within the 0-50 Myr period.b) same as a for 50-100 Ma; c) same as a for 100-150 Ma; d) same as a for 150-230 Ma; e) V2012 model for times between 0-50 Ma f) same as e for 50-100 Ma; g) same as e for 100-150 Ma h) same as e for 150-235 Ma. D R A F T July 18, 2019, 4:11pm D R A F T ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION X - 41 Figure 9. Caption on the next page D R A F T July 18, 2019, 4:11pm D R A F T X - 42 ULVROVA ET AL.: SUBDUCTION INITIATION AND CESSATION Figure 9. Cumulative distribution functions for distance between continents and points along subduction zones at di↵erent stages of their development for reconstructions from the Triassic to present (see supporting text and animations S2-S3). Each coloured line represents the distribution for a specific time in panels a, b, d and e. In panels c and f, the relatively small number of initiating and ceasing subduction segments are subdivided into broad time ranges encompassing the earlier and later stages of Pangea breakup. See text for further explanation. (a) CDF for random points falling within reconstructed extent of ocean basins; the grey background shows the envelope of these distributions based on random points, and is reproduced on the other panels for visual reference; (b) distance to continent for segments of active subduction zones for the kinematic reconstruction of Müller et al. [2016]; (c) distances to continent for initiating and ceasing subduction segments derived from Müller et al. [2016]; (d) distance to continent for remnants of past subduction mapped from seismic tomography [van der Meer et al., 2010, 2012]; (e) As (d), but only including ’primary’ subduction according to the definition of van der Meer et al. [2010]; (f) distance to continents for subduction zones in (e) at the beginning and end of their lifespans (assumed to approximate initiation and cessation) for the slab remnant reconstruction. D R A F T July 18, 2019, 4:11pm D R A F T