From strong plates to weak boundaries: strain localization in the lithospheric mantle with low- to high-temperature dislocation creep

Plate-like behavior in mantle convection models is commonly obtained using a strongly temperature-dependent viscosity together with a low yield stress that caps lithospheric strength. Yet, alternative mechanisms for limiting strength, including low-temperature plasticity, have been proposed. Here, we investigate how different rheological formulations promote lithospheric-scale strain localization, using 2-D thermo-mechanical simulations of upper mantle extension. We test a suite of olivine flow laws (diffusion creep and dislocation creep at low- and/or high-temperature), with and without a yield-stress cap, and introduce diagnostics that attribute plate weakening to either increasing strain rate (mechanical weakening) or increasing temperature (thermal weakening). Without a yield-stress cap, deformation remains distributed in all cases, indicating that weakening of the whole lithosphere is required for strain localization. In such simulations, a new extensional plate boundary develops in two stages: progressive narrowing of the deforming zone followed by rapid thinning of the weakened lithosphere. Across plate ages of 10–100 Myr, localization depends weakly on age but strongly on divergence rate, primarily through Stage-1 narrowing, and is controlled more by the stiffest lithospheric region than by the presence of a shallow weak layer. Including dislocation creep allows both mechanical and thermal weakening to operate within the deforming plate, enhancing feedbacks and accelerating localization relative to yield-stress plus diffusion creep. This implies that yield-stress/diffusion-only rheologies may overestimate break-up timescales in whole-mantle convection models. Low-temperature dislocation creep weakens the plate at 800–1000 K, providing a physically grounded mechanism for limiting strength in this temperature range. Finally, we propose that thermal weakening may in some cases promote rapid strength loss and accelerated extension during the late stages of natural rift evolution.

1 Introduction

The theory of plate tectonics describes the motion of rigid lithospheric plates separated by narrow, weak boundaries (e.g. Le Pichon, 1968). Observations from seismicity and geodetic strain-rate fields show that most deformation is concentrated within these plate boundaries, while most plate interiors deform only weakly (e.g. Kreemer et al., 2014). However, broad intraplate deformation zones also exist, for example between the India and Australia plates (e.g. Iaffaldano et al., 2018). In addition, plate reconstructions indicate that the number, size, and organization of plates have changed through geological time (e.g. Morra et al., 2013), implying that plate boundaries are created, reorganized, and abandoned. Such transient boundary formation is evident in incipient continental rifts (Müller et al., 2019), back-arc basins (Sdrolias and Müller, 2006), or subduction initiation settings (Lallemand and Arcay, 2021). New plate boundaries may develop by reactivating inherited weak zones (e.g. Foley and Bercovici, 2014; Heron et al., 2016; Mazzotti and Gueydan, 2018; Chen et al., 2020; Zhou et al., 2020; Fuchs and Becker, 2022; Brune et al., 2023; Gerya, 2024) or by localizing initially distributed deformation into a narrow shear zone (e.g. Lamb and Watts, 2010; Schmalholz et al., 2014; Asti et al., 2022; Whitney et al., 2023; Zwaan et al., 2024). Understanding when and how distributed deformation localizes remains central to explaining the dynamics of plate creation and the emergence of plate-like behavior in mantle convection.

A key difficulty is that the lithospheric mantle is expected to be strong in the creeping domain under low-temperature conditions. For olivine, the dominant mineral in the lithospheric mantle, high-temperature power-law dislocation creep predicts large flow stresses in the cold lithosphere (often exceeding hundreds MPa) (Kohlstedt et al., 1995; Karato, 2008; Burov, 2011). By contrast, estimates of available tectonic driving forces, of order 2 to 50 TN m⁻¹, suggest that plate interiors do not deform at deviatoric stresses below 200 MPa unless additional weakening processes operate (Parsons and Richter, 1980). Consistent with this, strongly temperature-dependent non-Newtonian viscosity alone is insufficient to promote lithospheric break-up (Bercovici, 1993), and instead tends to produce stagnant-lid convection (e.g. Solomatov and Moresi, 1997). Generating plate-like behavior therefore requires mechanisms that both reduce lithospheric strength and promote strain localization (e.g. Solomatov, 2004; Mulyukova and Bercovici, 2019).

Many candidate weakening processes have been explored at lithospheric scales, including shear heating (e.g. Yuen and Schubert, 1979; Brun and Cobbold, 1980; Fleitout and Froidevaux, 1980; Kaus and Podladchikov, 2006; Kiss et al., 2019), viscous anisotropy (e.g. Tommasi et al., 2009; Mameri et al., 2021; Duretz et al., 2025), strain softening or damage formulations (e.g. Frederiksen and Braun, 2001; Gueydan et al., 2014; Meyer et al., 2017; Tarayoun et al., 2019), and grain-size reduction leading to diffusion creep and/or grain-boundary sliding (e.g. Bercovici and Ricard, 2013; Gueydan and Précigout, 2014; Dannberg et al., 2017; Schierjott et al., 2020; Ruh et al., 2022). A widely used pragmatic approach in mantle convection modeling is to combine temperature-dependent creep with a yield-stress cap, which limits effective lithospheric strength and enables mobile-lid or plate-like regimes (e.g. Moresi and Solomatov, 1998; Tackley, 1998; Trompert and Hansen, 1998; Korenaga, 2010; Nakagawa and Iwamori, 2017; Coltice et al., 2019). However, the required stress limits are typically lower than ∼ 200–300 MPa (e.g. Richards et al., 2001; Crameri and Tackley, 2014; Mallard et al., 2016) and would correspond, assuming that yield stress represents brittle deformation, to effective friction coefficients (≲0.2, e.g. Moresi and Solomatov, 1998; Crameri et al., 2012; Nakagawa and Iwamori, 2017) substantially below laboratory values (∼ 0.6, e.g. Byerlee, 1978).

Low-temperature creep deformation has been proposed as an alternative route to limiting lithospheric strength without imposing an ad-hoc stress cap throughout the creeping domain (e.g. Jain et al., 2017). Experiments show that olivine can deform under high-stress, low-temperature conditions following an exponential (Peierls-type) flow law (Evans and Goetze, 1979; Frost and Ashby, 1982; Raterron et al., 2004; Mei et al., 2010; Demouchy et al., 2013; Hansen et al., 2019). Low-temperature dislocation creep has also been invoked to explain deformation of subducting slabs in the transition zone (Čížková et al., 2002; Garel et al., 2014). More recently, dislocation-dynamics models have provided a unified description of olivine dislocation creep that bridges classical high-temperature power-law creep and low-temperature Peierls creep, suggesting a continuous mechanical behavior across a wide range of temperatures and strain rates controlled by the interplay of glide and climb (Gouriet et al., 2019). This unified formulation has already been implemented in geodynamical simulations to investigate feedbacks that weaken the asthenosphere around sinking slabs (Garel et al., 2020).

Here, we evaluate the capacity of a unified low- to high-temperature (LT–HT) dislocation creep formulation (Gouriet et al., 2019; Garel et al., 2020) to promote lithospheric-scale strain localization in a controlled extensional setting. We compare its dynamical behavior with that obtained using commonly employed mantle rheologies (diffusion creep and/or high-temperature dislocation creep), with and without a yield-stress cap. Strain localization in geodynamical models is usually quantified using geometric measures such as “plateness” (Weinstein and Olson, 1992; Tackley, 2000) or relative shear-zone width “localization potential” (Montési, 2013), while bulk plate strength can be inferred from the tectonic force required to sustain imposed kinematics (Bialas et al., 2010; Brune et al., 2012). Recent studies have also introduced diagnostics of weakening to compare deformation mechanisms in subduction settings (Patočka et al., 2024), alongside theoretical analyses of weakening potentials (Montési and Zuber, 2002) or in simplified frameworks (Kaus and Podladchikov, 2006; Schmalholz and Fletcher, 2011). However, there remains a need for time-dependent diagnostics that track not only whether viscosity decreases, but also why it decreases as thermal and kinematic fields evolve self-consistently.

To address this, we introduce diagnostics that partition local viscosity change into contributions from increasing strain rate (mechanical weakening) and increasing temperature (thermal weakening) during localization. We apply these diagnostics to 2-D thermo-mechanical simulations of lithospheric extension across a range of initial plate ages and imposed extension rates. In this framework, deformation focusing and plate thinning are coupled to asthenospheric upwelling, enabling positive feedbacks that amplify strain rate and/or temperature increase depending on the rheological formulation. To isolate weakening processes within the lithospheric mantle, we model a single mantle material and intentionally neglect crustal layering and associated lithological and rheological contrasts. This simplified configuration is designed to clarify how low- to high-temperature dislocation creep, relative to yield-stress and diffusion-only parameterizations, controls the efficiency, depth distribution, and timescale of lithospheric-scale strain localization in the mantle.

In the remainder of the paper, we first describe the thermo-mechanical extension set-up, the rheological formulations tested, and the numerical implementation (Sect. 2). We then introduce post-processing diagnostics that quantify strain localization, plate thinning and bulk plate strength, and that partition viscosity change into strain-rate-driven and temperature-driven contributions (Sect. 3). We present results in four steps: (i) end-member deformation outcomes in a reference configuration, (ii) a two-stage localization trajectory and characteristic timescales in a representative localizing case, (iii) a comparison across rheological combinations to identify which dependencies control localization efficiency and the depth/temperature range of weakening, and (iv) sensitivity to initial plate age and imposed extension rate (Sect. 4). Finally, we discuss the implications for lithospheric weakening in natural rift systems, and for yield-stress-based parameterizations commonly used in whole-mantle convection models, and summarize the main conclusions (Sects. 5 and 6).

2 Methods: thermo-mechanical model

2.1 Extension Set-Up

We construct a 2-D thermo-mechanical model of upper-mantle extension in which lithospheric strength (viscosity) depends on both strain rate and temperature. The model domain is a 1200 km wide by 400 km deep rectangle composed entirely of mantle material (Table 1). This mantle-only configuration is designed to isolate strain localization within the lithospheric mantle and to facilitate comparison with large-scale plate-like convection models, which commonly neglect crustal structure and employ simplified rheologies.

Table 1 Parameter names and values for the model set-up.

Accordingly, we do not model crust-mantle interactions or full continental rift mechanics. This is a deliberate simplification, as crustal layering can reduce bulk lithospheric strength and introduce rheological contrasts that enhance localization in the mantle through crust-mantle mechanical coupling (e.g. Allemand and Brun, 1991; Burov and Watts, 2006; Rosenbaum et al., 2010). In many numerical rift models, strong coupling across rheological discontinuities favors a “narrow-rift” style of deformation (e.g. Huismans and Beaumont, 2005; Gueydan et al., 2008; Chenin et al., 2018; Duclaux et al., 2020; Zwaan et al., 2021; Wang et al., 2023). In our set-up, localization is instead controlled by mantle rheology and by the imposed yield-stress parameterization (Sect. 2.2), allowing us to isolate how these ingredients influence the onset and evolution of localized extension.

Figure 1 Simulation set-up in the 400 × 1200 km 2-D domain. Mechanical boundary conditions (in blue) are a free-surface top, and depth-dependent horizontal flow on the vertical sides, where $v_{\mathrm{1}/\mathrm{2}}$ is the plate half extension rate. The bottom boundary is no-slip and closed ($v_z=v_x=\mathrm{0}$), except over a 600 km-wide open segment between x=300 and 900 km where purely vertical flow is allowed (v_x=0). Thermal boundary conditions (in red) are constant temperatures at top (273 K) and bottom (1600 K), and insulating vertical boundaries.

The initial thermal structure is prescribed using the half-space cooling model for a laterally uniform plate age. Thermal boundary conditions are fixed temperatures of T_s=273 K at the surface and T_m=1600 K at the base, with thermally insulating vertical sides (Fig. 1). The initial deformation state is derived from a purely horizontal velocity field, producing a nearly uniform initial strain rate in the plate (${\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}$ = 3.73 × 10⁻¹⁶ s⁻¹ for $v_{\mathrm{1}/\mathrm{2}}$ = 1 cm yr⁻¹, Fig. 2).

Figure 2 Initial conditions showing (left) the temperature field calculated from the half-space cooling model (temperature profile T(z) for a 50 Myr-old plate) and (right) the second invariant strain rate ${\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}$ (here ∼ 3.73 × 10⁻¹⁶ s⁻¹ for a half extension rate $v_{\mathrm{1}/\mathrm{2}}$ = 1 cm yr⁻¹). The imposed initial velocity field is depicted by black arrows, with red lines indicating the 900 and 1500 K isotherms. This purely horizontal velocity field (v_z = 0) is calculated by assuming at each depth a linear increase of the horizontal velocity along x from the domain center ($v_x(x=\mathrm{600}\mathrm{km},z)=\mathrm{0}$) to the domain sides, where a Couette velocity is imposed as a boundary condition (cf. v_x(z) calculation for x=0 or 1200 km in the Supplement, Sect. S2.1).

Extension is imposed through symmetric, depth-dependent horizontal outflow velocities on the side boundaries (Fig. 1). The vertical profile v_x(z) is obtained from a 1-D Couette solution with surface velocity $v_{\mathrm{1}/\mathrm{2}}$ and zero velocity at the base, calculated using the depth- and temperature-dependent diffusion-creep viscosity of a 50 Myr-old plate (detailed in the Supplement, Sect. S2.1). This corresponds to an effective constant-velocity plate 89 km thick, defined as the depth above which horizontal velocity differs by less than 1 % from $v_{\mathrm{1}/\mathrm{2}}$ (Garel and Thoraval, 2021). Additional tests with alternative side-boundary velocity profiles produce no significant change in the results (detailed in the Supplement, Sect. S2.3). The upper boundary is treated as a free-surface (Kramer et al., 2012). At the base, no-slip and closed condition ($v_x=v_z=\mathrm{0}$) are imposed on the lateral segments, while vertical inflow with v_x=0 is allowed across a 600 km-wide central segment (see the Supplement, Sect. S2.2). This promotes localization at the domain center (see Fig. S3 in the Supplement, Sect. S2.3), while allowing basal inflow to adjust dynamically to the imposed lateral outflow (see Fig. S7 in the Supplement, Sect. S2.4) and evolving surface deformation.

We first analyze different rheological parameterizations in a reference configuration with $v_{\mathrm{1}/\mathrm{2}}$ = 1 cm yr⁻¹ and an initial plate age of 50 Myr (Sect. 4.1, 4.2 and 4.3). We then explore the influence of initial plate age (10–100 Myr) and half-extension rate (0.2–5 cm yr⁻¹) on localization behaviour (Sect. 4.4).

2.2 Rheological parameterization

2.2.1 Mechanisms of mantle deformation

Viscous deformation is represented by combinations of (i) olivine creep flow laws and (ii) a stress-limited viscosity (“yield stress”). The creep laws include diffusion creep and dislocation creep, based on experimental and numerical constraints on olivine deformation (e.g. Hirth and Kohlstedt, 1995a, b; Gouriet et al., 2019). We assume that multiple creep mechanisms may operate simultaneously (Sect. 2.2.2).

A yield-stress formulation is included in some simulations because such stress caps are widely used in geodynamical models to reproduce mobile-lid or plate-like behavior (e.g. Tackley, 2000; Mallard et al., 2016). Without this type of stress limitation, strongly temperature-dependent viscous rheologies typically produce stagnant-lid regimes (Solomatov, 1995). In that context, yield stress is commonly interpreted as a first-order proxy for pseudo-brittle deformation. Seismicity in oceanic lithosphere suggests that brittle failure may extend to temperatures of 600 °C (900–1000 K) (e.g. Engeln et al., 1986; Abercrombie and Ekström, 2001; McKenzie et al., 2005). We therefore test, in a subset of simulations, cases in which yielding is restricted to temperatures below 900–950 K. This prevents yielding from operating throughout the deeper creeping lithosphere while retaining a pseudo-brittle response in the shallow plate.

Elasticity is neglected in all simulations. Because our focus is on long-timescale viscous localization and plate-scale weakening, we expect this simplification to have limited influence on the localization trends examined here; we return to this point in the Discussion (Sect. 5.4).

2.2.2 Computation of mantle effective viscosity

Diffusion creep in olivine is represented by a temperature-dependent Newtonian viscosity:

$$\begin{array}{}\text{(1)}& {\mathit{\eta }}_{\mathrm{diff}}=A_{\mathrm{diff}}^{-\mathrm{1}}\cdot \exp \left({\displaystyle \frac{E_{\mathrm{diff}}+P\cdot V_{\mathrm{diff}}}{R\cdot (T+\mathit{\delta }T)}}\right),\end{array}$$

where A_diff is the pre-exponential constant, corresponding to a grain size of ∼ 3 mm (Garel et al., 2020), E_diff is the activation energy, V_diff is the activation volume, R the gas constant (Table 2), P is the lithostatic pressure (ρgz), T the temperature, and δT is a temperature correction applied only to viscosity calculations for diffusion (Eq. 1) and high-temperature dislocation creep (Eq. 2), to mimic the effect of a 0.5 K km⁻¹ adiabatic gradient within the incompressible approximation. We compare dynamic and lithostatic pressures in the Supplement (Sect. S3.1).

Table 2 Variables used in the rheological laws investigated in this study (Sect. 2.2).

Dislocation creep is non-Newtonian and strain-rate-dependent. We consider two alternative flow laws. First, we use conventional high-temperature (HT) power-law dislocation creep (Garel et al., 2020):

$$\begin{array}{}\text{(2)}& {\mathit{\eta }}_{\text{disl-HT}}=A_{\mathrm{disl}}^{-\frac{\mathrm{1}}{n}}\cdot {\dot{\mathit{\epsilon }}}_{\mathrm{disl}}^{\frac{\mathrm{1}-n}{n}}\cdot \exp \left({\displaystyle \frac{E_{\mathrm{disl}}+P\cdot V_{\mathrm{disl}}}{n\cdot R\cdot (T+\mathit{\delta }T)}}\right)\end{array}$$

where A_disl is the pre-exponential constant, E_disl is the activation energy, V_disl the activation volume, n the stress exponent and ${\dot{\mathit{\epsilon }}}_{\mathrm{disl}}$ the second invariant of the strain-rate tensor associated with the dislocation creep contribution. Second, we use a unified low- to high-temperature (LT–HT) dislocation creep law, derived from dislocation-dynamics models (Gouriet et al., 2019) and calibrated against macroscale effective mantle-viscosity constraints (Garel et al., 2020):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\eta }}_{\text{disl-LT–HT}}={\displaystyle \frac{{\mathit{\sigma }}_{\text{disl-LT–HT}}}{\mathrm{2}\cdot {\dot{\mathit{\epsilon }}}_{\mathrm{disl}}}}\text{with}\text{(4)}& {\displaystyle }{\mathit{\sigma }}_{\text{disl-LT–HT}}=A_{\mathrm{0}}\left(T\right)\left(\mathrm{1}+\tanh \left[A_{\mathrm{1}}\left(T\right)\left[{\log }_{\mathrm{10}}\left({\dot{\mathit{\epsilon }}}_{\mathrm{disl}}\right)-A_{\mathrm{2}}\left(T\right)\right]\right]\right)\end{array}$$

where A₀, A₁, A₂ are polynomial functions of temperature (Table 2). This formulation converges to HT dislocation creep (Eq. 2) at high temperature, and predicts similar low stresses at lower temperature compared to other parameterizations (e.g. Mei et al., 2010; Demouchy et al., 2013; Jain et al., 2017; Warren and Hansen, 2023) (Fig. S12 in the Supplement, Sect. S3.3). When both diffusion and dislocation creep are present in the rheological parameterization, the bulk creep viscosity is computed assuming that the two mechanisms act in parallel:

$$\begin{array}{}\text{(5)}& {\mathit{\eta }}_{\mathrm{creep}}={\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{diff}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{disl}}}}\right)}^{-\mathrm{1}}\end{array}$$

and that the second invariant of the total strain rate tensor (${\dot{\mathit{\epsilon }}}_{\mathrm{II}}$) is partitioned between them:

$$\begin{array}{}\text{(6)}& {\dot{\mathit{\epsilon }}}_{\mathrm{II}}={\dot{\mathit{\epsilon }}}_{\mathrm{diff}}+{\dot{\mathit{\epsilon }}}_{\mathrm{disl}}\end{array}$$

assuming a single deviatoric stress σ

$$\begin{array}{}\text{(7)}& \mathit{\sigma }=\mathrm{2}{\mathit{\eta }}_{\mathrm{creep}}{\dot{\mathit{\epsilon }}}_{\mathrm{II}}=\mathrm{2}{\mathit{\eta }}_{\mathrm{diff}}{\dot{\mathit{\epsilon }}}_{\mathrm{diff}}=\mathrm{2}{\mathit{\eta }}_{\mathrm{disl}}{\dot{\mathit{\epsilon }}}_{\mathrm{disl}}\end{array}$$

An iterative scheme is used to compute η_disl from the dislocation creep strain rate (detailed in the Supplement, Sect. S3.2) rather than the total strain rate, thereby avoiding artificial weakening. This allows us to quantify the fraction of total creep deformation accommodated by dislocation creep (Eq. S13 in the Supplement).

The pseudo-brittle contribution is represented by a non-Newtonian yield viscosity:

$$\begin{array}{}\text{(8)}& {\mathit{\eta }}_{\mathrm{yield}}={\displaystyle \frac{{\mathit{\sigma }}_y}{\mathrm{2}\cdot {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}}\end{array}$$

where σ_y is a constant yield stress (200–500 MPa depending on the simulation). The effective viscosity used in the momentum equation is then taken as the minimum of creep and yield viscosities:

$$\begin{array}{}\text{(9)}& \mathit{\eta }=min({\mathit{\eta }}_{\mathrm{creep}},{\mathit{\eta }}_{\mathrm{yield}}).\end{array}$$

Numerical stability is maintained by imposing viscosity cutoffs of 10¹⁸ and 10²⁵ Pa s. At each point, we also identify the dominant deformation mechanism: yielding versus creep, and within creep, diffusion versus dislocation, based on the creep mechanism contributing ≥ 50 % of the total strain rate.

2.2.3 Strategy for investigating rheological control on lithospheric extension

We investigate a suite of rheological combinations (Table 3):

• Diffusion creep only (D),

• Diffusion creep combined with either HT or a LT–HT dislocation creep (D−d_HT or D−d_LT–HT),

• Diffusion creep combined with a yield-stress rheology ($D-Y_{\mathrm{200}/\mathrm{500}\mathrm{MPa}}$),

• Diffusion creep, dislocation creep, and yield-stress combined ($D-d_{\text{HT or LT–HT}}-Y_{\mathrm{200}/\mathrm{500}\mathrm{MPa}}$).

The two constant yield stresses, 200 and 500 MPa, are chosen to span the range between values commonly used to generate plate-like behaviour in mantle convection models (e.g. Mallard et al., 2016), and the upper stress range over which the LT–HT dislocation creep formulation is calibrated (Gouriet et al., 2019; Garel et al., 2020). In a subset of simulations, yielding is restricted to the expected brittle domain (see Sect. 5.3), by allowing it only below 900 or 950 K ($Y_{\text{200 or 500\hspace{0.17em}MPa}}^{\mathrm{950}\mathrm{K}}$ or $Y_{\text{200 or 500\hspace{0.17em}MPa}}^{\mathrm{900}\mathrm{K}}$).

These rheological combinations differ in both their strain-rate and temperature dependencies, and in the depth and temperature at which deformation transitions from yield-dominated to creep-dominated behaviour. They therefore provide a controlled way to test how the depth distribution of weakening influences localization, asthenospheric upwelling, and the timing of lithospheric break-up under constant extension velocity.

2.3 Numerical methods

We solve the conservation equations of mass, momentum and energy for an incompressible Stokes fluid under the Boussinesq approximation using the finite-element, control-volume code Fluidity (e.g., Davies et al., 2011). This framework has been widely used in geodynamical applications and extensively benchmarked against analytical and numerical reference solutions (e.g. Garel et al., 2014; Le Voci et al., 2014; Kramer et al., 2021b; Duvernay et al., 2022). The equations are discretized on an unstructured Eulerian mesh that is dynamically adapted throughout each simulation. Refinement criteria are applied to temperature, velocity, strain-rate, and viscosity, ensuring adequate resolution of lithospheric thermal gradients, localization zones, and strong viscosity contrasts. Triangular element sizes range from 200–1000 m in the finest regions to 50 km in the coarsest regions, with the highest resolution concentrated in the lithosphere and regions of active deformation (Supplement, Sect. S1). A typical simulation contains approximately 50 000 nodes, although this varies as simulations evolve.

3 Post-processing diagnostics

For post-processing, physical fields are interpolated from the unstructured finite-element mesh onto a regular grid with 5 × 1 km² spacing. This allows consistent calculation of spatial averages and Eulerian time derivatives.

We define the lithosphere as the region colder than 1500 K. In this single-material model, that threshold provides a first-order proxy for the lithosphere-asthenosphere transition, which is otherwise gradual in temperature, deformation and velocity (Garel and Thoraval, 2021). We also pay particular attention to the deeper lithospheric interval between 800 and 1500 K, where the yield-creep transition and the majority of the weakening occur. Some diagnostics are vertically-averaged over the lithosphere and denoted 〈X〉_1500 K(x,t), using either an arithmetic or geometric mean depending on the range spanned by X. To compare simulations with different extension rates $v_{\mathrm{1}/\mathrm{2}}$, we define a bulk extensional strain $\mathit{\epsilon }\left(t\right)=\frac{\mathrm{2}\cdot t\cdot v_{\mathrm{1}/\mathrm{2}}}{L}$ where t is time and L is the domain width. All simulations are run to ∼ 40 % bulk strain, corresponding to ∼ 500 km cumulative horizontal extension.

3.1 Strain localization in the lithospheric mantle

We first quantify strain localization within the plate using the strain-rate amplification $A_{\mathrm{SR}}(x,z,t)$ relative to the initially uniform strain rate in the plate ${\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}$ (Fig. 2), and its geometric mean over lithospheric thickness at each horizontal position.

To track lateral focusing through time, we define the width of the deforming zone, w(t), as the region where $\langle A_{\mathrm{SR}}{\rangle }_{\mathrm{1500}\mathrm{K}}(x,t)>\mathrm{1}$. From this, we compute a narrowing rate. We also quantify the lateral viscosity contrast within the plate R_visc(t), as the ratio between the maximum vertically-averaged lithospheric viscosity along x and the minimum value in the central weak zone:

$$\begin{array}{}\text{(10)}& R_{\mathrm{visc}}\left(t\right)={\displaystyle \frac{\max \left({〈\mathit{\eta }〉}_{\mathrm{1500}\mathrm{K}}(x,t)\right)}{\min \left({〈\mathit{\eta }〉}_{\mathrm{1500}\mathrm{K}}(x,t)\right)}}\end{array}$$

Finally, we compute the plateness P(t) (Tackley, 2000; Crameri, 2018). Details of the calculation are given in the Supplement (Sect. S4.1).

3.2 Plate thinning, weakening, and bulk stiffness

To quantify lithospheric thinning, we define the plate thickness h_plate(t) as the minimum depth of the 1500 K isotherm along x. Its time evolution is used to estimate an upwelling rate. We track time-dependent weakening or hardening through the logarithmic rate of change in effective viscosity:

$$\begin{array}{}\text{(11)}& W(x,z,t)=-{\displaystyle \frac{\partial {\log }_{\mathrm{10}}\left(\mathit{\eta }\right)}{\partial t}}\approx -{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}{\log }_{\mathrm{10}}\left({\displaystyle \frac{{\mathit{\eta }}_{t+\frac{\mathrm{\Delta }t}{\mathrm{2}}}}{{\mathit{\eta }}_{t-\frac{\mathrm{\Delta }t}{\mathrm{2}}}}}\right)\end{array}$$

where Δt is chosen to be much smaller than the localization timescale (typically 0.05–1 Myr, depending on extension rate). Smaller values of Δt do not change the weakening amplitude. Positive values of W indicate weakening whereas negative values indicate hardening. We also compute the vertically averaged weakening rate 〈W〉_1500 K(x,t) using a geometric mean. A 10-fold decrease in viscosity in 1 Myr corresponds to a weakening rate W of about 3 × 10⁻¹⁴ s⁻¹.

As a bulk measure of lithospheric strength, we calculate the tectonic force TF(t) (in N m⁻¹) required to maintain the imposed extension velocity at the side boundaries (Bialas et al., 2010; Brune et al., 2012):

$$\begin{array}{}\text{(12)}& \mathrm{TF}\left(t\right)=H\cdot \sqrt{(S_{xx}-S_{zz}{)}^{\mathrm{2}}+S_{xz}^{\mathrm{2}}}\end{array}$$

with

$$\begin{array}{}\text{(13)}& S_{ij}\left(t\right)={\displaystyle \frac{\mathrm{1}}{H}}\underset{\mathrm{0}}{\overset{H}{\int }}{\mathit{\sigma }}_{ij}(z,t)\mathrm{d}z\end{array}$$

where S_ij is the stress tensor vertically averaged over the full domain height H at x=100 km, to minimize boundary effects (see the Supplement, Sect. S4.2).

3.3 Physical controls on plate weakening

To distinguish why viscosity decreases during localization, we partition weakening into contributions associated with temperature and strain-rate changes. This builds on the general idea that localization potential can be related to the sensitivity of viscosity to evolving state variables (Montési and Zuber, 2002). At each point on the interpolated grid, viscosity depends only on temperature and strain rate, so its time derivative may be written as:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial \mathit{\eta }(T,\dot{\mathit{\epsilon }})}{\partial t}}={\displaystyle \frac{\partial \mathit{\eta }}{\partial T}}{\displaystyle \frac{\partial T}{\partial t}}+{\displaystyle \frac{\partial \mathit{\eta }}{\partial \dot{\mathit{\epsilon }}}}{\displaystyle \frac{\partial \dot{\mathit{\epsilon }}}{\partial t}}\end{array}$$

We approximate this decomposition over a finite time interval Δt as:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}}={\left.{\displaystyle \frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}}\right|}_{\dot{\mathit{\epsilon }}}+{\left.{\displaystyle \frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}}\right|}_T\text{with}\\ & \begin{array}{l}{\left.\frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}\right|}_{\dot{\mathit{\epsilon }}}=\frac{\mathit{\eta }\left(T\left(x,z,t+\frac{\mathrm{\Delta }t}{\mathrm{2}}\right),\dot{\mathit{\epsilon }}(x,z,t)\right)-\mathit{\eta }\left(T\left(x,z,t-\frac{\mathrm{\Delta }t}{\mathrm{2}}\right),\dot{\mathit{\epsilon }}(x,z,t)\right)}{\mathrm{\Delta }t}\\ {\left.\frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}\right|}_T=\frac{\mathit{\eta }\left(T(x,z,t),\dot{\mathit{\epsilon }}\left(x,z,t+\frac{\mathrm{\Delta }t}{\mathrm{2}}\right)\right)-\mathit{\eta }\left(T(x,z,t),\dot{\mathit{\epsilon }}\left(x,z,t-\frac{\mathrm{\Delta }t}{\mathrm{2}}\right)\right)}{\mathrm{\Delta }t}\end{array}\end{array}\end{array}$$

The corresponding normalized contributions to weakening are then

$$\begin{array}{}\text{(16)}& F_{\mathrm{T}}={\displaystyle \frac{{\left.\frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}\right|}_{\dot{\mathit{\epsilon }}}}{\frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}}}\text{and}{F}_{\mathrm{SR}}={\displaystyle \frac{{\left.\frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}\right|}_T}{\frac{\mathrm{\Delta }\mathit{\eta }}{\mathrm{\Delta }t}}}\end{array}$$

where F_T measures the contribution from temperature evolution (“thermal” weakening) and F_SR, the contribution from strain-rate evolution (“mechanical” weakening). We use Δt≈0.25 Myr for this partitioning. Smaller values introduce short-period oscillations in F_T and F_SR without changing the broader trends. These diagnostics are evaluated only where the material is undergoing weakening (W>0). Mixed cases can nevertheless arise, for example when mechanical weakening due to strain-rate increase exceeds thermal hardening due to plate cooling. In such cases, F_T may become negative and F_SR may exceed 1 (as discussed in Sect. 4.4). For plotting purposes, the colour scale for F_SR is capped between 0 and 1.

4 Results

We present our results in four steps. First, we describe the range of deformation patterns obtained for a reference plate age (50 Myr) and half-extension rate ($v_{\mathrm{1}/\mathrm{2}}=\mathrm{1}$ cm yr⁻¹), from distributed deformation to localization into a narrow extensional boundary. Second, we use a representative localizing case to define characteristic times and a two-stage evolution of localization based on the co-evolution of deforming-zone width and plate thickness. Third, we compare rheological parameterizations to explain the differences in localization efficiency and in the depth/temperature range over which weakening occurs. Fourth, we explore sensitivity to initial plate age and imposed extension rate.

4.1 Diffuse vs. localized deformation during lithospheric extension

For an initially 50 Myr-old plate extended at $v_{\mathrm{1}/\mathrm{2}}=\mathrm{1}$ cm yr⁻¹, the rheological parameterization produces three distinct deformation outcomes by 40 % bulk strain (Table 3, Figs. 3, A1).

Figure 3 Evolution of the strain rate field (second invariant) during two end-member simulations, shown after 1, 6 and 12 Myr of extension: (a) Simulation D−d_HT (ref-0), with diffusion creep and high-temperature dislocation creep, Table 3 characterized by uniform and diffuse lithospheric thinning; (b) simulation $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ (ref-1), combining diffusion creep, low- and high-temperature dislocation creep, and a yield-stress set to 500 MPa, which successfully localizes plate deformation. For each snapshot, the horizontal velocity at the surface is plotted above. Velocity field is depicted by black arrows, and the 900 and 1500 K isotherms are indicated in red.

First, Simulation ref-1 serves as a reference for the strain localization scenario (Table 3, and Fig. 3b). Its rheology, labeled $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$, accounts for diffusion creep, low- and high-temperature dislocation creep (Eq. 3) and a yield stress of 500 MPa (Eq. 8). This simulation shows a progressive narrowing of the deformed domain around the box center, where most of the deformation is accommodated as shown by the significant increase in the horizontal velocity gradient at the surface (Fig. 3b, 1–6 Myr). In a second stage, significant asthenospheric upwelling develops beneath the most stretched lithosphere portion, resulting in a boundary separating two divergent plates that behaves as a horizontal velocity discontinuity (Fig. 3b, 12 Myr). Hence, in localizing cases, a narrow extensional plate boundary develops at the domain center, separating two rigid plates (Simulation ref-1 in Fig. 3b). These runs exhibit high plateness (P>0.95), a large lateral viscosity contrast (R_visc>100), and a deformed-zone width w<100 km, while plate thickness h_plate approaches a quasi-steady value after localization. The localized zone is underlain by focused asthenospheric upwelling.

Second, the non-localizing simulation ref-0 (rheology D−d_HT, Table 2 and Fig. 3a) exhibits a long-lasting and steady plate deformation. This diffuse deformation mode corresponds to a slow plate thinning (∼ 0.05 cm yr⁻¹), arising from the competition between thickening by diffusive cooling and lithosphere extension. In such cases where deformation is always distributed, the plate deforms broadly and remains laterally uniform, without any increase in the horizontal velocity gradient at the surface (e.g., Simulation ref-0 in Fig. 3a). These runs show low plateness (P<0.3) and weak lateral viscosity contrast (R_visc<10). Strain-rate amplification remains modest across the plate and deformation is expressed primarily as gradual, laterally uniform thinning.

Third, an intermediate regime occurs in which the deformation localization is very slow (Appendix A). Strain begins to focus but the area where deformation is enhanced remains relatively wide at 40 % strain (incomplete localization). These runs still exhibit narrowing of the deformed zone and plate thinning, with w≳150 km at the end of the simulations, comparable to the early-to-intermediate state of Simulation ref-1 (e.g., Fig. 3b at 6 Myr).

In the remainder of the Results, we use a representative localizing simulation (ref-1) to define a reference localization trajectory and characteristic timescales, and then compare rheological combinations and forcing parameters against that trajectory.

4.2 A two-stage scenario of strain localization and lithosphere weakening

We quantify the localization trajectory in Simulation ref-1 using the co-evolution of deformed-zone width w, the minimum plate thickness h_plate, and the weakening rate W (Sect. 3). The deforming region (defined by 〈A_SR〉_1500 K>1, Sect. 3.1) progressively narrows toward the domain center (Fig. 4a), while weakening intensifies within the lithosphere beneath the developing boundary (Fig. 4b and d).

Figure 4 Localization of lithospheric deformation and weakening in Simulation $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ (ref-1, Table 3). (a) Time-distance evolution of strain-rate amplification 〈A_SR〉_1500 K(x,t), calculated as the vertically-averaged ratio relative to the initial uniform strain rate in the plate (here ${\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}\sim \mathrm{3.73}$ × 10⁻¹⁶ s⁻¹, see Sect. 3.1). The thick black line features the isocontour 1 of this ratio. (b) Time-distance plot of the geometric mean of the vertically-averaged weakening rate 〈W〉_1500 K(x,t) (Sect. 3.2). Orange shades indicate where weakening (W>0) affects more than 60 % of lithospheric thickness, with plotted values representing the geometric mean 〈W_SR〉_1500 K(x,t) for depths where $W(x,z,t)>\mathrm{0}$. On the other hand, purple-blue shades indicate where weakening affects less than 40 % of lithospheric thickness, plotting the vertical geometric mean $-\langle -W_{\mathrm{SR}}{\rangle }_{\mathrm{1500}\mathrm{K}}(x,t)$ for depths where W≤0 (hardening or neutral). Areas where weakening affects between 40 % and 60 % of the lithospheric thickness are left white. (c, d) Depth-time evolution of the weakening rate W at x=300 km (c) and 600 km (d). The pink outline highlights zones of high weakening rate where $\frac{W}{{\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}}>\mathrm{30}$, and the red lines indicate isotherms from 800 to 1500 K. Regions with different dominant deformation mechanisms are delineated by dashed dark-blue contours. In all panels, vertical dashed lines indicate the transition time, t_tr and the plate boundary time, t_PB.

The localization trajectory is well described by two stages (Fig. 5a). Stage 1 is dominated by lateral focusing of deformation: w decreases rapidly whereas h_plate decreases only modestly. Stage 2 is dominated by rapid plate thinning driven by asthenospheric upwelling: h_plate decreases rapidly while the narrowing of w slows down. Stages 1 and 2 are detailed below. We define two characteristic times. The transition time t_tr marks a shift in trajectory in h_plate(t) : w(t) space (Appendix B); it coincides with the end of the approximately steady tectonic force plateau TF (Fig. 5c). The plate-boundary time t_PB marks the onset of a quasi-steady geometry, when both h_plate and w become approximately stationary and the tectonic force has dropped to its post-localization value (Fig. 5c).

Figure 5 Results for Simulation $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$, ref-1. (a) Temporal evolution (and corresponding cumulative strain in %) of the width of deformed zone w (black) and minimum lithosphere thickness (orange) as defined in Sect. 3.1. Average narrowing and upwelling rates are calculated for Stage 1 and Stage 2, respectively. Rates can be expressed in km Myr⁻¹ or in cm yr⁻¹ (where 1 cm yr⁻¹ = 10 km Myr⁻¹). Temporal evolution (and corresponding cumulative strain in %) of (b) the maximum weakening rate (blue) within the lithosphere, defined as the region between the 273 and 1500 K isotherms at x=600 km (Eq. 11), and plate lateral viscosity contrast (green) R_visc (Eq. 10), and (c) tectonic force (yellow) and plateness (pink). In all panels, vertical dashed lines indicate the transition time, t_tr and the plate boundary time, t_PB, also displayed on the curves by dots (for t_tr) and stars (for t_PB).

Table 3 List of numerical simulations with a half-extension rate $v_{\mathrm{1}/\mathrm{2}}$ of 1 cm yr⁻¹ and an initial plate age of 50 Myr. Reference simulations for scenarios of distributed deformation “distr. def.” (ref-0) or localized plate-boundary “loc. PB” (ref-1) are shown in Fig. 3a and b, respectively. The intermediate regime of incomplete localization “inc. loc.” is shown in Fig. D2c and f. The characteristic times of transition t_tr and plate-boundary formation t_PB are defined in Sect. 4.2.

Figure 6 Co-evolution of the minimum plate thickness h_plate(t) as a function of the width w(t) of the deformed zone (see Sect. 3.2) for experiments with the same initial age (50 Myr) and extension velocity $v_{\mathrm{1}/\mathrm{2}}$ = 1 cm yr⁻¹ for 8 different rheological combinations. The transition time t_tr between Stages 1 and 2 is indicated by a circle symbol and the localization time t_PB by a star symbol. Open and filled symbols correspond to rheologies with a yield stress set to 200 and 500 MPa, respectively. The square represent the initial time.

For Simulation ref-1, Stage 1 lasts 6.5 Myr and is characterized by rapid narrowing (mean 14 cm yr⁻¹) with limited thinning (100 to 82 km in Fig. 5a). During this stage, weakening at the domain center is distributed through much of the lithosphere thickness (Fig. 4d) and exceeds off-center weakening (Fig. 4c), producing a growing but moderate viscosity contrast (R_visc<50; Fig. 5b). In the plate interior outside the deforming zone, viscosity increases (hardening), particularly toward the end of Stage 1 (Fig. 4b).

Stage 2 lasts 5.5 Myr (from 6.5 to 12 Myr) and is characterized by rapid thinning (82 to 10 km) and a sharp reduction in tectonic force, whereas narrowing is slower (mean ∼ 4 cm yr⁻¹; Fig. 5a). Weakening peaks shortly before t_PB (Fig. 5b), while the lateral viscosity contrast continues to increase and reaches its maximum near t_PB (Fig. 5b). After t_PB, the system evolves toward a steady thermal structure around the mature boundary and the weakening rate in the narrow boundary region approaches zero (Fig. 4b).

4.3 Influence of rheological parameterization on localization scenario

Across all rheological combinations, simulations that localize deformation follow the same two-stage trajectory defined in Sect. 4.2 (Fig. 6). Stage 1 is dominated by progressive narrowing of the deforming zone with only modest thinning, whereas Stage 2 is dominated by rapid plate thinning associated with focused asthenospheric upwelling and only limited additional narrowing (Fig. 7a–b). The major rheological control is therefore whether deformation localization occurs and how quickly the system progresses through the two stages, rather than the geometrical trajectory or the maximum weakening rate (10⁻¹³ s⁻¹ for all simulations in Fig. 7f).

Figure 7 Temporal evolution (and corresponding cumulative strain in %) of diagnostics for simulations with 8 different rheological combinations, and the same initial plate age (50 Myr) and extension velocity ($v_{\mathrm{1}/\mathrm{2}}$ = 1 cm yr⁻¹). (a) Width of deformed zone w, (b) h_plate, minimum plate thickness along x, (c) tectonic force TF, (d) plateness P, (e) lateral viscosity contrast R_visc, and (f) maximum weakening rate W_max over the whole lithosphere thickness (T < 1500 K) at x = 600 km. In all panels, characteristic times are indicated by dots (t_tr) and stars (t_PB).

In simulations without a yield-stress cap, deformation remains distributed regardless of whether diffusion creep is combined with high-temperature dislocation creep or with low- to high-temperature dislocation creep (e.g., Simulation ref-0; Fig. 3a). Thus, in this lithospheric extension setup, dislocation and diffusion creep alone are insufficient to generate lithospheric-scale localization; additional weakening of the cold upper lithosphere (here provided by a yield-stress cap) is required.

Among simulations that do localize strain, the characteristic times vary substantially between rheologies (Table 3). Localization takes about twice as long for diffusion creep plus yield stress (e.g., D−Y_200 MPa and D−Y_500 MPa) and for cases in which yielding is restricted to temperatures below 950 K (e.g., $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$), than for rheologies that combine diffusion creep, dislocation creep, and yield stress (e.g., $D-d_{\mathrm{HT}}-Y_{\mathrm{200}/\mathrm{500}\mathrm{MPa}}$ and $D-d_{\text{LT-HT}}-Y_{\mathrm{200}/\mathrm{500}\mathrm{MPa}}$; Fig. 7a–b). Consistently, the tectonic force required to maintain the imposed extension velocity during Stage 1 is the largest for D−Y_500 MPa (Fig. 7c), indicating a stronger deforming zone at the plate center, as also suggested by the low lateral viscosity contrast R_visc (Fig. 7e).

Figure 8 Depth-time evolution of the relative contributions to weakening from strain rate (F_SR) and temperature (F_T), computed at box center x=600 km, for simulations (a) $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ (Ref-1); (b) $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$; (c) D−Y_500 MPa; and (d) $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$. The diagnostics are calculated only in regions in weakening (W>0), with hardening or stable regions (W≤0) shown in white. The color scale is capped between 0 and 1. Areas of intense weakening, where the weakening rate normalized by the initial strain rate in the plate $\frac{W}{{\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}}$ exceeds 30, are highlighted by pink contours. Boundaries between dominant deformation mechanisms are marked by dashed dark-blue lines, and isotherms from 800 to 1500 K are displayed in light blue.

We interpret these differences using the weakening-partition diagnostics (Sect. 3.3), which quantify the relative contributions of strain-rate increase (F_SR) and temperature increase (F_T) to viscosity reduction (Fig. 8). In all simulations, the shallow lithosphere is yield-dominated and weakening is purely strain-rate driven there (F_SR=1). Differences emerge in the deeper lithosphere, where deformation transitions to creep and where maximum weakening occurs. The temperature of the yield-to-creep transition depends strongly on the creep rheology: it is highest for diffusion creep alone (about 1300 K in D−Y_500 MPa), lower when high-temperature dislocation creep is included (about 1000 K in $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$), and the lowest when low- to high-temperature dislocation creep is included (about 800 K in $D-d_{\text{LT-HT}}-Y_{\mathrm{500}\mathrm{MPa}}$; Fig. 8). These shifts imply that a thicker portion of the lithosphere deforms in the creep regime when dislocation creep is included, extending thermo-mechanical weakening to lower temperatures.

The weakening-partition diagnostics further show that dislocation creep enables both mechanical and thermal weakening within the creeping lithosphere. In rheologies that include dislocation creep, weakening is predominantly strain-rate driven during Stage 1 (typically F_SR>60 % in the main weakening zone), while weakening becomes almost purely temperature-driven during Stage 2 (typically F_T>0.9; Fig. 8a, b). In contrast, for diffusion creep plus yield stress, weakening in the creeping domain is entirely temperature-driven (F_T≥1; Fig. 8c), limiting mechanical weakening during Stage 1 and delaying the onset of rapid thinning.

This interpretation is illustrated by Fig. 9, which shows theoretical viscosity profiles for two idealized end-member evolutions: a 10-fold increase in strain rate at constant plate age, used as a proxy for Stage 1, and plate thinning at constant strain rate, used as a proxy for Stage 2. For a given creep law, the effect of a Stage-1-type strain-rate increase can be estimated from the viscosity contrast between the yield-creep transition and the warm sub-lithospheric layer at about 75 km depth (Fig. 9b). By this measure, the mechanical weakening associated with LT–HT dislocation creep (6-fold) is intermediate between that associated with HT dislocation creep (5-fold) and that associated with the yield-stress-plus-diffusion-creep rheology (10 or 1-fold). However, the 10-fold strain rate increase during Stage 1 produces the lowest absolute viscosities below the yield-creep transition. The same approach can be used to estimate the amplitude of Stage-2-type thermal weakening from the viscosity change accompanying lithospheric thinning and warming (Fig. 9d). Although the viscosity reduction associated with lithospheric warming is largest for diffusion creep (>100) and smallest for LT–HT dislocation creep (<40), the absolute viscosity reached after a given amount of warming and thermal rejuvenation remains lower when dislocation creep is included (Fig. 9c).

Figure 9 Theoretical (a, c) viscosity profiles, and (b, d) viscosity ratio illustrating (a–b) the effect of a 10-fold increase in strain rate with constant plate age of 50 Myr (proxy of Stage 1 evolution if temperature increase is neglected), and (c–d) the effect of thermal-structure variations associated with decreasing plate age from 40 to 10 Myr at a constant strain rate of 3 × 10⁻¹⁵ s⁻¹ (proxy of plate thinning in Stage 2 assuming no mechanical contribution). Profiles are computed for 4 rheologies, with vertical black lines in panels (a) and (c) indicating yield viscosities corresponding to 200 or 500 MPa (Eq. 8).

Finally, in a subset of simulations, yielding is restricted to low temperature to disable pseudo-brittle mechanism above 950 K, leaving only creep mechanisms at higher temperatures (Sect. 2.2.1). This produces a marked delay even when dislocation creep is present (compare $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ and $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$; Table 3). Vertical profiles for the $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$ case show that imposing a thermal cutoff on the yield rheology creates a stiff layer less than 4 km thick at the base of the shallow yielding domain, between the colder lithosphere above and the warmer creeping mantle below (Appendix C). Although thin, this layer remains mechanically important because it forms a local viscosity peak that slows the overall decrease in viscosity throughout the whole lithosphere, which in turns delays deformation focusing. This highlights the importance of capping the lithosphere strength in the creep domain within the mid-lithosphere temperature range (roughly 800–1300 K in these setups): even a thin stiff layer at these depths (Fig. 9) is sufficient to impede Stage 1 narrowing and to delay both t_tr and t_PB (Figs. 6, 7; Table 3). This explains why lowering the yielding cutoff to 900 K inhibits localization within the simulated strain window (Table 3). Additional resolution tests show that, despite its striking sharpness, the presence of this high-stiffness layer is well resolved, and is thus a robust result (Sect. S1 in the Supplement).

4.4 Influence of initial plate age and imposed extension rate

We test the robustness of the two-stage localization pathway and its sensitivity to: (i) the initial thermo-mechanical state of the lithospheric mantle and (ii) the magnitude of far-field forcing. We vary initial plate age (from 10 to 100 Myr), and the imposed half-extension rate (from 0.2 to 5 cm yr⁻¹). We focus on two rheological combinations that exhibit markedly different rates of strain localization in the reference setup (Fig. 7): diffusion creep plus yield stress (D−Y_500 MPa), and diffusion creep plus low-to-high-temperature dislocation creep plus yield stress ($D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$) (Table 4).

Across most age-velocity combinations, the localization trajectory remains fundamentally characterized by a two-stage scenario: Stage 1 is characterized by progressive narrowing of the deforming zone, while plate thickness changes only modestly; Stage 2 is characterized by rapid plate thinning associated with asthenospheric upwelling, with comparatively minor additional narrowing (Appendix D). The basic control on characteristic times is the extension rate, whereas initial plate age exerts only a weak influence (Fig. 10a, b). This set of experiments further underlines that dislocation creep accelerates deformation localization relative to diffusion creep plus yield stress. For a given extension rate, the dislocation-creep rheology produces lower effective viscosity in the central deforming zone (typically by a factor of 5 to 10 at t_tr), which increases the lateral viscosity contrast between the weak zone and the plate interior and reduces the tectonic force required to maintain the imposed velocity (Fig. 10e, g). This enhanced strength contrast is expressed primarily through Stage 1: for all tested ages and velocities, the Stage 1 duration (measured by strain at t_tr) is shorter when using dislocation creep, while the differences in Stage 2 timing are smaller and diminish further as extension rate increases.

Table 4 List of simulations investigating different lithospheric ages at spreading onset and half-spreading rates $v_{\mathrm{1}/\mathrm{2}}$. Three rheological combinations are tested (first column). The last two columns list the average strain obtained at the characteristic times of transition t_tr and plate-boundary t_PB, respectively (Sect. 3). The initial strain rate in the plate (${\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}$) is spatially uniform and equal to: 7.47 × 10⁻¹⁷, 1.87 × 10⁻¹⁶, 3.73 × 10⁻¹⁶, 7.47 × 10⁻¹⁶ and 1.87 × 10⁻¹⁵ s⁻¹ respectively for $v_{\mathrm{1}/\mathrm{2}}$ = 0.2, 0.5, 1, 2 and 5 cm yr⁻¹.

Figure 10 Diagnostics as a function of half-extension velocity $v_{\mathrm{1}/\mathrm{2}}$ for two different rheologies, D−Y_500 MPa (red) and $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ (blue). Symbols represent initial plate ages ranging from 10 to 100 Myr. (a) Strain ε(t_tr), used as a proxy for the duration of Stage 1. (b) Difference ε(t_PB)−ε(t_tr) used as a proxy for duration of Stage 2. (c) Mean narrowing rate during Stage 1. (d) Mean upwelling rate during Stage 2. (e) Effective viscosity vertically-averaged over the lithosphere thickness (T < 1500 K) at t_tr and at horizontal position x=600 km. (f) Maximum lateral viscosity contrast (R_visc max.) during the two stages, with a maximum reached around t_PB. (g) Tectonic force TF at time t_tr. In all panels, solid lines connect the markers corresponding to the initial lithospheric age of 50 Myr, used as a reference age.

Extension rate strongly modulates Stage 1 duration. Faster extension increases strain rate within the deforming zone, lowering effective viscosity via both yielding and (where active) dislocation creep, thereby increasing the narrowing rate (Fig. 10c, e). As a result, the total strain required to reach the transition from Stage 1 to Stage 2 decreases systematically with extension rate for both rheologies (Fig. 10a). In contrast, the Stage 2 duration is comparatively insensitive to extension rate for $v_{\mathrm{1}/\mathrm{2}}\gtrsim \mathrm{1}$ cm yr⁻¹ (Fig. 10b), even though upwelling rates increase with extension rate (Fig. 10d). The maximum lateral viscosity contrast increases with extension rate up to ∼ 2 cm yr⁻¹ and then saturates (Fig. 10f), suggesting that once a sufficiently weak plate boundary is formed, it cannot weaken much further at higher extension rates.

Only the slowest forcing conditions fail to localize within the simulated strain window, and this occurs exclusively for the diffusion-plus-yield rheology. For old plates (50–100 Myr) extended at $v_{\mathrm{1}/\mathrm{2}}=\mathrm{0.2}$ cm yr⁻¹, localization remains incomplete by 40 % strain (Table 4). In these cases, thermal diffusion can outweigh extensional thinning during early deformation, leading to net lithosphere thickening and delayed weakening (Appendix D; Figs. D1 and D2). Even in cases that eventually localize at low velocities or for initially young plates, the resulting plate boundary remains relatively wide and thick, and the maximum viscosity contrast remains modest (R_visc≲200 at t_PB; Fig. 10f), indicating a weaker degree of focusing than in faster-extension cases.

The weakening-partition diagnostics clarify how forcing modulates feedbacks. At low extension rates (and in some young-plate cases, early in Stage 1), thermal hardening associated with cooling competes with mechanical weakening, such that net weakening in the creeping lithosphere is dominated by strain-rate effects (note that F_SR can exceed 1 where temperature decreases while strain rate increases, Sect. 3.3). For the dislocation-creep rheology, this strain-rate-driven weakening extends to higher temperatures within the creeping lithosphere (up to ∼ 1300 K), supporting continued narrowing despite thermal hardening. By contrast, for diffusion creep plus yield stress, weakening below the yielding layer remains predominantly thermal; when cooling dominates early, the system remains in a thick yielding-dominated configuration that inhibits rapid focusing, contributing to incomplete localization under the slowest forcing.

5 Discussion

Taken together, our results show that lithospheric-scale strain localization in this extensional setting is not controlled solely by the integrated strength of the plate. It also depends on where weakening develops within the lithosphere and on how rheology promotes deformation focusing. Indeed, rheology modulates the positive feedbacks linking strain narrowing, lithosphere thinning and asthenosphere upwelling, leading to further weakening and an increase in lateral viscosity contrasts, thereby driving the dynamics of plate boundary formation. Three main conclusions emerge. First, simulations showing deformation localization follow a robust two-stage evolution, with progressive narrowing followed by rapid plate thinning. Second, dislocation creep accelerates localization by extending thermo-mechanical weakening into a broader portion of the creeping lithosphere. Third, even a thin stiff layer in the mid-lithosphere can markedly delay plate-boundary formation, implying that simplified weak-plate parameterizations may miss an important control on localization efficiency. In the following, we discuss these points in turn and consider their implications for natural rifting and for large-scale geodynamical models.

5.1 A two-stage pathway to plate-boundary formation

A robust result of this study is that all simulations that successfully localize strain follow the same two-stage pathway to plate-boundary formation, irrespective of the rheological parameterization (Fig. 6). In Stage 1, deformation mainly focuses laterally: the width of the deforming zone decreases rapidly, whereas plate thickness changes only modestly. In Stage 2, the localization pattern shifts to rapid plate thinning driven by focused asthenospheric upwelling, while additional narrowing becomes limited. The transition between these two stages therefore marks the point at which a progressively focused deforming region evolves into a rapidly thinning lithosphere, at the end of which a mature plate boundary is achieved.

This two-stage behavior provides a simple physical framework for interpreting the localization process. Stage 1 is primarily a focusing problem, in which viscosity reduction within the deforming zone, allowed by non-Newtonian rheology, progressively increases the contrast with the surrounding plate interior. Stage 2 is primarily a thinning process, in which upwelling and heating beneath the localized zone further reduce temperature-dependent viscosity and accelerate lithospheric thinning. In that sense, the transition time t_tr represents a dynamical tipping point: beyond it, thermal weakening associated with asthenospheric upwelling becomes dominant and plate-boundary development proceeds rapidly.

A broadly comparable two-stage evolution has been inferred for some natural rift systems: the Atlantic and Australia-Antarctica rifts have probably undergone an early stage of deformation, at (total) extension rates slower than 1 cm yr⁻¹, lasting ∼ 25–50 Myr, followed by an abrupt acceleration of extension over only 2–10 Myr (Brune et al., 2016). This acceleration is interpreted as the consequence of “the rapid decrease of rift strength” in constant-force rifting models (Brune et al., 2016), and of “a positive feedback loop between extension velocity and rift strength loss” in whole-mantle convection models (Ulvrova et al., 2019). Our models do not reproduce such settings directly, particularly because extension is imposed kinematically rather than through a self-consistent far-field force balance.

We performed additional simulations with either an abrupt or a linear temporal increase of extension velocity imposed at side boundaries (Supplement, Sect. S5.1), in which the tectonic force TF remains quasi-constant during Stage 1 before abruptly dropping during Stage 2. As in the constant-extension velocity models, Stage-1 narrowing is dominated by mechanical weakening, while Stage 2 is associated with thermal weakening with rapid TF decrease and fast plate thinning (Figs. 7, 8). We speculate that imposing forces or stresses at side boundaries, instead of imposing velocities, would result in a decrease in plate strength at time t_tr, and an increase in extension velocity, as predicted by rifting models and reconstructed in the natural rifting cases (Brune et al., 2016). Thus, we propose that rifting acceleration in nature may correspond to a geodynamic tipping point with thermal weakening of the lithospheric mantle enhancing plate strength drop, thus further promoting plate extension.

5.2 Importance of enabling thermo-mechanical weakening through dislocation creep

The clearest rheological result of this study is that simulations including dislocation creep localize deformation substantially faster than those with yield stress plus diffusion creep alone (Table 3; Fig. 7a, b). This difference does not arise simply because the plate is weaker in some bulk sense. Rather, dislocation creep changes where and how weakening operates within the lithosphere and enhances the feedbacks leading to deformation localization: faster narrowing or thinning increases mechanical or thermal weakening, causing the decrease in viscosity at the domain center (associated with efficient focusing of deformation in Stage 1 and with rapid central upwelling in Stage 2), thereby increasing the lateral viscosity contrast, which in turn further enhances narrowing or thinning, and so on.

When dislocation creep is included, the transition from yield-dominated to creep-dominated deformation occurs at lower temperatures, so a larger fraction of the lithosphere participates in creep deformation (Fig. 8). In practice, this extends weakening into colder parts of the deforming plate, particularly through the approximate 800–1300 K interval. As a result, dislocation creep broadens the part of the lithosphere that can respond dynamically to both increasing strain rate and increasing temperature, thereby enabling positive feedbacks over a thicker region.

The weakening-partition diagnostics show that this has two important consequences. First, during Stage 1, dislocation creep allows significant strain-rate-driven weakening within the creeping lithosphere, which promotes faster narrowing of the deforming zone. Second, during Stage 2, it permits stronger temperature-driven weakening over a thicker part of the lithosphere, which accelerates plate thinning once focused asthenospheric upwelling is established. In contrast, in diffusion-creep-plus-yield-stress simulations, weakening below the shallow yielding layer is almost entirely thermal, limiting mechanical weakening during Stage 1 and delaying the onset of rapid thinning (Fig. 8).

This difference is also consistent with the analysis of the theoretical viscosity profiles shown in Fig. 9 (Sect. 4.3). Relative to diffusion creep plus yield stress, rheologies including dislocation creep produce the lowest viscosity associated with both mechanical weakening during a Stage 1-type strain-rate increase and stronger thermal weakening during a Stage 2-type thinning trajectory. In addition, dislocation creep extends both effects to lower temperatures and shallower depths. The main consequence is therefore not simply a lower viscosity, but a stronger coupling between deformation focusing, asthenospheric upwelling, and further weakening.

Our results further suggest that, if deformation is driven by far-field kinematics, as in our modeling conditions, overall lithospheric strength exerts only a second-order control on localization efficiency. For example, Simulation D−Y_200 MPa localizes more slowly than simulations $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ and $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$, even though the former has a lower tectonic force during Stage 1 (Fig. 7c). Extra simulations are performed with depth-dependent yield stress, keeping the peak strength of 500 MPa at the same depth as in the constant σ_Y simulations (Appendix E). These complementary simulations exhibit comparable localization characteristic times despite the lower strength at the plate surface, which supports the previous conclusion: reducing near-surface strength does not substantially change localization times if a deeper mechanically important layer remains. Localization is therefore not controlled solely by average plate weakness, but by rheology-dependent feedbacks and by the depth distribution of strength within the lithosphere. Nevertheless, we acknowledge that the situation would be different if lithospheric deformation were instead driven by far-field constant stresses or forces, especially in case of a low to moderate tectonic rate. In such a setting, imposing a yield stress of 500 MPa instead of 200 MPa might significantly modify the evolution of lithospheric spreading.

The aforementioned feedbacks between the temperature and strain rate fields enabled by the activation of dislocation creep are modeled in our experiments while shear heating is not included. Shear heating has been proposed to enhance shear localization by thermal runaway (e.g. Fleitout and Froidevaux, 1980; Kaus and Podladchikov, 2006; Thielmann and Kaus, 2012; Kiss et al., 2019). Using a posteriori calculations, we estimate a maximum heat dissipation rate for Simulation ref-1 ($D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$) larger than ∼ 10⁻⁶ W m⁻³ in a ∼ 15 km thick layer directly below the yield-creep transition (Supplement, Sect. S5.3). Considering that a heat dissipation rate of order 10⁻⁵ W m⁻³ can lead to a temperature increase close to ∼ 50 K (Kiss et al., 2020; Arcay et al., 2023), we therefore expect that the temperature rise due to shear heating would remain lower than 50 K in our simulations. Therefore, including shear heating could foster thermal weakening in our experiments, further enhancing plate weakening and strain localization, but to a low-to-moderate extent. This effect is predicted to be especially emphasized when using (LT-)dislocation creep for which the temperature-dependent creep layer is the thickest (Supplement, Sect. S5.3).

Taken together, these results show why yield-stress rheologies are not dynamically equivalent to rheologies that include dislocation creep at moderate to high temperatures. Dislocation creep accelerates localization because it expands the depth and temperature range over which both mechanical and thermal weakening can operate, thereby strengthening the feedbacks required to transform distributed extension into a localized plate boundary.

5.3 Implementing yield-stress or LT-dislocation creep as ductile strength-limiters

In these simulations, yield stress plays two distinct roles. First, it acts as a proxy for shallow pseudo-brittle weakening, enabling deformation of the very cold upper lithosphere. That role is essential in the present extension setting: without weakening of the shallow plate, deformation remains distributed and of low amplitude, even when dislocation creep is included. Some mechanism that weakens the cold upper lithosphere is therefore required in addition to any deeper strength-limiter.

Second, in rheologies that do not include low-temperature dislocation creep, yield stress also limits the effective strength in a part of the deeper lithosphere. This is evident from the comparison between simulations with and without an imposed thermal cutoff for yielding. When LT–HT dislocation creep is present, restricting yielding to temperatures below 950 K has no effect on localization timescales. By contrast, when only HT-dislocation creep is included, the same restriction markedly delays localization or prevents it within the simulated strain window (Table 3). In that sense, when low-temperature dislocation creep is absent, the yield stress partly substitutes for a strength-limiter within the approximate 800–1100 K interval, that is, across the temperature range over which deformation transitions from dislocation creep to yielding (Fig. 8). This interval overlaps at least partly with the ductile domain inferred for olivine at temperatures greater than about 800–900 K from laboratory and theoretical studies on melting and homologous temperatures (Hirschmann, 2000; Wang, 2016; Demouchy et al., 2023). It also overlaps with the mantle brittle-ductile transition expected up to approximately 600 °C (900–1000 K), as inferred from seismicity in the oceanic lithosphere (e.g. Engeln et al., 1986; Abercrombie and Ekström, 2001; McKenzie et al., 2005), and may extend to 900 °C (1200–1300 K) under hydrous mantle conditions (Kohli et al., 2021), making the physical interpretation of such stress caps intrinsically ambiguous because deformation around the brittle-ductile transition is likely governed by multiple poorly constrained processes (Kohlstedt et al., 1995; Meyer et al., 2019).

Our simulations further show that lowering stresses of the shallowest part of the lithosphere does neither accelerate nor modify the breakup process (Appendix E), whereas even a stiff layer less than 5 km thick within the mid-lithosphere can strongly delay Stage-1 narrowing and therefore delay plate-boundary formation (Sect. 4.3, Fig. C1; Appendix C). This layer forms a thin but mechanically important bottleneck between the shallow yielding domain and the warmer, weaker creeping mantle. As long as this bottleneck remains strong, it dampens the feedbacks between deformation focusing, plate thinning, and further weakening, with strain rate increase lessened in the whole lithospheric section, hence milder ensuing viscosity reduction (Fig. Cb, c). In that sense, the relevant control on localization is not simply the weakness of the shallowest lithosphere, but the peak strength of the lithosphere, and whether the stiff layer can itself weaken. This underlines the need to cap the lithosphere strength beyond the brittle realm, in a temperature range corresponding approximately to 800–1000 K, i.e. the thermal range within the mid-lithosphere where the activation of HT-dislocation creep would otherwise lead to very high viscosities.

Our results therefore support the view that yield stress in geodynamical models should not be interpreted only as a shallow brittle proxy. In practice, it may also act as a first-order proxy for deeper ductile strength limitation when low-temperature plasticity is not represented explicitly. Following Tackley (2000) and Van Heck and Tackley (2008), we therefore suggest distinguishing between a shallow pseudo-brittle stress cap and a deeper parameterization representing low-temperature ductile strength limitation, even if the latter is not dynamically equivalent to LT-dislocation creep itself (Sect. 5.2).

Simulations with (LT-)-HT dislocation creep exhibit significantly lower tectonic forces TF than simulations with only yielding and diffusion creep. In addition, using lower yield stress values also reduce the tectonic forces needed to maintain extension (Fig. 7c). Nevertheless, initial TF values of approximately 20–40 TN m⁻¹ remain on the high side compared to estimates proposed for natural rifting (∼ 1 TN m⁻¹, Brune et al., 2023), ridge push (∼ 2–3 TN m⁻¹, Parsons and Richter, 1980), but would be in agreement with estimation of slab pull (∼ 10–50 TN m⁻¹, e.g., Toth and Gurnis, 1998; Conrad and Lithgow-Bertelloni, 2002; Schellart, 2004; Wu et al., 2008). In large-scale mantle convection models, chosen yield stresses are low (≲ 200–300 MPa, e.g. Richards et al., 2001; Crameri and Tackley, 2014; Mallard et al., 2016; Coltice et al., 2019), or with a low yield-stress increase with pressure (often ≲ 0.2, e.g. Moresi and Solomatov, 1998; Korenaga, 2010; Crameri et al., 2012; Nakagawa and Iwamori, 2017) compared to values derived from friction coefficients inferred from laboratory experiments (∼ 0.6, e.g. Byerlee, 1978). These choices generally result in lithospheric strengths lower than laboratory-based estimates (> 500 MPa, e.g. Kohlstedt et al., 1995; Karato, 2008; Burov, 2011), but they allow models to reproduce plate-like behavior and reorganizations (Janin et al., 2025). Our results suggest that including low-(to high-)temperature dislocation creep provides an alternative mechanism to act as a strength-limiter in the lithospheric mantle, allowing a somewhat higher stress-cap, more consistent with deformation experiments, while still enabling plate weakening and strain localization. Other plausible mechanisms have also been proposed to generate plate-like behavior compatible with higher yield stress values, such as lateral strength heterogeneities associated with cratonic continental blocks (Rolf et al., 2012).

5.4 Implications for geodynamical models and limitations of the present study

A central implication of this study is that rheologies combining diffusion creep with a simple yield-stress cap (e.g. Moresi and Solomatov, 1998; Trompert and Hansen, 1998; Coltice et al., 2019) may overestimate the time required to localize strain and form new plate boundaries. In our simulations, adding dislocation creep substantially accelerates localization relative to rheology D−Y because it extends the thermal range of weakening in the creeping lithosphere and strengthens the feedbacks that link strain-rate increase, asthenospheric upwelling, and further viscosity reduction. This suggests that large-scale geodynamical models may localize too slowly if they omit this rheological contribution. Yet, including dislocation creep in whole-mantle convection simulations will have other dynamical consequences, such as a sub-plate asthenosphere weakening scaling with surface plate velocities (e.g. Patočka et al., 2024), or an increased mechanical decoupling at lithosphere-asthenosphere boundary that alters mantle flow and surface deformation pattern (Semple and Lenardic, 2021; Arnould et al., 2023).

More broadly, our results highlight a limitation of vertically-uniform stress-limiters. A low and constant yield-stress is not dynamically equivalent to a rheological structure in which shallow pseudo-brittle weakening coexists with a deeper creeping layer, since the latter can undergo strong thermo-mechanical weakening. The distinction matters because localization is sensitive not only to the integrated plate strength, but also to the depth distribution of strength and weakening. In particular, a thin stiff mid-lithosphere layer can act as a bottleneck that delays localization even when the plate is weak in a vertically averaged sense. Models developed to accurately estimate the absolute timescale of localization should distinguish a pseudo-brittle yield cap from a ductile strength limiter, either self-consistently arising from low-temperature dislocation creep (e.g. Mei et al., 2010; Jain et al., 2017; Gouriet et al., 2019; Demouchy et al., 2023; Warren and Hansen, 2023), or implemented with a so-called “ductile yield stress” as in Tackley (2000) and Van Heck and Tackley (2008).

A second important aspect of our study is the introduction of weakening-partition diagnostics, that provide a general framework for analyzing the rheological control on deformation localization. By separating strain-rate-driven and temperature-driven contributions to viscosity reduction, they allow for identifying when localization is controlled primarily either by mechanical weakening or by thermal feedbacks. This weakening-partition diagnostics should be transferable to other geodynamic settings, including subduction initiation (e.g. Gurnis et al., 2004; Billen and Hirth, 2005; Ueda et al., 2008; Zhong and Li, 2019; Zhang et al., 2021; Arcay et al., 2023) or plume-lithosphere interaction (e.g. Brune et al., 2013; Agrusta et al., 2015). It can also be extended to rheologies with additional state dependencies such as grain size, damage, or composition (e.g. Bercovici and Ricard, 2013; Fuchs and Becker, 2021; Dannberg et al., 2025).

The main limitation of the present study is the mantle-only configuration. By excluding crustal layering, we isolate mantle rheological controls but do not attempt to reproduce the full mechanical complexity of natural rift systems. Crust-mantle coupling, lithological discontinuities, and additional weakening processes such as strain softening, grain-size evolution, and damage are all expected to promote a more rapid localization in a “narrow-rift” setting featuring a strong lower crust (e.g. Allemand and Brun, 1991; Buck, 1991; Huismans and Beaumont, 2003; Gueydan et al., 2008; Gueydan and Précigout, 2014; Tetreault and Buiter, 2018; Heckenbach et al., 2021). Indeed, for rifting models under total extension around 1 cm yr⁻¹ (e.g. Brune et al., 2014; Gueydan and Précigout, 2014; Chenin et al., 2018) the corresponding transition time ( t_tr) and break-up time (t_PB) are approximately 5 and 9–16 Myr, respectively, compared to 14 and 26 Myr in our fastest model (rheology $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$). In that sense, the absolute times reported in this study are best interpreted as end-member values for a simplified mantle system, rather than as direct predictions for crust-bearing rifts. The interplay between low-temperature dislocation creep and additional weakening mechanisms will depend on the details of the rheological parameterizations and their evolution laws, but is beyond the scope of the present study (e.g. Precigout et al., 2007; Bercovici et al., 2015; Li and Gurnis, 2024). The presence of a crust also questions the existence of a brittle mantle (Burov and Watts, 2006), since the transition between pseudo-brittle yield-stress and LT-dislocation creep is expected around 800 K for a background strain rate of 10⁻¹⁶–10⁻¹⁵ s⁻¹ (Fig. 8), which corresponds to 27 km depth for a 50 Myr old plate, i.e. shallower than a continental Moho (∼ 30 km deep), but deeper than an oceanic one (∼ 4–10 km depth).

Neglecting elasticity is another limitation. This assumption may first lead to an overestimation (up to a factor of three) of shallow stresses in our models (Patočka et al., 2019). In our visco-plastic experiments (“plastic” being here a synonym for “brittle”), the central peak in shallow stresses (at 5 km depth, Supplement, Sect. S5.2) is located in the area where strain localizes, and where free-surface elevation is noticeable. If reducing shallow stresses by a factor 4 (as shown by additional tests with depth-dependent yield stress) indeed leads to smaller central depression (factor ∼ 3), it does not modify the two-stage weakening scenario or the timing of incipient plate-boundary formation (see Appendix E). This supports the conclusion that the localization mechanisms identified in our study are controlled primarily by the rheological structure of the mid- and lower lithosphere, rather than by shallow stress alone.

Neglecting elasticity by assuming a viscoplastic rheology also affects the prediction of brittle deformation and fault development. Previous studies comparing viscoplastic and visco-elasto-plastic formulations have shown that, although first-order stress patterns may be comparable, elasticity implementation will change fault geometry, spatial distribution, interaction and rotation, as well as the rift-topography evolution (e.g. Olive et al., 2016). In particular, our yield-stress formulation produces distributed viscoplastic flow over a finite volume rather than localized planar faulting, and therefore does not capture the detailed geometry of brittle shear-band development. Viscoplastic formulations may exhibit mesh sensitivity and numerical convergence issues (Duretz et al., 2020, 2021). Simulating a constant yield stress prevents the formation of shear bands in the pseudo-brittle layer (Watremez et al., 2013). For a moderately fast extension, the pseudo-brittle layer deforms by pure shear in a layered lithosphere if the viscosity underneath the upper brittle layer is high (∼ 10²³ Pa s, Huismans and Beaumont, 2005), which may explain the absence of shear bands in our models with a moderate yield stress increase with depth (viscosity ∼ 10²³ Pa s at the yield-creep transition for a 500 MPa yield stress at $\dot{\mathit{\epsilon }}={\mathrm{10}}^{-\mathrm{15}}$ s⁻¹, Figs. 9a, c and C1). More physically grounded approaches, including damage-based rheologies and visco-elasto-plastic formulations accounting for weakening by microcracks growth, may therefore improve the representation of shallow deformation (Petit et al., 2024), but remain computationnaly demanding. We acknowledge that treating the shallow brittle layer more realistically could significantly modify the deformation pattern at shallow depths, but expect the large-scale feedbacks of lithospheric weakening documented in this study to remain broadly relevant.

6 Conclusions

Our 2-D thermo-mechanical extension experiments show that lithospheric-scale strain localization is controlled by the interplay between rheology and self-reinforcing thermo-mechanical feedbacks. In all simulations that successfully localize deformation, plate-boundary formation follows a robust two-stage pathway: an initial stage of progressive narrowing of the deforming zone, followed by a second stage of rapid plate thinning driven by focused asthenospheric upwelling. This two-stage evolution provides a simple physical framework for understanding how distributed extension evolves into a localized plate boundary.

A central result is that dislocation creep substantially accelerates localization relative to yield stress plus diffusion creep alone. This effect arises not simply because viscosity is lower overall, but because dislocation creep extends weakening into a broader and colder part of the creeping lithosphere. In doing so, it enables both strain-rate-dominated weakening during Stage 1 and stronger temperature-driven weakening during Stage 2, thereby enhancing the feedbacks that promote localization. By contrast, rheologies based only on diffusion creep plus a yield-stress cap localize more slowly because weakening in the deeper lithosphere is more restricted.

Our results also show that localization depends on the depth distribution of weakening and strength. In particular, even a thin stiff layer within the mid-lithosphere can markedly delay Stage-1 narrowing and hence delay plate-boundary formation. This suggests that low-temperature dislocation creep acts as a ductile stress-limiter at temperatures of about 800–1100 K, similarly to the effect of common yield-stress parameterizations.

More broadly, this study introduces diagnostics that partition viscosity reduction into temperature-driven and strain-rate-driven contributions. These diagnostics directly link weakening to the evolving thermal and kinematic state of the system, clarifying why some rheological combinations localize efficiently whereas others do not. They also provide a transferable framework for analysing weakening processes in other geodynamic settings (e.g. subduction initiation and plume-lithosphere interaction) and in rheologies with additional state dependencies.

Taken together, these results identify dislocation creep as a key ingredient for efficient plate-boundary formation, and show that what matters most is not simply overall plate strength, but the spatial extent of thermo-mechanical weakening within the lithospheric mantle.

Appendix A: End-members of deformation pattern (distributed vs. localized)

The two end-member modes of lithospheric spreading, namely, distributed deformation and localized deformation, described in Sect. 4.1, are assigned at simulation end (40 % of strain, corresponding to 500 km total surface extension) based on two main distinct diagnostics: plateness and lateral viscosity contrasts (Fig. A1). In some cases we observe an intermediate case, that is, incomplete localization, for which the deformation pattern has not reached any quasi-steady state at 40 % of strain. For instance, the plate thickness h_plate has not reached any plateau after 40 % of strain (see Fig. 7b).

Figure A1 (a) Plateness P and (b) lateral viscosity contrast R_visc obtained after 40 % strain for various extension rates, initial plate age, and rheologies. The end-member simulations (distributed deformation, localized plate boundary) are depicted in red and blue, respectively. The incomplete localization case in yellow is intermediate between the two end-members.

Appendix B: Transition time t_tr calculation

For all “localized” simulations (loc. PB, Tables 3, 4), we estimate the transition time t_tr from the co-evolution of plate thickness h_plate(t) and width of the enhanced deformed zone w(t). The transition time t_tr is defined as the point furthest from the line connecting the lithospheric structures at t₀ and t_PB, to the evolution curve h_plate(t)−w(t) (Fig. B1). This point corresponds to the knee in the curve, i.e. a transition between the first stage of extension during which deformation narrowing dominates (strong decrease in w while h_plate only slightly lessens) and the second stage, in which the asthenospheric upwelling dominates (moderate w reduction but strong lithospheric thinning and h_plate reduction). Defining a transition with such a methodology was for example used to estimate the threshold in grain orientation spread to discriminate between recrystallized and relict grains of quartz from a given distribution of intracrystalline lattice orientation in a natural sample (e.g. Cross et al., 2017).

Figure B1 Definition of the transition time t_tr illustrated using the reference experiment, Simulation ref-1 (rheology $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$). The transition time t_tr is determined as the time at which the distance between the co-evolution of the deformed-zone width w and the plate thickness h_plate (blue line), and the straight line between the end-member states at t=0 and t=t_PB (black line) is maximum.

Appendix C: Rheological feedbacks in the presence of a stiff mid-lithosphere layer

Figure C1 Vertical profiles of (a) temperature, (b) strain rate and (c) effective viscosity at x = 600 km at three different times (0, 4 and 13 Myr) for two rheologies $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ (pink) and $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$ (green) differing only by the temperature cut-off imposed for the yielding rheology. The thin stiff layer forming for rheology $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$ is highlighted by the red hatched area Arrows depict the time evolution.

Simulation $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$ investigates the effect of a thermal limitation imposed on the yielding realm. In this experiment, the yield stress rheology is restricted to temperatures lower than 950 K (see Sect. 2.2.3) and is compared to Simulation $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ in which no thermal limitation is imposed (Table 3). Since brittle deformation in the lithospheric mantle has been argued to end at a temperature close to 900 K (e.g. Engeln et al., 1986; Abercrombie and Ekström, 2001; McKenzie et al., 2005), limiting yielding to a temperature around 950 K allows for modeling the lithosphere behavior if the brittle-ductile transition was controlled by a maximum temperature. Simulations performed without and with a temperature-capped yield rheology are compared using vertical profiles of temperature, strain rate and viscosity sampled at the center of the simulation box and computed at three different times (0, 4 and 13 Myr) at different stages of localization (Fig. C1).

At the start of simulations, the temperature of the yield-creep transition is close to 1000 K in Simulation $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ (Sect. 4.3) and is located at ∼ 41 km depth, while the  950 K isotherm is at 38 km depth.  As a consequence, if the yield stress is thermally-capped, dislocation creep is activated in a 3 km-thick layer below the yield boundary (38–41 km depth). In this layer, viscosities reach the highest values, locally exceeding 10²⁴ Pa s. After 4 Myr of lithospheric extension (end of Stage 1 in Simulation $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$, the strain rate has increased with respect to the initial state by a factor of ∼ 2.3 for the thermally-capped yield rheology, while this increase is doubled in Simulation $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}$. As a result, the viscosity in the yielding realm is decreased by a factor 2 for the temperature-capped yield, instead of 3 without the thermal limitation in yield.

The difference between the two simulations in the average lithospheric strain rate at the box center is amplified through time. As Stage 2 proceeds for Simulation $D-d_{\mathrm{HT}}-Y_{\mathrm{500}\mathrm{MPa}}^{\mathrm{950}\mathrm{K}}$, a strong lithospheric thinning at 13 Myr (1200 K increase at 5 km depth relative to the initial state) adds thermal weakening to mechanical weakening (200-fold increase in strain rate), which reduces viscosity by a factor ∼ 2400 with respect to the initial state. The corresponding viscosity weakening is limited to a ∼ 170-fold decrease for the thermally capped yield rheology.

Appendix D: Varying the extension rate and initial thermal plate age

Figure D1 (This figure uses the same representation as Fig. 8 in the main text.) Depth-time evolution of the relative contributions to weakening F_T and F_SR calculated at box center x = 600 km with rheology $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$. Rows correspond to varying initial plate age and columns to varying half-extension velocity $v_{\mathrm{1}/\mathrm{2}}$. The central panel represents the reference setup (initial plate age of 50 Myr, half-extension rate of 1 cm yr⁻¹). In each panel, time varies from 0 to t_PB (which differs from one panel to another), and the transition time t_tr is indicated by a vertical dashed line. The diagnostics are calculated only in regions experiencing weakening (W>0), hardening or stable regions (W≤0) being displayed in white. The color-scale is saturated between 0 and 1. Zones of intense weakening, where the weakening rate normalized by the initial strain rate in the plate $\frac{W}{{\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}}$ exceeds 30, are outlined in pink. Boundaries between dominant deformation mechanisms are outlined with dashed dark-blue lines as in Fig. 4, and isotherms from 800 to 1500 K are displayed in light blue.

Figures D1 and D2 display the thermomechanical evolution modeled at the box center obtained for rheologies $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$ and D−Y_500 MPa, respectively, when the initial plate age and half-extension rate are varied (Sect. 4.4). At low velocity, the early extension stage is characterized by the cooling of the lithosphere (thermal hardening). As a consequence, the material weakening results from mechanical weakening only even in the creeping part of the lithosphere (dark green, Fig. D1a, b, d, g). In particular, Fig. D2c and f illustrate a scenario of incomplete localization (Sect. 4.1, Table 4), where complete localization of deformation is not achieved at the end of the simulation when strain reaches 40 %.

Figure D2 (This figure uses the same representation as Fig. 8 in the main text.) Depth-time evolution of the relative contributions to weakening F_T and F_SR calculated at box center x = 600 km with rheology D−Y_500 MPa. Panels are organized with rows corresponding to varying initial plate age and columns corresponding to varying half-extension velocity $v_{\mathrm{1}/\mathrm{2}}$. The central panel represents the reference setup (initial plate age of 50 Myr, half-extension rate of 1 cm yr⁻¹). In each panel, time varies from 0 to t_PB (which differs between panels), and the transition time t_tr is indicated by a vertical dashed line, except in panels (c)–(f) where localization is not achieved and the results are shown for the total duration of the simulation. The diagnostics are calculated only in regions in weakening (W>0), with hardening or stable regions (W≤0) shown in white. The colorscale is capped between 0 and 1. Zones of intense weakening, where the weakening rate normalized by the initial strain rate in the plate $\frac{W}{{\dot{\mathit{\epsilon }}}_{t_{\mathrm{0}}}}$ exceeds 30, are highlighted in pink. Boundaries between dominant deformation mechanisms are marked by dashed dark-blue lines as in Fig. 4, and isotherms from 800 to 1500 K are displayed in light blue.

Figure D3 provides a global illustration of the two-stage evolution in simulations for various initial plate ages and half-extension velocities (Sect. 4.4). The curves for the same initial plate age are almost superimposed, except for low velocity $v_{\mathrm{1}/\mathrm{2}}=\mathrm{0.2}$ cm yr⁻¹ which exhibits thickening during the early Stage 1. The “incomplete localization” scenario of the simulation with rheology D−Y_500 MPa for an initial 50 Myr-old plate and $v_{\mathrm{1}/\mathrm{2}}$ = 0.2 cm yr⁻¹ (Sect. 4.1, Table 4) is illustrated by the pink curve in Fig. D3b: at the end of the simulation (strain of 40 %), the plate is still thicker than 700 km and the deformed zone wider than 150 km.

Figure D3 Co-evolution of the minimum plate thickness h_plate(t) as a function of the width of the most deformed zone for two rheological combinations: (a) $D-d_{\text{LT–HT}}-Y_{\mathrm{500}\mathrm{MPa}}$. (b) D−Y_500 MPa. Curves are shown for varying initial plate age (indicated by line-styles and half-spreading velocity (indicated by colors). The transition times t_tr between Stages 1 and 2 are marked by dots.

Appendix E: Simulations with a depth-dependent yield stress

We perform two simulations with a depth-dependent yield stress for the following two creep rheologies: (i) diffusion creep (D), (ii) diffusion creep plus low- and high-temperature dislocation creep (D−d_LT–HT). We define a depth-dependent yield stress:

$$\begin{array}{}\text{(E1)}& {\mathit{\sigma }}_y\left(z\right)={\mathit{\sigma }}_{\mathrm{0}}+z\cdot \mathrm{\Delta }{\mathit{\sigma }}_z\end{array}$$

with

$$\begin{array}{}\text{(E2)}& \mathrm{\Delta }{\mathit{\sigma }}_z=\mathit{\rho }\cdot g\cdot \mathit{\mu }\end{array}$$

where z is depth, σ₀ = 50 MPa is the yield stress at the surface (cohesion), Δσ_z (Pa m⁻¹) is the yield stress gradient and μ is the yield stress increase with pressure related to the friction coefficient, f_s. The parameter μ corresponds to the ratio between the horizontal tectonic deviatoric stress and the lithostatic pressure (neglecting pore fluid pressure), while the friction coefficient f_s is the ratio between the shear stress and the normal stress acting on a plane fault (e.g., Doin and Henry, 2001; Turcotte and Schubert, 2002). The yield stress gradients are chosen to ensure that the transition of deformation mechanism from yield to creep corresponds to a stress of approximately 500 MPa and occurs at the same depth as in the case of a constant yield stress with depth. This transition depth depends on the creep law considered. Consequently, for the following creep rheologies: diffusion creep (D) and diffusion creep plus low- to high-temperature dislocation creep (D−d_LT–HT), one may get: Δσ_z(D) = 6.181 MPa km⁻¹, μ(D) ≈ 0.19 and Δσ_z(D−d_LT–HT) = 16.423 MPa km⁻¹, μ(D−d_LT–HT) ≈ 0.51, respectively (Fig. E1), keeping unchanged the depth of the yield-creep transition with respect to a constant yield stress ensures to maintain the rheological structure of the underneath creeping layer in the lithosphere.

Figure E1 Strength profile computed for a uniform strain rate $\dot{\mathit{\epsilon }}$ = 3 × 10⁻¹⁵ s⁻¹ and a temperature profile corresponding to a 50 Myr-old lithosphere. The red and dark blue curves are respectively: diffusion, and diffusion plus low- and high-temperature dislocation creep rheologies, in combination with a depth-dependent yield stress (dashed lines), or with a constant yield stress of 200 or 500 MPa (orange dotted and light blue solid lines, respectively). The yield-creep transition for the depth-dependent yield stress is reached at a maximum strength of 500 MPa, which is consistent to the yield-creep transition depth for constant yield stress of 500 MPa combinations.

Figure E2 Temporal evolution (and corresponding cumulative strain in %) of diagnostics for two sets of simulations corresponding to two creep viscosities (diffusion creep and diffusion creep with LT–HT dislocation creep). In each set, the yield stress is either constant, and set to 200 or 500 MPa, or increasing with depth. The color and style coding used to caption simulations is the same as in Fig. E1. Note that this figure is a similar representation of diagnostics as Fig. 7) in the main text. Simulations have the same initial plate age (50 Myr) and extension velocity ($v_{\mathrm{1}/\mathrm{2}}$ = 1 cm yr⁻¹). (a) Width of the enhanced deformed zone w, (b) plate thickness h_plate at x = 600 km, (c) tectonic force TF, (d) plateness P, (e) lateral viscosity contrast R_visc, and (f) maximum weakening rate W_max over the whole lithosphere thickness (T < 1500 K) at x = 600 km.

The two chosen stress-increase with depth are higher than in commonly used depth-dependent yield stress formulations in large mantle convection simulations (usually lower than 1.1 MPa km⁻¹, e.g. Rolf and Tackley, 2011; Coltice et al., 2019; Arnould et al., 2023), and corresponds to relatively high friction coefficients (typically lower than 0.2 in convection experiments, e.g. Moresi and Solomatov, 1998; Korenaga, 2010; Crameri et al., 2012; Nakagawa and Iwamori, 2017) and closer to experiment values (∼ 0.6, e.g. Byerlee, 1978).

Figure E2 shows that the transition time t_tr, the plate boundary time t_PB and the overall scenario of localization obtained for a depth-dependent yield stress are very close to the reference cases performed with a constant yield stress, for the two investigated creep rheologies. Tectonic force and the lateral viscosity contrasts are lower and higher, respectively, when comparing a depth-dependent yield stress to a constant yield stress (Fig. E2c and e, at initial state and at time t_tr, respectively). Moreover, the free-surface topography is not significantly modified when the yield stress is depth-dependent instead of being constant, except in the noticeable narrow area where the plate boundary forms at the very end of the simulation (Supplement, Sect. S5.2).