Progressive weakening within the overriding plate during dual inward dipping subduction

The evolution of dual inward dipping subduction (DIDS) is crucial to understand multiple slab interaction. Yet, how DIDS influences the thermo-mechanical behaviour of the overriding plate remains unclear, as previous DIDS investigations all applied a compositional or Newtonian rheology that excludes temperature dependency. Here we apply a composite rheology, including temperature dependent creep deformation mechanisms, in 2-D thermo-mechanical numerical modelling to investigate how DIDS modifies the rheological structure of the overriding plate. Three variables on plate sizes are investigated to understand what may control the maximum degree of plate weakening. We find that reducing the initial length or initial thickness of the overriding plate and increasing the initial thickness of the subducting plate can enhance the viscosity reduction within the overriding plate. The progressive weakening can result in a variety of stretching states ranging from 1) little or no litho-sphere thinning and extension ( < 5% accumulation of strain), to 2) limited thermal lithosphere thinning ( < 30% accumulation of strain), and 3) localised rifting followed by spreading extension. Compared with single sided subduction, DIDS further reduces the magnitude of viscosity in the overriding plate. It does this by creating a dynamic fixed boundary condition for the overriding plate and forming a stronger upwelling mantle flow, both of which promote the progressive weakening in the overriding plate. The result implies that these generic DIDS effects are important aspects to consider to understand extension developed in natural DIDS cases. We also demonstrate that both temperature dependent creep rheologies and yielding deformation mechanism contribute significantly to the continuous viscosity reduction. The finding may also have a broader implication for more general processes that involve plate scale weakening, strain localisation and the formation of new plate boundaries.


Introduction
Subduction can pose a fundamental tectonic overprint on the overriding plate by generating a volcanic arc (Perfit et al., 1980;Straub et al., 2020), back-arc basin (Uyeda, 1981), orogeny (Faccenna et al., 2021), or even continental breakup (Dal Zilio et al., 2018).Most subduction zones involve only one subducting slab.Here we consider multiple subducting slabs, in particular dual inward dipping subduction (DIDS).DIDS occurs when the overriding plate is decoupled with two subducting slabs dipping towards each other.It is one of the four most commonly described subduction zones with multiple slabs, i.e., inward-dipping, same-dip, outward-dipping and oppositely dipping adjacent subduction zones (Holt et al., 2017;Király et al., 2021).
Seismic tomography shows that DIDS exists at the Caribbean plate between the Cocos slab and Lesser-Antilles subduction zone (Van Benthem et al., 2013), and in South-East Asia between the Philippine and the Sumatra subduction (Hall and Spakman, 2015;Huang et al., 2015;Maruyama et al., 2007).In combination with seismic tomography, recent plate reconstructions have made it more evident that DIDS could have developed in some regions in the past (Faccenna et al., 2010;Hall and Spakman, 2015) constrained by suture zone petrology demonstrating the existence of paleo-subduction.One extinct DIDS example is the North China Craton.Suture zone studies reveal that multiple inward dipping subduction may have surrounded the North China Craton from Early Paleozoic to Tertiary (Santosh, 2010;Windley et al., 2010).
Global Strain Rate Map (Kreemer et al., 2014) shows that high strain rate regions in the overriding plate is often wider and spatially more complex in multiple slab subduction zones, including the DIDS cases.Besides, crustal scale seismic surveys in the Caribbean Sea have found a good amount of extensional basins spreading widely across the overriding plate since the establishment of dual inward dipping subduction at ~70 Ma (e.g.Boschman et al., 2014;Braszus et al., 2021).Despite these observations, the role dual inward dipping subduction may play to form the high strain rate region, extensional basins or spreading centres remain unclear.
Numerical investigations have been conducted to understand the dynamics of dual inward dipping subduction.Research shows that the initial slab dip of the subducting plate affects the upper mantle dynamic pressure between the convergent slabs and stress state within the overriding plate (Holt et al., 2017).Varying the distance between the trenches, convergence rate, and asymmetry of subducting plates can alter the topography of the overriding plate (Dasgupta and Mandal, 2018).The size of the overriding plate, the viscosity ratio of overriding plate over asthenosphere and lower mantle over upper mantle are tested to investigate their impact upon the slab geometry and the magnitude of mantle upwelling flow underlying the overriding plate (Lyu et al., 2019).
These pioneering investigations show that dual inward dipping subduction can generate a variety of upper mantle flow patterns which regulate the stress state and topography of the overriding plate.However, these models all applied a constant viscosity or Newtonian rheology that excludes temperature dependency for both plates and convective mantle flow during simulation.Mineral deformation experiments indicate that viscosity varies as a function of multiple parameters, e.g., temperature, pressure, stress, strain rate etc. (Bürgmann and Dresen, 2008;Burov, 2011;Hirth and Kohlstedt, 2003;Karato, 2010;Lynch and Morgan, 1987).Thus, previous dual inward dipping subduction models with simplified rheology were unable to fully reflect the weakening process within the overriding plate, e.g., high strain rate spreading centres in the back-arc region.
Single sided subduction models considering temperature dependent composite rheology, e.g., dislocation creep, diffusion creep, yielding etc., has improved our understanding of subduction's impact upon the overriding plate (Alsaif et al., 2020;Bessat et al., 2020;Čížková and Bina, 2013;Garel et al., 2020Garel et al., , 2014;;Schliffke et al., 2022;Suchoy et al., 2021).However, it remains much less explored in terms of evaluating different deformation mechanism's contribution to the weakening process observed in the overriding plate, or how different deformation mechanisms interplay with each other during subduction.
In this research, a series of 2-D thermo-mechanical models incorporating temperature dependent rheology laws are run to investigate how dual inward dipping subduction may differ from single sided subduction and previous DIDS models in deforming the overriding plate.We also identify and quantify different dominant deformation mechanisms' contribution to induce progressive weakening in the overriding plate and investigate the interplay among different deformation mechanisms.

Methods
We ran the thermally-driven dual inward dipping subduction models using the code Fluidity (Davies et al., 2011;Kramer et al., 2012), a finiteelement control-volume computational modelling framework.Our model incorporated an adaptive mesh that can accurately capture the dynamic changes in velocity, temperature, and viscosity etc., with a maximum resolution of 0.4 km.

Governing equations and rheology setup
Under the Boussinesq approximation (McKenzie et al., 1974), the equations governing thermally driven subduction process are derived from conservation of mass, momentum, and energy, for an incompressible Stokes flow in which u, g, σ, T, κ are the velocity, gravity, stress, temperature, and thermal diffusivity, respectively (Table 1).In particular, the full stress tensor σ ij consists of deviatoric and lithostatic components via where τ ij represents the deviatoric stress tensor, p the dynamic pressure, and δ ij the Kronecker delta function.
The deviatoric stress tensor and strain rate tensor εij are related according to with μ the viscosity.The density difference due to temperature is defined as  (Hirth andKohlstedt, 2003, 1995a;Ranalli, 1995).The UM and LM stands for "upper mantle" and "lower mantle," respectively.b The activation parameters and stress-dependent exponent used for dislocation creep are in agreement with dry olivine deformation experiments (Hirth and Kohlstedt, 1995b).c The parameterisation (based on Kameyama et al., 1999) makes Peierls creep tend to be weaker than yielding in the upper mantle, thus enabling trench retreat and creating richer slab morphology in the upper mantle (Garel et al., 2014).d A very high maximum yield strength value is used here to ensure that yielding only dominates at the depth of crustal scale.A friction coefficient of 0.2 is following numerical models (Garel et al., 2014;Gülcher et al., 2020), and it is intermediate between lower values of previous subduction models (Crameri et al., 2012;Di Giuseppe et al., 2008) and the actual friction coefficient of the Byerlee law (Byerlee, 1978).
Z. Lei and J.H. Davies where α is the coefficient of thermal expansion, ρ s is the reference density at the surface temperature T s (Table 1).
The key rheology difference of the model setup with previous dual inward dipping subduction models (Dasgupta and Mandal, 2018;Holt et al., 2017;Lyu et al., 2019) is that the magnitude of viscosity here considers temperature dependent creep deformation mechanisms (Fig. S1).The governing rheological laws are identical throughout the model domain, though the rheology parameters we use differ to match different deformation mechanisms potentially dominating at different depths in the Earth.In detail, a uniform composite viscosity is used to take account of four deformation mechanisms under different temperature-pressure conditions: diffusion creep, dislocation creep, Peierls mechanism, and yielding (Garel et al., 2014).The effective composite viscosity in the computational domain is given by μ = ( 1 where μ diff , μ disl , μ y define the creep viscosity following in which A is a prefactor, n the stress component, E the activation energy, P the lithostatic pressure, V the activation volume, R the gas constant, T r the temperature obtained by adding an adiabatic gradient of 0.5 K/km in the upper mantle and 0.3 K/km in the lower mantle to the Boussinesq solution (Fowler, 2005), εII the second invariant of the strain rate tensor.Note that in the lower mantle only diffusion creep applies and the lower mantle is 30 times more viscous than the upper mantle.
While the fourth deformation mechanism, yielding, is defined by a brittle-failure type yield-stress law as with μ y the yielding viscosity and τ y the yield strength.τ y is determined by with τ 0 the surface yield strength, f c the friction coefficient, P the lithostatic pressure, and τ y,max the maximum yield strength (Table 1).

Model setup
The computational domain is 10,000 km by 2900 km, with x (width) coordinates and z (depth) coordinates extending from the surface to the bottom of the lower mantle (Fig. 1).Such a wide domain reduces the influence of side and bottom boundary conditions (Chertova et al., 2012).The thermal boundary conditions at the surface and bottom are defined by two isothermal values: T = T s and T = T m for surface and base of lower mantle respectively, while the sidewalls are insulating.As for mechanical boundary conditions, a free-surface is applied at the top boundary to facilitate trench mobility, and create more accurate topography and lithosphere stress state, while the other boundaries are free-slip.
To simplify the model, a laterally symmetric dual inward dipping subduction is applied.The model is strictly symmetric along the vertical middle line of the domain (5000 km away from the side boundaries) in all aspects, e.g., the geometry and rheology properties.Age 0 SP and Age 0 OP represent the initial ages of subducting plate and overriding plate at the trench, where the two plates meet at the surface.Laterally on the surface, the age of the subducting plates increases linearly with their distance away from the mid-ocean ridge on either side.While vertically, the age of the plate at surface defines the initial thermal structure through a half-space cooling model (Turcotte and Schubert, 2014), with x the distance away from the mid-ocean ridge, erf the error function, z the depth, κ the thermal diffusivity.All parameters are listed in Table 1.The whole overriding plate is set up with a constant age.Thus, the thermal structure within the overriding plate is laterally homogeneous.The bottom of the thermal lithosphere is defined as the isotherm of 1300 K, where the temperature gradient starts to drop quickly (Garel and Thoraval, 2021).The initial thickness of the subducting plate (H 0 SP ) and overriding plate (H 0 OP ) can be calculated using Z. Lei and J.H. Davies where H 0 Plate is the initial thickness of the plate thermal lithosphere and erfinv is the inverse error function.
The free surface boundary condition together with the mid-ocean ridge setup allows the subducting slabs, the overriding plate and therefore the trench to move freely as subduction evolves.To initiate self-driven subduction without implementing external forces, the subducting plate is set up with a bend into the mantle and an 8 km thick low-viscosity decoupling layer on the top.This weak layer has the same rheology laws as the rest of the domain, other than its maximum viscosity is 10 20 Pa s, and its friction coefficient is 0.02 (i.e., an order of magnitude lower).Such unique rheological properties of the weak layer are graded out below a depth of 200 km to boost simulation efficiency, since the main objective of the weak layer is to decouple the subducting plate from the overriding plate at their interface.The initial slab bending radius is 250 km and initially the slab bends over 77 • from the trench (Fig. 1).
The resolution of the adaptive mesh ranges from 0.4 km to 200 km.The initial maximum resolution is in the weak zone (Fig. 1, b), and the minimum resolution of 200 km is in the lower mantle.The meshes are refined by the spatial gradients of the velocity, temperature, viscosity, and material types (weak layer versus normal material).

Model variables
Three variables are investigated here: the initial length of the overriding plate (L 0 OP ), the initial thickness of the subducting plate (H 0 SP ) and overriding plate (H 0 OP ) (Table 2).These are parameters also varied in previous research and therefore will allow easier comparison.H 0 SP and H 0 OP are dependent on plate age and calculated using Eq. ( 12).The magnitude of L 0 OP that has been tested in previous models ranges from 500 km to 4000 km (Dasgupta and Mandal, 2018;Holt et al., 2017;Lyu et al., 2019), and the result shows that L 0 OP >2500 km has little impact on the result (Lyu et al., 2019).Here L 0 OP is tested in the range from 1000 km to 2100 km.The values of H 0 SP and H 0 OP that has been tested before ranges from 75 to 125 km and 75-150 km separately and those models suggest that H 0 SP is more important in deciding the magnitude of upwelling mantle flow than H 0 OP (Lyu et al., 2019).So the range of H 0 SP is extended to 94-141 km (90-200 Ma, Table 2) while the range of H 0 OP is narrowed down to 67-100 km (45-100 Ma).

Varying viscosity in an evolving model: An example
The thermal-mechanical model setup of this research enables selfconsistent subduction, subduction that is driven just by the model's internal buoyancy and not by velocity boundary conditions.Similar to the self-consistent single subduction numerical and analogue models, subduction initiates as negative buoyancy pulls the slab to sink into the deeper mantle, followed by a second stage when slab starts to interact with the lower mantle (e.g., Capitanio et al., 2010;Gerya et al., 2008;Schellart and Moresi, 2013).We next describe in detail the dynamic evolution of dual inward dipping subduction for the model 'L 0 OP = 1700 km'.

Subduction through the upper mantle
When slabs subduct through the upper mantle, symmetric subduction develops about the midline of the overriding plate (~5000 km away from the side boundaries).As more slab is pulled into the mantle, the negative buoyancy grows gradually.It takes ~5.8 Myr before the slab starts to interact with the lower mantle (Fig. 2).
A pair of convective mantle flows with opposite sense of rotation are generated as the subducting slabs bend and sink in the upper mantle.The size of the convective cell grows with time and forms a crescent shape as wide as ~500 km before the slab reaches the depth of lower mantle.The convective cell is composed of a narrow downwelling flow coupling close to the sinking slab and a wide upwelling flow further away.The upwelling flow fades gradually as its distance away from the subducting slab increases.In the model 'L 0 OP = 1700 km', the two sets of wedge flow have little interaction and can be considered as two separate units.This is because the length of the overriding plate is 1700 km, which is greater than two times the width (~500 km) of a convection cell.
The overriding plate exhibits a widespread tensile deviatoric stress field as a result of continuous subduction and the induced convective mantle wedge flows.Only a limited area close to the interface with the bending slabs develops compression (Fig. 2, a).The widespread positive normal deviatoric stress field (τ xx ), up to ~100 MPa, implies that the overriding plate holds an overall stretching tendency.Within the overriding plate, the governing deformation mechanism is spatially layered (Fig. 2, b).At depths shallower than 30 km within the overriding plate, yielding (pseudo-plastic) deformation dominates.Underlying the yielding layer lies ~10 km thick Peierls creep layer.While for depths from ~40 km to the bottom of the thermal lithosphere deformation is dominated by dislocation creep.High strain rate areas are observed within the overriding plate where isoviscos contour necks (Fig. 2, c).The thermal thickness of the overriding plate, defined by the 1300 K isotherm contour, remains nearly constant throughout the simulation.
As subduction initiates, it creates rheology heterogeneities within what initially was a laterally homogeneous overriding plate, allowing part of it to become weaker than the other parts.To visualize the variation in lithosphere viscosity, several levels of isoviscous contours are plotted, e.g., 10 24 , 10 23 , 10 22 , 10 21 , 10 20 Pa⋅s (Fig. 2, light yellow contours).As viscosity is a direct indicator of lithosphere's resistance to deformation at a given rate, here, we define the progressive weakening in the overriding plate as continuous viscosity reduction, i.e., necking of isoviscous contours.The weakening level is defined as the maximum order of viscosity magnitude drop.That is, weakening level 'I', 'II', 'III', 'IV' represents that the isoviscous contour 10 24 , 10 23 , 10 22 , 10 21 Pa⋅s is necked within the overriding plate respectively.The minimum viscosity achieved in the overriding plate for model 'L 0 OP = 1700 km' throughout Models are named with the variable tested, e.g., H 0 SP = 94 km and H 0 SP = 122 km corresponds to the initial subducting plate thickness of 94 km and 122 km separately, while the initial overriding plate length and thickness in both models remain the same as 500 km and 67 km.
the 10 Myr simulation is 10 21 -10 22 Pa⋅s, i.e., weakening level 'III'.The screenshots show that the homogeneous overriding plate is gradually segmented into three strong cores connected with two low viscosity necking centres (Fig. 2, a-b), where minimum viscosity develops and high strain rate is likely to localise (Fig. 2, c).It is noted that the two necking centres and high strain rate regions are spatially coupled with the underneath convective mantle wedge flow induced by the two slabs (Fig. 2, c), suggesting a possible causal connection between plate weakening and the slab induced mantle convection.

Subduction into the lower mantle
Due to the viscosity jump at the transition zone, subduction Z. Lei and J.H. Davies approaching the lower mantle would experience a short period of deceleration, where slab sinking rate and mantle wedge flow rate both reduce.Meanwhile, the necking of the isoviscous contours is reversed by a cooling (plate thickening) and strengthening (ε II reduction and viscosity increase) process within the overriding plate (Figs. 2,.Besides, the high tensile τ xx reduces and it becomes more heterogenous within the overriding plate.At the end of the 10 Myr simulation, the dip between the top of bending slab and the transition zone is ~45 • .

Length of the overriding plate
The first series of models investigate shortening the initial length of the overriding plate (L 0 OP ) or the initial distance between the two trenches at the surface from 2100 km to 1000 km, while keeping the initial thickness of the subducting and overriding plate as constants.As L 0 OP shortens, the two symmetric subducting slabs get closer, and the slab induced two separate mantle wedge flow cells start to combine with each other gradually, forming a joint and stronger diverging mantle flow underneath the overriding plate (Fig. 3, a-c, 5.6 Myr).Above the convective upper mantle, the two separate necking centres (Fig. 3, black squares) in the overriding plate get closer and merge into a single one Necking centres with minimum viscosity in the overriding plate are marked as black squares, which shows that the two separate necking regions tend to get closer and merge into a single one as the length of the overriding plate shortens from 1700 km to 1100 km.A detailed explanation of the contours and other symbols can be found in the caption of Fig. 2.
Z. Lei and J.H. Davies when L 0 OP is less than ~1200 km.As L 0 OP is reduced, it also takes less time to lower each level of viscosity within the overriding plate (Fig. 3, a-c).More importantly, the progressive weakening process can go further and neck the 10 21 Pa s isoviscous contour (weakening level 'IV') when L 0 OP is less than ~1300 km, initiating significant lithosphere thinning and even rifting and spreading extension within the overriding plate (Fig. 3, b-c).In this paper, we define rifting as a process where the plate's thermal thickness reduces to ~0 km, forming a mid-ocean ridge like structure and dividing the overriding plate into two separable plates.Spreading extension corresponds to complete lithosphere separation after rifting and generation of a new oceanic floor.It is noted that the rifting and spreading are characterized as thermal structures in this research, ant it does not include melting behaviour.
The extension deformation is favoured by the overall tensile deviatoric stress in the overriding plate before the slab reaches the lower mantle (Fig. 3).While the extension usually only lasts <1 Myr, it induces substantial modification to the dual inward dipping subduction system.For example, significant slab rollback starts to develop, creating a flattening slab geometry in the upper mantle and steepening dip angle (45 • to 75 • , Fig. 3) at the transition zone depth by the end of the 10 Myr simulation.
To take a closer look at the extension behaviour and its connection with viscosity reduction within the overriding plate, we plot the evolving magnitude of viscosity and second invariant of strain rate at 5 km depth (Fig. 4).The filled region in the figure represents the overriding plate, therefore its widening represents extension within it.The stretching factor β (McKenzie, 1978), defined as the final length of the overriding plate divided by its initial length, is used to quantify the bulk extension that develops within the overriding plate.As L 0 OP shortens, the stretching factor (β) of the overriding plate increases from 1.07 (Fig. 4, a, c) to 1.65 (Fig. 4, b, d) over the 10 Myr simulations, and the corresponding total extension increases from ~100 km to ~600 km.Meanwhile, the highest weakening level achieved within the overriding plate increases from level 'III' (L 0 OP = 1700 km) to level 'IV' (L 0 OP ≤ 1300 km).The contour maps show that notable extension during 5-6 Myr ties temporally and spatially with the necking of isoviscous contour ≤ 10 22 Pa⋅s (Fig. 4, a, b).During the same time, the everincreasing strain rate tends to narrow in width and localise around the middle of the overriding plate where minimum viscosity is located (Fig. 4, d), indicating the presence of strain localisation in the OP.The strong spatial correlation between areas of high strain rate and low viscosity is not surprising according to how viscosity is defined with Eq. (8,9), while we recognise the continuous necking of viscosity magnitude below ∼ 10 22 Pa⋅s acts as a good indicator for identifying the occurrence of strain localisation.
The edge of the filled contour in the lateral direction represents the interface between the overriding plate and subducting plate, i.e., the trench.It is noted that stagnant trenches are observed throughout the simulation, except for a short period (~1 Myr) of trench retreat coupling with OP extension when plate weakening reaches level 'III' (viscosity Z. Lei and J.H. Davies reduced to less than 10 22 Pa⋅s) in the overriding plate (Fig. 4, a-b).The stagnant trenches indicate that dual inward dipping subduction selfconsistently forms a "fixed" side boundary condition for the overriding plate.

Thickness of the overriding plate
The second series of models increase the initial thermal thickness (defined by the 1300 K contour) of the overriding plate (H 0 OP ) from 67 km to 100 km, while keeping the initial subducting plate thickness and the initial length of the overriding plate as constants.As H 0 OP increases, the time it takes for the slab to sink to the depth of 660 km gradually increases from 6.0 Myr to 8.8 Myr (Fig. 5).Meanwhile, the maximum weakening level developed within the overriding plate drops from 'IV' (H 0 OP = 67 km) to 'III' (H 0 OP = 74 km) and less than 'I' (H 0 OP = 100 km).The time it takes to lower each order of viscosity magnitude also increases, indicating that thicker overriding plate is more resistant to deformation during subduction.Besides, the overall tensile deviatoric stress state in the overriding plate (Fig. 5, a-b) is replaced by a more heterogeneous stress state when slabs sink in the upper mantle (Fig. 5,  c).
As H 0 OP thickens, the stretching factor (β) of the overriding plate decreases from 1.65 (Fig. 6, a) to 1.01 (Fig. 6, b) over the 10 Myr simulation.The corresponding total extension decreases from ~600 km to ~10 km.It is noted that model 'H 0 OP = 100 km' holds stagnant trenches throughout the 10 Myr simulation, indicating the presence of a seemingly immovable boundary condition ("fixed") imposed on the overriding plate during dual inward dipping subduction.Z. Lei and J.H. Davies Both the section view and contour map show that the highest viscosity reduction occurs along the midline of the overriding plate in models with varying H 0 OP .To investigate the progressive weakening in the necking centre, diagnostics are tracked along the vertical slice in the middle of the overriding plate (5000 km away from both side boundaries).Three diagnostics are integrated along the vertical slice and then divided by the thickness of the plate (Eq.( 13)), in which D represent the diagnostic.The averaged results include magnitude of viscosity (μ), second invariant of strain rate (ε II ), deviatoric normal stress component (τ xx ).The real-time lithosphere thickness (H OP ) in the middle of the overriding plate and the magnitude of subduction rate at trench (v SP ) are also recorded for further analysis (Fig. 7).Broadly, the evolution of mean viscosity (μ) in the overriding plate correlates well with the variation of the other diagnostics.Take the model 'H 0 OP = 67 km' for example (blue line or points in Fig. 7).There is gradual decrease of μ, and gradual increase of τ xx , εII and v SP during the simulation between 1 Myr to 2.5 Myr, while H OP remains nearly constant.During 2.5 Myr to 5 Myr, μ as a function of time remains a consistent negative slope as before, while τ xx starts to decrease gently after peaking at ~4 Myr.Meanwhile, H OP as a function of time starts to decrease with a gradually steepening negative slope.During the rifting and spreading extension between 5 Myr to 6 Myr, all diagnostics are varying more rapidly, with μ, H OP , τ xx dropping sharply and εII climbing steeply (log scale).Then the slab reaches the depth of 660 km (Fig. 7, marked by downward triangles), and starts to interact with the lower mantle.The weakening process is replaced by a gradual cooling and strengthening process in the overriding plate, where μ and H OP both increase while τ xx and εII decrease in magnitude.
As the thickness of the overriding plate increases from 67 km to 100 km, the minimum μ observed in the necking centre increases from ∼ 2 × 10 18 Pa⋅s to ∼ 2 × 10 23 Pa⋅s, which suggests that a thicker overriding plate is more resistant to deformation induced by the dual inward dipping subduction.In detail, Fig. 7-a (grey dashed line) shows that if μ is above ∼ 2 × 10 22 Pa⋅s, there is little lithospheric thinning (< 5 km) in the necking centre (Fig. 7, a, c).In the timesteps when μ is in the range of 10 21 − 2 × 10 22 Pa⋅s, thinning (< 25 km) starts to build up, but it is not weak enough to develop rifting extension (Fig. 7, c, green line).Only when the μ drops below 10 21 Pa⋅s does spreading extension (Fig. 7, c, blue and orange lines) develop within the overriding plate.It is noted in ).These results confirm that the evolution of viscosity magnitude is a good indicator for the various progressive weakening developed in the overriding plate during subduction.In Fig. 7-d, each scatter point represents 0.2 Myr simulation.The temporal path shows that as the initial thickness of the overriding plate reduces, εII in the overriding plate becomes increasingly sensitive to increasing subduction rate before the slab sinks to the depth of 660 km.Then the temporal path inflects upwards to higher εII for models that induce spreading extension (model 'H 0 SP = 67 km' and 'H 0 SP = 70 km'), and downwards for models that fail to spread.The upward inflection suggests that once extension initiates, it does not take much subduction rate to maintain the high strain rate developed in the overriding plate.As slabs slowly sink to the lower mantle, the magnitude of subduction rate and εII both gradually decrease, and the temporal path returns to the starting point for all models.
In the models with varying H 0 OP , the maximum subduction rate ranges from 20 to 35 cm/yr, while the corresponding maximum εII in the overriding plate ranges from 10 − 15 to 10 − 12 s − 1 .The general positive correlation between these two variables suggests that maximum subduction rate may play an important role in weakening the overriding plate.Over the 10 Myr simulation, only ~15% scatter points fall in the range where subduction velocity is over 10 cm/yr, implying that high strain rate deformation correlated with high subduction rate only lasts for a short period.

Thickness of the subducting plate
The third series of models investigate increasing the initial thermal thickness (defined by the 1300 K contour) of the subducting plate (H 0 SP ) from 94 km to 141 km, while keeping the initial overriding plate thickness and the initial length of the overriding plate as constants.As H 0 SP increases, the time it takes for the slab to sink to the depth of 660 km gradually reduces from 6.6 Myr to 6.0 Myr.Meanwhile, the maximum weakening level developed within the overriding plate increases from 'II' (H 0 SP = 100 km) to 'III' (H 0 SP = 122 km) and 'IV' (H 0 SP = 141 km).Besides, the time it takes to lower each order of viscosity magnitude reduces, indicating a faster progressive weakening as H 0 SP thickens (Fig. 8, b-c).In contrast to models with varying H 0 OP , varying H 0 SP does not affect the overall tensile deviatoric stress state developed in the overriding plate before slabs reach the lower mantle.
As H 0 SP increases, the stretching factor (β) of the overriding plate increases from 0.98 (Fig. 9, a) to 1.65 (Fig. 9, b) over the 10 Myr simulations, and the corresponding total extension increases from − 15 km (H 0 SP ≤ 100 km) to 600 km (H 0 SP = 141 km).In line with models testing different L 0 OP and H 0 OP , stagnant trenches and "fixed" boundary effect for the overriding plate are also observed here.The only difference is that cases with thin H 0 SP ( ≤ 100 km) may develop a small amount of trench advance, i.e., shortening instead of stretching the overriding plate.

Regime of stretching state
A variety of deformation patterns and stretching state within the overriding plate have been observed when varying L 0 OP , H 0 OP and H 0 SP .Several diagnostics are reported together to quantify the deformation developed within the overriding plate during the 10 Myr simulation (Table 3).The detail of each diagnostic is described as follows.
As introduced in Section 3.1.1,weakening levels 'I', 'II', 'III', 'IV' are determined by the minimum viscosity contour which is necked in the overriding plate during subduction.The higher the weakening level, the stronger the localised rheology modification observed within the overriding plate.All three groups of dual inward dipping subduction models manage to yield a variety of weakening levels in the overriding plate (Fig. 10, a-c).
t rift indicates the timestep when the overriding plate develops rifting extension (weakening level 'IV'), and a void value means that model fails to generate rifting extension within the overriding plate.t rift increases with thicker or longer overriding plate, and it decreases with thicker subducting plate.
t 660 equals how much time the subducting plate (defined by its 1300 K isotherm) takes to sink to the depth of 660 km.It is most sensitive to the variation of H 0 OP , while varying L 0 OP and H 0 SP generates less than ~1 Myr difference of t 660 compared with a ~ 3 Myr difference when modifying H 0 OP .v sink is the average sinking speed before slab reaches the lower mantle, and it equals 460 km (the vertical distance from the initial slab tip depth to the depth of 660 km) divided by t 660 .v sink may range from 5.2 to 8.5 cm/yr, and it correlates positively with the weakening level in models that vary H 0 SP and H 0 OP .In contrast, increasing L 0 OP can generate higher v sink , which, however, fails to induce greater weakening level in the overriding plate.
The stretching factor β (McKenzie, 1978), defined as the final length of the overriding plate divided by its initial length, is used to quantify the overall extension that develops within the overriding plate in 10 Myr simulation.β ranges from 0.95 to 1.04 for models that only develop weakening level 'I' or 'II'.In models that develop weakening level 'III', β ranges from 1 to 1.19.β can go up from 1.25 to 1.65 in models that develop weakening level 'IV'.
While the stretching factor (β) quantifies the bulk extension in the OP, it fails to reflect the fact that strain is likely to concentrate around the middle of the overriding plate with high strain rate instead of evenly distributed across the OP.To quantify the deformation accumulated along the midline of the overriding plate, we integrate the average second invariant of strain rate (ε II based on Eq. ( 13)) with time throughout the 10 Myr simulation.All three groups of models generate a wide range of accumulation of strain at the end of the simulation (Fig. 10, d-f).For all models that develop rifting or spreading extension, i.e., weakening level 'IV', the accumulation of strain is >100%.Accumulation of strain falls in the range of 5% to 100% for models that only reach weakening level 'III'.Models achieving weakening level 'II' and 'I', yield <5% strain accumulation, where the deformation is hardly observable in the overriding plate.The magnitude of accumulation of strain is significantly higher than the strain (β-1) derived from the streching factor β in models that reach weakening level 'III' and level 'IV', confirming that strain localisation characterizes the OP extension during dual inward dipping subduction.
All the diagnostics tie closely with the weakening level developed in the OP.By combining the qualitative and quantitative diagnostics presented in the results, we classify three stretching states in the overriding plate during the dual inward dipping subduction: 1) little or no lithosphere thinning and extension, discriminated by low weakening level up to level 'II', β up to 1.04, accumulation of strain up to 5% and almost no thermal lithosphere thinning in the necking centre; 2) limited lithosphere thinning and extension, identified by medium weakening level up to level 'III', β up to 1.19, accumulation of strain up to 30%, and limited thermal lithosphere thinning, e.g., ~15 km thinning for model 'H 0 OP = 74 km' (Fig. 7, c); 3) rifting and spreading extension, characterized by high weakening level up to level 'IV', β >1.25, accumulation of strain over 100%, and total thinning of the thermal lithosphere during rifting extension.

Discussion
The results show that dual inward dipping subduction can induce progressive weakening within a homogeneous overriding plate.With appropriate conditions, tested in this research, e.g., thick enough H 0 SP , thin enough H 0 OP , short enough L 0 OP , different levels of stretching state, ranging from no thinning nor extension to rifting and spreading extension, can develop within the overriding plate.The role that dual inward dipping subduction plays during the progressive weakening and the rheological origin of the weakening process are worth discussion.

Fixed overriding plate boundary condition
The stagnant tendency of the trenches in the viscosity contour maps (Fig. 4, Fig. 6, Fig. 9) show that dual inward dipping subduction can selfconsistently form a fixed lateral boundary condition for the overriding plate.This is due to the symmetric model setup, where subducting plates on both sides are prone to (at least initially) advance or retreat simultaneously.It creates roughly equal, symmetric and competing force from both ends of the overriding plate during subduction.As a result, the mobility of the overriding plate is inhibited, as indicated by the low velocity (<1 cm/yr, white contour) region within the overriding plate before extension develops (Fig. 11, a-b).It would be as if the mechanical boundary condition on the overriding plate was fixed.As the overriding plate keeps weakening during dual inward dipping subduction, divergent velocity can build up within the overriding plate, indicating initiation of extension (Fig. 11, a).Then, a short period (~1 Myr) of fast trench retreat is accommodated by ~600 km extension within the overriding plate.
Previous studies on single-sided subduction cases imply that the mobility of the overriding plate plays an important role in producing extension, especially in the back-arc region of the overriding plate.A mobile overriding plate can move as a whole to inhibit the build-up of deviatoric stress within the plate (Capitanio et al., 2010;Chen et al., 2016;Garel et al., 2014;Holt et al., 2015;Nakakuki and Mura, 2013), while the immobile overriding plate can facilitate strain localisation which accounts for the increased degree of deformation in the overriding plate compared with mobile plates (Arcay et al., 2008;Capitanio et al., 2010;Chen et al., 2016;Erdős et al., 2021;Nakakuki and Mura, 2013;Sternai et al., 2014;Yang et al., 2019).
To investigate the role of fixed trailing boundary condition in promoting extension, we consider a single sided subduction (SSS) model  a The accumulation of strain is calculated along the middle vertical slice (5000 km away from side boundaries) within the overriding plate.For models L 0 OP ≥ 1300 km, the necking centres are away from this middle vertical slice (Fig. 3, b).So, the accumulation of strain could be underestimated for these models.Considering that only model 'L 0 OP = 1300 km' achieved weakening level 'IV', the corrected accumulation of strain along its necking centre is around 600%.The underestimation for other models is moderate and will not change the conclusion of this research.with a free and mobile overriding plate, based on previous SSS research (Garel et al., 2014).The SSS model has the same parameters as the dual inward dipping subduction model 'H 0 OP = 67 km' in every aspect, e.g., rheology, initial subduction plate thickness (141 km) and initial overriding plate thickness (67 km) at the trench, except that there is only one subducting plate and the overriding plate holds a mobile side boundary.During the 10 Myr simulation, the SSS model only reached weakening level 'I', much lower than weakening level 'IV' in its reference model 'H 0 OP = 67 km'.The results also show that much less extension, evidenced by the absence of divergent velocity, is observed in the overriding plate of SSS model (Fig. 11, c) relative to that in the model 'H 0 OP = 67 km' (Fig. 11, a).Thus, the lack of mobility of the overriding plate plays a key role in promoting the degree of weakening during dual inward dipping subduction (DIDS).
Previous research also find that symmetric DIDS configuration can limit the trench mobility (Dasgupta and Mandal, 2018;Holt et al., 2017;Lyu et al., 2019), yet its role in promoting extension in the overriding plate has not been addressed.As mentioned before, previous DIDS models all apply a rheology law that excludes thermal effects for the plates and mantle during simulation.This is likely to create a strong overriding plate where strain or extension can hardly develop in hotter regions during subduction (Fig. 11, b).As it takes extension to accommodate the space trench retreat creates within a fixed OP, trenches tend to remain stagnant across all previous DIDS models throughout the simulation.
To validate this hypothesis, we run an additional model that reduces the thermal dependency for the rheology law applied in the overriding plate.In model 'μ OP disl ' , we removed dislocation creep, a thermal dependent deformation mechanism (Eq.( 8)), from the composite rheology for the overriding plate while keeping every other aspect the same as model 'H 0 OP = 67 km'.The result shows that ~60 km of total extension with maximum weakening level 'III' is observed in the overriding plate (Fig. 11, d), only 10% of the total extension (trench retreat) relative to the reference model 'H 0 OP = 67 km'.According to how bulk viscosity is calculated with Eq. ( 7), it will become even more difficult to deform the overriding plate, if other thermal dependent deformation mechanisms, i. e., Peierls creep and diffusion creep, are removed from the composite rheology for the overriding plate.Bearing this in mind, the result of additional model 'μ OP disl ' reveals why previous DIDS models (Dasgupta and Mandal, 2018;Holt et al., 2017;Lyu et al., 2019), only have nearly stagnant trench and fail to induce spreading extension in the overriding plate, even though the fixed overriding plate boundary effect is also observed.The result also highlights that incorporating thermal dependent rheology can significantly reduce the strength of the overriding plate during dual inward dipping subduction.A detailed evaluation on each deformation mechanism's role in weakening the overriding plate will be discussed in Section 4.2.
To summarise, dual inward dipping subduction can self-consistently form a fixed boundary condition for the overriding plate.This creates an environment promoting the development of strain localisation in the stagnant overriding plate relative to single sided subduction with a free mobile boundary condition.In addition, the incorporation of thermal dependent creep rheology allows higher weaken level to develop in the OP and it facilitates the formation of spreading extension within the OP during dual inward dipping subduction.

Stronger mantle wedge flow
Our results show that varying the size of the overriding plate can affect the degree of extension developed within the overriding plate.Previous research on single sided subduction shows that sinking slab can induce vigorous mantle return flow, which has been suggested to account for extensional deformation within the overriding plate, e.g., back-arc extension, supercontinent breakup (Chen et al., 2016;Dal Zilio et al., 2018;Erdős et al., 2021;Gerardi and Ribe, 2018;Husson, 2012;Sleep and Toksöz, 1971;Sternai et al., 2014).Multiple parameters can influence the magnitude of slab induced mantle wedge flow: a) agedependent thickness of the subducting plate (H 0 SP ) is associated with the magnitude of slab net negative buoyancy (e.g., Garel et al., 2014)  that correlates positively with the intensity of mantle wedge flow (Capitanio and Faccenda, 2012); b) lowering the thickness of the overriding plate (H 0 OP ) can not only increase the net negative buoyancy by increasing the hanging slab area in the upper mantle, but also reduce the shear force along the interface between two coupled plates (Hertgen et al., 2020), thus generating faster subduction rate and inducing stronger mantle wedge flow.These two parameters are tested in this research with a dual inward dipping subduction (DIDS) setup, and they also show effective control on the subduction rate (t 660 and v sink in Table 3) which correlates positively with the maximum weakening level (Table 3) and the second invariant of strain rate in the overriding plate (Fig. 7, d).Thickening H 0 SP or thinning H 0 OP can increase the velocity of descending slab is also reported in previous DIDS models (Lyu et al., 2019).It is also noted that higher slab sinking rate can induce stronger mantle wedge flow, which accounts for growing tensile deviatoric stress in the overriding plate.DIDS models in (Dasgupta and Mandal, 2018) use a velocity boundary condition for the subducting plate, so subduction rate is prescribed and constant ranging from 1 to 5 cm/yr for five models.The results also show that higher subduction rate can induce stronger mantle wedge flow which creates increasing tensile stress field in the central part of the OP.Nevertheless, the potential of slab induced mantle wedge flow in weakening (e.g., thinning and stretching) the overriding plate, is not addressed in previous DIDS models, as no significant deformation that relates to spreading extension is observed in the OP throughout the simulation.This could be ascribed to the lack of thermal dependency for the rheology law applied, as indicated in Fig. 11  (a, d).The topic will be further discussed in Section 4.2.
Apart from varying H 0 SP and H 0 OP , we find that increasing the initial length of the overriding plate (L 0 OP ) can also increase the subduction rate (Table 3).However, it fails to induce higher weakening level within the  13)) and the velocity gradient of the divergent mantle flow underlying the overriding plate midline (∇v x ).The maximum velocity gradient of divergent mantle flow in the range of t 660 for each model is marked by a vertical dashed line.
Z. Lei and J.H. Davies overriding plate.Actually, a negative correlation between subduction rate (v sink ) and the OP weakening level is observed instead.As H 0 OP remains constant in models with varying L 0 OP , the OP's strength or its resistance to deformation induced by the underlying mantle flow remains unchanged.In this case, it suggests that L 0 OP may dominate over subduction rate in affecting the intensity of slab induced mantle wedge flow or its ability to weaken the OP.To evaluate how L 0 OP modifies the intensity of slab induced mantle wedge flow, we quantify the vertical (v y ) and horizontal (v x ) component of velocity field in the mantle along a lateral slice at the depth of 75 km (Fig. 12), which is ~8 km below the overriding plate (bottom defined by 1300 K isotherm).It shows that the gradient of diverging mantle flow (slope of v x , ∇v x ) increases gradually underlying the OP midline as L 0 OP shortens (Fig. 12, c-d), suggesting a growing shear stress applied at the bottom of the stagnant OP.It is also noted that models with L 0 OP ≥ 1500km tend to form two separate mantle return flow cells, above which two individual necking centres develop symmetrically along the midline of the overriding plate (Fig. 12, a).In contrast, slab induced convective mantle flows tend to unite with each other and act upon a single necking centre in models with L 0 OP ≤ 1200 km, generating higher magnitude of upwelling component of mantle flow (Fig. 12, a).Resulting from the stronger diverging mantle flow underlying the OP, higher magnitude of strain rate and weakening level develop in models with shorter L 0 OP (Fig. 12, d).Considering the slowing down of the slab sinking rate as L 0 OP shortens, previous research suggests that the stronger mantle wedge flow can generate greater velocity gradients in the direction perpendicular to the slab surface, creating a strong shear field that can flatten the slab dip (Dasgupta and Mandal, 2018;Enns et al., 2005;Holt et al., 2015).The flattened slab dip is also observed in Fig. 3 over the 10 Myr simulation.Another possible explanation is that the sinking rolling back subducting slabs require flow from the deeper mantle to enter the mantle wedge for conservation of mass.This becomes harder as the gap between the slabs becomes narrower.Therefore, this reduces how quickly the slabs with shorter L 0 OP can sink.Here, we propose that the distance between slabs can play a dominant role over the slab sinking rate in regulating the mantle wedge flow intensity, which affects the mantle flow's potential to weaken the overriding plate and modulate the slab motion.
Our result is consistent with previous DIDS models, which also find that reducing the L 0 OP can slow down the sinking rate.This effect is supported by evidence such as increased t 660 (Lyu et al., 2019) and gentler slab dip in the upper mantle (Dasgupta and Mandal, 2018).Besides, the effect of creating a stronger mantle flow as L 0 OP reduces is addressed to account for growing tensile deviatoric stress in the OP (Dasgupta and Mandal, 2018;Lyu et al., 2019).Yet, no significant extension in the OP correlating with the strengthened mantle flow was observed.This may be because previous DIDS models do not consider thermal dependency for the rheology, making the OP too strong to be stretched as suggested in Fig. 11 (a, d).A detailed evaluation on the role of temperature dependent rheology in assisting plate weakening in the overriding plate will be discussed in Section 4.2.
A similar effect is also reported in other multiple slab subduction models, where mantle flow induced by two neighbouring 3D slabs with opposite dips tends to merge with each other and forms a stronger one as the lateral distance between the two slabs reduces (Király et al., 2016).It is also noted that the integrated mantle flow can regulate the slab motion and trench retreat velocity (Király et al., 2016).
In brief, previous research on subduction shows that the slab induced return flow can apply basal traction, which can propagate upwards and create high enough tensile deviatoric stress that overcome the strength of the overriding plate (e.g.Capitanio et al., 2010;Dal Zilio et al., 2018;Holt et al., 2015;Husson, 2012;Sternai et al., 2014;Yang et al., 2019).A united, thus stronger return flow is likely to strengthen the mantle flow's ability to weaken the overriding plate during dual inward dipping subduction (Fig. 12, d).

Overriding plate weakening mechanism
As introduced in the methods, we applied composite rheology which incorporates four deformation mechanisms everywhere in the simulation domain.Here, the dominant deformation mechanism (DDM) is defined as the rheology law that yields the minimum magnitude of viscosity at a certain point, similar with deformation mechanism partitioning in previous research (Bessat et al., 2020;Garel et al., 2020;Ruh et al., 2022).We try to understand the temporal and spatial evolution of the DDM within the overriding plate, especially around the midline of the OP where strain localisation tends to develop.Then we evaluate the contribution of each deformation mechanism in promoting viscosity reduction within the overriding plate, especially the temperature dependent creep rheologies.

Dominant deformation mechanism analysis
Model 'L 0 OP = 1700 km' with limited extension has shown that the DDM is stratified with yielding, Peierls creep and dislocation creep as temperature increases with depth in the overriding plate (Fig. 2, b).Here, we further investigate how the DDM evolves in models that develop rifting and spreading extension within the overriding plate, e.g., model 'H 0 OP = 67 km'.Therein, the temporal phases show that the DDM is also spatially layered (Fig. 13, a-g), with yielding initially dominating from the surface to the depth of ~35 km, underlain by Peierls creep dominating for the next ~10 km and then dislocation creep dominating for ~25 km (Fig. 13, b-d).Among all the DDM at different depths throughout the simulation, the DDM of yielding layer is always the thickest and the DDM of dislocation creep layer comes as the second.To be noted, the DDM of diffusion creep with limited area is observed around the bottom of the overriding plate during the initial plate weakening (Fig. 13, b), and it is completely replaced by dislocation creep after 3 Myr.During the thinning process of the overriding plate, the deformation mechanism of Peierls creep gives way to yielding and dislocation creep as DDM (Fig. 13, d-g).The replacement and interplay among different DDM will be discussed in the next subsection.
While we do not implement a multi-material approach to define the rheology of different layers (e.g., crust) in the lithosphere, the uniform compositional rheology law self-consistently generates the layered structure in Fig. 13.In detail, yielding only dominates over other creep mechanisms in the cold regions, corresponding to the crustal depth range.The temperature-dependent creep mechanisms dominate over yielding in the hot bottom region, equivalent to the depth range of mantle lithosphere.The continuous necking process shows that the viscosity reduction initiates from the surface (yielding) and the bottom of the plate (dislocation or diffusion creep).Then the viscosity contour necks in the middle depth (Peierls creep) of the plate as seen in Fig. 13,  (b-e).This suggests that both thermal dependent rheology and nonthermal dependent rheology (yielding) contribute to the progressive weakening in the overriding plate (OP).To confirm whether rifting and spreading extension still occur in the OP when thermal component of rheology is reduced or removed, we run two additional models.These models are identical with model 'H 0 OP = 67 km', except that one excludes dislocation creep for the overriding plate (model 'μ OP disl ') and the other only keeps the yielding deformation mechanism for the overriding plate (model 'μ OP Y only ').The result shows that the maximum weakening level additional model 'μ OP disl ' achieves is level 'III', i.e., it fails to neck the isoviscous contour of 10 21 Pa⋅s or generate spreading extension before the slab reaches the 660 km at ~7 Myr (Fig. 13, h-l).It indicates that dislocation creep's role in inducing OP spreading extension is irreplaceable.Based on this result and how the bulk viscosity is calculated (Eq.( 7)), removing any further thermal dependent creep deformation mechanisms, i.e., Peierls creep and diffusion creep, will only create an even stronger overriding plate, i.e., no ridge formation or spreading extension will ever occur.Here, we demonstrate this prediction by presenting model 'μ OP Y only ' where the model has a strong lithospheric bottom, and it fails to even neck the isoviscous contour of 10 24 Pa⋅s (weakening level 'I') in the OP after 6 Myr simulation (Fig. 13, m).The result of additional model 'μ OP disl ' and model 'μ OP Y only ' also reveal why previous DIDS models, without considering temperature dependent rheology, fail to induce any significant extension or trench retreat in the overriding plate (Dasgupta and Mandal, 2018;Holt et al., 2017;Lyu et al., 2019).
Another end member of rheology setup would be only incorporating temperature dependent creep rheology and excluding yielding.Such rheology setup is commonly used in large-scale mantle convection Z. Lei and J.H. Davies models, which focus on studying deformation in the upper and lower mantle.Composite rheology incorporating dislocation creep and diffusion creep is usually considered, and it will generate a strong layer (μ ≥ 10 24 Pa⋅s) with extremely low strain rate (≤ 10 − 18 s − 1 ) in cold regions (less than ~600 K) near the surface (Schulz et al., 2019).Similar stronger surface layer are widely observed when only dislocation and diffusion creep are considered for the composite rheology (Asaadi et al., 2011;Becker, 2006;Becker et al., 2008).We do not run models with only temperature dependent creep rheology laws here, but we briefly investigate the case by plotting the viscosity contour map as a function of temperature, second invariant of strain rate, and depth (lithostatic pressure) in the supporting material (Fig. S1, g-i).It also suggests a strong layer in the cold region (<600 K) near the surface, and it takes several orders of magnitude higher strain rate to reduce viscosity to lower than 10 21 Pa⋅s relative to cases that consider both yielding and temperature dependent rheology (Fig. S1, d-f).As model 'μ OP Y only ' has demonstrated that it takes viscosity reduction from both the surface and the bottom to develop continuous necking, a lack of viscosity reduction around the surface will inhibit the development of plate weakening in the OP.Thus, we conclude that it takes both yielding and temperature dependent rheology laws to promote continuous viscosity reduction and induce spreading extension in the overriding plate.

Weakening contribution analysis
The previous section has shown that the DDM may vary at different depth range within the overriding plate.To evaluate the contribution of each DDM to inducing rifting and spreading extension throughout the simulation, we slice the overriding plate vertically through its midline where the most intensive necking and strain localisation take place.Then we group the points along the midline by the type of DDM at each timestep.For each group with the same DDM, we calculate at each timestep the arithmetic average of the second invariant of strain rate and temperature state, which are then plotted on the viscosity contour map computed through Eq. (7-9).As the magnitude of viscosity is not sensitive to lithostatic pressure in the depth range of 70 km (Fig. S1, e-f), we use the depth of 50 km to create the background viscosity contour map (Fig. 14, a).
One way to evaluate the contribution of deformation mechanisms to plate weakening is to quantify how much (order of) viscosity reduction each DDM achieves.For model 'H 0 OP = 70 km', both yielding and dislocation creep reduce the viscosity to magnitude lower than 10 21 Pa⋅s (Fig. 14, a), which is the critical magnitude to initiate rifting and spreading extension (Fig. 7, a).This indicates that yielding and dislocation creep play an active role throughout the plate weakening process, which includes lithospheric thinning, rifting and spreading extension within the OP.Peierls creep can reduce viscosity to the range of 10 21 -10 22 Pa⋅s (Fig. 7, a), suggesting that it also contributes substantially to induce plate thinning.The absence of Peierls creep in the necking centre when rifting and spreading extension develops in the OP implies that it is secondary deformation mechanism relative to yielding and dislocation creep as strain localizes in this stage.Diffusion creep induces the least viscosity reduction to ~5 × 10 22 Pa⋅s, which suggests that it only softens the plate's bottom layer for further deformation.For models that do not develop rifting or spreading extension, the temporal paths of the DDM are similar with model 'H 0 OP = 70 km' except that the minimum viscosity is never less than 10 21 Pa⋅s.
To be noted, we observe an abrupt accelerating viscosity reduction in the range of 10 20 -10 22 Pa⋅s for the DDM of yielding and dislocation creep (Fig. 14, a).Concurrently, the strain rate also jumps to a faster rate, displayed by the less dense scatter of points (each represents a timestep).As demonstrated in Fig. 4 (b, d), the necking of isoviscous contour of 10 22 Pa⋅s also indicates the development of strain localisation with narrowing width of high strain rate region in the OP.The accelerating viscosity reduction implies that the overriding plate along its midline (necking centre) falls into positive feedback between plate weakening and strain localisation.A similar feedback process is widely observed in simulations that involve plate-scale weakening and generation of new plate boundaries (Fuchs andBecker, 2022, 2019;Gueydan and Précigout, 2014;Wenker and Beaumont, 2018).In the case of uniaxial stretching, the plate strength is proportional to μ × H OP (Ribe, 2001), both of which in our models are reducing during the plate thinning process.Considering that the plate strength measures the very resistance to the underlying mantle flow, the reduction of viscosity and plate thickness will incur further plate weakening, which in turn allows higher strain rate to develop.As multiple rheology laws (yielding, Peierls creep and dislocation creep) are strain rate dependent here, the everincreasing strain rate can induce further viscosity reduction, which marks another round of feedback weakening.It has been proposed that the power law exponent over 1 for the creep mechanisms can promote the development of strain localisation (Wenker and Beaumont, 2018), which can abruptly accelerate plate weakening by inducing nonlinear decay of the plate strength (Brune et al., 2016).The interpretation is consistent with our results.Here, we use power law exponent of 3.5 and 20 for dislocation creep and Peierls creep mechanism (Garel et al., 2014;Maunder et al., 2020), which yields a similar variation of effective viscosity with published dislocation and Peierls creep laws (Hirth and Kohlstedt, 2003;Kameyama et al., 1999;Katayama and Karato, 2008).
While the rheology law (Eq.(8-9)) of the four deformation mechanisms shows that the magnitude of viscosity is dependent on evolving temperature, strain rate, and lithostatic pressure, the diagram indicates that the viscosity reduction is most sensitive to the ever-increasing strain rate, which results from both yielding and temperature dependent creep deformation mechanisms (Fig. 14, a).Specifically, the incorporation of temperature dependent rheology enables the thermally activated strain rate weakening, which is key to weaken the lithosphere hotter than ~900 K and induce spreading extension within the OP.This is a major improvement in rheology and its impact on promoting extension within the OP relative to the previous DIDS models that do not consider thermal effects (Dasgupta and Mandal, 2018;Holt et al., 2017;Lyu et al., 2019).The results also demonstrate that the slab induced mantle flow weakens the mantle lithosphere mainly by inducing thermally activated strain rate weakening rather than by heating up the lithosphere via conduction.The dominant role of strain rate induced weakening over heat conduction (a much slower process) is also reported in the interaction between upwelling plumes and overlying lithosphere (Burov and Guillou-Frottier, 2005).That is to say, the rheology and buoyancy parameters will be more important than the heat conduction parameters in producing different levels of rheological weakening within the overlying plate.

Limitations
Relative to the previous dual inward dipping subduction models, a major improvement of this work is incorporating a composite rheology that considers temperature dependent creep deformation mechanisms, which promotes the continuous viscosity reduction and the development of strain localisation in the overriding plate.Though, we recognise that many other processes, which we do not consider, can also affect the rheological evolution of the lithosphere, e.g., migration of hydrous fluids, partial melting, and grain size evolution (Arcay et al., 2008;Braun et al., 1999;England and Katz, 2010;Fuchs and Becker, 2019;Hicks et al., 2023;Montési and Hirth, 2003).For example, grain size reduction is likely to take place when strain builds up and it may make diffusion creep become the dominant deformation mechanism, replacing dislocation creep, in the mantle lithosphere (Gueydan et al., 2014;Ruh et al., 2022).Taking all these parameters into consideration is likely to generate a more realistic simulation but at the cost of making the computation much more expensive for plate scale simulation (e.g., Dannberg et al., 2017).
Subduction can generate convective mantle flow that includes both poloidal and toroidal components.The 2D models tested here neglect the effects of toroidal flow and the third dimension.This could amplify the magnitude of poloidal flow and its weakening effect applied within the overriding plate.Considering that the poloidal component dominates over the toroidal component when slab subducts through the upper mantle (Funiciello et al., 2004) and the magnitude toroidal component correlates positively with trench retreat rate (e.g., Li et al., 2014;Schellart and Moresi, 2013), the lack of toroidal flow would have limited impact on the progressive weakening presented here as trench retreat is insignificant before spreading extension initiates.

Implication for natural DIDS case
Bearing all the limits in mind, we cautiously compare our model predictions with observations in the Caribbean Sea plate.While we do not aim to reproduce the current tectonic framework in the Caribbean Plate, we still find some significant implications in applying our research to understand specific features in this region, e.g., widespread extensional basins, active (back-arc) spreading centres etc.
The Caribbean plate has established dual inward dipping subduction (DIDS) since 90-70 Ma (Boschman et al., 2014;Braszus et al., 2021;Riel et al., 2023), with the Farallon plate (subsequently, Cocos and Nazca plates) subducting at the Central America Trench in the west and Proto-Caribbean plate (followed by Atlantic plate) subducting at the Lesser Antilles Trench in the east.One intriguing observation from plate reconstruction is the apparent increase in the distance between the two trenches from ~1100 km to over 2000 km since the establishment of the dual inward dipping subduction (Barrera-Lopez et al., 2022;Boschman et al., 2014;Braszus et al., 2021;Riel et al., 2023;Romito and Mann, 2021), implying that the Caribbean plate has undergone extension.The extension includes the formation of multiple extensional (back-arc) basins throughout the Caribbean Sea plate, e.g., Tobago Basin, Grenada Basin, Venezuela Basin and Colombia Basin since ~60 Ma (Allen et al., 2019;Braszus et al., 2021;Romito and Mann, 2021;Steel and Davison, 2021).Our generic models demonstrate that DIDS can self-consistently form fixed boundary conditions, which can promote plate weakening and strain localisation within the OP.Furthermore, DIDS can lead to the formation of a united, thus stronger mantle flow that can strengthen its ability to weaken the overriding plate.These two effects have become active since the formation of dual inward dipping subduction framework at 90-70 Ma, and they can be important aspects to consider to understand extensional deformation history in this region.
On the other hand, we acknowledge that the development of the Caribbean Large Igneous Province (CLIP) at ~140-70 Ma (Hoernle et al., 2004) may have also contributed significantly to the widespread extensional units in this region (Pindell et al., 2006).Further discussion on either subduction or CLIP dominates the regional tectonic framework goes beyond the scope of this research, while we recognise that the establishment of DIDS since 90-70 Ma can form an environment that promote the upwelling flow in the upper mantle, as indicated by a recent numerical study that relates the CLIP event with subduction initiation when DIDS is established in this region (Riel et al., 2023).
We note that there is uncertainty in plate reconstructions on plate sizes and limited evidence on the timing of extension.Therefore, we must consider this comparison as somewhat more generic and qualitative rather than specific and quantitative.

Conclusion
Relative to the previous dual inward dipping subduction (DIDS) models, the 2-D thermo-mechanical models here demonstrate that DIDS, after considering temperature dependent rheology, can induce progressive weakening in an initially homogeneous overriding plate by lowering its viscosity and forming high strain rate necking centres within it.Three variables on plate sizes are investigated to understand what may control the maximum degree of plate weakening.It shows that the initial length (L 0 OP ) and thickness (H 0 OP ) of the overriding plate are Z.Lei and J.H. Davies negatively correlated with the maximum degree of weakening (Fig. 15).While the initial thickness of the subducting plate (H 0 SP ) positively relates to the maximum weakening level.The progressive weakening can result in a variety of irreversible stretching states ranging from 1) little or no lithosphere thinning and extension (<5% accumulation of strain), to 2) limited thermal lithosphere thinning (<30% accumulation of strain), and 3) localised rifting followed by spreading extension.
Relative to single-sided subduction with a mobile overriding plate, DIDS can reduce the magnitude of viscosity to a lower level within the overriding plate.It achieves this by effectively creating a dynamic fixed boundary condition for the middle (overriding) plate (Fig. 15).This inhibits the mobility of the plate and helps promote localised strain to accommodate the slab rollback tendency on both sides.Besides, when the initial length of the overriding plate is short enough (L 0 OP ≤ 1300 km), dual inward dipping subduction can form a united, thus stronger upwelling mantle flow which can reduce the viscosity of the overriding plate to a lower magnitude than models with longer L 0 OP .While the DIDS effect of self-consistently forming a fixed boundary condition and generating a stronger upwelling mantle flow are also reported in previous DIDS models, their role in promoting extension in the overriding plate has been neglected.This research implies that the generic connections between the DIDS plate sizes and plate weakening, along with the DIDS impact on limiting the trench mobility and generating stronger mantle flow are important aspects to consider to understand extensional deformation history in natural DIDS cases.
We also demonstrate that the incorporation of temperature dependent creep rheologies is critical to enable thermal-activated weakening in the lithosphere hotter than ~900 K.In addition to the strain rate weakening induced by yielding in the cold proportion of the lithosphere, a continuous viscosity reduction followed by the development of strain localisation is observed in the overriding plate as rifting and spreading extension develops.Both the temperature dependent creep deformation mechanisms and yielding deformation mechanism contribute significantly to the continuous viscosity reduction, which is also likely to be promoted by the positive feedback between plate weakening and strain localisation.Removing the thermal dependent creep deformation mechanisms for the overriding plate will create a strong mantle lithosphere that fails to develop rifting or spreading extension within it.This reveals why previous DIDS models, without considering creep rheology, fail to induce significant extension within the overriding plate.The finding may also have a broader rheological implication for more general processes that involve plate scale weakening, strain localisation and the formation of new plate boundaries.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Dual inward dipping model geometry and initial setup illustrated with the initial temperature field as the background.(a) The whole computational domain.Age 0 SP and Age 0 OP represent the initial ages of subducting plate and overriding plate at trench.The viscosity jump (Δμ) between upper and lower mantle at 660 km transition zone is set up with a fixed value of 30.To be noted, L 0 OP represents the initial length of the overriding plate measured by the lateral distance between the two trenches.(b) Enlarged area of the trench zone where the bending slab meets the flat overriding plate.The 1100 K and 1300 K isotherms are marked as white contours.The unstructured initial mesh after refinement is displayed as dark blue triangular meshes.(c) Vertical profile of viscosity against depth within the overriding plate.The plot line is ~300 km away from the initial trench.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. The simulation screenshots of model 'L 0 OP = 1700 km'.(a) Deviatoric normal stress component (τ xx ) where positive value represents tensile and negative value denotes compressive.(b) Temporal evolution of the dominating deformation mechanism, which is defined as the rheology law that yields the minimum magnitude of viscosity in a specific region.(c) Second invariant of strain rate (ε II ).The progressive weakening process within the overriding plate is demonstrated by the necking of the isoviscous contours.The 5 groups of yellow contours encompassing the plates in each screenshot are isoviscous contours of 10 20 , 10 21 , 10 22 , 10 23 , 10 24 Pa⋅s from outward to inward.The bold 'X' underlying the overriding plate denotes the absence of viscosity reduction for the isoviscous contour, whose value is stated on the far right-hand side of Fig. 2, i.e., 10 20 , 10 21 Pa⋅s for model 'L 0 OP = 1700 km'.Also, the first screenshot with bold 'X' notes the maximum viscosity reduction that model 'L 0 OP = 1700 km' can achieve is 10 21 -10 22 Pa⋅s throughout the 10 Myr simulation.The two sets of green solid lines are 700 K and 1300 K isotherm contours to image the thermal geometry of the plate.The transition zone at the depth of 660 km is marked by the horizontal white dashed line.The white right-angle scale bar lying above the right end of the transition zone represents 200 km in both directions.The bottom left corner caption shows the elapsed simulation time and bottom right corner is the name of the model.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Progressive weakening of the overriding plate during dual inward dipping subduction with decreasing length of the overriding plate, (a) model 'L 0 OP = 1700 km', (b) model 'L 0 OP = 1300 km', (c) model 'L 0 OP = 1100 km'.Section view is illustrated by the deviatoric normal stress field (τ xx ).Necking centres with minimum viscosity in the overriding plate are marked as black squares, which shows that the two separate necking regions tend to get closer and merge into a single one as the length of the overriding plate shortens from 1700 km to 1100 km.A detailed explanation of the contours and other symbols can be found in the caption of Fig.2.

Fig. 4 .
Fig. 4. Temporal evolution of viscosity magnitude and second invariant of strain rate along the depth of 5 km within the overriding plate for models with different initial length of the overriding plate.(a, c) Model 'L 0 OP = 1700 km ' .(b, d) Model 'L 0 OP = 1000 km ' .The edge of the filled contour in the lateral direction represents the interface between the overriding plate and subducting plate, i.e., the trench.The white arrows indicate the necking process or strain localisation within the overriding plate.β is stretching factor, calculated as the final OP length divided by the initial length at the depth of 5 km.

Fig. 5 .
Fig. 5. Progressive weakening of the overriding plate during dual inward dipping subduction with increasing thickness of the overriding plate.(a) Model 'H 0 OP = 67 km', (b) model 'H 0 OP = 74 km' and (c) model 'H 0 OP = 100 km'.Section view illustrates the deviatoric normal stress field (τ xx ).A detailed explanation of the contours and symbols can be found in the caption of Fig. 2.

Fig. 6 .
Fig. 6.Temporal evolution of viscosity along a horizontal line at the depth of 5 km in the overriding plate for models with different initial thickness of the overriding plate.(a) Model 'H 0 OP = 67 km'.(b) Model 'H 0 OP = 100 km'.The edge of the contour filling in the lateral direction represents the interface between the overriding plate and subducting plate.

Fig. 7 .
Fig. 7. Temporal evolution of diagnostics for models with varying initial overriding plate thickness.(a) viscosity (μ), (b) deviatoric normal stress component (τ xx ), (c) real-time lithosphere thickness (H OP ), (d) second invariant of strain rate (ε II ) against subduction rate (v SP ) at trench.In d, each scatter point represents 0.2 Myr over 10 Myr simulation.The downward triangles mark the timestep when the slab reaches the depth of 660 km.

Fig. 7
Fig.7-b, τ xx reduction may develop when spreading extension occurs in the overriding plate (Fig.7, b, blue and orange lines), while the other cases witness τ xx reduction only after slab reaches the depth of 660 km (Fig.7, b).These results confirm that the evolution of viscosity magnitude is a good indicator for the various progressive weakening developed in the overriding plate during subduction.In Fig.7-d, each scatter point represents 0.2 Myr simulation.The temporal path shows that as the initial thickness of the overriding plate reduces, εII in the overriding plate becomes increasingly sensitive to increasing subduction rate before the slab sinks to the depth of 660 km.Then the temporal path inflects upwards to higher εII for models that induce spreading extension (model 'H 0 SP = 67 km' and 'H 0 SP = 70 km'), and downwards for models that fail to spread.The upward inflection suggests that once extension initiates, it does not take much subduction rate to maintain the high strain rate developed in the overriding plate.As slabs slowly sink to the lower mantle, the magnitude of subduction rate and εII both gradually decrease, and the temporal path returns to the starting point for all models.In the models with varying H 0 OP , the maximum subduction rate ranges from 20 to 35 cm/yr, while the corresponding maximum εII in the overriding plate ranges from 10 − 15 to 10 − 12 s − 1 .The general positive correlation between these two variables suggests that maximum subduction rate may play an important role in weakening the overriding plate.Over the 10 Myr simulation, only ~15% scatter points fall in the range where subduction velocity is over 10 cm/yr, implying that high strain rate deformation correlated with high subduction rate only lasts for a short period.

Fig. 8 .
Fig. 8. Progressive weakening of the overriding plate during dual inward dipping subduction with increasing age of the subducting plate.(a) Model 'H 0 SP = 100 km'.(b) Model 'H 0 SP = 122 km'.(c) Model 'H 0 SP = 141 km'.A detailed explanation of the contours and symbols can be found in the caption of Fig. 2.

Fig. 9 .
Fig. 9. Temporal evolution of viscosity along a horizontal line at the depth of 5 km in the overriding plate with different initial thickness of the subducting plate.(a) Model 'H 0 SP = 100 km'.(b) Model 'H 0 SP = 141 km'.The edge of the contour filling in the lateral direction represents the interface between the overriding plate and subducting plate.

Fig. 10 .
Fig. 10.Key diagnostics used to characterise the rheology modification within the overriding plate.(a-c) Weakening level developed within the overriding plate.(df) Accumulation of strain in the middle of the overriding plate (5000 km away from the side boundaries).

Fig. 11 .
Fig. 11.Temporal evolution of horizontal velocity component along a horizontal slice, x coordinate from 3900 km to 6100 km, at the depth of 20 km from the surface.(a) model 'H 0 OP = 67 km', (b) model 'H 0 OP 100 km', (c) single sided subduction with a mobile overriding plate referring to (a), (d) model 'μ OP disl ' identical with (a) except that it excludes dislocation creep for the overriding plate.The contour filling represents the variation of horizontal component of velocity vector throughout the 10 Myr simulation.Positive value means right-ward motion and negative value is left-ward motion.The white area represents that the plate is nearly stagnant.And the edge of the white area marks the interface between the subducting plate and overriding plate or rifting and spreading centre within the overriding plate.SP and OP are short for subducting plate and overriding plate separately.

Fig. 12 .
Fig. 12. Intensity of mantle flow and its correlation with plate weakening in the overriding plate.(a-b) Vertical and horizontal velocity component of mantle flow along a slice at the depth of 75 km, ~8 km below the overriding plate, after 4.4 Myr simulation.Positive value means rightward or upward motion, and negative value represents leftward or downward motion.In (a), the horizontal coordinate of the necking centres in the overriding plate is plotted as black square to visualize its spatial correlation with upwelling component of mantle wedge flow.(d) Linear correlation in log scale between the progressive weakening in the overriding plate (ε II , defined by Eq. (13)) and the velocity gradient of the divergent mantle flow underlying the overriding plate midline (∇v x ).The maximum velocity gradient of divergent mantle flow in the range of t 660 for each model is marked by a vertical dashed line.

Fig. 13 .
Fig. 13.The temporal and spatial evolution of the dominant deformation mechanism in the overriding plate during progressive weakening, represented by the continuous necking of viscosity contours.The model in the left column is model 'H 0 OP = 67 km' (a-g).The two models in the right column are identical with the left one except that one excludes dislocation creep for the overriding plate (model 'μ OP disl ', h-l) and the other only keeps yielding deformation mechanism for the overriding plate (model 'μ OP Y only ', m).The dashed zoom-in block in (a) and (h) shows the location of screenshots in (b-g) and (i-l) separately.The progressive weakening process within the overriding plate is demonstrated by the necking of the isoviscous contours.The 5 groups of yellow contours encompassing the plates in each screenshot are isoviscous contours of 10 20 , 10 21 , 10 22 , 10 23 , 10 24 Pa⋅s from outward to inward.The bottom of the overriding plate is defined by the green isothermal contour of 1300 K.The bottom left corner caption shows the elapsed simulation time and bottom right corner is the name of the model.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .
Fig. 14.Temporal path of each dominant deformation mechanism over the background viscosity contour map along the midline of the overriding plate, i. e., the necking centre in model 'H 0 OP = 70 km'.The background viscosity contour map is calculated based on Eq. (7-9) with prescribed lithostatic pressure at the depth of 50 km.The phase diagram is divided by the white dashed lines into five domains based on the calculation of which deformation mechanism yields the minimum viscosity at the given second invariant of strain rate, temperature.Each scatter points represents a timestep.

Fig. 15 .
Fig. 15.Synoptic comparison of different model setup's role in affecting the necking behaviour developed within the overriding plate.(a) Single sided subduction.(b) Dual inward dipping subduction.(c) Thickness of the subducting plate or overriding plate (H 0 SP , H 0 OP ).(d) Length of the overriding plate (L 0 OP ).

Table 1
Key parameters used in this research.
a The rheology parameter of diffusion creep is guided by dry olivine deformation experiments

Table 2
List of model setup.

Table 3
Summary of diagnostics for all models.For further description of the diagnostics please see the main text.