Evolution of Subduction Cusps From the Perspective of Trench Migration and Slab Morphology

Introduction

During the evolution of plate tectonics, trenches located at the junction of subducting and overriding plates can develop various kinds of geometries (Schellart et al., 2007; Müller et al., 2016). When two trenches intersect with each other to form a corner, we define it as a subduction cusp. For example, at 40 Ma, the trenches along the Kurile Islands and northeast Japan intersected and formed a subduction cusp (Figure 1A) (Vaes et al., 2019). At 35 Ma, the Pacific plate subducted beneath the Eurasian and Philippine Sea plates, and the trenches along the Eastern Japan and Izu-Bonin arc were generally perpendicular to each other, forming a subduction cusp (Figure 1D) (Hall, 2002; Ma et al., 2019).

FIGURE 1

FIGURE 1. (A), (B), and (C) show the plate boundary configurations near the Kurile Islands at (A) 40 Ma, (B) 25 Ma, and (C) 0 Ma (after Vaes et al., 2019). (D), (E), and (F) show plate reconstruction results near the Philippine Sea plate at (D) 35 Ma, (E) 30 Ma, and (F) 25 Ma (after Müller et al., 2016). Red lines with filled triangles and arrows represent trenches and strike-slip faults, respectively. The Kurile and Izu-Bonin cusps are represented by blue angles. Abbreviations: PAC = Pacific plate; EUR = Eurasian plate; PSP = Philippine Sea plate; Sa = Sakhalin; and Ho = Hokkaido.

Subduction cusps can be generated under various tectonic settings. For example, a series of numerical models and analog models have shown that a subduction cusp can be formed by aseismic ridge or plateau subduction (Martinod et al., 2005; Martinod et al., 2013; Zeumann and Hampel, 2016). It is possible that the cusps linking the Kurile Islands and northeast Japan trenches and the Izu-Bonin arc and East Japan trenches (Figure 1) were generated by the subduction of oceanic plateaus that have since been entirely consumed (Rosenbaum and Mo, 2011).

Once a subduction cusp is formed, its evolutionary pathway is of particular interest in terms of trench migration and slab morphology. The evolution of the Kurile cusp and Izu-Bonin cusp has some common features regarding their trench migration. At 40 Ma, the angle between the trenches along northeast Japan and the Kurile Islands was ∼90° (Figure 1A). Then, the trench along the Kurile Islands began to retreat and rotate counterclockwise, forming a dextral strike-slip fault between Sakhalin and central Hokkaido (Figure 1B). After that, the trench along northeast Japan began to retreat at approximately 25 Ma. The retreat of these two trenches facilitated the opening of the Kurile and Japan Sea basins. The trench retreated faster at the subduction cusp, causing an increasing angle between these two trenches (Figures 1A–C) (Vaes et al., 2019). Similarly, at 35 Ma, the Pacific plate subducted beneath the Philippine Sea plate and East Asian margin. The trenches along these two overriding plates were almost perpendicular to each other (Figure 1D). The Philippine Sea plate moved northward and rotated clockwise, with the subduction of the Pacific plate, thereby forcing the triple junction (between the Pacific plate, Philippine Sea plate and Eurasian plate) to move northeastward, and the angle between the Izu-Bonin arc and Japan trench gradually became larger (Figures 1D–F) (Hall, 2002; Ma et al., 2019).

The Kurile cusp and Izu-Bonin cusp also show similar characteristics in slab morphology. From seismic tomography results, we can find that the slab dip angle varies along these two subduction zones and reaches its minimum beneath the cuspate area (Miller and Kennett, 2006; Zhang et al., 2019). Along the Kurile trench, the subducting slab steepens gradually from the cuspate area to the northeast Kurile and Kamchatka (Miller and Kennett, 2006), and the slab dip angle along the Izu-Bonin subduction zone increases southward (Zhang et al., 2019). Such variations in slab dip angle have formed a unique slab morphology beneath the Kurile Islands and Izu-Bonin arc. Part of the subducting slab penetrates into the lower mantle through the mantle transition zone, while other parts of the slab remain stagnant in the mantle transition zone (Torii and Yoshioka, 2007; Zhao et al., 2012; Zhang et al., 2019).

The 3-D evolution of subduction zones has been studied before using numerical models. For example, Bengtson and van Keken (2012) and Kneller and van Keken (2008) used a 3-D subduction model to investigate the influence of slab geometries on mantle flow, shear wave anisotropy, and temperature structure. Morishige and Honda (2013) focused on the influence of rheology on mantle flow, slab morphology and seismic anisotropy in a subduction model with a triple junction. However, few studies have investigated the evolution of subduction cusps and related trench migration patterns and slab morphologies.

In this paper, we explore the evolution of subduction cusps using a 3-D dynamic subduction model. In particular, we systematically investigate the controlling effects of the overriding plate strength, initial slab-pull force, and initial cusp angle on trench migration and slab morphology in subduction cusps.

Materials and Methods

To investigate the trench migration process and slab morphology of subduction cusps, a 3-D dynamic subduction model is developed. Leng and Gurnis (2015) modified CitcomCU (Zhong, 2006) by adding viscoelastic and material tracking methods. Our calculation method follows Leng and Gurnis (2015) to solve the mass, momentum, and energy conservation equations. Here, only the basic model setup and material rheology are introduced. For more details on the methodology, please refer to Leng and Gurnis (2011) and Leng and Gurnis (2015).

The model domain is 660 km deep, 2,970 km wide and 2,970 km in length. The “particle in cell” method makes it possible to calculate a 3-D numerical model with a relatively low resolution. We use a uniform grid of $256\times 256\times 64$ ($X$, $Y$, and $Z$ directions) in our model, leading to a resolution of $11.6\times 11.6\times 10.3$ km ($X$, $Y$, and $Z$ directions). We set 27 particles in each element. Our model includes a square subducting oceanic plate and an overriding continental plate (Figure 2A). The subducting plate (1947 $\times$ 1947 km) is located 33 km away from the lateral boundaries of the model to ensure that the subducting plate can be detached from the lateral model boundaries and subduct freely. The overriding plate is attached to the lateral boundaries (Figure 2A) and remains stable during the subduction process. Thus, trench retreat is mainly caused by the extension of the weak back-arc region (Arc1 and Arc2 in Figure 2A). The overriding plate thickness is 60 km and consists of a 20-km-thick continental crust and a 40-km-thick continental lithospheric mantle. A 7-km-thick oceanic crust overlies the subducting plate (Figure 2B). To induce subduction at the initial time step, the tip of the oceanic lithosphere is bent to reach a certain depth (165 km in the reference model) with a dip angle of 30°.

FIGURE 2

FIGURE 2. (A) Top view of our geodynamic model, including a subducting oceanic plate and an overriding continental plate. We select profiles $AA\text{'}$ ($Y$=1,155 km), $BB\text{'}$ ($X$=1815 km), and $CC\text{'}$ (diagonal direction) to show slab morphology and quantify the distance of trench migration. θ represents the initial angle of the subduction cusp. (B) Schematic diagram of the compositional field near the trench. On the left is the layered overriding plate (OP), which consists of 20-km-thick continental crust (light gray) and 40-km-thick continental lithospheric mantle (deep gray). On the right is the subducting plate (SP), where deep blue, light blue, and yellow colors represent oceanic crust and two layers of oceanic lithospheric mantle with different rheological properties. The red zone is a weak zone at the interface between the overriding and subducting plates. (C) The initial temperature field of our model along profile $AA\text{'}$. We use a linear temperature profile for the overriding continental plate and a half-space cooling model for the subducting oceanic plate.

We impose isothermal boundary conditions with a fixed temperature at 0°C at the top boundary and 1,400°C at the bottom ignoring adiabatic heating. A half-space cooling model with an age of 60 Ma is imposed as the initial temperature field of the subducting plate, and a linear temperature field is employed on the overriding continental plate (Figure 2C). The influences of the subducting plate age and overriding plate thickness on subduction dynamics have been discussed previously (Holt et al., 2015; Agrusta et al., 2017). A thick overriding plate prevents trench retreat and has the tendency to develop a large slab dip angle (Holt et al., 2015), whereas old plate subduction promotes trench retreat and usually develops a small slab dip angle, leading to a stagnant slab in the mantle transition zone (Agrusta et al., 2017). Therefore, we fix the overriding plate thickness to 60 km and the subducting plate age to 60 Ma in our model. We apply free-slip boundary conditions to all model boundaries.

An incompressible Maxwell body is used to describe the viscoelasticity of the material in our model. The following equation is employed to calculate the contribution of viscous and elastic components to the strain rate (Moresi et al., 2002):

$\begin{array}{c}{\dot{\epsilon }}_{ij}=\frac{1}{2G}{\dot{\tau }}_{ij}+\frac{1}{2\eta }{\tau }_{ij}\end{array}(1)$

where ${\dot{\epsilon }}_{ij}$, $G$, $\eta$, ${\tau }_{ij}$, and ${\dot{\tau }}_{ij}$ represent the strain rate, shear modulus, dynamic viscosity, deviatoric stress, and time rate of deviatoric stress, respectively. The non-Newtonian viscosity is temperature- and stress-dependent and can be expressed by:

$\begin{array}{c}\eta ={\eta }_0{\left(\frac{{\dot{\epsilon }}_{II}}{{\dot{\epsilon }}_0}\right)}^{\left(\frac{1}{n}-1\right)}exp\left[\frac{E}{nR}\left(\frac{1}{T}-\frac{1}{T_0}\right)\right]\end{array}(2)$

where ${\eta }_0$, ${\dot{\epsilon }}_{II}$, ${\dot{\epsilon }}_0$, $n$, $E$, $R$, $T$, and $T_0$ represent the reference viscosity, second invariant of the deviatoric strain rate tensor, reference strain rate, strain exponent, activation energy, gas constant, absolute temperature, and reference temperature, respectively. We set an upper limit, ${\eta }_{max}$, and a lower limit, ${\eta }_{min}$, for viscosity to guarantee the convergence and efficiency of our model. The values of the model constraints can be found in Table 1.

TABLE 1

TABLE 1. Model constants.

The yielding stress, ${\tau }_y$, is expressed as:

$\begin{array}{c}{\tau }_y=\mu P+C\end{array}(3)$

where $\mu$, $P$, and $C$ represent the coefficients of friction, pressure, and cohesion, respectively. In our model, $\mu$ and $C$ decrease with increasing accumulated plastic strains ${\epsilon }_p$ following these equations:

$\begin{array}{c}\mu ={\mu }_0\left[1-min\left(1,\frac{{\epsilon }_p}{{\epsilon }_f}\right)\right]\end{array}(4)$

$\begin{array}{c}C=C_f+\left(C_0-C_f\right)\left[1-min\left(1,\frac{{\epsilon }_p}{{\epsilon }_f}\right)\right]\end{array}(5)$

where ${\mu }_0$, $C_0$, $C_f$, and ${\epsilon }_f$ represent the initial coefficients of friction, initial cohesion, minimum cohesion, and reference plastic strain, respectively. In our model, ${\mu }_0$, $C_0$, and ${\epsilon }_f$ are set as 0.5, 40 MPa, and 0.3, respectively. A 22.8-km-thick pure viscous weak zone is imposed at the interface between the subducting and overriding plates (Figure 2B), and we set the viscosity of the weak zone as ${10}^{20}$ Pa s for plate decoupling. To prevent plate rupture during the subduction process, we also set a 20-km-thick pure viscous layer at a depth of 30 km in the subducting oceanic plate (Figure 2B) (Stegman et al., 2010; Schellart and Moresi, 2013). By systematically comparing their model results with interseismic deformation, Itoh et al. (2019) suggested that the Southern Kurile arc and its back-arc areas are weak. Here, we set the initial value of ${\epsilon }_p$ to 0.2 to reflect the weak rheology in the back-arc area.

Results

The Reference Model

We first run a reference case_ref (Table 2), with two weak arc regions (Arc1 and Arc2) in the overriding continental plate. The weak region is 297 km wide and 165 km away from the trench (Figure 2A). The tip of the subducting slab is bent to reach a depth of 165 km, which can provide enough initial slab-pull force to initiate subduction.

TABLE 2

TABLE 2. Parameter values for cases in this study. ${\epsilon }_{p1}$ and ${\epsilon }_{p2}$ represent the initial plastic strains of Arc1 and Arc2, respectively. $d_y$ represents the maximum depth that the front of the subducting plate can reach initially and the unit of $d_y$ is km. θ is the initial angle of subduction cusp (Figure 2A). All parameters in case_r320 are same to case_ref, except a higher model resolution, which is 320 $\times$ 320 $\times$ 80 ($X$, $Y$, and $Z$ directions) in case_r320.

At the beginning of the subduction process, the tip of the subducting slab changes its dip angle from 30° to nearly 90°, and the subducting plate begins to move slowly trenchward due to a relatively small slab-pull force (Figure 3A). Then, as the subduction process continues, the subducting plate begins to accelerate, reaching the bottom boundary of the model domain after 8.83 million years (Myr) of subduction (Figure 3B). Finally, the subducting slab stagnates on the bottom boundary after 13.25 Myr (Figure 3C). The subducting slab develops different dip angles in different directions. Here, we select profiles $AA\text{'}$, $BB\text{'}$, and $CC\text{'}$ in Figure 2A to show the variation in slab dip angle in the $X$, $Y$, and diagonal directions (Figure 4). At 13.25 Myr, the slab stagnates on the bottom boundary and develops a similar morphology along profiles $AA\text{'}$ and $BB\text{'}$ due to the symmetric model set up in the $X$ and $Y$ directions. However, in the diagonal direction along profile $CC\text{'}$, the slab dip angle is much smaller than that of profiles $AA\text{'}$ and $BB\text{'}$ (Figure 4C). Such a small slab dip angle can resist mantle flow reaching the corner of mantle wedge, thus causing a cold mantle wedge (Figure 4C).

FIGURE 3

FIGURE 3. The evolution of slab morphology at 5.52 (A), 8.83 (B), and 13.25 (C) Myr in the reference model. The color indicates the depth of the slab. We put some black tracking lines on the slab surface to better display the three-dimensional structure of the slab.

FIGURE 4

FIGURE 4. Temperature field along the profiles of $AA\text{'}$ (A), $BB\text{'}$ (B), and $CC\text{'}$ (C) at 13.25 Myr in the reference model. The locations of the $AA\text{'}$, $BB\text{'}$, and $CC\text{'}$ profiles are shown in Figure 2A.

We placed tracking tracers at the end of the subducting plate (66 km away from the plate edge) to obtain its velocity. Figure 5A shows that the subducting velocity (red lines in Figure 5A) is almost the same in the $X$ and $Y$ directions. The slab-pull force becomes larger because of a longer slab tip as the subduction process proceeds, and the velocity of the subducting plate shows a significant increase (Figure 5A). After approximately 8 Myr of subduction, the subduction velocities $v_x$ and $v_y$ increase from 0 to 8 cm/yr (Figure 5A). Then, the tip of the subducting slab reaches the bottom boundary, which resists the subduction process. After ∼4 Myr of fluctuation, the velocity of the subducting plate begins to decrease (Figure 5A). Figure 5B shows the corresponding trench locations at different times. It can be noted that the trench retreats faster in the diagonal direction than in the $X$ and $Y$ directions, leading to an increase in the cuspate corner angle from 90° to ∼150° (Figure 5B). However, the trench retreat speed and the angular speed of opening of subduction cusp doesn’t have a proportional relationship to subducting velocity (Figure 5B), because trench retreat can be influenced by many factors (Holt et al., 2015; Agrusta et al., 2017).

FIGURE 5

FIGURE 5. (A) The velocities of the subducting plate in the $X$ ($v_x$, solid line) and $Y$ ($v_y$, dashed line) directions. Red and orange colors show the results of case_ref and case_320, respectively. (B), (C), and (D) show the trench migration history in case_ref, case_ss, and case_sw, respectively. (E) The migration distances of the trench along profiles $AA\text{'}$ (solid line), $BB\text{'}$ (dashed line), and $CC\text{'}$ (dotted line) in Figure 2A. Orange, red, and blue colors show the results of case_ref, case_ss, and case_sw, respectively.

Resolution test has been made to determine the proper resolution required to resolve the subduction zone accurately (case_r320 in Table 2). The results show that the velocities of subducting plate are very comparable for a resolution of 256 $\times$ 256 $\times$ 64 (red lines in Figure 5A) and 320 $\times$ 320 $\times$ 80 (orange lines in Figure 5A). Thus, we choose a resolution of 256 $\times$ 256 $\times$ 64 for other models presented in this paper, considering both calculation accuracy and efficiency.

Influence of the Overriding Plate Strength

The overriding plate is relatively weak due to the existence of weak arc regions in the reference model. Previous 2-D studies have shown that the overriding plate strength has great influences on subduction dynamics (Garel et al., 2014; Holt et al., 2015; Yang et al., 2018). Here, we investigate the influence of the asymmetry of the weak arc region on the evolution of subduction cusps (case_sw and case_ss in Table 2). In our reference model case_ref, the overriding plate is weak in both the Arc1 and Arc2 regions (Figure 2A). In case_sw, only the Arc2 region is weak, whereas in case_ss, both weak arc regions are removed.

Figure 5C–E shows the influence of the overriding plate strength on the trench retreat process. Figures 5C, D show the evolution of trench geometries in case_ss and case_sw, respectively. In case_ss, the trench retreats slightly in the diagonal direction, changing the sharp subduction cusp to a gentle curvature locally. Nevertheless, the trench positions are generally stable in both the $X$ and $Y$ directions, and the intersection angle of these two trenches remains ∼90° throughout the whole subduction process (Figure 5C). In case_sw, the trench position remains stable in the $X$ direction generally; however, the trench retreats significantly in the $Y$ direction, causing an increasing cusp angle during the subduction process (Figure 5D). Such an asymmetric trench geometry is completely different from the symmetric subduction regime in case_ref and case_ss (Figures 5B,C). It can be observed that the trench retreats symmetrically in both case_ref and case_ss, with similar retreating distances in the $X$ and $Y$ directions but a larger retreating distance in the diagonal direction (Figure 5E). Comparing case_ss to case_ref, we can find that the trench retreat distance is smaller in case_ss because of the strong overriding plate, which makes it difficult for the trench to retreat (Figure 5E). In case_sw, the overriding plate is strong in Arc1 but weak in the Arc2 region; thus, the trench retreats asymmetrically. A larger retreat distance appears in the $Y$ direction, while in the $X$ direction, the trench retreat distance is similar to the results in case_ss (Figure 5E). Together, the results of these three cases (case_ref, case_ss, and case_sw) suggest that the weakening of the overriding plate promotes trench retreat, which is consistent with the results of previous studies (Holt et al., 2015; Agrusta et al., 2017). Thus, the asymmetric distribution of weak regions in the overriding continental plate can cause asymmetric trench migration processes. It can also be noted that the trench retreats more significantly in the diagonal direction regardless of the overriding plate strength, which typically causes the subduction cusp to become smooth and disappear during the subduction process (Figures 5B–E).

Figure 6 shows the morphology of the subducting slab in the $X$, $Y$, and diagonal directions in case_ss and case_sw. In case_ss, the slab subducts nearly vertically with a slab dip angle of ∼80° in the $X$ and $Y$ directions (Figures 6A, C). In case_sw, the subducting slab subducts asymmetrically; the slab dip angle is smaller in the $Y$ direction with the weak Arc2 region, while in the $X$ direction, the slab subducts nearly vertically (Figures 6B, D). A weak overriding plate is easier to deform and thus promotes trench retreat, which forms a smaller dip angle in the subducting slab (Holt et al., 2015). In the diagonal direction, the slab dip angle is smallest regardless of how we change the overriding plate strength (Figures 6E, F).

FIGURE 6

FIGURE 6. Temperature field along the profiles of $AA\text{'}$, $BB\text{'}$, and $CC\text{'}$ at 12.15 Myr in case_ss (A–E) and case_sw (B–F). The locations of the $AA\text{'}$, $BB\text{'}$, and $CC\text{'}$ profiles are shown in Figure 2A.

Influence of the Initial Slab-Pull Force

The slab-pull force is considered to be the main driving force behind subduction (Schellart, 2004). Here, we change the initial length of the bending tip of the subducting slab to control the initial slab-pull force (case_d215 and case_d265 in Table 2). A longer bending tip of the subducting slab provides a larger slab-pull force. In the reference case_ref, the tip of the subducting slab reaches 165 km depth beneath both the Arc1 and Arc2 regions. For the cases of case_d215 and case_d265, we keep the slab tip at 165 km depth beneath the Arc1 region but increase the slab tip to 215 and 265 km depths beneath the Arc2 region.

Figure 7A shows the influence of the initial slab-pull force on the velocities of the subducting plate in the $X$ and $Y$ directions. In case_ref, the velocities of the subducting plate are similar in the $X$ and $Y$ directions, and the slab subducts symmetrically (Figure 7A). However, when the initial slab-pull force becomes larger beneath the Arc2 region (Figure 2A) for case_d215 and case_d265, the subducting velocity in the $Y$ direction becomes significantly larger than that in the $X$ direction (Figure 7A), causing the subducting plate to rotate counterclockwise, and the subduction becomes asymmetric.

FIGURE 7

FIGURE 7. (A) The velocities of the subducting plate in case_ref (orange), case_d215 (red), and case_d265 (blue). The solid line and dashed line represent the velocities of the subducting plate in the $X$ and $Y$ directions, respectively. (B) The migration distances of the trench along profiles $AA\text{'}$ (solid line), $BB\text{'}$ (dashed line), and $CC\text{'}$ (dotted line) in Figure 2A. Orange, red, and blue colors show the results of case_ref, case_d215, and case_d265, respectively.

Figure 7B shows the influence of the initial slab-pull force on trench migration. When the initial slab-pull force in the $Y$ direction becomes larger, the trench migration distance becomes larger in the $X$ direction. Meanwhile, in the $Y$ direction, the trench migration distance becomes smaller (Figure 7B) and (Figures 8A,B). Previous studies have pointed out that when the subducting velocity increases, it becomes more difficult for the trench to retreat (Schellart, 2005). Our results confirm this point in subduction cusp evolution (Figure 7B) and (Figures 8A,B). Regardless of how we change the initial slab-pull force, the trench retreats more significantly in the diagonal direction, causing a smaller slab dip angle (Figure 7B) and (Figures 8A,B).

FIGURE 8

FIGURE 8. Trench migration history in case_d215 (A), case_d265 (B), case_a120 (C), and case_a180 (D).

Influence of the Initial Cusp Angle

In the previous cases, the slab dip angle remains smallest in the diagonal direction, even for case_ss, in which the trench retreats a very limited distance (Figures 5C) and (Figure 6E). Therefore, the overriding plate strength and the initial slab-pull force cannot be the key controlling factors driving the slab dip angle in the diagonal direction. On the other hand, the initial cusp angle may play important roles in the resulting slab morphology. We run case_a120 and case_a180 with different initial cusp angles of 120° and 180° to test this point (Table 2).

With increasing initial cusp angles, the subducting plate becomes smaller. The subduction process continues for only ∼12 Myr in case_a180 (Figure 8D). Comparing case_ref, case_a120, and case_a180 (Figures 5B) and (Figures 8C,D), it can be observed that the trench only retreats insignificantly when the initial cusp angle is large. When the initial cusp angle becomes smaller than 120°, the trench retreat distance in the diagonal direction increases drastically (Figures 5B) and (Figures 8C,D). Therefore, a smaller cusp angle leads to a faster trench retreat in the diagonal direction, which eventually smooths and destroys the cusp. Figure 9 shows the slab dip angle in the diagonal direction for case_ref, case_a120, and case_a180. With increases in the initial cusp angle from 90° to 120° and 180°, the slab dip angle in the diagonal direction increases from ∼30° to ∼45° and ∼90°. Thus, a larger cuspate corner angle causes a larger slab dip angle.

FIGURE 9

FIGURE 9. Temperature field along profile $CC\text{'}$ (Figure 2A) in case_ref (at 12.15 Myr) (A), case_a120 (at 11.04 Myr) (B), and case_a180 (at 6.63 Myr) (C).

Discussion

Evolution of the Kurile and Izu-Bonin Cusps

Plate reconstruction results suggest that the Pacific plate began to subduct beneath Kurile and Japan at ∼60 Ma, following the subduction of the Izanagi plate (Figure 1A) (Müller et al., 2016; Vaes et al., 2019). The cuspate corner angle at the Kurile Islands was approximately 90° at 40 Ma and it became increasingly larger during the subduction of the Pacific plate (Figures 1A–C), which is consistent with our model results (Figures 5B–D). The trench along the Kurile Islands retreated faster near the cuspate corner, developing a wedge-shaped extensional basin in the back-arc area (Schellart et al., 2003), which is also in high agreement with the results of our model (Figures 8A,B). The morphology of the subducting slab beneath the Kurile Islands can be obtained from seismic tomography results. Our results show that the slab dip angle is smallest in the diagonal direction and becomes increasingly larger toward two lateral edges when the initial cusp angle is relatively small (Figure 3). Seismic tomography shows similar results beneath the Kurile Islands in which the slab dip angle is small beneath the cuspate corner area and becomes larger northeastward (Miller and Kennett, 2006).

One interesting observation for the Kurile cusp is that the retreat distance of the Kurile trench is larger than that along the northeast Japan trench, suggesting asymmetric subduction. Both the overriding plate strength and initial slab-pull force may cause asymmetric cusp subduction (Figure 5D) and (Figure 8A). Many geological observations suggest that the overriding Okhotsk and northeast Japan plates both experienced extensive tectonic events before the subduction of the Pacific plate (Itoh et al., 2019; Vaes et al., 2019); therefore, the strength of both plates were likely weak. As a result, the asymmetric subduction of the Kurile cusp may have been caused by different initial slab-pull forces. According to the results of our model, we suggest that the initial slab-pull force beneath northeast Japan is larger than that beneath the Kurile Islands. After the Izanagi plate subducted beneath the Eurasian plate, the Pacific plate began to subduct, and the two plates broke up at the mid-ocean ridges (Thorkelson, 1996; Seton et al., 2015). If the plates broke earlier beneath northeast Japan, the early subduction of the Pacific plate would have caused a larger slab-pull force beneath northeast Japan than beneath the Kurile Islands.

At the Izu-Bonin cusp, the trench migration history and slab morphology are similar to those of the Kurile cusp. The angle of the cuspate corner changed from ∼90° to ∼120° between 35 Ma and 25 Ma (Figures 1D–F). The Izu-Bonin trench retreated faster near the cuspate corner and formed a wedge-shaped extensional basin in the back-arc region on the Philippine Sea plate (Hall, 2002; Ma et al., 2019). Seismic tomography results show that the slab dip angle becomes larger southward (Zhang et al., 2019). In our model, the overriding plate is a continental plate that is different from the Izu-Bonin cusp. Nevertheless, the Philippine Sea plate is a young oceanic plate, which is presumably weaker than a weak continental plate. Thus, based on our model results, it is reasonable to assume that the Izu-Bonin trench retreats dominantly due to a weak overriding plate.

Evolution of Other Cusps

Aside from the Kurile and Izu-Bonin cusps, we can find another cusp located at Solomon Sea in plate tectonic history. The Pacific plate was subducting beneath Philippine Sea plate and Solomon Sea plate, and the cusp angle became smaller continuously between 35 Ma and 20 Ma (Zahirovic et al., 2014; Seton et al., 2016). There are many other cusps on the Earth’s surface, such as the Kamchatka cusp and Alaskan cusp, whose cusp angles have remained small during the last several million years (Müller et al., 2016). These cusps have strong correlations to the subduction of aseismic ridges or buoyant blocks (Rosenbaum and Mo, 2011). At the Kamchatka cusp, the Hawaii-Emperor seamount trail subducts beneath the overriding continental plate, possibly causing a subduction cusp at the collision area. The angle of the cuspate corner remains small at the Gulf of Alaska, where the collision of the Yakutat block and Kodiak and Cobb seamount trails occurs (Mazzotti and Hyndman, 2002; Rosenbaum and Mo, 2011). Comparing the Kamchatka and Alaskan cusps with the Kurile and Izu-Bonin cusps, it should be noted that the ongoing subduction of buoyant seamount trails beneath the Kamchatka and Alaskan cusps corresponds to small cusp angles, whereas the Kurile and Izu-Bonin cusps have a tendency to become smooth and disappear. Thus, we can speculate that the subduction process of aseismic ridges or buoyant blocks causes the formation of subduction cusps. Once the subduction of aseismic ridges or buoyant blocks terminates, the subduction cusp tends to be smooth and disappear, as shown in our model results.

3-D Effects of the Trench Retreat and Slab Dip Angle

In previous 2-D studies, it has been proposed that trench retreat can strongly influence the dip angle of subducting slabs (Agrusta et al., 2017). When the trench retreats fast, the subducting slab has a small dip angle. When the trench retreats slowly, it is difficult for the subducting slab to incline, and it prefers to subduct vertically with a large slab dip angle. Here, in our 3-D subduction cusp model, we can see from Figure 5C and Figure 6 that the trench retreating distance is limited in case_ss, but the slab dip angle is still small in the diagonal direction, showing that the 3-D effects of a subduction cusp are important for slab geometry.

Effects of the Lower Mantle

Here, we did not include the lower mantle in our model for computational efficiency. Since both the viscosity jump and endothermic phase transition at a depth of 660 km could make the slab stagnate in the mantle transition zone (Torii and Yoshioka, 2007; Yoshida, 2013; Agrusta et al., 2017), such a simplification is probably reasonable. Nevertheless, a larger box with a lower mantle included is helpful for future studies of subduction cusp evolution. Moreover, our model results show that the slab dip angle in the diagonal direction is the smallest compared with that in the $X$ and $Y$ directions. Thus, the subducting slab in the diagonal direction more easily stagnates in the mantle transition zone, whereas the subducting slab with a larger dip angle in the $X$ and $Y$ directions tends to roll back or subduct vertically into the lower mantle. Seismic tomography results have shown that the subducting Pacific slab stagnates on the mantle transition zone beneath Izu-Bonin cusp, and the south part of it rolls back beneath Izu-Bonin-Mariana trench, indicating a tear of Pacific plate (Zhang et al., 2019). The special slab morphology from subduction cusp evolution may be employed to better explain seismic tomography results, especially at Kurile cusp and Izu-Bonin cusp (Miller and Kennett, 2006; Zhang et al., 2019).

Conclusion

In summary, using a 3-D dynamic subduction model, we investigate the influence of overriding plate strength, initial slab-pull force, and initial cusp angle on slab morphology and trench migration of subduction cusps. Following conclusions are obtained:

Subduction cusps have the tendency to become smooth and disappear during the subduction process.

The slab dip angle is the smallest in the diagonal direction of subduction cusps, and a larger cuspate corner angle leads to a larger slab dip angle.

The asymmetric distribution of the overriding plate strength and initial slab-pull force determines the asymmetric evolutionary pathway of subduction cusps.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the author, without undue reservation.

Author Contributions

HZ and WL contributed to conception and design of this study. XS found plate reconstruction data and plot them. HZ wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work is supported by the National Natural Science Foundation of China (41774105, 41820104004, 41688103) and the Fundamental Research Funds for the Central Universities (WK2080000144).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Agrusta, R., Goes, S., and van Hunen, J. (2017). Subducting-slab Transition-Zone Interaction: Stagnation, Penetration and Mode Switches. Earth Planet. Sci. Lett. 464, 10–23. doi:10.1016/j.epsl.2017.02.005

CrossRef Full Text | Google Scholar

Bengtson, A. K., and van Keken, P. E. (2012). Three-dimensional thermal Structure of Subduction Zones: Effects of Obliquity and Curvature. Solid Earth 3 (2), 365–373. doi:10.5194/se-3-365-2012

CrossRef Full Text | Google Scholar

Garel, F., Goes, S., Davies, D. R., Davies, J. H., Kramer, S. C., and Wilson, C. R. (2014). Interaction of Subducted Slabs with the Mantle Transition‐zone: A Regime Diagram from 2‐D Thermo‐mechanical Models with a mobile Trench and an Overriding Plate. Geochem. Geophys. Geosyst. 15 (5), 1739–1765. doi:10.1002/2014gc005257

CrossRef Full Text | Google Scholar

Hall, R. (2002). Cenozoic Geological and Plate Tectonic Evolution of SE Asia and the SW Pacific: Computer-Based Reconstructions, Model and Animations. J. Asian Earth Sci. 20 (4), 353–431. Article Pii s1367-9120(01)00069-4. doi:10.1016/s1367-9120(01)00069-4

CrossRef Full Text | Google Scholar

Holt, A. F., Becker, T. W., and Buffett, B. A. (2015). Trench Migration and Overriding Plate Stress in Dynamic Subduction Models. Geophys. J. Int. 201 (1), 172–192. doi:10.1093/gji/ggv011

CrossRef Full Text | Google Scholar

Itoh, Y., Wang, K., Nishimura, T., and He, J. (2019). Compliant Volcanic Arc and Backarc Crust in Southern Kurile Suggested by Interseismic Geodetic Deformation. Geophys. Res. Lett. 46 (21), 11790–11798. doi:10.1029/2019gl084656

CrossRef Full Text | Google Scholar

Kneller, E. A., and van Keken, P. E. (2008). Effect of Three-Dimensional Slab Geometry on Deformation in the Mantle Wedge: Implications for Shear Wave Anisotropy. Geochem. Geophys. Geosyst. 9, a–n. Article Q01003. doi:10.1029/2007gc001677

CrossRef Full Text | Google Scholar

Leng, W., and Gurnis, M. (2011). Dynamics of Subduction Initiation with Different Evolutionary Pathways. Geochem. Geophys. Geosyst. 12, a–n. Article Q12018. doi:10.1029/2011gc003877

CrossRef Full Text | Google Scholar

Leng, W., and Gurnis, M. (2015). Subduction Initiation at Relic Arcs. Geophys. Res. Lett. 42 (17), 7014–7021. doi:10.1002/2015gl064985

CrossRef Full Text | Google Scholar

Ma, P., Liu, S., Gurnis, M., and Zhang, B. (2019). Slab Horizontal Subduction and Slab Tearing beneath East Asia. Geophys. Res. Lett. 46 (10), 5161–5169. doi:10.1029/2018gl081703

CrossRef Full Text | Google Scholar

Martinod, J., Funiciello, F., Faccenna, C., Labanieh, S., and Regard, V. (2005). Dynamical Effects of Subducting Ridges: Insights from 3-D Laboratory Models. Geophys. J. Int. 163 (3), 1137–1150. doi:10.1111/j.1365-246X.2005.02797.x

CrossRef Full Text | Google Scholar

Martinod, J., Guillaume, B., Espurt, N., Faccenna, C., Funiciello, F., and Regard, V. (2013). Effect of Aseismic ridge Subduction on Slab Geometry and Overriding Plate Deformation: Insights from Analogue Modeling. Tectonophysics 588, 39–55. doi:10.1016/j.tecto.2012.12.010

CrossRef Full Text | Google Scholar

Mazzotti, S., and Hyndman, R. D. (2002). Yakutat Collision and Strain Transfer across the Northern Canadian Cordillera. Geol 30 (6), 4952–5498. doi:10.1130/0091-7613(2002)030<0495:Ycasta>2.0.Co;2

CrossRef Full Text | Google Scholar

Miller, M. S., and Kennett, B. L. N. (2006). Evolution of Mantle Structure beneath the Northwest Pacific: Evidence from Seismic Tomography and Paleogeographic Reconstructions. Tectonics 25 (4), a–n. Article Tc4002. doi:10.1029/2005tc001909

CrossRef Full Text | Google Scholar

Moresi, L., Dufour, F., and Mühlhaus, H.-B. (2002). Mantle Convection Modeling with Viscoelastic/brittle Lithosphere: Numerical Methodology and Plate Tectonic Modeling. Pure Appl. Geophys. 159 (10), 2335–2356. doi:10.1007/s00024-002-8738-3

CrossRef Full Text | Google Scholar

Morishige, M., and Honda, S. (2013). Mantle Flow and Deformation of Subducting Slab at a Plate junction. Earth Planet. Sci. Lett. 365, 132–142. doi:10.1016/j.epsl.2013.01.033

CrossRef Full Text | Google Scholar

Müller, R. D., Seton, M., Zahirovic, S., Williams, S. E., Matthews, K. J., Wright, N. M., et al. (2016). “Ocean Basin Evolution and Global-Scale Plate Reorganization Events since Pangea Breakup,” in Annual Review of Earth and Planetary Sciences. Editors R. Jeanloz, and K. H. Freeman, 44, 107–138. doi:10.1146/annurev-earth-060115-012211

CrossRef Full Text | Google Scholar

Rosenbaum, G., and Mo, W. (2011). Tectonic and Magmatic Responses to the Subduction of High Bathymetric Relief. Gondwana Res. 19 (3), 571–582. doi:10.1016/j.gr.2010.10.007

CrossRef Full Text | Google Scholar

Schellart, W. P., Freeman, J., Stegman, D. R., Moresi, L., and May, D. (2007). Evolution and Diversity of Subduction Zones Controlled by Slab Width. Nature 446 (7133), 308–311. doi:10.1038/nature05615

PubMed Abstract | CrossRef Full Text | Google Scholar

Schellart, W. P. (2005). Influence of the Subducting Plate Velocity on the Geometry of the Slab and Migration of the Subduction Hinge. Earth Planet. Sci. Lett. 231 (3-4), 197–219. doi:10.1016/j.epsl.2004.12.019

CrossRef Full Text | Google Scholar

Schellart, W. P., Jessell, M. W., and Lister, G. S. (2003). Asymmetric Deformation in the Backarc Region of the Kuril Arc, Northwest Pacific: New Insights from Analogue Modeling. Tectonics 22 (5), a–n. doi:10.1029/2002tc001473

CrossRef Full Text | Google Scholar

Schellart, W. P., and Moresi, L. (2013). A New Driving Mechanism for Backarc Extension and Backarc Shortening through Slab Sinking Induced Toroidal and Poloidal Mantle Flow: Results from Dynamic Subduction Models with an Overriding Plate. J. Geophys. Res. Solid Earth 118 (6), 3221–3248. doi:10.1002/jgrb.50173

CrossRef Full Text | Google Scholar

Schellart, W. P. (2004). Quantifying the Net Slab Pull Force as a Driving Mechanism for Plate Tectonics. Geophys. Res. Lett. 31 (7), a–n. Article L07611. doi:10.1029/2004gl019528

CrossRef Full Text | Google Scholar

Seton, M., Flament, N., Whittaker, J., Müller, R. D., Gurnis, M., and Bower, D. J. (2015). Ridge Subduction Sparked Reorganization of the Pacific Plate‐mantle System 60-50 Million Years Ago. Geophys. Res. Lett. 42 (6), 1732–1740. doi:10.1002/2015gl063057

CrossRef Full Text | Google Scholar

Seton, M., Mortimer, N., Williams, S., Quilty, P., Gans, P., Meffre, S., et al. (2016). Melanesian Back-Arc basin and Arc Development: Constraints from the Eastern Coral Sea. Gondwana Res. 39, 77–95. doi:10.1016/j.gr.2016.06.011

CrossRef Full Text | Google Scholar

Stegman, D. R., Schellart, W. P., and Freeman, J. (2010). Competing Influences of Plate Width and Far-Field Boundary Conditions on Trench Migration and Morphology of Subducted Slabs in the Upper Mantle. Tectonophysics 483 (1-2), 46–57. doi:10.1016/j.tecto.2009.08.026

CrossRef Full Text | Google Scholar

Thorkelson, D. J. (1996). Subduction of Diverging Plates and the Principles of Slab Window Formation. Tectonophysics 255 (1-2), 47–63. doi:10.1016/0040-1951(95)00106-9

CrossRef Full Text | Google Scholar

Torii, Y., and Yoshioka, S. (2007). Physical Conditions Producing Slab Stagnation: Constraints of the Clapeyron Slope, Mantle Viscosity, Trench Retreat, and Dip Angles. Tectonophysics 445 (3-4), 200–209. doi:10.1016/j.tecto.2007.08.003

CrossRef Full Text | Google Scholar

Vaes, B., Hinsbergen, D. J. J., and Boschman, L. M. (2019). Reconstruction of Subduction and Back‐Arc Spreading in the NW Pacific and Aleutian Basin: Clues to Causes of Cretaceous and Eocene Plate Reorganizations. Tectonics 38 (4), 1367–1413. doi:10.1029/2018tc005164

CrossRef Full Text | Google Scholar

Yang, T., Moresi, L., Zhao, D., Sandiford, D., and Whittaker, J. (2018). Cenozoic Lithospheric Deformation in Northeast Asia and the Rapidly-Aging Pacific Plate. Earth Planet. Sci. Lett. 492, 1–11. doi:10.1016/j.epsl.2018.03.057

CrossRef Full Text | Google Scholar

Yoshida, M. (2013). The Role of Harzburgite Layers in the Morphology of Subducting Plates and the Behavior of Oceanic Crustal Layers. Geophys. Res. Lett. 40 (20), 5387–5392. doi:10.1002/2013gl057578

CrossRef Full Text | Google Scholar

Zahirovic, S., Seton, M., and Müller, R. D. (2014). The Cretaceous and Cenozoic Tectonic Evolution of Southeast Asia. Solid Earth 5 (1), 227–273. doi:10.5194/se-5-227-2014

CrossRef Full Text | Google Scholar

Zeumann, S., and Hampel, A. (2016). Three-dimensional Finite-Element Models on the Deformation of Forearcs Caused by Aseismic ridge Subduction: The Role of ridge Shape, Friction Coefficient of the Plate Interface and Mechanical Properties of the Forearc. Tectonophysics 684, 76–91. doi:10.1016/j.tecto.2015.12.022

CrossRef Full Text | Google Scholar

Zhang, H., Wang, F., Myhill, R., and Guo, H. (2019). Slab Morphology and Deformation beneath Izu-Bonin. Nat. Commun. 10, 1310. Article. doi:10.1038/s41467-019-09279-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, D., Yanada, T., Hasegawa, A., Umino, N., and Wei, W. (2012). Imaging the Subducting Slabs and Mantle Upwelling under the Japan Islands. Geophys. J. Int. 190 (2), 816–828. doi:10.1111/j.1365-246X.2012.05550.x

CrossRef Full Text | Google Scholar

Zhong, S. (2006). Constraints on Thermochemical Convection of the Mantle from Plume Heat Flux, Plume Excess Temperature, and Upper Mantle Temperature. J. Geophys. Res. 111 (B4), B04409. Article. doi:10.1029/2005jb003972

CrossRef Full Text | Google Scholar