Where does subduction initiate and cease ? A global 1 scale perspective

Abstract The thermo-mechanical evolution of the Earth's mantle is largely controlled by the dynamics of subduction zones, which connect the surface tectonic plates with the interior. However, little is known about the systematics of where subduction initiates and ceases within the framework of global plate motions and evolving continental configurations. Here, we investigate where new subduction zones preferentially form, and where they endure and cease using statistical analysis of large-scale simulations of mantle convection that feature self-consistent plate–like lithospheric behaviour and continental drift in the spherical annulus geometry. We juxtapose the results of numerical modelling with subduction histories retrieved from plate tectonic reconstruction models and from seismic tomography. Numerical models show that subduction initiation is largely controlled by the strength of the lithosphere and by the length of continental margins (for 2D models, the number of continental margins). Strong lithosphere favours subduction inception in the vicinity of the continents while for weak lithosphere the distribution of subduction initiation follows a random process distribution. Reconstructions suggest that subduction initiation and cessation on Earth is also not randomly distributed within the oceans, and more subduction zones cease in the vicinity of continental margins compared to subduction initiation. Our model results also suggest that intra-oceanic subduction initiation is more prevalent during times of supercontinent assembly (e.g. Pangea) compared to more recent continental dispersal, consistent with recent interpretations of relict slabs in seismic tomography.

where 0 is the surface yield stress, d is the depth and 0 Y is the yield stress depth 124 derivative. If the stress reaches Y , we calculate the e↵ective viscosity ⌘ e↵ on the grid billions of years.

170
The system is driven by heating from radioactive elements and from heat conducted is no-slip above the air layer.

177
The convective vigour of the system is measured by the Rayleigh number Ra where ⇢ 0 is the reference density, g the gravitational acceleration, ↵ the thermal expan-178 sivity, T is the superadiabatic temperature drop across the mantle with the thickness 179 D,  is the thermal di↵usivity and ⌘ 0 is the reference viscosity. For all experiments, we 180 keep Ra = 10 6 (calculated with ⌘ 0 = 6 ⇥ 10 22 Pa s). This is 10 100⇥ lower than the

Model analysis
To detect subduction inception and follow its evolution, we analyze the divergence of 209 the surface velocity. As soon as a negative peak of the surface velocity divergence appears

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. The

219
long duration of each model allows us to collect up to several hundreds of subduction 220 initiations, typically several tens of them.

221
To investigate systematic patterns within the distributions retrieved from the models,

Results
The organization of the system dictates its dynamic evolution: sinking slabs drive that there is a strong subadiabatic gradient (Figure 1c,f,i). This is partly because the 249 system is internally heated but more importantly negative gradients arise due to pressure 250 dependent viscosity [Christensen, 1985].

251
Subduction zones (and in this study we refer to subduction zones as convective down-252 welling currents in the numerical simulations) next to continental margins are one-sided 253 while a distinctive asymmetry is observed for intra-oceanic subduction zones (cf. Figure 1 254 and animation S1 in the supporting material). The degree of asymmetry in the latter case 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 as 295 intra-oceanic reside in the oceans where they either merge with another subduction zone 296 or migrate in the oceans before reaching a continental margin where they endure for some 297 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 ( Figure   299 2). Rarely, a subduction zone is formed and ceases as intra-oceanic without colliding with 300 another convergence zone ( Figure 5). In the weak lithosphere case, termination is random  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 continent 320 to 67% for three continents (cf. Figure 4b and 6). Two fundamental results are consistent 321 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 more 323 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.