Variable plate kinematics promotes changes in back-arc deformation regime along the north-eastern Eurasia plate boundary

The stretching of the lithosphere leading to back-arc basins formation generally develops behind arc-trench systems and is considered the consequence of slab retreat relative to the upper plate. Here, we examine the deformation regime evolution within the overriding plate due to subduction processes, using thermo-mechanical numerical simulations. We explore the north-eastern Eurasia plate boundary and the mechanisms of subducting Pacific plate since 57 Ma. During this time interval, several extensional basins formed along the Eurasia margin, such as the East China Sea, the Japan Sea, and the Kuril basin. Here, we increased the simulation complexity, with the inclusion of (i) the kinematic variability of the Pacific plate over the geological past with respect to a fixed Eurasia, incorporating time-dependent (i.e., temporally evolving) velocities computed from plate motion reconstructions; (ii) a Low-Velocity Zone within the asthenosphere, and (iii) a horizontal eastward mantle flow. Our results show a crucial role of the mantle flow for the development of lithospheric extension and back-arc basin opening, and a main kinematic control of the subduction trench position, which advances and retreats, into distance intervals in the order of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}∼ 100 km, and providing stages of compression and extension in a back-arc basin.

Table S1.Kinematic evolution of the Pacific plate with respect to Eurasia fixed.List of time at which the velocity was changed in the models 1,2 .

Figure S1
. Model setups.Panel A represents setup for Model 1.The numerical domain includes timedependent convergence velocities 1,2 , a LVZ between 100 and 200 km depth, and no mantle flow.Panel B represents setup for Model 2, which includes time-dependent convergence velocities 1,2 , a LVZ between 100 and 200 km depth, and a horizontal mantle flow applied with a constant velocity of 3 cm/a, which corresponds to half of the average convergence velocity of the Pacific plate [3][4][5] .Boundary conditions for top and bottom domain boundaries are free slip, whereas left and right domain boundaries are set to be periodic, to allow a throughgoing mantle flow, as in 6 .
In red the direction of application of the subducting plate velocities, whereas in blue the direction of application of the mantle flow.

Analysis of the horizontal deviatoric stress
As a comparison with Figure 5 in the main text, we computed the horizontal deviatoric stress (Sxx, Fig. S4) for Model 2, in a point situated at 350 km from the trench, and 25 km depth in the lithosphere, within the back-arc basin.
Figure S4 shows that for the first 2 Ma compression occurred.At 55 Ma compression started to decrease and converted to extension at 49 Ma.From 47 Ma extension prevailed, which decreased at 39 Ma, whereas at 30 Ma some weak inversions started to occur.These latter became more effective in the last phases of the model run, from 10 to 0 Ma.This evolution of Sxx agrees with the general behavior of Model 2 (Figs. 4 and 5c in the main text), whereas the observed differences derive from local state of stress within the lithosphere, in the basin: in fact, for consistency we kept the measurement point at the same distance from the trench for all timesteps.However, the state of stress within the lithosphere at subduction zones changes with time and location, and is also related to other features such as, for example, the slab dynamics and its interactions with mantle transition zones 10,11 .In this zoom on the trench area, the black squares represent the area of lithospheric weakness which forms where the maximum thinning occurs, at 32 Ma, and from which the recent (i.e., 6.0 Ma) incipient subduction with opposite subduction polarity originates.

Figure S2 .
Figure S2.Results of models with LVZ and no mantle flow.These model results show the last 57 Ma subduction evolution of the Pacific plate under Eurasia.It includes time-dependent plate kinematics for the Pacific plate relative to Eurasia by updated model by9 , which include corrections to the Pacific rotations prior to 83 Ma 1,2 .These results show no back-arc basin opening throughout the entire subduction duration, which is inconsistent with the geodynamic history of the eastern margin of the Eurasia plate.From velocities, at depths, no counterflow within the mantle is generated to induce slab retreat, and eventually leading to back-arc basin opening.

Figure S3 .
Figure S3.Results of models with LVZ and mantle flow.An extensional phase can be detected, beginning at ∼ 46.5 Ma (panel a), with a developed back-arc already at ∼ 32 Ma (panel b).At 32 Ma, a compressional phase starts (panel c), ended at ∼ 23 Ma, when a new extensional phase begins.At ∼ 15 Ma (panel d) a new inversion of the extensional trend starts.A very short opening trend starts at ∼ 7 Ma (panel e), whereas at ∼ 6 Ma a new inversion starts (panel f).Panel gshows the Present-day state of the north-eastern Eurasia margin.Velocities within the mantle show a slab which is entirely influenced by the eastward mantle flow in the initial opening phases of the back-arc basin, whereas it influences the slab mainly below 200 km in the subsequence phases.In fact, where the LVZ decollement level is located (100-200 km), the velocity field has no unique direction and shows, in the subduction wedge, several phases of westward and eastward directions.Moreover, a stagnation is observable at about the 660 km discontinuity.The black squares point to the area of lithospheric weakness which forms where the maximum thinning occurs, at 32 Ma, and from which the recent (i.e., 6.0 Ma) incipient subduction with opposite subduction polarity originates.The red squares point to the extended area in the continental upper plate.

Figure S4 .
Figure S4.Evolution of the horizontal deviatoric stress (Sxx, black solid line, with circles) for Model 2.Sxx>0 represents horizontal extension (h-extension), whereas Sxx<0 represents horizontal compression (h-compression).The zero in the figure corresponds to the dashed black line, whereas red and grey areas correspond to trench advance and trench retreat phases, respectively, of Fig.5in the main text.The white area is for the stable trench location (stable TL, Fig.5, main text).

Figure S5 .
Figure S5.Log10 of the second invariant of the strain rate resulted from our numerical Model 2. Panel a shows the domain between 1000 and 1500 km.Every other panel shows the domain between 1500 and 2000 km.Basically, each panel shows the strain rate in the middle of extension (B, D, F) and compression (C, E, G) phases.Higher values of the second invariant in the trench area during the phases of the Pacific plate velocity increase (i.e., from 32 Ma) are shown.In this zoom on the trench area, the black squares represent the area of lithospheric weakness which forms where the maximum thinning occurs, at 32 Ma, and from which the recent (i.e., 6.0 Ma) incipient subduction with opposite subduction polarity originates.

Table S2 .
Rheological parameters.Rheological 7 and thermal 8 parameters of materials used for the experiments.

evaluated since
the short time of subduction evolution (~6-10 Ma).