A possible mechanism for spontaneous cyclic back-arc spreading

Back-arc spreading is a non-steady-state process exemplified by the repeated cycles of spreading of the South Fiji and the Lau Basins behind the Tonga arc, and the Parece Vela Basin and the Mariana Trough behind the Mariana arc. Spreading in these regions starts with rifting within the volcanic arc before shifting to the back-arc region where it develops into a phase of well-defined spreading. 2D thermo-mechanical subduction modeling incorporating phase transitions at depths of 410 km and 660 km suggests the presence of a low-viscosity and low-density mantle wedge is an important condition for arc rifting to occur. Back-arc spreading starts when a nearly vertical slab impinges upon the 660 km discontinuity causing downdip compressive stress that is transmitted up the slab resulting in extensional within-arc stress. Trench retreat during a phase of back-arc spreading causes a decrease in slab dip angle and buckling of the slab. Back-arc spreading ceases during this buckling phase. Rifting starts once more when the nearly vertically dipping ‘heel’ of the buckled slab again impinges upon the 660-km boundary. The second phase of rifting initially focuses within the arc but subsequently shifts to the back-arc region leading to renewed back-arc spreading. Our modeling predicts that subduction of thick (old age) and weak (small yield stress) slabs, which have intermediate resistance to slab bending, leads to cyclic back-arc spreading. In contrast, continuous back-arc spreading is predicted for thick and strong slabs with a large resistance to bending, and no back-arc spreading is predicted for slabs with a small resistance to bending (thin slabs). Geological processes such as toroidal mantle flow around the lateral edges of a slab, collisions with buoyant lithosphere and interactions with third plates may have important roles in the development of cyclic back-arc spreading in specific cases. However, the presence of a common timescale of ~ 20 Myr suggests there a general underlying control on back-arc basin formation that is common to many if not all subduction zones. The new model presented here can account for the main features of cyclic back-arc spreading seen in the Tonga-Kermadec and the Calabrian arcs.


Introduction
The formation of back-arc basins due to back-arc spreading ranks with arc volcanism and the generation of large subduction megathrust earthquakes as a first-order characteristic of subduction zones. In general, backarc spreading is not a steady-state process (Sdrolias and Müller 2006;Clark et al. 2008) and there are well-documented examples where repeated phases of back-arc spreading are recognized, such as the South Fiji and the Lau Basins behind the Tonga arc, and the Parece Vela Basin and the Mariana Trough behind the Mariana arc (Sdrolias and Müller 2006). These examples of cyclic back-arc spreading show some key features in common, suggesting that they are controlled by the same underlying physical process. In particular, we emphasize that (i) spreading events start with rifting focused within the arc before shifting to the back-arc domain; and (ii) each Page 2 of 19 Ishii and Wallis Progress in Earth and Planetary Science (2022) 9:27 spreading event is punctuated by a pause in the spreading (Clark et al. 2008). The main driving force for subduction of oceanic plates is the negative buoyancy of the subducting slab (Forsyth and Uyeda 1975). Numerical mechanical models of subduction zones show that the negative buoyancy of the subducting slab tends to cause sinking of the slab away from the subduction zone-commonly referred to as slab rollback-and this tends to result in an extensional stress state of the overriding plate. In more detail, this stress state depends on several conditions (Nakakuki and Mura 2013;Schellart and Moresi 2013;Holt et al. 2015): (1) the extensional stress is larger for an overriding plate that is fixed to the sidewall of the model than for an overriding plate with an unpinned trailing edge; (2) mantle corner flow induced by the subducting slab drags the overriding plate toward the trench resulting in compressive stress in the fore-arc area; and (3) toroidal mantle flow around the lateral edges of a slab enhances the slab rollback leading to larger extensional stress in the back-arc lithosphere for a shorter width of slab. To cause rifting within the arc, this extensional stress state needs to act on locally weakened lithosphere. The weakening of arc lithosphere is thought to be achieved by thinning (due to mechanical or thermal erosion) and/or injection of melt and fluid (Arcay et al. 2005;Gerya and Meilick 2011;Nakao, et al. 2016).
The boundary between the upper and lower mantle-the 660-km boundary-also affects the subducting process due to a resistance caused by the endothermic phase transition and an associated increase in viscosity (Faccenda and Dal Zilio 2017;Li et al. 2019 and references therein). Recent high-resolution seismic tomography reveals a significant variation of slab geometries in this domain characterized by two endmembers: (i) stagnation and flattening at the 660-km boundary and (ii) penetration into the lower mantle (Li et al. 2008;Fukao and Obayashi 2013;Goes et al. 2017;Li et al. 2019). The positive correlation recognized between the length of stagnant slab at 600 km depth and the distance of trench retreat (Goes et al. 2017) indicates a close mechanical relationship between the trench motion and the interaction between the subducting slab and the 660-km boundary. Numerical models have shown that the interaction can cause a cyclic change in trench motion and/or in the deformation of the overriding plate (Čížková and Bina 2013;Nakakuki and Mura 2013;Schellart and Moresi 2013). However, previous models have not satisfactorily accounted for the cyclicity of back-arc spreading, which is the main focus of this study.
In this study of cyclic back-arc spreading, we focus on the conditions necessary for the back-arc spreading to occur, and the timing of both the start and finish of each phase of spreading. Our approach is based on two-dimensional thermo-mechanical subduction modeling with 410 km and 660 km phase transitions. In this model, a mantle wedge with a relatively low viscosity and low density is assumed to induce the initial stages of arc rifting. In addition, to evaluate the effect of mechanical interaction between the slab and the 660-km boundary, the age of the subducting plate, the maximum yield strength of the slab and the strength of the arc lithosphere are all treated as variables. Earlier studies have emphasized the importance of flow around the edges of subducting slabs. Our new approach shows that incorporating other recognized features of subduction-in particular the influence of buckling due to the interaction of slabs with the 660 km discontinuity-can account for many of the features of back-arc spreading, including its cyclicity, without appealing to more complex 3D effects.

Model
The governing equations for the flow and temperature fields are the conservation equations of mass, momentum and energy, with the extended Boussinesq approximation. Conservation of mass is approximated by the incompressible continuity equation: where v x and v z are the horizontal and vertical components of the velocity vector, respectively. The two-dimensional Stokes equations for creeping flow take the form where σ xx , σ zz and σ xz are the components of the deviatoric stress tensor, P is the pressure, g is the acceleration due to gravity, and density ρ is given by the expression (Gerya et al. 2021) where ρ 0 is the reference density at the reference pressure P 0 and temperature T 0 . α and β are the thermal expansion and compressibility coefficients, respectively. These values are shown in Table 1.
The energy equation takes the form where q x and q z are horizontal and vertical heat fluxes, and ε xx , ε zz and ε xz are components of the strain rate tensor. The terms H L , H r , H a and H s are the values of latent heat for phase transitions, radiogenic heating, adiabatic heating and viscous dissipation, respectively, and all have units of W/m 3 . Q L is the latent heat per unit mass, γ is the Clapeyron slope, and ∆ρ is the density difference between the phases (Christensen 1998;Faccenda and Dal Zilio 2017). We assume ∆ρ 410 = 124 kg/m 3 and ∆ρ 660 = 193 kg/ m 3 , to account for the olivine phase transitions (olivine to wadsleyite, and ringwoodite to bridgmanite and ferropericlase) and ~ 60% modal content of olivine in peridotite (Faccenda and Dal Zilio 2017, Fig. 1a). These changes in density are derived from the values of ρ 0 given in Table 1. The Clapeyron slopes of the 410-km and 660km boundaries are assumed to be γ 410 = 2 MPa/K and γ 660 = − 2 MPa/K, respectively. The isobaric heat capacity C P is assumed to be 1200 J/kg/K. The values of thermal conductivity k and H r are listed in Table 1. To incorporate the effects of the fluid circulation, the upward migration of melt and the associated rapid transport of heat around the back-arc spreading center, the present model assumes a high effective thermal conductivity (k eff ) for the upper mantle of the overriding plate within the domain z < 20 km. The k eff is assumed to vary horizontally as a function of temperature at z = 20 km (T z = 20 km ).
where T min = 600 °C and T max = 1200 °C are minimum and maximum cutoff temperature, respectively. The values of the thermal conductivities, k 0 and k 1 , are 3.1 W/ mK and 10 W/mK, respectively. Therefore, k eff falls in the range 3.1-13.1 W/mK. The constitutive relationship is given by where ε ij is the deviatoric strain rate tensor and the effective viscosity, η eff , is determined from the minimum equivalent viscosity among the viscous flow, the frictional yield stress and maximum yield stress: The values used for the rheological parameters are listed in Table 2. For viscous flow, we assume that the power-law rheology that is associated with dislocation  (2002) *4 Peacock and Wang (1999) *5 Linearly increase with the depth of 15-60 km to incorporate the effect of eclogitization *6 Eq. 6 is applied for z > 20 km to incorporate the effect of heat transport due to fluid and/or melt migration * 7 Upper mantle for subducting plate (purple in Fig. 1b  where A dis , E dis , V dis and n are a material constant, the activation energy, the activation volume and the stress exponent for dislocation creep, respectively. A N , E N and V N represent the values used for the Newtonian approximation. R is the universal gas constant. The frictional stress and the equivalent viscosity are taken to be (9) where c and φ are the cohesion and effective coefficient of friction, respectively, and ε II is the second invariant of the strain rate tensor. To incorporate the effect of strain weakening, values of c and φ decrease linearly from c inital and φ inital to c weakened and φ weakened as strain values increase from 0 to 0.5. For strain larger than 0.5, c weakened and φ weakened are used. A weak frictional strength of the oceanic crust facilitates decoupling between the subducting oceanic plate and the overriding plate. The model treats the maximum yield stress as a variable parameter: Finally, the lower and upper bounds for the effective viscosity are set to be 10 17 and 10 25 Pa s, respectively.
The values of activation energy (E N ) for the crust, the upper mantle and the mantle transition zone are inferred from experimentally derived rheological parameters ( Table 2). The rheological parameters (A N in Eq. 9) for crustal domains are determined so that viscosity values are nearly equal to the values of dislocation creep at the relevant strain rate (~ 10 -14 /s) and temperature (300-400 °C for quartz and 400-600 °C for plagioclase) assuming V N = 0 cm 3 /mol. Other values of the rheological parameters, including those for the lower mantle, are set to be consistent with the viscosity stratification estimated from geophysical observations (King 2016 and references therein). These values result in the rheological layering shown in Fig. 1a, that includes the viscosity jump of ~ 35 at the 660-km boundary and a viscosity peak in the middle of the lower mantle.
Dehydration of the subducting oceanic plate causes hydration and partial melting of the mantle wedge, both of which affect the physical properties of this domain. The mantle wedge has been estimated to have both a lower viscosity (by at least a factor 10) and a lower density (density difference of ~ 20 kg/m 3 ) than the asthenosphere based on comparisons between observations and dynamic modeling of the Tonga-Kermadec subduction zone (Billen and Gurnis 2001;Billen et al. 2003). In addition, our model uses a low bulk effective activation energy for the mantle wedge, because the upward migration of melt and fluid should result in a relatively homogeneous distribution of temperature and therefore also effective viscosity, which is strongly controlled by temperature.
The analysis uses the code I2VIS (Gerya and Yuen 2003), which is based on finite differences with a marker-in-cell technique and a fully staggered rectangular Eulerian grid. The numerical model domain (7000 km wide × 2900 km deep) includes non-uniform σ yield = 200, 500MPa, and (11) η yield = σ yield 2ε II .  645 × 222 rectangular grids with a resolution varying from 4 to 50 km. The upper central area (x = 3400-5200 km, z = 0-640 km) has the highest resolution of 4 × 4 km. The model includes ~ 5.1 × 10 6 markers. Two plates are defined in the initial state: a subducting oceanic plate and an overriding continental plate (Fig. 1b). The oceanic plate includes oceanic crust with a thickness of 8 km, and the continental plate includes upper and lower continental crustal domains both with thicknesses of 15 km. The initial temperature distribution of oceanic plate follows a half-space cooling model, using a thermal diffusivity of κ = 0.78 × 10 -6 m 2 /s. The initial temperature distribution of continental plate follows a 1D steady-state conduction model incorporating radiogenic heat production that yields a surface heat flow of 60 mW/m 2 . In the domain below the plate, the temperature follows an adiabatic profile with a potential temperature of 1300 °C and a thermal gradient that decreases exponentially with depth from 0.45 °C/km (z = 0 km) to 0.27 °C/km (z = 2900 km) (Fig. 1a). The thermal age of the oceanic plate is set to be 40, 60 or 100 Ma at the trench (A tr ) for the start of the calculation. A nearly constant thermal age of the oceanic plate at the trench is maintained throughout the calculation, by defining an appropriate decrease in the age of the oceanic plate toward the right side of the model in the initial state (Fig. 1b). The following thermal boundary conditions are used for the calculation domain: isothermal at the top (T = 0 °C) and bottom (T = 2236 °C) and zero heat flux through the two insulating sidewalls. The velocity boundary conditions are free-slip everywhere (Fig. 1c).
During the initial stage of the simulation, constant velocities are applied to the subducting plate (V sp = 5 cm/ yr) (Fig. 2a). This velocity is applied to grids with x = 6200 km and z = 0-100 km. In order to make the plates movable, the viscosity of the rectangular area in the upper right corner is reduced to 10 19 Pa s (Fig. 1c). In addition, the material and temperature distributions of the upper right corner domain (800 km wide × 300 km depth) are reset every simulation step to keep the distribution of temperature and material in the oceanic plate that will move toward the trench constant (Fig. 1c). The distinction between oceanic crust and mantle is maintained until reaching a depth of 150 km where both domains are treated as upper mantle (Fig. 2b). Our simulation also incorporates several changes that occur when the slab tip first reaches a depth of 250 km. These changes help model the effects of fluid migration and melting. The four changes, or new model conditions, are as follows: (1) the mantle wedge domain is replaced by materials with a lower viscosity and density; (2) the frictional properties for mantle and crust above the mantle wedge are decreased to c = 1 MPa and φ = φ AL , irrespective of a strain magnitude; (3) the depth at which the oceanic crust is treated as the same as the upper mantle is changed to 100 km; and (4) the subducting plate becomes free to move in response to the distribution of density and viscosity (Fig. 2c).
Viscous coupling between the subducting plate and the overriding mantle drives flow in the mantle wedge domain. Observation of surface heat flow combined with subduction zone thermal modeling studies has shown that at relatively shallow depths the slab and the wedge mantle are decoupled over geological timescales, and the transition from a decoupled to a coupled interface is considered to occur at a depth of ~ 80 km irrespective of the subduction zone (Wada and Wang 2009). Subduction of the oceanic plate leads to the dehydration of both the oceanic crust and slab mantle. The released fluid migrates upward and causes hydration and flux melting of the overriding wedge mantle. These processes result in vigorous arc volcanism and lead to the strong contrast between the hot sub-arc mantle and the cold fore-arc nose (Abers et al. 2017;Wada and Wang 2009). Attempting to incorporate the dehydration, fluid migration and melting directly into models of subduction zones (cf. Nakao et al. 2016;Li et al. 2019) is a complex task and subject to considerable uncertainties. Here, we choose a simpler approach to modeling these processes based on conditions 1-4 listed above. Condition (1) is implemented to incorporate the effect of hydration, partial melting of wedge mantle and upward migration of fluid and melt. Condition (2) is implemented to incorporate a weakening of the arc lithosphere due to additions of fluid and melt, where φ AL is treated as a variable parameter. Conditions (1) and (3) are implemented to reproduce a nearly constant decoupling-coupling transition depth along the subduction interface and the associated flow pattern in the mantle wedge. Note that a large effective thermal conductivity (Eq. 6) and φ AL is not applied to the upper mantle for the subducting plate (indicated by purple in Fig. 1b and 2), and the depth of the slab tip and the location of the slab surface are determined from the distribution of the upper mantle domain of the subducting plate.
Low-viscosity and low-density properties used to model the mantle wedge (Tables 1 and 2) are applied to a mantle domain bounded by the slab surface at its base and a vertical projection of the segment of slab located at depths 80-200 km. The upper boundary of the mantle domain is defined by a depth that varies between a minimum of 30 km below the arc that gradually increases to 130 km in the far limits of the back-arc region ( Fig. 2c and d). To maintain this definition, the mantle wedge domain is adjusted every calculation step according to the change in slab geometry ( Fig. 2c and d). The horizontal range of the mantle wedge used in our model is supported by Fig. 2 Calculation procedure. a Initial conditions. b Distinction between oceanic crust and upper mantle is removed at > 150 km depth. c When the slab tip reaches a depth of 250 km, (1) area of mantle wedge (bounded by the slab surface, a horizontal range between the location of the slab surface at depths of 80-200 km and depths z > 30 km) is replaced by materials with lower viscosity and density, and (2) the frictional properties for mantle and crust above the mantle wedge are decreased to c = 1 MPa and φ = φ AL . After this time, (3) the depth where oceanic crust is no longer distinguished from mantle is changed to 100 km, and (4) the subducting plate is free to move in response to the underlying dynamics. d Geometry of the mantle wedge changes according to the calculated slab geometry. A red circle indicates the position of the marker used to measure trench velocity (V tr ) Page 7 of 19 Ishii and Wallis Progress in Earth and Planetary Science (2022) 9:27 observations of the locations of arc volcanoes and the underlying slabs. The slab surface depths below volcanic fronts are mainly in the range 70-170 km (England and Katz 2010). In addition, arc volcanoes in NE Japan and northern Chile are distributed over a wide horizontal distance corresponding to range of underlying slab surface depths of 100-200 km (Yoshida et al. 2014;Tassara et al. 2006). This range is consistent with the estimated depth range for slab dehydration: Dehydration of oceanic crust starts at the transition between the decoupled and coupled interface domains (~ 80 km depth), and the maximum dehydration depth of the slab mantle is ~ 200 km (van Keken et al. 2011). The depth for the upper boundary of the mantle wedge used in our model is based on geophysical observations. For instance, in NE Japan seismic tomography reveals the presence of a low-velocity zone in the middle of the mantle wedge that deepens toward back-arc side from the uppermost mantle below the arc (Nakajima et al. 2001;Huang et al. 2011). The low-velocity zone is interpreted to be a domain with a high melt content and may reflect the flow pattern in the mantle wedge (Nakajima et al. 2001;Kogiso et al. 2009).
Our two-dimensional model only considers the interaction between a single subducting oceanic plate and a corresponding overriding continental plate and ignores any three-dimensional effects such as toroidal flow around the slab edge (e.g., Schellart and Moresi 2013) and resistance from side plates (e.g., Clark et al. 2008). The continental plate has a low velocity, and thus its motion is assumed to be controlled by global dynamics rather than interactions with the subducting plate. Therefore, the overriding plate is fixed to the sidewall of the model. In contrast, the subducting oceanic plate is assumed to be free to move in response to the underlying dynamics including the effects of negative buoyancy of the slab. Figure 3 shows the time evolution for A tr = 100 Ma, σ yield = 200 MPa and φ AL = 0.1 (case 1) as a reference case. The velocities of the overriding plate (V op ) and subducting plate (V sp ) are determined at locations x = 3000 km and x = 6200 km, respectively. The trench velocity (V tr ) is measured by a marker in the fore-arc domain (red circle in Fig. 2). Because the overriding plate is fixed to the left margin, a positive V op value indicates stretching of the overriding plate, and the difference between V tr and V op indicates the velocity of the intra-or back-arc spreading.

A typical case for back-arc spreading
After the initial subduction with an imposed velocity (V sp = 5 cm/yr), the slab accelerates due to the additional negative buoyancy from the upward deflection of the exothermic phase transition when the slab tip passes through the 410-km boundary (Fig. 3h). As the slab tip approaches the 660-km boundary, the slab slows down due to the combined resistance of the endothermic phase transition at the 660-km boundary and an associated increase in viscosity ( Fig. 3a and h). A torque resulting from the resistance at the 660-km boundary and the negative buoyancy of the shallower part of slab causes buckling of the slab (Fig. 3b). Irrespective of the buckling, downdip compressional stresses due to the resistance at the 660-km boundary are transmitted up the slab, indicated by the increased extensional stress within the arc area after the slab tip reaches the 660-km boundary ( Fig. 3h and i). This is the main cause of the rifting located within the arc (Fig. 3b) that gradually transitions to back-arc spreading (Fig. 3c).
During back-arc spreading, a large fraction of the subduction rate is accommodated by trench retreat (V sp < V tr ). The decrease in slab dip angle that results from the trench retreat ( Fig. 3c and h) causes a second phase of buckling of the slab (Fig. 3d). After this type of buckle has formed, the shallower portion of the slab resumes a largely free descent, and back-arc spreading stops (Fig. 3h). The present model assumes a high effective thermal conductivity for the back-arc lithosphere (Eq. 6) as a way to incorporate the effects of the fluid circulation, the upward migration of melt and the associated rapid transport of heat. The rapid cooling in the region around the locus of back-arc spreading, or back-arc ridge, restores the strength of the back-arc lithosphere during the period of no back-arc spreading. Therefore, when the steeply dipping 'heel' of the buckled slab reaches the 660-km boundary, rifting restarts within the arc instead of causing renewed spreading in the same back-arc ridge region (Fig. 3e). The stress state of the arc area is neutral during the period of back-arc spreading and is extensional during the period of no back-arc spreading and the stage of arc rifting ( Fig. 3h and i). When the backarc spreading stops, the velocity of the subducting slab increases rapidly, resulting in an increase in trenchward flow velocity of the asthenosphere below the overriding plate ( Fig. 3d and g). This flow drags the plate causing a short period of weekly compressional stress state of the arc (Fig. 3h and i). In our modeling, the same sequence of events repeats two more times (Figs. 3f, g, 4a and d) and cyclic backarc spreading is accompanied by oscillatory changes in velocities of the subducting plate (V sp ) and trench (V tr ) (Fig. 3h). The slab dip angle also shows cyclical changes: increasing during the period of no back-arc spreading and decreasing during the period of back-arc spreading (Fig. 3h). During the cyclic back-arc spreading associated with trench retreat, the folded slab near the 660-km boundary sinks very slowly. After the third back-arc spreading, the 'heel' of the buckled slab penetrates the 660-km boundary without causing a fourth phase of back-arc spreading ( Fig. 4a and d).

Effects of slab age
For A tr = 60 Ma, σ yield = 200 MPa and φ AL = 0.1 (case 2), the model shows two cycles of back-arc spreading, and each spreading starts when the 'heel' of the buckled slab reaches the 660-km boundary ( Fig. 4b and e). The amount of back-arc spreading and trench retreat is smaller than that for case 1. For A tr = 40 Ma, σ yield = 200 MPa and φ AL = 0.1 (case 3), the arc lithosphere only shows slight extension and lacks any clear back-arc spreading ( Fig. 4c and f ). Because subduction starts from a trench that is essentially fixed in location, the geometry of the folded slab in the lower mantle shows an overall near vertical dip.

Effects of maximum yield stress
For A tr = 100 Ma, σ yield = 500 MPa and φ AL = 0.1 (case 4), continuous and steady back-arc spreading is predicted. The slab dip angle in the upper mantle shows little change, and the slab flattens at the 660-km boundary ( Fig. 5a and  d). For A tr = 60 Ma, σ yield = 500 MPa and φ AL = 0.1 (case 5), the model shows three cycles of the back-arc spreading ( Fig. 5b and e). The velocity of back-arc spreading is smaller, and the period of back-arc spreading is longer than the reference case. For A tr = 40 Ma, σ yield = 500 MPa and φ AL = 0.1 (case 6), the model shows two cycles of the back-arc spreading, although the amount of back-arc spreading is very small (Fig. 5c and f ).

Effects of strength of arc lithosphere
For A tr = 100 Ma, σ yield = 200 MPa and φ AL = 0.2 (case 7), back-arc spreading starts later than the reference case (φ AL = 0.1), and the model shows three cycles of back-arc spreading although the extent of the third spreading is small (Fig. 6a and c). For A tr = 100 Ma, σ yield = 200 MPa and φ AL = 0.05 (case 8), back-arc spreading starts earlier than the reference case, and the model shows continuous spreading ( Fig. 6b and d).

Effects of Clapeyron slope of phase transitions
To examine the effect of the Clapeyron slope, results for different values of the Clapeyron slope are shown in Additional file 1: Fig. S1. The effect is relatively small. A case with a larger γ 410 results in tighter folded slab, and a case with a smaller absolute value of γ 660 results in earlier slab penetration into the lower mantle.

Effects of viscosity and density of mantle wedge
To examine the effects of imposed values of viscosity and density for the mantle wedge, we performed calculations with different values of parameters (A N in Eq. 9 and ρ 0 in Eq. 4) for the mantle wedge. A N = 1 × 10 -8 /s/MPa and A N = 6 × 10 -8 /s/MPa correspond to 2 and 1/3 times more viscous than the reference model, respectively (Additional file 1: Fig. S2a and c). Densities of ρ 0 = 3512 kg/m 3 and ρ 0 = 3472 kg/m 3 correspond to density contrast of ∆ ∆ρ 0wedge = 0 kg/m 3 and ∆ρ 0wedge = 40 kg/m 3 , respectively (Additional file 1: Fig. S2b and d). The reductions in both viscosity and density have similar effects. The results for cases of a high-viscosity wedge and ∆ρ 0wedge = 0 kg/m 3 show three cycles of back-arc spreading; the third cycle for high-viscosity wedge (Additional file 1: Fig. S2a) and the second and third cycles for the high density wedge (Additional file 1: Fig. S2b) all start at the site of the ridge of the previous spreading. For cases with a low-viscosity wedge and ∆ρ 0wedge = 40 kg/m 3 , back-arc spreading starts earlier than the reference case, and two continuous cycles of spreading occur without a pause followed by another two cycles of spreading (Additional file 1: Fig. S2c and d). The slabs at 660-km boundary formed by the continuous spreading are flatter than that in the reference case, causing a slower rate of foundering of the flattened slab.

Tests of model robustness
The results for resolution tests are shown in Additional file 1: Fig. S3. The model with low resolutions (6 × 6 km in upper central area) shows earlier start of back-arc spreading and thus flatter slab around the 660-km boundary than the reference case (Additional file 1: Fig. S3a). As a result, the slab sinks slower than the reference case and the back-arc spreading repeats four times. The model with high resolutions (3 km width × 2 km depth in upper central area) shows a very similar evolution to the reference case except for the later start of back-arc spreading and earlier slab penetration into the lower mantle than the reference case (Additional file 1: Fig. S3a). The present model uses the standard formulation of I2VIS, in which densities are defined as a function of temperature, pressure and composition (e.g., Gerya et al. 2021). However, a numerical model with the extended Boussinesq approximation generally assumes pressure-independent density (Tosi et al. 2010). To examine the effect of pressure dependency of density in the model, we performed calculation with β = 0 in Eq. 4. In this calculation, the physical parameters are changed to keep the viscosity stratification and the density difference of the phase transitions (∆ρ 410, ∆ρ 660 ) unchanged (Additional file 1: Fig. S4a). The result shows a very similar evolution to the reference case (Additional file 1: Fig. S4b and c).
The constant velocity imposed on the subducting plate (V sp = 5 cm/yr) during the early stage of subduction may cause extra viscous dissipation (H s ). To examine this possible artificial effect, we performed a calculation without viscous dissipation until the slab tip reaches a depth of 250 km. The result shows an almost identical evolution to the reference case (Additional file 1: Fig. S5), indicating that this effect is negligible.

Mechanism of back-arc spreading
Subduction of the oceanic plate is mainly driven by negative buoyancy of the subducting slab which is controlled by the age of the slab: Older and therefore cooler subducting slabs have a larger negative buoyancy. Old oceanic plates are also thicker than young plates and thus exhibit a greater resistance to bending. For an old slab of A tr = 100 Ma (case 1), after the slab tip reaches the 660-km boundary, downdip compressional stresses due to the resistance at the 660-km boundary are transmitted up the slab and the extensional stress of the arc area increases and the back-arc spreading starts (Fig. 3b, h and i). In contrast, for a young subducting plate (A tr = 40 Ma, case 3), the downdip compressional stresses due to the resistance at the 660-km boundary are largely absorbed  Ishii and Wallis Progress in Earth and Planetary Science (2022) 9:27 by the buckling of the slab that has a smaller resistance to bending (Fig. 4f ). The low viscosity and low density of the mantle wedge facilitate upward flow in the wedge corner and assist rifting within the arc (Fig. 3b and e). At the beginning of arc rifting, the width of the mantle wedge is narrow due to the steep dip of the slab (Fig. 3b and e). The transition from arc rifting to back-arc spreading is accompanied by a decrease in slab dip angle, an increase in the width of the mantle wedge and a backward migration of the locus of upward flow with respect to the trench position ( Fig. 3c and f ). Models without a mantle wedge but with weak arc lithosphere (c = 1 MPa and φ AL = 0.1 after the slab tip reaches 250 km depth) show neither an upward flow beneath the arc nor back-arc spreading (Additional file 1: Fig. S6). A model lacking a viscosity jump at the 660-km boundary also does not exhibit backarc spreading (Additional file 1: Fig. S7). These results indicate that both the existence of the low-density and low-viscosity mantle wedge and the viscosity jump at the 660-km boundary play significant roles in the development of back-arc spreading. In addition, models with higher values of viscosity and density for mantle wedge show enhanced arc rifting and reducing these values suppresses arc rifting (Additional file 1: Fig. S2).
In summary, the negative buoyancy of the oceanic slab drives the subduction. Upon encountering the resistance of the 660-km boundary, oceanic plates with different Page 13 of 19 Ishii and Wallis Progress in Earth and Planetary Science (2022) 9:27 ages continue the subduction process in different ways: Subduction of young oceanic plates is associated with buckling of the slab without associated back-arc spreading, while subduction of old oceanic plates is associated with slab rollback accompanied by back-arc spreading. This difference in behavior is mainly due to the difference in the resistance to bending of the slab. However, case 4 (σ yield = 500 MPa and φ AL = 0.1) has a larger σ yield and φ AL than case 8 (σ yield = 200 MPa and φ AL = 0.05), and both cases show continuous back-arc spreading (Figs. 5a  and 6b). In addition, case 2 (A tr = 60 Ma and φ AL = 0.1) has a smaller A tr and φ AL than case 7 (A tr = 100 Ma and φ AL = 0.2), and both cases show two cycles of back-arc spreading with relatively shorter time duration, and additional small spreading or stretching (Figs. 4b and  6a). These results indicate that the balance between the strength of arc lithosphere (φ AL ) and the resistance of slab bending (A tr and σ yield ) control the subduction processes. The back-arc spreading starts when the nearly vertical slab impinges upon the 660-km boundary. The back-arc spreading stops when the slab dip angle decreases to a minimum which marks the timing when the slab starts to buckle (Fig. 3). After this stage, the slab in the upper mantle is in a state of largely free descent and does not 'feel' the resistance at the 660-km boundary.
Our model assumes a high effective thermal conductivity for the back-arc lithosphere to incorporate the effect of the heat transport by the fluid circulation and the upward migration of melt (Eq. 6). Models that lack this domain of high effective thermal conductivity, but otherwise have the same conditions as the reference case, show repeated back-arc spreading with the same locus of spreading and lack any jumps in the location of the spreading center (Additional file 1: Fig. S8). For this reason, the quiescent interval between the back-arc spreading events is shorter than the reference case (Fig. 3). These observations indicate that rapid cooling of back-arc ridge area and restoration of the lithosphere strength contribute to the next phase of back-arc spreading starting from arc rifting and that the interval between the back-arc spreading events is the time required for the descending buckled slab to reach the 660-km boundary and the rifting within the arc to begin.

Comparison with previous models
Schellart and Moresi (2013) used 2D and 3D numerical models to suggest that where the overriding plate is free to move, trench retreat occurs cyclically. During periods of subduction hinge retreat, the slab becomes passively draped along the 660 km discontinuity, which these workers assume to be impenetrable. This stage is followed by a phase where the slab is tightly folded during which the trench remains largely immobile. In contrast to the results of Schellart and Moresi (2013), in our study the slab-flattening occurs due to buckling and descent of the slab during a period of no back-arc spreading (Fig. 3d and  e). While the model of Schellart and Moresi (2013) used a constant and relatively low viscosity for the lithosphere, our models used a temperature-dependent viscosity with a relatively low yield stress, which results in a highly heterogeneous distribution of effective viscosity within the slab. The constant and relatively low viscosity of the slab in the model of Schellart and Moresi (2013) results in a smaller resistance to bending of the slab than ours. The 2D numerical model of Nakakuki and Mura (2013) takes a similar approach to ours in using a temperaturedependent viscosity and a penetrable 660-km boundary. Their model also incorporates a weak rheology around the arc area that plays a role in back-arc spreading for a case of a fixed overriding plate. However, this model did not exhibit cyclic back-arc spreading. We were not able to unambiguously identify the main reasons for the differences between the present model and the models of Nakakuki and Mura (2013). We think it likely that the cause lies in the differences in material properties around the arc area. The 2D numerical model of Nakao et al. (2016) incorporates the processes of dehydration of the subducting slab and hydration of the mantle above the slab along with the associated changes in viscosity and density of rocks, and used a temperature-dependent viscosity and a penetrable 660-km boundary. This model shows different subduction behavior depending on the degree to which the viscosity and density are reduced due to hydration. For large reductions in viscosity and relatively small reductions in density related to hydration, the model exhibits back-arc spreading. Their model does not show cyclic back-arc spreading. This may in part be due to the short computation time of 20-40 Myr. Billen (2017) shows arc rifting developed in a 2D numerical model, and in this model arc rifting is preceded by the thickening of a low-viscosity boundary layer between the two plates and intrusion of mantle rocks into this low-viscosity domain. The development of a low-viscosity interplate layer may correspond to the formation of an accretionary prism in natural subduction zones. The validity of applying these model concepts to a natural subduction system requires further investigation.

Model limitations
Two significant limitations of the present two-dimensional model are that it does not take into account possible effects of toroidal mantle flow around the slab edge and it only considers the interaction of two plates ignoring possible effects of other adjacent plates (Schellart and Moresi 2013 and references therein). These limitations mean care is needed when comparing the model to natural subduction zones. For instance, the Tonga arc shows rapid trench retreat (Schellart et al. 2008) associated with toroidal mantle flow around the slab edge (Long and Silver 2008;Lynner, et al. 2017) and the Cenozoic subduction process around the Japanese islands involves multiple plates including the Pacific, the Philippine Sea and the Eurasian plates (Holt et al. 2018). In addition, the formation of the Shikoku and Parece Vela Basin and the Mariana Trough behind the Izu-Bonin-Mariana arc is accompanied by the clockwise rotation and westward subduction of the Philippine Sea plate (Holt et al. 2018;Wu et al. 2016). In some cases, collision of buoyant lithosphere may also affect back-arc spreading (Wallace et al. 2009).
However, since some common features are recognized between the present model and the natural subduction system as discussed below, we believe the present simple model expresses fundamental processes related to the phenomenon of back-arc spreading.

Geological implications 4.4.1 Conditions for cyclic back-arc spreading
Cyclic back-arc spreading in nature generally shows similar timescales of ~ 20 Myr. The South Fiji and the Lau Basins that developed behind the Tonga-Kermadec arc started spreading at ~ 27 Ma and ~ 7 Ma, respectively (van de Lagemaat et al. 2018). The Parece Vela Basin and the Mariana Trough that developed behind the Mariana arc started spreading at ~ 30 Ma and ~ 7 Ma, respectively (Wu et al. 2016). The Liguro-Provençal and the Tyrrhenian Basins that developed behind the Calabrian arc started spreading at ~ 30 Ma and ~ 12 Ma, respectively (Faccenna et al. 2001;Chen et al. 2015). The presence of a common timescale for cyclic back-arc spreading strongly suggests that the cyclicity is an intrinsic part of the backarc spreading process.
In the present model, the period of cyclicity is ~ 15 Myr (Fig. 3h). A simple geometrical consideration indicates that this period corresponds to the time required for a length of slab to subduct during a cycle of backarc spreading (Fig. 7). Assuming this length is ~ 900 km and the average subduction velocity (V sp + V tr ) is 6 cm/ yr, the period will be 15 Myr (Fig. 7a). In natural subduction zones, the effect of toroidal mantle flow around the slab edge will enhance trench rollback and thus lead to a longer distance of back-arc spreading and a longer period for the cyclicity than predicted in the present model, which does not consider toroidal flow (Fig. 7b).
In the present model, the cyclic back-arc spreading occurs for cases of intermediate resistance to bending of the slab. A large resistance to slab bending with respect to the strength of the arc lithosphere results in continuous back-arc spreading (cases 4 and 8), and a small resistance to bending results in subduction without back-arc spreading (case 3). In natural subduction zones, cyclic back-arc spreading occurs at subduction zones where old oceanic plate is subducting, such as the Mariana arc and the Tonga-Kermadec arc of the western Pacific. The condition of σ yield = 500 MPa may be too large for the present Earth, because there is no known example of continuous back-arc spreading lasting more than 30 Myr. This inference about slab strength is consistent with the conclusion of Kaneko et al. (2019), who estimated the slab strength from the ratio of sinking velocity of slab in the lower mantle to the velocity of subducting surface plate motion.

Possible changes in slab geometries during cyclic back-arc spreading
The spreading of the South Fiji Basin behind the Tonga-Kermadec arc continued until ~ 15 Ma, and this was followed by a period of quiescence until spreading of the Lau Basin started at around 5-7 Ma. This more recent spreading is ongoing and currently propagating southwards (van de Lagemaat et al. 2018). Active rifting is now focused in the Havre Trough, which is located in the southern part of the back-arc basin. The observed changes in spreading activity can be related to an increase in the slab dip angle from north to south (Fig. 8a). The relation between the slab dip angle and the kinematics of back-arc spreading can be accounted for if back-arc Fig. 7 Cartoons illustrating how one particular segment of the slab subducts during a cycle of back-arc spreading. The distance between the two red circles corresponds to a fixed length, and their locations highlight the movement of the slab. a In the present model, the length is ~ 900 km for the case shown in Fig. 3. b In natural subduction zones, the effect of toroidal mantle flow around the slab edge will enhance trench rollback and the subducted length may be greater than shown in Fig. 7a Page 15 of 19 Ishii and Wallis Progress in Earth and Planetary Science (2022) 9:27 spreading of the Lau Basin-Havre Trough started when the dip angle of the slab was nearly vertical, and the rapid trench retreat of the Tonga arc decreased the slab dip angle in this area. In contrast, the slab beneath the Kermadec arc remains nearly vertical. The southward propagation of the back-arc spreading in this area may be  Ishii and Wallis Progress in Earth and Planetary Science (2022) 9:27 related to the contemporaneous southward migration of the position of the collision with the Louisville Seamount Chain, which is now located at the junction between the Tonga and Kermadec arcs (Ruellan et al. 2003). Toroidal mantle flow around the northern edge of the slab may also play a role. Based on the present results, the flattened slab beneath the Tonga-Kermadec arc is considered to have formed after the cessation of the spreading of the South Fiji Basin and before the start of the spreading of the Lau Basin. While the flattened slab beneath the Tonga arc remains at the 660-km boundary during the spreading of the Lau Basin, the flattened slab beneath the Kermadec arc has penetrated the boundary due to subduction from a nearly immobile trench (Fig. 8c).
In the western Mediterranean, the slab below the Calabrian arc flattens above the 660-km boundary until it finally penetrates into the lower mantle below the western margin of the Liguro-Provençal Basin (Fig. 8b). The model of Schellart and Moresi (2013) implies that a pause in the cyclic spreading accompanies a folding of the slab above the 660-km boundary. However, seismic tomography does not show such folding (Fig. 8b). In light of our results, we suggest the onset of spreading of the Tyrrhenian Basin started at ~ 12 Ma, 4 million years after the cessation of the spreading of the associated Liguro-Provençal Basin when the 'heel' of the buckled slab reached the 660-km boundary (Fig. 8d).

Relative scarcity of the cyclic back-arc spreading
Although many subduction zones show back-arc spreading, cyclic spreading is more limited. Differences in tectonic setting may control whether or not back-arc spreading becomes cyclic. One example may be pinning of the subduction zone by the collision of the Louisville Seamount Chain at the boundary between the Havre Trough and the Lau Basin, as mentioned above.
The Izu-Bonin-Mariana arc shows an increase in slab dip to the south (Gudmundsson and Sambridge 1998;Li et al. 2008). A northward increase in trench retreat associated with the clockwise rotation of the Philippine Sea plate has been suggested as the primary cause of this observed change in the slab dip angle (van der Hilst and Seno 1993; Wu et al. 2016). Spreading of the Mariana Trough behind the Mariana arc started at ~ 7 Ma shortly after cessation of the spreading that formed the Parece Vela Basin at ~ 8 Ma, but no similar spreading started behind the Izu-Bonin arc (Tani et al. 2011;Wu et al. 2016). After the cessation of the clockwise rotation of the Philippine Sea plate, the slab started to fall vertically with respect to the trench position. Our modeling suggests back-arc spreading of the Izu-Bonin arc should start from the south and propagate northward when the 'heel' of the buckled slab reaches the 660-km boundary.
The NE Japan subduction zone shows a nearly constant low-angle slab dip, and no spreading has occurred since the cessation of spreading of the Japan Sea at ~ 15 Ma. Using an analytical model of global mantle flow, Holt and Royden (2020) showed that the global pattern of slab dips including the NE Japan subduction zone can be explained by a force balance between slab buoyancy and dynamic pressure associated with global mantle flow. The dynamic pressure may have prevented the low-angle slab from falling and back-arc spreading restarting. In this case, the back-arc spreading is also controlled by the global distribution of subduction zones and associated global mantle flow pattern.

Evolution of subduction zones
The present model assumes that the subduction starts from subduction beneath the overriding continental plate. After the start of back-arc spreading, the oceanic plate subducts beneath the newly formed back-arc oceanic lithosphere. Calculated surface heat flow around the back-arc basin (~80 mW/m 2 , Additional file 1: Fig. S9) is consistent with the observations (Currie and Hyndman 2006). However, the present model does not include the formation of new crust and crust is thin or absent in areas that have undergone spreading. The fore-arc domain preserves the original thickness of crust (Additional file 1: Fig. S9).
The processes that lead to subduction initiation are under debate (Lallemand and Arcay 2021), but the subduction of the Calabrian arc and the Tonga-Kermadec arc is thought to have started when an oceanic domain foundered below a neighboring continental domain van de Lagemaat et al. 2018 (Ishizuka et al. 2013;2022): (1) subduction associated with the formation of the present Izu-Bonin-Mariana arc started with subduction beneath Mesozoic remnant arc terranes, (2) the subduction initiation is almost contemporaneous with spreading of the West Philippine Basin within the remnant arc terranes, and (3) the spreading of the West Philippine Basin was associated with magmatism caused by a mantle plume.
Page 17 of 19 Ishii and Wallis Progress in Earth and Planetary Science (2022) 9:27 Although the present model does not incorporate these specific processes during the subduction initiation, the subduction beneath the newly formed back-arc oceanic lithosphere (the West Philippine and the Shikoku-Parece Vela Basins and the Mariana Trough) during the cyclic back-arc spreading is consistent with the present model.

Conclusions
(1) Cyclic back-arc spreading is reproduced by using a two-dimensional thermo-mechanical model in which back-arc spreading starts when the nearly vertical slab impinges upon the 660-km boundary. The decrease in slab dip angle due to the trench retreat during the back-arc spreading causes buckling of the slab and the back-arc spreading stops. When the 'heel' of a buckled slab with a nearly vertical dip angle impinges upon the 660-km boundary, rifting restarts within the arc before shifting to the back-arc domain and causing a renewed phase of back-arc spreading. The relationship between slab buckling and back-arc spreading indicates that the extensional stress within the arc caused by the downdip compressive stress transmitted up the slab is a necessary condition for the start of each phase of the back-arc spreading.
(2) In the present model, cyclic back-arc spreading occurs for cases of intermediate resistance to bending of the slab. A large resistance to slab bending results in a continuous back-arc spreading, and a small resistance to slab bending results in subduction without back-arc spreading. These results indicate that in addition to the negative buoyancy of the subducting slab the resistance to bending of the slab is also a major factor that controls back-arc spreading.
(3) Cyclic back-arc spreading in nature shows a common timescale of ~ 20 Myr, suggesting that the cyclicity is an intrinsic part of the back-arc spreading mechanism. The present model shows a slightly shorter period of cyclicity (~ 15 Myr). This may result from two-dimensionality of the present model that ignores the effects of toroidal mantle flow around slab edges. (4) The relation between the processes of cyclic backarc spreading and the slab geometries in the Tonga-Kermadec and the Calabrian arcs can be explained on the basis of the present model.
Additional file 1: Figure S1. Effects of changes in the values of the Clapeyron slope. a The change in velocities of the subducting plate (V sp ), the overriding plate (V op ) and the trench (V tr ) combined with the average dip angle of the slab surface (Dip) measured in the depth range of 200-300 km. b Slab geometries. The distribution of effective viscosity is shown by different colors. Blue lines indicate the 410 km and 660 km phase boundaries. Figure S2. Effects of reducing viscosity (a, c) and density (b, d) of the mantle wedge. a Result for a case of A N = 1 × 10 -8 /s/MPa for the mantle wedge, otherwise the same conditions as the reference case. b Result for a case where ∆ρ 0wedge = 0 kg/m 3 for the mantle wedge, otherwise the same conditions as the reference case. c Result for a case where A N =6 × 10 -8 /s/MPa for the mantle wedge, otherwise the same conditions as the reference case. d Result for a case where ∆ρ 0wedge = 40 kg/m 3 for the mantle wedge, otherwise the same conditions as the reference case. Upper panel shows changes in velocities of the subducting plate (V sp ), the overriding plate (V op ) and the trench (V tr ). The average dip angle of the slab surface (Dip) measured in the depth range of 200-300 km is also shown. Lower panel shows the slab geometry. The distribution of effective viscosity is shown by different colors and arrows indicate the velocity field. Figure  S3. The effect of changing resolution. a High resolution. b Low resolution. Upper panel shows the change in velocities of the subducting plate (V sp ), the overriding plate (V op ) and the trench (V tr ), in addition to the average dip angle of slab surface (Dip) measured in the depth range of 200-300 km. Lower panels show the distribution of effective viscosity indicated by different colors and the velocity field shown by arrows. Figure S4. Results including pressure-independent density. a Initial distribution of temperature, viscosity and density. b The change in velocities of the subducting plate (V sp ), the overriding plate (V op ) and the trench (V tr ). The average dip angle of slab surface (Dip) measured in the depth range of 200-300 km is also shown. c The distribution of effective viscosity is indicated by different colors and the velocity field is shown by arrows. Figure S5. Result for a case without viscous dissipation until the slab tip reaches a depth of 250 km, otherwise the same conditions as the reference case. a The change in velocities of the subducting plate (V sp ), the overriding plate (V op ) and the trench (V tr ). The average dip angle of slab surface (Dip) measured in the depth range of 200-300 km is also shown. b The distribution of effective viscosity is indicated by different colors and the velocity field is shown by arrows. Figure S6. Result for the case without mantle wedge, otherwise the same conditions as the reference case. a The change in velocities of the subducting plate (V sp ), the overriding plate (V op ) and the trench (V tr ). The average dip angle of slab surface (Dip) measured in the depth range of 200-300 km is also shown. b The distribution of effective viscosity is shown by different colors. Right panels show enlargements of the arc area. Arrows indicate the velocity field. Figure S7. a The slab geometry for the case without a viscosity jump at the 660-km boundary, otherwise the same conditions as the reference case. The distribution of effective viscosity is shown by different colors. b Viscosity distribution. Figure S8. Result for the case without a high effective thermal conductivity for back-arc lithosphere, otherwise the same conditions as the reference case. a Time evolution of slab geometry and velocity field (arrows). The distribution of effective viscosity is shown by different colors. b The changes in velocities of the subducting plate (V sp ), the overriding plate (V op ) and trench (V tr ). Arrows indicate the time for Fig. S8a. Figure S9. Surface heat flow (upper panel) and the distribution of the continental crust (lower panel) at t = 33.6 Myr for the reference case (Fig. 3f ).