Oceanic crust recycling controlled by weakening at slab edges

Retreating subduction zones such as the Lesser Antilles, Gibraltar and Scotia have been migrating towards the Atlantic Ocean by cutting their way through the oceanic crust. This spontaneously retreating subduction is enabled by the development of faults at the edges of the slab, but the physical mechanisms controlling fault propagation and direction remain unknown. Here, using 3D numerical subduction models we show that oceanic lithosphere recycling is mainly controlled by the intensity of strain-induced weakening of fractures forming at the edges of the slab. Intense strain-induced weakening causes predominantly brittle fault propagation and slab narrowing until detachment. Without weakening, preponderantly ductile slab edge propagation occurs, which causes slab widening. This rheological control is not affected by the proximity of non-weakened passive continental margins. Natural examples suggest that slab edges follow convergent paths that could be controlled by fractures weakening due to deep water penetration into the oceanic lithosphere. Retreating subduction zones are enabled by the development of faults at the edges of the slab, but the physical mechanisms controlling fault propagation remain debated. Here, the authors find that oceanic crust recycling is controlled by weakening of fractures forming at the edges of slabs.

T he evolution of boundaries in the Earth's plates mosaic remains one of the key contemporary research questions 1-3 . It has recently been shown that the emergence of modernstyle subduction and its catalysts are closely tied to geological processes at the surface of the Earth 1 and that Archean subduction was likely of predominantly retreating character 1,2,4 . Modern retreating intra-oceanic subduction zones can efficiently consume ocean floor and bring material into the mantle by cutting their ways through lithospheric lid interiors (e.g. Scotia, Gibraltar and Caribbean subduction zones, all of which are or have been migrating into the Atlantic [5][6][7][8][9][10][11][12][13][14]. Slab retreat is enabled by fractures propagating at the edges of the slab 15 , at the so-called STEPfaults (Subduction-Transform Edge Propagator) 16,17 . The fracture trajectories and dynamics are thus crucial to determine how the slab migrates. Different scenarios can be envisaged ( Fig. 1): STEP-faults could either exhibit divergent fracturing paths, leading to the consumption of increasingly large portions of the oceanic floor; or parallel fracturing paths, consuming a band of oceanic floor of constant-width; or convergent fracturing paths, resulting in a narrowing of the retreating slab until its detachment. The aim of this study is to investigate which of these scenarios is typical on present-day Earth and if the dominant scenario of lid recycling was different back in geological times.
Fault propagation theory identifies three modes of fracturing 18 : opening (mode I), sliding (mode II) and out-of-plane tearing (mode III). The STEP-fault migration will be strongly controlled by the predominant fracturing mode at the slab edges, in turn controlling the fault migration characteristics and the resulting lid recycling pattern (narrowing, stable or widening). To determine which fracturing mode is favoured in nature and what is the resulting slab retreat trajectory, we performed 3D high-resolution visco-plastic thermomechanical numerical simulations 19 of spontaneously retreating oceanic subduction. Our models start with subduction initiation (Methods, Supplementary Fig. 1a, b) and evolve towards self-sustained slab retreat (Fig. 2).
Our numerical results suggest that both the deformation mechanism and the direction of STEP faults propagation are mainly controlled by the intensity of strain-induced weakening (SIW) of fractures forming at the edges of the slab. Intense SIW causes predominantly brittle STEP-fault propagation associated with the narrowing of the slab until its detachment. In contrast, without such fractures weakening, a preponderantly ductile slab edges propagation occurs that causes slab widening. Comparison of numerical results to existing natural examples suggests that in nature, STEPfaults follow convergent paths that could be controlled by intense fractures weakening due to deep water penetration into the relatively cold oceanic lithosphere of the present-day Earth.

Results
General slab evolution in oceanic domains. Our numerical experiments systematically demonstrate that STEP-faults can spontaneously propagate into the oceanic plate interior starting from pre-existing lithospheric-scale perturbations (Fig. 2, Supplementary Table 2, Supplementary movie 1). The geometries of tearing paths are mainly controlled by the intensity of SIW of brittle/plastic fractures forming inside the lithosphere (Methods). For models with an intense SIW, narrowing slabs (Fig. 2a) and converging tear paths (Fig. 2c, Supplementary Fig. 2a) controlled by out-of-plane (mode III) lithospheric-scale shear are characteristic. In these models, the slab eventually gets so narrow (~120-150 km) that it breaks off and subduction ceases ( Supplementary Fig. 2a). In contrast, widening slabs (Fig. 2b) and diverging tear paths (Fig. 2d, Supplementary Fig. 2b) controlled by the combination of out-ofplane shear (mode III) and extensional necking (mode I) of the lithosphere develop in models without SIW. In these models, the slab retreat slows with time ( Supplementary Fig. 2b), but subduction typically remains stable and does not stop within the model domain. We also observe the formation of linear topographic highs on the overriding plate (Fig. 2e, Supplementary Fig. 4) that are parallel to the STEP faults similarly to what is observed in nature (e.g. island chains in Fig. 1). In contrast, models without SIW (Fig. 2f) show the development of a wide straight arc with a deep (>9 km) trench that looks unlike the Scotia subduction zone or other modern narrow retreating subduction zones examples (cf. Figs. 1 and 2f).
Retreating slabs along continental margins. We further tested how the two retreating subduction scenarios can be affected by the presence of continental margins. Continental margins are found in proximity to the Caribbean, Scotia (Fig. 1) and also Gibraltar subduction zones, and represent lithospheric-scale heterogeneities [20][21][22] . These heterogeneities may potentially affect the plate tearing paths and deformation mechanisms. In particular, in the Caribbean, the Lesser Antilles slab has retreated towards the ocean between two continental blocks 23,24 (North and South America) with diverging continental margins. If the slab were to retreat further, STEP-faults could potentially follow these margins, thereby allowing the Caribbean subduction to widen and invade the Atlantic 5,25 . We investigated the potential for this future scenario with a simplified oceanic-continental initial plate configuration in which the slab retreats between two continents with equidistant non-weakened (i.e., without any preexisting weak zones) passive margins in the area where subduction is initiating, the distance between the margins increasing in the direction of slab retreat (Methods, Supplementary Fig. 1b).
Characteristic STEP-fault trajectories are not critically affected by the presence of diverging continental margins both for cases with (Fig. 3a, c) and without (Fig. 3b, d) SIW: the slab retreat trajectories are similar to what we observed in the full oceanic domain and the continental margins are not followed by the tearing paths ( Fig. 3) unless pre-existing rheologically weak zones are present along the margins 8,9 . This affirmation is valid both for faster and slower retreating slabs ( Supplementary  Fig. 3) and confirms what was suggested by Nijholt and Govers 17 . We notice however that the slab narrowing occurs at a lower rate when the slab retreat is slower and the overriding plate becomes respectively thicker (Supplementary Figs. 2 and 3, and Table 2, slower slabs are resulting from a higher mantle activation volume.).
SIW of brittle/plastic fractures (Methods section) has a strong influence on the subduction zone's retreat since it changes the predominant rheological mechanism of the lithospheric-scale deformation. An intense SIW favours deep penetration of brittle/ plastic faults into the lithosphere, whereas absence of the weakening increases the relative influence of ductile lithospheric deformation mechanisms (Methods section).  SIW (b, d). e, f: Topography of the retreating subduction zones for the models with (e) and without (f) SIW shown in a, c, and b, d, respectively. Note the similarity of topography in e to the Scotia region (Fig. 1).
topographic highs we notice in the track of the subducting slab can also be compared to the structures described in the North and South Scotia ridges 12,26 . Our models for a full oceanic domain show that the STEP-faults are best defined at shallow depth, between 60 and 80 km below the surface. This observation is in accordance with the focal depth of earthquakes measured at the edges of the South Sandwich trench 26,29 . In the Caribbean, seismic data analysis has been used to evidence a detachment of the oceanic Caribbean plate from the continental South American plate via STEP-fault propagation 30 , which is coherent with our model results. Moreover, a narrowing slab trajectory has been suggested by studies using wide angle seismics 14 .
All of our experiments involving SIW show a narrowing of the slab due to convergent STEP-faults, including when the lithosphere gets younger in the direction of retreat ( Supplementary  Fig. 4). The slower the retreat the more equidistantly the STEP faults propagate.
It is however hard to assess the exact influence of the South Sandwich fracture zone or small continental blocks around the Caribbean or Scotia (e.g. South Georgia) on the slab retreat. From our modelling results, the tendency of a retreating subduction zone to follow a fracture zone greatly depends on the rheology of the fracture zone: large weak discontinuities can be followed by STEP-faults 8,16 whereas small and less weak structures will be predominantly ignored.

Discussion
The impact of SIW intensity has strong implications for both modern and Archean Earth, and other terrestrial planets (e.g. Venus). Modern Earth is thought to have a colder and thus more rigid crust and mantle than Archean Earth [31][32][33] . Colder conditions favour brittle deformation 34 and deep water penetration along faults 8,20 , causing the development of weak hydrous minerals in the deforming lithosphere 35,36 , thereby creating favourable conditions for an intense SIW (Methods) and converging lid consumption patterns observed today, but this may not be valid during the onset of plate tectonics in the Archean.
Spontaneously retreating subduction zones that cut through the plate interior do not necessarily require modern-style plate tectonics. They could also develop within a single/stagnant lid as the result of other (e.g. plume-induced 2,37,38 and/or meteorite impact-induced 39,40 ) processes. This peculiar type of subduction may therefore control partial or complete recycling of the lid, including global or regional overturn events, on early Earth and other terrestrial planets 2,38,40,41 . If, on Modern Earth, converging STEP-faults and narrowing slabs seem to be favoured; with warmer and/or dry lithospheric conditions, where there is no significant brittle/plastic strain weakening due to limited or absent hydration along fractures, such as on early Earth or on Venus 42 , the simulations without SIW which favour ductile lithospheric deformation may be more representative. There, the retreating subduction zones should widen with time, ultimately leading to global lid recycling as suggested for Venus 42 . On early Earth, upon gradual cooling of the mantle, the propagation of STEP-faults during episodic retreating subduction may also have enabled the creation of new plate boundaries and contributed to the gradual emergence of a global plate mosaic that is crucial for establishing modern-style plate tectonics 42 .

Methods
Numerical approach. We model subduction initiation and slab retreat using the 3D thermo-mechanical code I3ELVIS 2 .
Mass, momentum and energy conservation are solved on a non-deforming Eulerian staggered grid and the physical properties of rocks (temperature, rock type, etc.) are advected via Lagrangian markers. Detailed code description is given in 19 .
Model rheology. The code accounts for visco-plastic deformation using an effective viscosity formulation for the material. The ductile deformation comprises diffusion and dislocation creep as described in Gerya et al. 2 . The ductile creep viscosity η ductile is calculated taking into account dislocation and diffusion creep as: where η diff , the effective viscosity for diffusion creep, is computed as and η disl , the effective viscosity for dislocation creep is computed as  where A D , E, V and n are respectively the pre-exponential factor, the activation energy and the activation volume; they are experimentally determined constants. R is the noble gas constant, P is the pressure, T the temperature, _ ε II the second invariant of the deviatoric strain rate tensor and σ cr the diffusion-dislocation transition stress (Supplementary Table 1). Partial melting, melt extraction and melt-induced weakening are modeled following Gerya et al. 2 .
Rheological layers are defined as follows (Supplementary Table 1). The ductile rheology of the oceanic crust is defined using the wet quartzite (upper crust) and anorthosite (lower crust) flow laws 43 , the continental crust ductile rheology is defined using the wet quartzite flow law 43 , the lithospheric and asthenospheric mantle rheology is defined based on dry olivine flow law 43 whereas the shear zone material is assigned a wet olivine rheology 43 . The boundary between the lithospheric and asthenospheric mantle is purely thermal and is set at 1273 K.
The ductile rheology is combined with brittle/plastic rheology based on a Drucker-Prager yielding criterion: where C eff is the material cohesion, γ eff is effective friction coefficient, λ melt the long-term melt-induced weakening factor 2 , _ ε II is the second invariant of the strain rate tensor, P melt is the melt pressure and P is the total pressure.
It should be noted that our rheological model neglects the elastic deformation. As the dominant large strain associated with the STEP faults propagation is viscoplastic, this simplification should not affect the model behaviour in a critical manner.
Strain-induced weakening. One of the critical factors in our numerical experiments is the SIW of brittle fractures forming in the lithosphere 35,44 . This weakening is related, in particular, to deep water percolation along deforming oceanic fault zones leading to their intense hydration 32,35 which strongly decreases the strength of fractured fault rocks 35,36,45 . In our numerical model, we use a simplified strain weakening model that prescribes a linear decrease of the internal friction coefficient with increasing cumulative brittle/plastic strain ε c 44,46 .
when ε 1 <ε c <ε 2 Where γ eff is the internal friction coefficient, ε 1 and ε 2 are two threshold strains, b 0 and b 1 are the initial and final internal friction coefficient values, _ ε ij plastic ð Þ is the brittle/plastic strain rate tensor 47 .
Numerical model setup. We model two different settings using regional-scale models: a purely oceanic domain and a mixed oceanic/continental domain. The models are 1212 × 1212 km wide in both horizontal directions and 392 km deep. We use a grid resolution of 2 and 3 km in the vertical and horizontal directions, respectively, which is sufficient to resolve plate boundaries and the lithological and rheological structures of both oceanic and continental domains. In both cases, we model slab evolution from subduction initiation to slab retreat. We initiate subduction using a strong lateral lithospheric age contrast (and hence density contrast) in the oceanic lithosphere. The lateral boundaries of the box are mechanically defined as free slip conditions, and no heat flux across these boundaries is allowed. The top boundary is free slip, and fixed temperature. No external forcing is applied, and the model is evolving entirely due to its internal dynamics. A new upper crust is gradually forming by cooling along with slab retreat. In addition, the uppermost part of the model interior (20 km) is composed of a sticky air layer 38 , which is a low density, low viscosity layer that emulates a free surface condition on top of the lithosphere. The lower boundary of the model is set as an open boundary both mechanically and for the temperature. Temperature is set at 273 K at the surface of the lithosphere and in the sticky air layer. The temperature in the lithosphere is calculated according to a cooling half-space model, where we set an initial age of 40 Ma for the lithosphere at the beginning of the simulation. In the asthenosphere, a positive gradient of 0.5 K/km is defined. Our limited model domain -which is virtually extended by the open boundarythus does not implicitly include the 660km upper-lower mantle transition that can potentially affect the retreating subduction dynamics, but is not subjected to undesirable interaction with a shallow (392 km) model boundary either. Since in most of our numerical experiments the modeled amount of slab subduction is <600 km this simplification is an approvable compromise. Due to the computational restrictions, our study focuses on highresolution 3D modeling of relatively narrow retreating subduction zones limited by STEP faults that evolve in the uppermost 660 km.
We investigated numerically two model settings (with and without continents) shown in Supplementary Fig. 1a, b. In both model setups, we start the simulation from a homogeneous 40-Ma old lithosphere throughout, with the exception of a rectangular 200 × 460 km domain of younger lithospheric age at the left model boundary (1-10 ka, Supplementary Fig. 1b, circled in red). No slab is prescribed. In this "lithospheric window" the cooling age of the lithosphere decreases from 10 to 1 ka from the border of the box towards the area where we want to trigger slab formation (the 1-ka part of the young window being located next to the area where the slab forms, Supplementary Fig. 1a, b). In order to enable the retreating subduction initiation, we prescribe a "H" shaped weak (γ eff = 0) shear zone pattern at the lithospheric window that extends by~200 km into the old lithosphere ( Supplementary Fig. 1a). This highly simplified model geometry is used to prescribe both the initial slab width and the direction of free subduction retreat without a need to impede internal dynamics by forcing convergence along the boundaries, nor setting an initial slab dip, the model is hence fully dynamic.
In the mixed oceanic/continental simulation, we added two continental margins at the lateral model boundaries (Supplementary Fig. 1b). These are equidistant around the area where the slab forms and the distance between the margins increases in the direction of slab retreat. The continents are located~500 km apart of each other and their margins stay parallel for~400 km (Supplementary Fig. 1b).

Data availability
Due to the very large size of the data, the model results are available from the corresponding author upon request.