2-D Numerical Simulations on Formation and Descent of Stagnant Slabs: Important Roles of Trench Migration and Its Temporal Change

We conducted numerical simulations of thermal convection of highly viscous fluids in a 2-D spherical annulus in order to study what mechanisms control the dynamic behaviors of subducting slabs such as the formation of “stagnant slabs” in the mantle transition zone (MTZ) and their descent into the lower mantle. Two series of experiments are carried out, by applying the different histories of migrating motion of “trench” where the slab of cold fluids descends from the top surface, together with systematically varying the velocities of subducting slabs and trench migration, the Clapeyron slope at around 660 km depth, and the viscosity jump between the upper and lower mantle. In the first series of experiments where the migration rate of trench is kept constant with time, our model successfully reproduces the diverse morphology of subducting slabs depending on the delicate combinations of control parameters, including five fundamental types named penetration, accumulation, entrainment, long-term stagnation and short-term stagnation (SS). In addition to the above fundamental types, we found two distinctive types of slab behaviors (named pS and cS) where the slabs are stagnated at the base of MTZ after they experience the penetration or collapse into the lower mantle and yield the snapshot behaviors very close to those of SS at some time instances. In the second series of experiments where we imposed a step-like change in the trench retreat velocity with time, we found that a deceleration of trench retreat enhances the collapse of the horizontally-lying or stagnant slabs into the lower mantle by, for example, inducing the avalanche of stagnant slabs from its hinge. In particular, we obtained in some cases the flattening of slabs well below the 660 km discontinuity without any discontinuous changes in mantle properties at the depth. Our results suggest that the formation and descent of stagnant slabs are strongly related to the trench retreat, particularly through its temporal changes. In other words, the variations in the shapes of subducting slabs in nature are most likely to reflect the difference in the history of trench migration.


INTRODUCTION
Seismic tomography studies (e.g., Fukao and Obayashi, 2013) have revealed that the slabs of subducted plates show diverse behaviors in the mantle beneath various subduction zones. In several regions such as Japan and Izu-Bonin arc, for example, so-called "stagnant slabs" are observed at the base of the mantle transition zone (MTZ), without penetrating into the lower mantle owing to the horizontal flattening at around the 660 km discontinuity. In other subduction zones where (parts of) the slabs have entered into the lower mantle, not all slabs can reach the deep lower mantle readily and smoothly as observed in, for example, beneath Central America arc: It is now recognized that in some regions subducted slabs seem to be stagnated at about 1,000 km depths (e.g., Java arc), experiencing a stagnation at the MTZ and penetration into the lower mantle simultaneously (e.g., Tonga arc), or stagnating at the MTZ before collapsing from its hinge (e.g., Mariana arc).
Since the numerical studies of mantle convection first arose (e.g., McKenzie et al., 1973), the dynamic behaviors of subducted plates in the mantle have been intensively studied, particularly focusing on the mechanisms of the formation of stagnant slabs in the MTZ. Christensen (1996) developed a 2-D numerical model in Cartesian geometry which incorporates the effects of an endothermic phase boundary and viscosity jump at the 660 km discontinuity together with imposed motions of subducting plate and trench at the surface, and reported that the trench retreat plays an important role in flattening of subducting slabs at the base of MTZ. Torii and Yoshioka (2007) systematically studied the influence of these physical parameters on the formation of stagnant slabs, and showed that stagnant slabs can be generated in realistic parameter ranges. These numerical models are further improved by many workers through incorporating other complicated factors such as the effects of the age of the subducting and overriding plate (Garel et al., 2014), the viscosity reduction in slabs due to the grain-size reduction (Čížková et al., 2002;Nakakuki et al., 2010), the presence of metastable phase transitions in the cold slabs (Tetzlaff and Schmeling, 2000;Agrusta et al., 2014).
In addition to the dynamics of the slabs at the base of MTZ, a growing number of numerical studies have been conducted on the avalanche of the subducting slabs into the lower mantle. Yoshioka and Naganoda (2010) carried out numerical experiments in a 2-D box model with imposed surface velocity,  (1) with r η = 1, and blue dashed lines indicate those with r η = 1.8, 3, 5.6, 10, 18, and 30, respectively. and pointed out that the trench advance and trench jumps cause the eventual fall of stagnant slabs into the lower mantle as proposed by van der Hilst and Seno (1993). Agrusta et al. (2017) used a 2-D self-consistent subduction models to investigate the transition between the stagnation and penetration of slabs, and have revealed that these transitions are strongly influenced by many agents such as the temporal changes in trench retreat velocity, the Clapeyron slope at the 660 km discontinuity, the subducting plate age, and so on. On the other hand, recent numerical models of the subduction in the Mariana arc (e.g., Čížková and Bina, 2015;Faccenna et al., 2018) suggested the double subduction system as a possible mechanism that triggers the change from trench retreat to trench advance in this region.
In this study, we will further develop numerical model of subducted slabs in the convecting mantle using 2-D model in spherical annulus, in order to study what mechanisms control the behaviors of subducted slabs in the mantle at various depths, such as the formation of stagnant slabs in the MTZ and the course of its avalanche into the lower mantle. Compared with earlier works, we will perform numerical calculations by varying the values of control parameters systematically over broad ranges. As our initial attempt to study the effects of various parameters, we will follow the strategy by Torii and Yoshioka (2007) and focus on the effects of trench migration rate and its temporal change. By carefully comparing the shapes of slabs obtained in our simulation models with those estimated from seismic tomography, we aim at deepening the understanding on the dynamics and temporal evolutions of subducting slabs in the mantle.

Model Setup
Our numerical model is schematically shown in Figure 1A. We consider thermal convection of mantle materials which is driven by a subducting and retreating motion of a cold slab in one eighth (azimuth θ min ≤ θ ≤ θ max where θ max − θ min = 45 • ) of a twodimensional spherical annulus whose outer and inner radii are r max = 6, 400 km and r min = 4, 400 km, respectively. We take into account the effects of the phase transition from olivine to wadsleyite at around 410 km depth and that from ringwoodite to bridgmanite and ferro-periclase at around 660 km depth. The density contrasts associated with these phase transitions are taken to be 8.3 and 7.8%, respectively (Dziewonski and Anderson, 1981). The Clapeyron slope of the phase transition at 410 km depth is taken to be γ 410 = +3 MPa/K, while that of 660 km depth γ 660 is varied from −3 MPa/K or −1 MPa/K, based on the result of laboratory experiments (e.g., Ito and Takahashi, 1989;Katsura et al., 2003;Fei et al., 2004). The meanings and values of parameters used in this study are given in Table 1.
The viscosity η of mantle material is assumed to depend exponentially on the temperature T and pressure p (or the depth from the top surface), and given by (1) where T surf = 0 • C is the temperature at the top surface, r (r min ≤ r ≤ r max ) is the distance from the center of the Earth, η surf is the viscosity at the top surface (r = r max and T = T surf ), T is the temperature scale which is arbitrarily chosen to be 2, 000 • C, and E T = ln(10 6 ) and E p = ln(10 4 ) express temperature and depth dependence of the viscosity. In addition, the effect of the viscosity jump between the upper and lower mantle is introduced by f (Ŵ 660 ) through the "phase function" Ŵ 660 of the endothermic phase transition at around 660 km depth (Christensen and Yuen, 1985); f (Ŵ 660 ) is taken to be 1 and r η in the upper and lower mantle, respectively. Here, we assume the value of r η to be 1 ≤ r η ≤ 30 (see Figure 1B) as had been employed in earlier studies (e.g., Torii and Yoshioka, 2007;Yoshioka and Naganoda, 2010;Kameyama and Nishioka, 2012). In this study, the kinematic motions of subducting and overriding plates are imposed near the top surface, as in our earlier work (Kameyama and Nishioka, 2012). First, as in the works by Yoshioka and coworkers (Yoshioka and Sanshadokoro, 2002;Torii and Yoshioka, 2007;Yoshioka and Naganoda, 2010), the subducting plate is modeled by a downward flow of cold and viscous fluid at a uniform velocity v pl along with a channel of 80 km thick which guides the cold fluid from the top surface into the MTZ. The channel is bent downward with the dip angle 45 • at the trench which is located θ tr = 5.625 • away from the left-side of the top surface (θ = θ max ). In addition, the length of the channel is assumed to increase with time, in order to model the descent of cold slab. The maximum depth of the upper plane of the 80 km thick channel is taken to be 400 km. Second, we take into account the effect of trench retreat by the similar method as van Hunen et al. (2000). Here the overriding plate is assumed to migrate together with the trench and channel at a rate of v tr in the horizontal direction with respect to the deep mantle. That is, the effect of trench migration is numerically implemented through the coordinate transformation between the In the present numerical model, the rate of plate convergence v c at the trench is given by the sum of the subducting plate velocity (v pl ) and the overriding plate velocity (v tr ).
When non-dimensionalized using a velocity scale of v scale ≡ κ/(r max − r min ) (where κ = 10 −6 m 2 /s is thermal diffusivity), the dimensionless rate of plate convergence v ′ c = v c /v scale is kept to be 5, 000 throughout this study, while that of trench migration v ′ tr = v tr /v scale are changed from 0 to 3, 000. That is, v c = 7.88 cm/yr and v tr ≤ 4.73 cm/yr are employed. The initial distribution of temperature T is given by where T pot = 1280 • C is the potential temperature of the mantle, α is thermal expansivity, and C p is specific heat. In the righthand side, the first and second terms represent the effects of half-space cooling (Turcotte and Schubert, 2002) and adiabatic compression ( Figure 1B), respectively. In this equation, τ (θ ) is the age of plate as a function of the horizontal position θ . At the trench (θ = θ max − θ tr ), τ (θ ) is taken to be τ pl (=101Ma). On the other hand, τ (θ ) decreases with θ in the oceanic region (θ ≥ θ max − θ tr ), while τ (θ ) = τ pl is fixed in the continental region (θ ≤ θ max − θ tr ). The boundary conditions for temperature T are (i) fixed temperature T = T surf = 0 • C at the top surface (r = r max ), (ii) fixed temperature at the left boundary (θ = θ max ) in order to maintain the initial thermal structure of subducting oceanic plate given by Equation (3), and (iii) adiabatic conditions at the bottom (r = r min ) and right (θ = θ min ) boundaries. On the other hand, the boundary conditions for velocity v are given by (i) v r = 0 and fixed v θ at the top surface, (ii) v r = 0 and ∂v θ ∂θ = 0 at the both side boundaries (θ = θ min , θ max ) and (iii) v θ = 0 and ∂v r ∂r = 0 at the bottom boundary, where v θ and v r are θ -and r-components of the velocity, respectively. That is, no slip condition is assumed at the top surface, while at the other boundaries the flow is taken to be perpendicular to the boundaries and zero normal stress. In addition, the condition is imposed at the right-hand side wall so as to keep the net mass flux zero across the boundary.

Numerical Techniques
We calculate the flow field, the distributions of temperature T and their temporal variations by solving the governing equations for conservations of mass, momentum and the thermal energy under the extended Boussinesq approximation (Christensen and Yuen, 1985). The conservation of the thermal energy includes the effects of adiabatic heating and latent heat due to the phase transitions. However, the effect of viscous dissipation (or frictional heating) is ignored in this study, in order to avoid a spuriously high dissipation near the trench where the channel is sharply bent in an artificial manner. The equations for conservation of mass and momentum are written in terms of the stream function ψ defined by instead of directly solving for the velocity v and pressure p. On the other hand, the time integration of energy equation is carried out by the Crank-Nicolson scheme. The discretization of the basic equations are carried out based on the finite volume method. The region is uniformly divided into 256 and 640 meshes in r-and θ -directions, respectively. This pattern of mesh division yields almost square numerical meshes near the top surface. The discretized equations for temperature T and flow fields (velocity v and pressure p) are solved by our thermal convection simulation code for 2-D spherical annulus (Kameyama and Harada, 2017). The numerical validity of our code has been already verified by comparing with previous study (Zhong and Gurnis, 1993). The time increment δt is taken to be the half of the value determined from the numerical stability of energy equation, yielding δt ∼ 1.6 × 10 4 years. We confirmed that δt defined thus is small enough, by comparing the numerical results for several cases with the ones done using smaller time increments by factors of 0.5 and 0.25. The time integrations are carried out for more than 100 million years (My) for all the cases presented here, and for up to 150 My in some cases where the slab stagnation is observed (see below for details).

RESULTS
Two series of experiments are presented in this study, using different histories of the migrating motion of trench. In the first series of experiments, we carried out calculations where the rate v tr of trench retreat is assumed to be constant with time. In the second series of experiments, on the other hand, we take into account the temporal changes in v tr , by imposing its step-like change with time.

Case
Type Case Type Case Type

Slab Behaviors Under Constant Rate of Trench Retreat
In the first series of experiments, we investigate the subduction dynamics obtained from a series of 84 numerical experiments where we systematically varied the values of viscosity jump r η between the upper and lower mantle, Clapeyron slope γ 660 at 660 km discontinuity, and the rate v tr of trench retreat (see Table 2). The velocity v pl of subducting plate is accordingly varied in order to fix the rate of plate convergence v c at the trench. We here note that, rather than v pl , v c controls the amount of the slab of cold fluid which is "injected" into the mantle during the period of calculation. In other words, for all the cases with fixed v c presented here, the effects of negative buoyancy of cold slabs are equal at the same elapsed times from the onset of subduction.

The Effects of Control Parameters on Slab Dynamics
We first discuss the effects of the control parameters in this study (r η , γ 660 , and v tr ) on the behaviors and shapes of subducting slab around the MTZ and in the lower mantle. In Figure 2 we show the snapshots of the distributions of temperature T and viscosity η at the elapsed time of about 50 million years (My) after the onset of subduction, obtained for four selected cases from the first series of calculations. Shown in Figure 2A are the snapshots of the case where we employed r η = 1, γ 660 = −1 MPa/K, and v tr = 0 cm/yr. As can be clearly seen from the figure, a slab of cold and viscous fluid is subducting through the channel down to the top of MTZ (see the phase transition at around 410 km depth indicated by the upper red lines). The slab further sinks across the phase transition at around 660 km depth (indicated by the lower red lines), and finally penetrates into the lower mantle directly and smoothly. In order to see the influence of the trench retreat, we show in Figure 2B the snapshot of the case where a non-zero rate of trench retreat is given (v tr = 3.15 cm/yr) while the values of other parameters are the same as those in Figure 2A. By comparing Figure 2B and Figure 2A we can clearly see that the trench retreat enhances the flattening and stagnation of subducted slabs. Indeed, a horizontally-lying portion of cold slab can be observed in the MTZ in Figure 2B, which is in a stark contrast to the case in Figure 2A. However, as can be seen from Figure 2B, the entire portion of subducted slab is not stagnated around the 660 km depth; the tip of cold slab penetrates and sinks in the lower mantle. This indicates that the combined effect of v tr = 3.15 cm/yr and γ 660 = −1 MPa/K (and r η =1) employed in this case is not strong enough for the entire cold slab to be stagnated.
We next present the snapshots of Figures 2C,D in order to see the influences of γ 660 and r η under a non-zero v tr . Shown in Figures 2C,D are the cases where a steeper γ 660 = −3 MPa/K and a larger jump r η = 10 in viscosity is employed at the 660 km discontinuity than those in Figure 2B, respectively, while the values of other parameters are kept unchanged. We can clearly see that, from the comparison with Figure 2B, the entire portion of subducted slab stagnates at the 660 km discontinuity, either because of a positive buoyancy in the cold slab induced by the phase transition with strongly negative Clapeyron slope γ 660 (in Figure 2C) or because of a strong viscous resistance coming from highly viscous lower mantle with large r η (in Figure 2D). These results are consistent with the earlier findings by Torii and Yoshioka (2007) that classified the manner of slab stagnation into buoyancy-dominated and viscosity-dominated regimes.
Our results with 2-D spherical annulus showed that all of the control parameters in this study (the viscosity jump r η and Clapeyron slope γ 660 at the 660 km discontinuity, and the rate v tr of trench retreat) enhances the slab stagnation at the base of MTZ. Among them, r η and γ 660 tend to prevent subducting slabs from penetrating into the lower mantle, whereas v tr flattens the subducted slabs at around the MTZ. These results are consistent with those of earlier numerical studies in the 2-D box model (e.g, Torii and Yoshioka, 2007;Kameyama and Nishioka, 2012). The close similarity between the results of our experiments and earlier ones also suggests that the effect of the Earth's curvature is minor on the dynamic behavior of slabs in the shallow mantle.

Five Fundamental Types of Slab Behaviors
Our systematic calculations of 84 cases showed that there are five fundamental types for the behaviors of subducted slabs at around the MTZ, depending on the combinations of three control parameters (r η , γ 660 , and v tr ). Figure 3 shows the temporal development of the five types of subducting slab behaviors. In Table 2 and Figure 4, on the other hand, we summarize the classification of subducting slab behaviors for all 84 cases.
Among the five types of slab behaviors, first three are the ones where subducted slabs are not stagnated at the base of MTZ owing to extreme choices of control parameters other than γ 660 . Figure 3A shows the "penetration" (P)-type behavior of slabs. The P-type is characterized by the penetration of cold slabs into the lower mantle immediately and vertically. This type of behavior is observed for 15 cases out of 84 calculations, and occurs for v tr = 0 and small r η (r η ≤ 10). On the other hand, Figure 3B shows the "accumulation" (A)-type behavior of slabs. In the A-type, subducted slabs do not penetrate into the lower mantle but accumulate at around the 660 km discontinuity,  , red open circles, and orange filled circles, respectively. In addition, the types of stagnation after penetration (pS) and stagnation after collapse (cS) are represented by small letters of "pS" and "cS" in the right side of filled orange circles, respectively. In the schematic illustration of the regimes, the black triangles indicate the initial trench location, and the dashed lines indicate the 660 km discontinuity.
owing to the large viscosity jump between the upper and lower mantle. This type of behavior is observed for 6 cases out of 84 calculations, and occurs for v tr = 0 and large r η (r η ≥ 18). In Figure 3C we show the third type of slab behaviors, called "entrainment" (E)-type, observed for the cases when both v tr and r η are large. In the E-type, subducting slabs are entrained and pulled up by a local circulations in the shallow mantle driven by a trench retreat and is broken into pieces, since the effects of the imposed v tr and r η are too strong. This type is observed in 18 cases of our calculations.
The fourth and fifth types of slab behaviors are the ones where the stagnant slabs are observed in the presence of trench retreat (v tr = 0) in addition to the effects of γ 660 and r η . Shown in Figure 3D is an example of "long-term stagnation" (LS)-type, while in Figure 3E is that of "short-term stagnation" (SS)-type. For both types, the subducting slab is bent and flattened at the 660 km discontinuity, forming a horizontally-lying stagnant slab at the base of MTZ in the earliest stage of subduction. In the LS-type, subducted slabs are stagnated around the 660 km discontinuity for a long time interval (at least more than 80 My) without penetrating into the lower mantle. The horizontallylying portion of cold slab becomes longer with time along with the trench retreat. In the SS-type, in contrast, the stagnant slab starts to collapse into the lower mantle at its tip after the first several ten million years, because the negative buoyancy caused by the mass of the stagnant slab dominates the positive buoyancy owing to the negative Clapeyron slope γ 660 and the resistance due to the viscosity jump between the upper and lower mantle. The LS-and SS-type of behaviors are observed for 20 and 6 cases out of 84 calculations, respectively.
We however note that not all the patterns of slab behaviors obtained here fall in any of five fundamental types: In several cases we observed distinctive temporal evolutions whose snapshots are very similar to those of SS-type at some time instances. In the next subsubsection, we will study the slab behaviors of such cases in more detail.

Distinctive Slab Behaviors: Stagnation After Penetration or Collapse
In the cases of short-term stagnation (SS) described in section 3.1.2, the cold slabs sink into the lower mantle after they lie horizontally at the base of MTZ. In some other cases, in contrast, we observed that a horizontally-lying (or "stagnated") portion of cold slab develops after a part of slab sinks into the lower mantle. By the letters "pS" or "cS" attached with the filled orange circles in Figure 4, we indicate the values of parameters which yield a "stagnation" of slab after its penetration (hereafter called "pS") or collapse ("cS") into the lower mantle. In Figure 5, on the other hand, we show the temporal evolution of the distributions of temperature T and horizontal stress σ θ θ around the subducting slabs for the cases with pS and cS, together with that of SS for comparison. As can be seen from the snapshots at around t = 10 My, the behaviors of subducting slabs are very similar for all the three types before the tips of subducting slabs impinge the 660 km discontinuity. In addition, the behaviors of slabs are quite similar for sufficiently large t (e.g., t > 60 My) between the types. However, during the period of temporal evolution except for very early and very late stages, the behaviors of slabs significantly differ between the types. Figure 5A shows the temporal evolution of the slabs for the case with stagnation after penetration (pS) where we assumed γ 660 = −1 MPa/K, r η = 1 and v tr = 3.15 cm/yr. From the snapshot at t ≃ 15 My, we can see the region with σ θ θ > 0 (i.e., vertical compression) in the lower part of the slab just above the 660 km discontinuity both for the cases with pS and SS. However, the degree of vertical compression in this region is weaker than that for the case with SS (see Figure 5C), because of the weaker resistance coming from the gentler Clapeyron slope γ 660 for given v pl and v tr . This results in a slab penetration to a greater depth by about 20 My than that for SS. At t = 19.98 My, the stress state near the tip of the subducting slab in the lower mantle is strong horizontal compression (vertical extension) in the upper part, indicating that the cold slab tends to sink owing to its negative buoyancy. During the subsequent period of t < 39.95 My, however, the horizontal compression becomes significantly weakened in the tip of slab in the lower mantle because the temperature difference becomes smaller between the slab and the surrounding mantle, leading to a slowing down of the descent of slab in the lower mantle. The deceleration of slab can be also seen from the stress state at the position of slab penetration at 660 km discontinuity; the horizontal compression gradually decreases with time there. On the other hand, the position of slab is almost unchanged at the base of the MTZ for t > 29.96 My, because the slab penetration also acts as a strong anchor of cold slab in the uppermost part of the lower mantle. Taken together with the effect of trench retreat, a subducting slab is flattened to form a horizontally-lying portion in the MTZ whose horizontal dimension increases with time. As can be seen from the snapshot of t > 40 My, the slab subduction for the pS-type finally yields a structure quite similar with those of SS for t > 40 My.
On the other hand, Figure 5B shows the temporal evolution of the slabs for the case with stagnation after collapse (cS) where we assumed γ 660 = −3 MPa/K, r η = 1 and v tr = 1.58 cm/yr. Comparing the snapshots for t less than about 20 My in Figures 5B,C, the behavior of slab for the case with cS is quite similar with that for SS. Indeed, for both cases, horizontally-lying portions of cold slabs are about to be formed at the base of MTZ. In addition, there are regions with vertical compression (σ θ θ > 0) in the lower part of the slab just above the 660 km discontinuity. However, the degree of vertical compression in this region is stronger for cS than that for SS (see Figure 5C), because of the slower trench retreat and faster imposed subduction together with the steeper Clapeyron slope γ 660 . During the subsequent period of t ≤ 30.01 My, the cold slab sinks into the lower mantle to a greater depth than for the case in Figure 5C. On the other hand, the imposed subduction also induces a horizontal compression near the tip of cold slab, which leads a buckling of slab and forms a concave structure of collapsed slab in the uppermost lower mantle at t = 30.01 My. During the subsequent period of t ≥ 39.93 My, the collapsed slab further sinks in the lower mantle owing to its negative buoyancy. In addition, as in the cases of pS, the descent of slab in the lower mantle strongly anchors the cold slab in the uppermost part of the lower mantle, and a horizontally-lying portion develops at the base of MTZ along with the trench retreat. As can be seen from the snapshot of t = 60.01 My, the slab subduction for the cS-type finally yields a structure quite similar with those of SS.

Stagnant Slab Behaviors With Temporal Changes in Trench Retreat Velocity
In the first series of experiments where the rate of trench retreat v tr is kept unchanged with time, we found that the trench retreat greatly helps to form and sustain the horizontally-lying or stagnant slabs (see section 3.1). In this section, we carry out the second series of experiments where the effect of trench retreat weakens with time, and study how the behaviors of subducting slabs are affected by the temporal changes in v tr particularly for the cases where horizontally-lying or stagnant slabs are formed under the continuous retreating motion of trench. In this study, for simplicity, we assumed that the trench retreat occurs at a constant non-zero rate for the elapsed time 0 ≤ t ≤ t stop while it does not for t > t stop . Here we choose t stop to be around 40 My, since the horizontally-lying slabs are well-developed or ready to be developed by t = t stop for all the cases with LS, SS, cS, and pS.
We conducted 45 numerical simulations by imposing sudden stop of the trench retreat at around 40 My after the onset of slab subduction, using the values of control parameters (r η , γ 660 , and non-zero v tr for t ≤ t stop ) which yield horizontally-lying slabs (LS, SS, cS, and pS in section 3.1). Among them, we observed the collapse of the stagnant slabs into the lower mantle only for the cases with r η ≤ 10. In particular, the avalanche styles of stagnant slabs observed in our simulations can be classified into four styles, depending not only on the types of slab behaviors described in section 3.1 but also on the slab shapes at the time instance of the sudden stop of the trench retreat (t = t stop ). In Table 3 and Figure 6, we summarize the classification of slab behaviors for all   Figure 6, we also plot 6 cases indicated by "did not fall" in Table 3, where the slabs did not fall into the lower mantle.) In Table 3, the symbol "type(40)" is meant for the type of the slab behavior for the cases with t stop = 40 My, while "type(∞)" stands for the type presented in section 3.1 with the continuous trench retreat (i.e., t stop = ∞). Figure 7A shows the slab behavior obtained for the case where we employed r η = 3, γ 660 = −2 MPa/K and v tr = 3.15 cm/yr is applied only for t ≤ 40 My after the onset of subduction (For comparison, we also show in small panels attached with large ones the snapshots for the same time instances of t > 40 My for the corresponding case without a sudden stop of trench retreat). The choice of control parameters yields the LS-type without the sudden stop of trench retreat; the stagnant slab is formed at the base of MTZ by t = 39.96 My after the onset of subduction, and it hardly falls into the lower mantle for more than 80 My. In the case where we impose a sudden stop of trench retreat at 40 My after the onset of subduction, in contrast, we observe an avalanche of stagnant slab. As can be seen from the snapshot at t = 60.08 My, a cold plate is bent downward at the hinge (the left-hand end) of the stagnant slab along with the imposed downward motion of subducting plate, and starts to penetrate into the lower mantle. At t = 80.00 My, a large portion of subducted slab collapses into the lower mantle, owing to its negative buoyancy. The type I of slab behavior is observed for 15 cases where type(∞) is LS. We also obtained the same type for 5 cases where type(∞) is SS. This reflects the fact that the shapes of slabs for SS are quite similar with those for LS at t = t stop before the collapse of slabs for SS. This result is also consistent with those of Yoshioka and Naganoda (2010) and Agrusta et al. (2017) where the collapse of stagnant slab enhanced by the change from trench retreat to trench advance. Figure 7B shows the slab behavior for the case where the type(∞) is pS with r η = 3, γ 660 = −2 MPa/K and v tr = 1.58 cm/yr. At t = 39.93 My just before the sudden stop of trench retreat, the subducting slab penetrates into the lower mantle almost smoothly, in spite of a small kink at the 660 km discontinuity. The descent of cold slab continues after the sudden stop of trench retreat along with the imposed downward motion of subducting plate. During the subsequent period of t ≤ 79.98 My, the kink of subducting slabs gradually sinks down to about 1,500 km depth. The tip of subducting slab, on the other hand, reaches down to about 1,500 km depth by t = 60.03 My, and is almost anchored at this depth during t ≤ 79.98 My. The difference in the rate of descent in the slab results in a horizontally-lying portion of cold slab in the lower mantle at around 1,500 km depth at this time instance, as if it is stagnated at this depth.

Type II: Stagnation Below 660 km
This type of slab behaviors can be seen in total of 8 cases where type(∞) is pS while a horizontally-lying portion of cold slab is not well-developed yet by t = t stop . We also note that in two cases with r η = 10 a thickening of cold slabs, similar to the accumulation in section 3.1.2, is observed at the 660 km discontinuity in addition to the trapping of slabs in the lower mantle.  (40); triangles, squares, diamonds, and inverted triangles indicate Types I, II, III, and IV, respectively. The cases with "did not fall" are also indicated by the letter "N" in the figure. On the other hand, the colors of symbols stand for the slab morphologies of Type(∞); red open symbols indicate LS, while yellow filled ones are other types whose Type(∞) are indicated by small letters ("SS", "pS," and "cS") attached to the symbols. In the schematic illustration of the regimes, slab morphologies before (black contour) and after (gray fill) changing v tr for Type I, II, III, and IV are shown. The black triangles and the dashed lines indicate the initial trench location and the 660 km discontinuity, respectively.

Type III: Buckling After Penetration in the Lower Mantle
In Figure 7C we show the slab behaviors for the case where the type(∞) is cS with r η = 1, γ 660 = −3 MPa/K and v tr = 1.58 cm/yr (see also Figure 5B in section 3.1). During the period of t ≤ 39.93 My where the trench retreat continuously occurs, the subducting slab collapses into the lower mantle at its hinge owing to its negative buoyancy. In addition, during the subsequent period of t < 59.98 My where the trench retreat is turned off, the behavior of subducting slab is quite similar with those for the case with cS where the trench retreat continuously occurs (see also Figure 5B). However, because of the absence of trench retreat and enhanced rate of imposed subduction, the subducted slab is intensively buckled in the uppermost part of the lower mantle.
This type of slab behavior is observed in a total of 9 cases where type(∞) is pS or cS. This is because in these cases the horizontally-lying portion of cold slab is not significantly developed by the time of the sudden stop of trench retreat (t = t stop ). On the other hand, for the cases with 3 ≤ r η ≤ 5.6 and γ 660 = −3 MPa/K, we observed that the slabs accumulate at around the MTZ, owing to the resistance force from large r η and steep γ 660 . Figure 7D shows the slab behaviors in the case with type(∞) is pS where we employed r η = 1, γ 660 = −1 MPa/K and v tr = 4.73 cm/yr. At t = 40.01 My just after the sudden stop of trench retreat, a well-developed SS-type slab takes place, which is characterized by a very long horizontally-lying portion and its tip sunk in the lower mantle. However, as can be seen from the snapshot of t = 59.92 My after the sudden stop of trench retreat, the enhanced motion of imposed subduction newly induces a penetration of cold slab at the hinge of stagnant slab in addition to that at its tip. During the subsequent period, the cold subducting plate is preferentially supplied to the newly penetrating portion, while in the old penetrating portion the temperature difference becomes smaller between the slab and the surrounding mantle.

Type IV: Avalanche From the Hinge After Penetration From the Head of a Slab
This type can be observed in 2 cases out of 45 calculations for the cases where type(∞) is pS or SS, particularly when a horizontally-lying stagnant portion of cold slab is sufficiently long at around the 660 km discontinuity by the time of the sudden stop of trench retreat.

DISCUSSION AND CONCLUDING REMARKS
In this paper we conducted numerical experiments of mantle convection with imposed kinematic motion of cold slabs and trench migration, aiming to understand the cause of the diverse morphology of subducting slabs observed by the seismic tomography. In particular, we studied the effects of control parameters (the viscosity jump r η and the Clapeyron slope γ 660 at the 660 km discontinuity, and the rate v tr of trench retreat) on the formation of "stagnant slabs" in the MTZ and their avalanche into the lower mantle. In the first series of our experiments where v tr is kept constant with time, we confirmed that the trench retreat greatly helps to form horizontally-lying or stagnant slabs at the 660 km discontinuity. In particular, we found distinctive temporal evolutions of slabs in some ranges of control parameters, in addition to the fundamental five types named penetration (P), accumulation (A), entrainment (E), longterm stagnation (LS), and short-term stagnation (SS). In these distinctive types (named pS and cS), the subducted slabs tend to stagnate at the 660 km discontinuity after they experience the penetration or the collapse into the lower mantle, yielding the shapes of slabs quite similar to those of SS at some time instances. In the second series of our experiments, we investigated the effects of the sudden stop of trench retreat on the horizontallylying slabs (LS, SS, cS, and pS types in the first series of experiments). We found that the styles of avalanche of stagnant slabs vary depending not only on the types of slab behaviors described in section 3.1 but also on the slab shapes at the time instance of the sudden stop of the trench retreat. In these experiments, we also observed wider variations of slab behaviors than in the cases with continuous trench retreat, such as the slabs which look as if they are stagnated below the 660 km discontinuity. The results of our experiments suggest that the history of trench migration (the retreat velocity and its temporal change) plays the greatest role in the formation and avalanche of stagnant slabs, among the control parameters which may affect the dynamics of subducting slabs.
The fundamental five types of slab behaviors, particularly the LS and SS, obtained in the first series of our experiments can be well-compared with those of natural ones observed by the seismic tomography. For example, the "stagnant slab" behaviors observed beneath the Northeast Japan or Izu-Bonin arcs (e.g., Zhao, 2009;Fukao and Obayashi, 2013) can be well-explained either by the LS or SS types in this study: These slabs seem to be LS, if the western tip of the flat slab in the MTZ observed beneath the East Asia arc is the tip of the Pacific slab whose subduction started about 60 Ma (e.g., Seton et al., 2012Seton et al., , 2015 and the Pacific slab has not yet penetrated into the lower mantle (Honda, 2016). On the other hand, these slabs seem to be SS, if the western tip of the stagnant Pacific slab represents the subduction which occurred about 20-30 My as suggested in Liu et al. (2017). Similarly, the slab morphology observed in the Tonga subduction zone, characterized by a horizontally-flattening portion of the slab in MTZ and the remnant slabs in the lower mantle (van der Hilst, 1995;Richards et al., 2011), is also comparable with that of SS type in this study, although the seismological interpretation is still controversial for Tonga slab (Fukao and Obayashi, 2013).
On the other hand, the occurrence of pS and cS types (stagnation after penetration or collapse) in addition to the fundamental five types in the first series of our experiments strongly indicate that the morphology of slabs similar to SS can also come out from the slabs which start to stagnate horizontally in the MTZ after their penetration or collapse into the lower mantle. In other words, it is difficult to distinguish the SS type from the pS or cS type solely from the comparison of slab shapes at some time instances given by the snapshot of slab behavior observed by seismic tomography. We can certainly speculate that, if simply based on the comparison of control parameters relevant for each type with those of the subduction zones listed above, the pS and/or cS types are more likely to occur than SS in natural subduction zones. However, more careful inspection should be necessary to accurately determine the types of slab behaviors for particular subduction zones by, for example, taking into account the detailed history of trench migration.
In addition, the occurrence of four types of slab behaviors in the second series of our experiments suggests that the temporal variations in the rate of trench retreat can significantly induce the complexities in the collapse of stagnant slabs into the lower mantle which can be hardly observed for the cases with constant rate of trench retreat. For example, the collapse of the stagnant slabs at their hinges obtained for the types I and IV can be well-compared with that observed in the Mariana arc (Fukao and Obayashi, 2013). In particular, the occurrence of horizontally-lying portion of slabs significantly below the 660 km discontinuity for the type II is very similar to that observed in the Western Java arc which looks to be "trapped"at about 1,000 km depth (Fukao and Obayashi, 2013). The occurrence of type II simply due to the sudden stop of trench retreat may therefore suggests that a "trapped" slab at some depths in the lower mantle does not necessarily imply the existence of discontinuous changes in physical properties such as viscosity (Marquardt and Miyagi, 2015;Rudolph et al., 2015) at that depth.
Although we do believe that our model successfully reproduces the slab behaviors such as the stagnant slabs at the MTZ and their avalanches into the lower mantle, we need further numerical calculations without several simplifications which we made in the present study. For example, the slab behaviors in the lower mantle, such as the flattened slab well below the 660 km discontinuity, are probably affected by the limited model thickness of 2,000 km as well as the imposed boundary condition at the bottom boundary, which may therefore need more careful validation using improved model with a full mantle thickness. In addition, the assumption of the imposed motions of subducting and overriding plates in this study are most likely to modulate the flow in the mantle wedge and, in turn, the behaviors of subducting slabs. Implicit in this assumption is a further presumption that the trench migration occurs independently of the deformation of slabs or other model parameters such as the Clapeyron slope and the viscosity jump at around 660 km depth, which may not necessarily hold in natural subduction zones. Besides, we employed the prescribed kinematic motions of subducting plate and trench (v pl and v tr ) using a highly simplified history of trench migration (constant through time or sudden stop at 40 My after the initiation of the subduction), which may not necessarily occur in natural subduction zone. In particular, the assumption of the constant rate of plate convergence at the trench (v pl + v tr ) throughout the calculations may not be fully consistent with the dynamics in the underlying mantle. Hence, it is quite important to conduct further numerical simulations by using models with the selfconsistent motions of subducting and overriding plate, including the complexities in geometry of the subduction zones such as the presence of back-arc opening (e.g., Nakakuki and Mura, 2013) and the double subduction system (e.g., Čížková and Bina, 2015;Faccenna et al., 2018), in order to fully understand the diversity of slab behaviors in the natural subduction zones.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.