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

The geometries of trenches vary worldwide due to continuous plate boundary reorganization. When two trenches intersect to generate a corner, a subduction cusp is formed. Although subduction cusps are frequently observed throughout historical plate movement reconstructions, few studies have been conducted to explore the controlling factors of trench migration and slab morphology along subduction cusps. Here, we use a 3-D dynamic subduction model to explore the influence of the overriding plate strength, initial slab-pull force, and initial cusp angle on the evolution of subduction cusps. Our numerical model results suggest the following: 1) subduction cusps have a tendency to become smooth and disappear during the subduction process; 2) the slab dip angle is smallest in the diagonal direction of the subduction cusp, and a larger cuspate corner angle leads to a larger slab dip angle; 3) the asymmetric distribution of the overriding plate strength and initial slab-pull force determine the asymmetric evolutionary pathway of subduction cusps. Our results provide new insights for reconstructing the evolution of subduction cusps from seismological and geological observations.


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).
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 strikeslip 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 . 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  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.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 783409 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 × 256 × 64 (X, Y, and Z directions) in our model, leading to a resolution of 11.6 × 11.6 × 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 × 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 40km-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 FIGURE 2 | (A) Top view of our geodynamic model, including a subducting oceanic plate and an overriding continental plate. We select profiles AA' (Y 1,155 km), BB' (X 1815 km), and CC' (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'. We use a linear temperature profile for the overriding continental plate and a half-space cooling model for the subducting oceanic plate.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 783409 lithosphere is bent to reach a certain depth (165 km in the reference model) with a dip angle of 30°. 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): where _ ε ij , G, η, τ ij , and _ τ 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: where η 0 , _ ε II , _ ε 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, η max , and a lower limit, η min , for viscosity to guarantee the convergence and efficiency of our model. The values of the model constraints can be found in Table 1.
The yielding stress, τ y , is expressed as: where μ, P, and C represent the coefficients of friction, pressure, and cohesion, respectively. In our model, μ and C decrease with increasing accumulated plastic strains ε p following these equations: where μ 0 , C 0 , C f , and ε f represent the initial coefficients of friction, initial cohesion, minimum cohesion, and reference plastic strain, respectively. In our model, μ 0 , C 0 , and ε 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 ε p to 0.2 to reflect the weak rheology in the back-arc area.

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  2 | Parameter values for cases in this study. ε p1 and ε 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 × 320 × 80 (X, Y , and Z directions) in case_r320. 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 AA', BB', and CC' 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' and BB' due to the symmetric model set up in the X and Y directions. However, in the diagonal direction along profile CC', the slab dip angle is much smaller than that of profiles AA' and BB' (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. Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 783409 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.

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°t hroughout 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  Figure 2A. Orange, red, and blue colors show the results of case_ref, case_ss, and case_sw, respectively.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 783409 7 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).

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

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  Figure 2A.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 783409 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.

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 . 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  Figure 2A. Orange, red, and blue colors show the results of case_ref, case_d215, and case_d265, respectively.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 783409 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

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 Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 783409 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 . 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.