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.

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 A A ' , B B ' , and C C ' 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 A A ' and B B ' due to the symmetric model set up in the X and Y directions. However, in the diagonal direction along profile C C ' , the slab dip angle is much smaller than that of profiles A A ' and B B ' (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).

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).

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 × 256 × 64 (red lines in Figure 5A) and 320 × 320 × 80 (orange lines in Figure 5A). Thus, we choose a resolution of 256 × 256 × 64 for other models presented in this paper, considering both calculation accuracy and efficiency.

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).

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 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).

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.

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.

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.

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.

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).