Dynamic slab segmentation due to brittle–ductile damage in the outer rise

Subduction is the major plate driving force, and the strength of the subducting plate controls many aspects of the thermochemical evolution of Earth. Each subducting plate experiences intense normal faulting1–9 during bending that accommodates the transition from horizontal to downwards motion at the outer rise at trenches. Here we investigate the consequences of this bending-induced plate damage using numerical subduction models in which both brittle and ductile deformation, including grain damage, are tracked and coupled self-consistently. Pervasive slab weakening and pronounced segmentation can occur at the outer-rise region owing to the strong feedback between brittle and ductile damage localization. This slab-damage phenomenon explains the subduction dichotomy of strong plates and weak slabs10, the development of large-offset normal faults6,7 near trenches, the occurrence of segmented seismic velocity anomalies11 and distinct interfaces imaged within subducted slabs12,13, and the appearance of deep, localized intraplate areas of reduced effective viscosity14 observed at trenches. Furthermore, brittle–viscously damaged slabs show a tendency for detachment at elevated mantle temperatures. Given Earth’s planetary cooling history15, this implies that intermittent subduction with frequent slab break-off episodes16 may have been characteristic for Earth until more recent times than previously suggested17. Numerical subduction models used to determine the consequences of bending-induced plate damage show that slab weakening and segmentation can occur at the outer-rise region of the subducting plate.

Subduction is the major plate driving force, and the strength of the subducting plate controls many aspects of the thermochemical evolution of Earth. Each subducting plate experiences intense normal faulting [1][2][3][4][5][6][7][8][9] during bending that accommodates the transition from horizontal to downwards motion at the outer rise at trenches. Here we investigate the consequences of this bending-induced plate damage using numerical subduction models in which both brittle and ductile deformation, including grain damage, are tracked and coupled self-consistently. Pervasive slab weakening and pronounced segmentation can occur at the outer-rise region owing to the strong feedback between brittle and ductile damage localization. This slab-damage phenomenon explains the subduction dichotomy of strong plates and weak slabs 10 , the development of large-offset normal faults 6,7 near trenches, the occurrence of segmented seismic velocity anomalies 11 and distinct interfaces imaged within subducted slabs 12,13 , and the appearance of deep, localized intraplate areas of reduced effective viscosity 14 observed at trenches. Furthermore, brittle-viscously damaged slabs show a tendency for detachment at elevated mantle temperatures. Given Earth's planetary cooling history 15 , this implies that intermittent subduction with frequent slab break-off episodes 16 may have been characteristic for Earth until more recent times than previously suggested 17 .
The subduction of negatively buoyant oceanic lithosphere is a key driver of terrestrial tectonics. Subduction results from buoyancy forces that bend and pull the lithosphere into the interior of the Earth's mantle where mechanical properties of subducted lithospheric slabs are strongly modified by various physical-chemical processes 10,18,19 . One of the large-scale effects of this modification is the pronounced mechanical dichotomy of stronger lithospheric plates at the surface and weaker slabs in Earth's interior, which has been proposed on the basis of various geological-geophysical observations combined with some analytical and numerical modelling 10,[19][20][21][22][23] . Whereas strong plates at the surface are a prerequisite for terrestrial-style, one-sided subduction 10,24 , weakened lithospheric subducted slabs, which bend, stretch, pile, segment or lie flat at the top of the lower mantle, are needed to reproduce the spectrum of slab morphologies 25,26 , observed bending curvatures 23 , lateral changes in Earth's geoid 18,27 and variations of intraslab seismicity with depth 22 . Slab weakening is in apparent contradiction with experimentally calibrated, thermally activated rheological laws for the lithospheric mantle 28 , which predict high effective viscosity (>10 24 Pa s) of subducted slabs in the upper mantle 18 , implying a slab/ mantle viscosity contrast of ≫1,000. By contrast, a number of natural observables suggest that this contrast should be much lower (around 100-500) 20,22,23,25,29 , for example, to reconcile the dynamics with seismic tomography models and slab seismicity distribution with depth that indicate strong bending, stretching and even disruption of subducted slabs in the upper mantle 22,23,25,26,30 .
The apparent rheological paradox of strong plates and weak slabs 10,18,25 can be resolved by assuming some additional lithospheric weakening processes that are intrinsically related to the transition from horizontal plate motion to its subduction into the mantle. In this respect, plate bending at the outer rise is a primary candidate for changing the mechanical properties of the lithosphere 4,19 . Each subducting plate, irrespective of its age and earlier history, unavoidably experiences a transition from its horizontal to vertical motion through bending when passing through the outer-rise region that is present at every subduction trench. Plate bending is not fully elastic and is associated with a number of irreversible physical-chemical processes that can markedly weaken the mechanical properties of the subducting plate 1,4,[6][7][8]19 .
It is well understood that the colder, brittle top region of subducting plates is pervasively damaged by extensional outer-rise normal faulting 1,4,7,9,19 , possibly associated with downwards water penetration and mantle lithosphere serpentinization [2][3][4][5]8 . By contrast, the deeper and warmer portions of the plate deform in compression by both viscous creep 14 and compressional faulting 31 . This deformation style indicates the existence of a barrier to the penetration of water and deeper chemical alteration of the plate. As the result, the deeper portion of the plate can be predominantly weakened by ductile damage processes such as grain-size reduction assisted by Zener pinning [32][33][34][35][36] . Relating these brittle and ductile damage processes to the plate bending at the outer rise can advance our understanding of whether and how the pervasive weakening of subducting plates may occur in nature and what Article consequences that deformation memory will have for the dynamics and stability of subduction.
Here we investigate the consequences of bending-induced plate damage by using two-dimensional numerical thermomechanical subduction models in which both brittle-plastic and ductile deformation as well as grain-size evolution are coupled self-consistently (Methods). Our results indicate that brittle-ductile slab weakening and plate segmentation should occur at the outer-rise region owing to the strong feedback between brittle and ductile damage localization, and that this process can explain a range of observations from well imaged subduction zones. Figure 1 shows an example of a typical model evolution resulting in subducting slab damage and segmentation. The numerical experiment starts from subduction initiation at a transform fault 24,37 , which separates two oceanic plates of different ages, leading to gravitational instability. The older plate starts to subduct in a retreating manner whereas the younger, now overriding, plate is subjected to extension and horizontal motion towards the retreating trench. After the initial period of retreat associated with gradual steepening of the slab, the slab angle stabilizes and the advancing horizontal motion of the subducting plate begins. This mode of subduction continues until the slab reaches the mantle transition zone and starts to flatten due to the negative Clapeyron slope of the spinel-perovskite phase transition (Methods), which initiates a new episode of trench retreat. The viscosity of the subducting plate that results from this dynamically self-consistent subduction scenario shows a pattern of 150-200-km-wide segments separated by narrow low-viscosity zones. This damage pattern forms as the result of linked, localized brittle-plastic deformation and grain-size reduction at the outer rise (Fig. 2a). As a consequence, the subducting slab deforms easily in a chain-like fashion in response to its interaction with the mantle transition zone. The resulting slab morphologies are markedly different from those of viscoplastic slabs without damage 21,25 and reflect the influence of the strong coupling between the brittle and ductile strain-localization mechanisms during plate bending.

Slab segmentation at the outer rise
Model sensitivity studies show that the style of slab deformation depends strongly on both the ductile damage and the strain-induced weakening of faults as well as on the age of the subducting plate and the mantle potential temperature. In particular, deactivation of ductile damage (grain-size reduction) and/or strain-induced weakening produces smoothly bent slabs (compare Extended Data Fig. 2a, e and Extended Data Figs. 2b, f, 2c, g, 2d, h) comparable to previous subduction models 16,21,24,25 . A lack of ductile damage also makes subduction initiation more difficult (compare Extended Data Fig. 3a, e and Extended Data Fig. 3c, g) owing to the increased bending resistance of the subducting plate, which should become older and denser to allow for spontaneous subduction initiation (compare Extended Data Fig. 3a, e and Extended Data Fig. 4c, g). The combined effects of fault weakening and ductile damage on slab segmentation are notably distinct from the individual forms of weakening. In the absence of normal fault weakening, slab segments are less pronounced and grain-size reduction inside the subducting plate is distributed more evenly owing to the weaker feedback between faulting and grain-size reduction (compare Extended Data Fig. 2a, e and Extended Data Figs. 2c, g, 4c, g, 4d, h). This, in turn, causes more distributed rather than strongly localized, segmented weakening of slabs. A similar effect is achieved by decreasing the rate of fault-weakening with strain (compare Extended Data Fig. 3b, f and Extended Data Fig. 3c, g). In the absence of ductile damage, slab segments remain but the displacement along individual normal faults is reduced (Extended Data Fig. 2b, f). As a result, the number of segments increases together with a decrease in the characteristic segment width (compare Extended Data Fig. 2a, e and Extended Data Fig. 2b, f). The total amount of deformation in each segment also decreases.
Slab segmentation is thus primarily driven by normal fault weakening whereas the ductile damage makes this process more intense, localized and laterally extensive. Ductile creep of mature faults has no influence on segmentation and mainly affects the deformation of segmented slabs in the mantle: weaker (serpentine type) rheology facilitates slab bending and break-off whereas stronger (dry olivine type) rheology produces less deformed slabs (Extended Data Table 2).
The age of the subducting plate also controls the characteristic width of slab segments: older and thus thicker subducting plates show wider slab segments than those of younger plates (compare Extended Data Fig. 3c, g and Extended Data Fig. 3d, h). An increase in the mantle potential temperature, by contrast, promotes slab segmentation and bending by steepening of the slab angles owing to the increased negative slab buoyancy and the reduced amount of viscous resistance of the asthenosphere to slab penetration (compare Fig. 1 and Extended Data Figs. 2a, 4a). This also leads to notable acceleration of subduction and more frequent slab break-off 16 (compare Fig. 1 and Extended Data Figs. 2a, e, 4a, e).
Many outer-rise fault systems are orientated along pre-existing faulting structures within the oceanic plate 38,39 . We therefore tested the influence of pre-existing faults by running several models in which these faults have been initially prescribed with different spacing (5 km, 10 km and 20 km) and polarity (towards/outwards the trench) (Extended Data Fig. 5). Similar to the reference model, slab segmentation also developed in the models with pre-existing faults, which have a minor effect on the characteristic segment width (Extended Data Fig. 5). As expected, pre-existing faults are easily reactivated by the plate bending and control the predominant location and polarity of the much deeper outer-rise normal faults (Extended Data Fig. 5).

Weakening of outer-rise faults
Our numerical models suggest that strain weakening of outer-rise normal faults is a key process controlling subducting slab segmentation, and is thus discussed further. Strain weakening of faults is a common assumption of geodynamic models 40 and is crucial for reproducing a number of strain-localization phenomena in both oceanic and continental lithosphere such as large-offset normal faults 40 , oceanic transform faults 41 , and oceanic and continental core complexes 42 . The physics of this process is incompletely understood and may include (but is not limited to) pressurized fluid percolation 24 , growth of hydrous minerals 43 , structural softening 44 , shear heating 45 , coseismic weakening and grain-size reduction 46 , and intergranular cavitation 47 . Intense hydration of outer-rise normal faults has been suggested on the basis of seismological, geophysical and theoretical arguments [2][3][4][5]8,31 .
In particular, numerical models 4 find that the dynamic pressure associated with plate bending may be large enough to overcome the confining lithostatic pressure and cause downwards water suction along outer-rise normal faults. By contrast, analytical models 9 suggests that such large dynamic pressure cannot be achieved, and lowered seismic velocities within the oceanic lithosphere under the outer rise may instead be explained by thermal cracking. Putting aside such uncertainties regarding the extent of outer-rise fault hydration, it is clear that these structures systematically reveal lowered friction coefficients (≤0.3) 31,48 compared with dry oceanic lithosphere (0.6-0.85) 49 , and only such overall weakening is required for our simplified strain-weakening models (Methods).

Natural evidence for slab segmentation
Our model predicts that, owing to the interplay between the fault weakening and grain-size reduction in the outer rise, subducting slabs are pervasively damaged and can deform easily in a chain-like fashion in the mantle. This new, emergent mechanism provides a self-consistent explanation for the previously proposed subduction dichotomy of strong plates and weak slabs 10,[19][20][21][22][23]30 , which controls the observed slab morphologies 25,26 in the deep mantle. One of the testable consequences of slab segmentation is the punctuated development of large-offset (throw ≥400 m) normal faults located above the localized regions of intense grain-size reduction (Fig. 2a, Extended Data Fig. 6). Mirroring the deep slab morphology, this more localized faulting pattern contrasts with broadly distributed, moderate-offset (throw ≤300 m) normal faults as produced by models without grain-size reduction 4 and/ or without intense fault weakening (Fig. 3a, Extended Data Fig. 2b-d, f-h). We predict large-offset normal faults to be transient structures that episodically form close (<30 km) to the trench and subsequently subduct under the forearc (Extended Data Fig. 6). These structures should thus only be observed at those trenches in nature where the slab-segmentation process is in its mature stage. Fault throw patterns may then vary between different subduction zones and even along the same trench 1 owing to lateral changes in slab-segmentation maturity presumably caused by regional variations in subducting-plate velocity, crustal thickness and cooling age 1 . This behaviour is indeed what is observed. Large-offset normal faults are only found very close (<30 km) to some subduction trenches or under their frontal prisms, such as in the Japan Trench 6,7,50,51 (Fig. 2c, d) and in the northwest portion of the Middle America Trench 1,51 , but are absent in other trenches, for example, in most of the Mariana Trench 51 and in the southeast portion of Middle America Trench 1,51 (Fig. 3c, d). The deepest part of the slab in the bending region in Japan has likewise been inferred to display a localized, intraplate area of reduced effective viscosity based on post-seismic deformation 14 . Such deep intraplate weakening is predicted by our model, where this zone corresponds to areas of grain-size reduction that develop in either a strongly localized (mature segment, Fig. 2a) or a more widely distributed (immature segment, Fig. 3a) manner. We analysed (Methods) large-scale fault throw patterns for the models with mature and immature slab segmentation and found significant differences (compare Fig. 2b and Fig. 3b). The model stages of immature slab segmentation show more broadly distributed grain-size reduction and moderate-offset faults with throw that gradually increases towards the trench and then show a plateau at 150-250 m within 40-50 km of Only faults with throw >20 m are considered. c, The seismic profile with the forearc slope morphology, the regional basal unconformity, and subducting Pacific Plate and normal faults indicated 50 . The largest offset is characteristic for the normal fault buried under the frontal prism. d, The distribution of fault throws at the Japan Trench measured by ref. 7 . The data from 13 seismic profiles 7 across the trench are shown by different symbols. The zones of development of large-offset (throw ≥400 m) and moderate-offset (throw ≤300 m) normal faults in b and d are highlighted by pink and blue colours, respectively.
Article the trench (Fig. 3b). By contrast, mature slab segmentation shows an abrupt increase and high variability of throws in a narrow zone close to the trench and under the frontal prism (Fig. 2b). These two contrasting cases find their analogy in nature: whereas an abrupt appearance of large-offset normal faults (throw 400-550 m) is documented within the 25-km-wide zone at the Japan Trench (Fig. 2c, d), much more uniformly distributed moderate-offset normal faults (throw <200 m) are found at the southeast part of the Middle America Trench (Fig. 3c, d). Although slab segmentation may explain the development of the contrasting fault throw patterns near trenches, the absence of (exposed) large-offset normal faults near some trenches does not imply the absence of slab segmentation as those are transient features (Fig. 2a, c, Extended Data Fig. 6). Broadly variable outer-rise normal faulting patterns worldwide show no correlation with major subduction parameters such as oceanic plate age, subduction velocity and slab pull 52 . Our model provides a self-consistent explanation for this lack of correlation as the spatial distribution and the offsets of outer-rise faults change with time and are mainly dependent on the local maturity of slab segmentation process rather than on the subduction parameters. Further support for the occurrence of slab segmentation in nature can also be found in high-resolution seismic tomography (Fig. 4, Extended Data Figs. 7, 8). Recent full-waveform imaging yields tomographic models that show unprecedented resolution, such as in the case of Japan 11 . Both the seismically active part of the slab transiting the upper mantle, and the flat-lying, apparently stalled, part of the slab show segmentation in terms of their seismic velocity. The latter is sensitive to grain size 53 and, for Japan, the inferred along-slab positions of the segments' boundaries (Methods) cluster into seven trench-parallel zones of lowered seismic velocity within the slab (Fig. 4d), which bound six individual, 200-425-km-wide segments. The segment width variation along the slab is comparable to our numerical modelling predictions (compare Fig. 4a, b and 4c, d, Extended Data Figs. 3d, h, 4c, g). Compared with our reference model, the wider slab segments inferred for Japan may be due to older seafloor ages (compare Extended Data Fig. 3c, g and Figs. 3d, h, 4c, g). Depth-dependent sensitivity of seismic waves to temperature and imaging imperfections complicate further interpretation, and other explanations for regional slab disruption have been discussed (for example, ref. 54 ). However, our grain-size evolution models provide a general explanation that does not rely on a specific convective setting.
Slab segmentation also allows deeper penetration of normal faults into the slab (up to 35-40 km, compare Fig. 2a and Fig. 3a). Global slab seismicity in the outer rise is often characterized by the presence of deep (20-50 km) compressional, thrust-faulting earthquakes. These earthquakes often occur below the elastic core of the subducting plate at temperatures up to 600 °C (ref. 31 ), that is, close to the brittle-ductile transition in the mantle lithosphere. The flow-to-friction transition near the base of the seismogenic zone may be characterized by a runaway transition from dislocation and diffusion creep to dilatant deformation, involving incompletely accommodated grain-boundary sliding 47 . Localization of the compressional ductile deformation observed in our numerical experiments (Fig. 2a) and controlled by grain-size reduction may thus create more favourable conditions for deep thrust-faulting; however, is not an ultimate cause for it. The majority of deep compressive earthquakes occur very close (<30 km) to trenches 31 where normal faults penetrate deeper into the lithosphere (Figs. 2a, 3a).
Lastly, the formation of a grain-size reduction zone in the bottom part of the slab may contribute to the development of distinct intraslab interfaces that are imaged within both relatively shallow descending slabs 12 (Fig. 5a) and within stagnant slabs in the mantle transition zone 13 (Fig. 5c). The upper intraslab interface with positive polarity (Fig. 5a, c) is interpreted as the subducted oceanic Moho 12,13 . By contrast, the lower interface (Fig. 5a, c) is characterized by the large shear-wave velocity reduction and probably requires the presence of melts at the bottom of the lithosphere 12,13 . This appears consistent with our models as the intense grain-size reduction zone that forms in the bottom of the slab (Fig. 5b, d) should both reduce shear wave velocity 53 and decrease permeability, and hence act as a barrier for percolation of melt 55 , thus enabling its accumulation at the bottom of subducted lithosphere.

Occurrence of intermittent subduction
Besides explaining various subduction zone features, brittle-ductile slab segmentation may have further consequences, including for the stability of the modern style of subduction. In the Archean eon, higher mantle potential temperatures 17 may have significantly reduced both plate strength and mantle resistance to slab penetration, and segmentation may have then induced frequent slab break-off. This could have caused punctuated, episodic (intermittent) subduction 16,17 . Indeed, our numerical experiments performed even at modestly higher temperatures of 100-150 K compared with the present-day show a strong tendency for slab disruption (Extended Data Fig. 2a, e) and break-off

Article
(Extended Data Fig. 4a, e) in models with grain-size reduction compared with those without (Extended Data Figs. 2b,f, 4b,f). This may imply a sensitive dependence of subduction on mantle temperature, and thus possibly a larger role for an intermittent style of subduction in plate tectonics during our planet's cooling 15 compared with what is expected from mantle convection with memory-free rheologies 17 . More generally, our models show that even oceanic lithosphere may be pervasively affected by deformation history, and the associated damage memory affects not only surface deformation at trenches but also how slabs deform and stir their surroundings as they descend into the lower mantle.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-03937-x.

Modelling approach
The thermomechanical two-dimensional numerical code I2VIS is used for the modelling of subduction initiation. It is based on a combination of a finite-difference method, applied on a staggered Eulerian grid, and a marker-in-cell technique 56,57 . The momentum, mass and energy conservation equations are solved in an Eulerian frame, and physical properties are transported by a Lagrangian. Non-Newtonian, visco-plastic rheologies and variable thermal conductivity are used in the model [58][59][60][61][62][63] (Extended Data Table 1), which accounts for major phase transitions in the oceanic crust and mantle as well as adiabatic, radiogenic and frictional internal heating sources. The full details of this method, allowing for its reproduction, are provided elsewhere 56,57 . Numerical C codes used for the modelling are provided with the paper.

Numerical model design
The initial model setup (Extended Data Fig. 1) corresponds to the one used for spontaneous subduction initiation at an oceanic transform fault 24,37 . The computational domain is equivalent to 3,000 × 3,000 km (Extended Data Fig. 1) and is resolved with an irregular rectangular grid of 1,261 × 511 nodes in the horizontal and vertical directions, respectively, and contains 19 million randomly distributed markers. The irregular grid resolution varies from 10 × 10 km at the model boundaries to 1 × 1 km in 1,000-km-wide and 200-km-deep subduction zone and outer-rise faulting area. All sides of the model have free-slip mechanical boundary conditions. The free-surface boundary condition atop the crust is implemented by using a 12-km-thick 'sticky' air/water layer 64 with low density (1 kg m −3 above 9 km, 1,000 kg m −3 below 9 km) and viscosity (10 17 Pa s). The initial thermal structure and thickness of the plate (Extended Data Fig. 1) is defined by prescribing a laterally uniform cooling age and respective geotherm 60 with 273 K at the surface and a mantle potential temperature of 1,523-1,823 K, varied in different experiments (Extended Data Table 2). We explored two types of model setup with different initial conditions. In models with free subducting plate, within 500 km at the right model boundary, the subducting-plate age gradually decreases towards 1,000 yr, which corresponds to a weak mid-ocean ridge located at the boundary (Extended Data Fig. 1a).
In models with subducting plate attached to the right model boundary, the subducting-plate age remains unchanged towards the boundary (Extended Data Fig. 1b). An adiabatic gradient of 0.5 K km −1 is initially prescribed in the asthenospheric mantle (Extended Data Fig. 1). Within 500 km at the lower boundary, the temperature increases linearly by 744 K to mimic the hot boundary layer at the core-mantle boundary. Temperature-dependent thermal conductivity is used for the mantle and the crust (Extended Data Table 1). The thermal boundary conditions are 273 K at the top, 3,567-3,717 K (depending on the mantle potential temperature) at the bottom and zero heat flux on two other sides of the model. To ensure efficient heat transfer from the surface of the crust, the temperature of the 'sticky' air/water is kept constant at 273 K. Gravitational acceleration of 9.81 m s −2 has been used in the model. The surface of the lithosphere evolves by erosion and sedimentation according to the following Eulerian transport equation 65 when z > 9 km. The maximal surface slope for the accumulated sedimentary prism is limited by 17°. Surface processes, however, have a relatively minor role in subduction dynamics and slab morphology in our numerical experiments as follows from test runs without surface processes (Extended Data Table 2).

Density model
We use the extended Boussinesq approximation with the incompressible continuity equation and variable density in the momentum and energy conservation equations. The density of rocks varies with pressure (P) and temperature (T) according to the equation where ρ 0 is the standard density at P 0 = 1 MPa and T 0 = 298 K, and α = 2 × 10 −5 Κ −1 and β = 4.5 × 10 −12 Pa −1 are the coefficients of thermal expansion and compressibility, respectively (Extended Data Table 1). Our models take into account the phase transformations of olivine into wadsleite and ringwoodite 67 and into bridgmanite in the mantle 68 . Eclogitization of subducted basaltic and gabbroic crust is taken into account by linearly increasing the density of the crust with pressure from 0% to 16% in the P-T region between the experimentally determined garnet-in and plagioclase-out phase transitions in basalt 69 . Stishovite and perovskite 68 transitions in the crust are also taken into account for density changes. The physical parameters for each experiment are presented in Extended Data Table 2.

Viscoplastic rheological model
The viscous and brittle (plastic) properties (Extended Data Table 1) are implemented via evaluation of the effective viscosity of the material. For the ductile rheology, the contributions from different flow laws such as dislocation and diffusion creep are taken into account by composite rheology for η ductile η η η where η diff and η disl are effective viscosities for diffusion and dislocation creep, respectively. For the crust, constant grain size is assumed and η diff and η disl are computed as is the square root of the second invariant of the strain rate tensor, σ cr is the assumed diffusion-dislocation transition stress, and A, E, V and n are the experimentally determined pre-exponential factor, activation energy, activation volume and stress exponent of the viscous creep, respectively (Extended Data Table 1).
For the mantle, the ductile creep model also takes into account grain-size reduction and growth processes assisted by Zener pinning, and η diff and η disl are calculated according to refs. 32,33,[70][71][72] . The rheology follows a composite law as in equation (1), wherein where h is the mean grain size and m is the grain-size exponent. The interplay between diffusion and dislocation creep is controlled by a grain-size-evolution equation dependent on the mechanical work and temperature. The grain-size evolution model relies on several assumptions. (1) The mantle peridotite is assumed to be composed of two well mixed phases: olivine and pyroxene with a fixed volume fraction of 60% and 40%, respectively. These phases are considered to have the same density and rheology. (2) In both phases, the relative motion is considered to be negligible and therefore their velocity is the same. (3) It is assumed that the grain-size distribution is close to a self-similar log-normal distribution. Therefore, it always retains the same shape and its mean variance and amplitude are fully characterized by a unique grain size. We make the further assumption that the system is in a state known as the pinned-state limit 32,72 wherein the grain-size evolution is controlled by the pinning of phases by each other (that is, Zener pinning is dominant) 32 . In these conditions, the grain size is controlled by the roughness r of the interface between the two phases. A relation between the mean grain size h (sufficient to fully describe the system) and the roughness r is given by h = r h g , where h g ≈ π/2 for the phase volume fraction in our model 72 . The roughness evolution is described by the following equations 32,33,70,73 r t ηG qr f r γ η I 0

2.9
where G I is the interface coarsening, G g is the grain growth rate, G fac = 100 is the grain growth rate factor, q = 4 is the roughness coarsening exponent, p = 2 is the grain-size coarsening exponent, γ I is the surface tension, A g = 2 × 10 (4-6p) the is pre-exponential factor, E g = 3 × 10 5 is the grain-growth activation energy, V g = V diff is the grain-growth activation volume, f I is the fraction of mechanical work Ψ converted to interface damage resulting in grain-size reduction; f 0 = 0.001 is the interface damage at 1,000 K, and η = 3φ ol φ px is interface area density depending on the volume fractions of olivine (φ ol = 0.6) and pyroxene (φ px = 0.4) in the mantle. The ductile rheology is combined with a brittle (plastic) rheology to yield an effective viscous-plastic rheology using the following upper limit for the ductile viscosityη ductile II μ μ γμ = − γ 0 for γ ≤ γ 0 and μ = μ 1 for γ > γ 0 , where μ is the internal friction coefficient, (μ 0 and μ 1 are the initial and final internal friction coefficient, respectively, Extended Data Table 1), μ γ = (μ 0 − μ 1 )/γ 0 is the rate of faults weakening with integrated plastic strain γ (γ 0 is the upper strain limit for the fracture-related weakening), C is the rock compressive strength at P = 0 (Extended Data Table 1), t is time and ε ij(plastic) is the plastic strain rate tensor. It is also assumed that the mantle inside outer-rise normal faults that reached the upper strain limit (γ 0 ) is serpentinized and has the respective rheology (Extended Data  Table 1).

Computing fault throw from numerical models
Fault throw was evaluated on the basis of integrated plastic strain γ stored on markers. To have a conservative estimate of the plastic strain related to normal faulting, we first defined a reference surface (white solid line in Figs. 2a, 3a) located inside the rheologically strong gabbroic layer at 3.5-km vertical distance above the oceanic Moho. Then, strain on markers found within 250-m vertical distance around this reference surface has been averaged at the nearest vertical grid lines intercepting the reference surface. Faults were then identified as the local maxima of the plastic strain γ, for which an affective fault throw magnitude M was estimated as where D = 1 km is a characteristic fault width in our models (as faults always localize within one grid cell) and α = 56° is the characteristic normal fault dip angle in the models. We used a throw of 20 m as the sensitivity threshold for this numerical approach. The Matlab code throw.m used for this calculation is provided with the paper.

Defining of slab-segment width from seismic tomography
We analysed five regularly spaced along-dip seismic tomography profiles of the Japan slab 11 (Extended Data Fig. 4) for both v p (pressure wave) and v s (shear wave) velocity variations (Extended Data Fig. 5). First, we inferred middle-slab lines along each profile that were drawn along the most intense positive seismic velocity anomalies within the slab. Then we identified segments boundaries along the middle slab lines, separately for v p and v s velocity, by assuming that they correspond to zones of clear intraslab seismic velocity decrease located in between individual positive seismic velocity anomalies. Segment boundary depths and along-slab positions were then defined at their locations within the middle-slab line and plotted on the same plot for different profiles and seismic velocities (Fig. 5d). By doing that we assumed that, owing to the variable quality of seismic data, not all segment boundaries can be always visible in each profile (Extended Data Fig. 8) whereas they should statistically cluster into similar positions measured along the slab (Fig. 5d). Error bars for this visual inspection correspond to around ±50 km for the depth of the middle-slab line and around ±100 km for defining thesegment boundary position along the line (Fig. 5d).

Data availability
All input files used in the numerical modelling are available at https://doi. org/10.17605/OSF.IO/bnvth. Source data are provided with this paper.

Code availability
The C and Matlab codes used for numerical experiments and visualization are available at https://doi.org/10.17605/OSF.IO/bnvth.  Table 2). b, f, Model with 10 km spaced faults dipping toward the trench (model xbeqm, Extended Data Table 2). c, g, Model with 5 km spaced faults dipping toward the trench (model xbeqn, Extended Data Table 2). d, h, Model with 10 km spaced faults dipping outward the trench (model xbeqo, Extended Data Table 2). Pre-existing faults are prescribed as 1 km wide and 14 km deep zones of weak basaltic crust and serpentinized mantle within stronger gabbroic crust and lithospheric mantle, respectively (Extended Data Table 1). Initial fault dip is 63°. Fig. 6 | Gradual development of large-offset normal faults in the reference model (Fig. 1). a-d Fig. 7). The distribution of v p (left column, panels a-e) and v s (right column, panels f-j) seismic velocity anomaly is based on the tomography model of Tao et al. 13 . Positions of segment boundaries (red triangles) defined along the middle-slab line (red solid lines) are inferred on the basis of visual inspection (Methods).