The Mode of Trench-Parallel Subduction of the Middle Ocean Ridge

Introduction

A mid-ocean ridge (MOR) is a divergent boundary of plate tectonics where new plates are continuously generated (Turcotte and Schubert, 2002). Although MORs are generally located at the centers of ocean basins, such as the present mid-Atlantic ridge, they may subduct into the deep mantle when tectonic environments cause MORs to move to ocean trenches (Georgieva et al., 2019). Typical examples include the subducted MOR between the Izanagi plate and the Pacific plate, which subducted at 55–45 Ma (Müller et al., 2016) and the subducted MOR between the Farallon plate and the Pacific plate, which subducted at 30–0 Ma (Ferrari et al., 2018). According to the contact angle between the MOR and the trench, MOR subduction can be divided into high-angle (trench-vertical) subduction and low-angle (trench-parallel) subduction (Wu and Wu, 2019).

Different methods of MOR subduction may cause different evolutionary pathways and surface expression. On the one hand, for trench-vertical MOR subduction, large-scale MORB-like igneous rocks erupt at the intersections between MORs and trenches, which have a small influence range but a long duration. Symmetrical igneous rock sequences, such as adakite-granitoid-adakite sequences, along the MOR are formed on the overriding plate (Ling et al., 2009; Isozaki et al., 2010). Analog model experiments further show that ridge subduction has a strong impact on the accumulation and migration of the accretionary wedge and the shape of the subduction zone edge (Wang et al., 2019). On the other hand, for trench-parallel MOR subduction, large-scale MORB-like igneous rocks erupt at the trench, which has a large influence range but a relatively short duration (Gutiérrez et al., 2005; Müller et al., 2016; Liu et al., 2020). For instance, early Cenozoic MORB bedrocks have also been found in Hokkado, Japan (Maeda and Kagami, 1996; Nanayama et al., 2019). The seafloor age can also affect the local morphology of the slab (Hu et al., 2016).

Many numerical studies have focused on the subduction evolution of a continuous oceanic plate (e.g., Gurnis and Hager, 1988; Christensen, 1996; Torii and Yoshioka, 2007; Čížková and Bina, 2013; Garel et al., 2014; King et al., 2015; Yang et al., 2018), whereas only a few studies have been conducted on MOR subduction at plate boundaries. Trench-vertical aseismic ridge subduction can affect the slab geometry and the deformation of the overriding continental plate based on analog model results (Martinod et al., 2013). Specifically, aseismic ridge subduction can provide extra buoyancy to reduce the dip angle of the slab, such as flat slab subduction beneath the South American plate (van Hunen et al., 2002; Gutiérrez et al., 2005; Rosenbaum et al., 2005; Espurt et al., 2008; Hu et al., 2016). Furthermore, when an active ridge subducts into the mantle, a slab window forming strong magmatism can be generated (Li et al., 2020). Groome and Thorkelson, (2009) showed that the slab window migrates with the subduction of an active ridge. Moreover, the spreading velocity and cooling age of ridges have a strong influence on ridge-inversed subduction (Qing et al., 2021). Burkett and Billen (2009) suggested that the change in buoyancy of the weakened zone of lithosphere close to the MOR causes the slab detachment from the subducting plates. The morphology of the MOR is affected by the MOR spreading speed (Püthe and Gerya, 2014). Overall, since much attention has been given to MOR subduction, the specific mode of trench-parallel subduction of the MOR needs to be thoroughly analyzed.

MOR subduction along the western boundary of the Pacific plate in the early Cenozoic and along the eastern boundary of the Pacific plate at present is primarily viewed as trench-parallel (Seton et al., 2015; Müller et al., 2016; Tang et al., 2018; Liu et al., 2020). In this study, we focus on the different modes of trench-parallel MOR subduction and the related key controlling parameters. Specifically, the pulling force from the slab connecting to the foregoing oceanic plate and the motion of the following oceanic plate may play important roles in MOR subduction. The initial dip angle of the slab and the viscosity structure of the mantle also affect the plate movement. Therefore, we investigate the effects of foregoing slab length, the following plate speed, the initial dip angle of the slab and the viscosity structure of the mantle on the modes of MOR subduction and analyze the surface expression, including the topography and volcanic activity, on the overriding plate.

Methods

We construct a 2-D Cartesian viscoelastoplastic model to study the evolution of plate movement during MOR subduction based on the finite element method. We use the Ellipsis approach (Moresi et al., 2003) to solve the governing equations of mass, momentum and energy. The movement and deformation of the lithosphere and mantle are tracked by Lagrangian particles (Moresi et al., 2003; Leng and Gurnis, 2015; Chen et al., 2020).

An incompressible Maxwell material body is used to describe the viscoelasticity of the model composition. The deformation rate of the material composition is described by the viscous part and elastic part:

where ${\dot{\epsilon }}_{ij}$ is the strain rate, $G$ is the shear modulus, ${\dot{\tau }}_{ij}$ is the time variation rate of deviatoric stress, $\eta$ is the viscosity, and ${\tau }_{ij}$ is the deviatoric stress. $i$ and $j$ represent spatial indices.

The viscosity equation we use is non-Newtonian and is temperature- and depth-dependent (Karato and Wu, 1993):

where ${\eta }_0$ is the reference viscosity, ${\dot{\epsilon }}_{{\rm I}{\rm I}}$ is the second invariant of the deviatoric strain rate tensor, ${\dot{\epsilon }}_0$ is the reference strain rate, $n$ is the strain exponent, $E$ is the activation energy, $R$ is the gas constant, $T$ is the absolute temperature (in Kelvin), $T_0$ is the reference temperature, ${\eta }_z$ is the pre-exponential factor of depth-dependent viscosity, $d$ is the model box thickness, and $z$ is the depth (Table 1). The depth-dependent mantle viscosity structure in the upper mantle is variable with uncertainties (Steinberger and Calderwood, 2006; Čížková et al., 2012). Therefore, we consider two possible viscosity structures: Temperature-dependent viscosity (ViscT) versus Temperature- and depth-dependent viscosity (ViscT&D). For the ViscT viscosity structure, ${\eta }_z=0$. For the ViscT&D viscosity structure, ${\eta }_z=\ln \left(50\right)=3.912$.

TABLE 1. Physical parameters in the model.

The plasticity of the material is applied with the Drucker-Prager criterion through the yield stress:

where ${\tau }_y$ is the yield stress, $\mu$ is the coefficient of friction, $P$ is the pressure and $C$ is the cohesion, which is the strength of the material at zero hydrostatic pressure. As the plastic strain increases, the mantle material becomes weaker (Buck and Poliakov, 1998; Poliakov and Buck, 1998). Therefore, both $\mu$ and $C$ linearly decrease when the accumulated plastic strain increases:

where ${\mu }_0$ is the initial coefficient of friction, ${\epsilon }_p$ is the accumulated plastic strain, ${\epsilon }_f$ is the reference plastic strain, $C_0$ is the initial cohesion and $C_f$ is the minimum cohesion (Table 1). The material no longer weakens when ${\epsilon }_p$ becomes larger than the reference plastic strain. To decouple the overriding plate and the subducting plate, we placed a long strip of weak zones along the interface between the continental lithosphere and the subducted oceanic plate (28 km in width, Figure 1). The plastic parameters of the weak zone, ${\mu }_0$, ${\epsilon }_f$, $C_0$ and $C_f$, are set to 0.0, 0.01, 0.6 and 0.6 MPa, respectively (Muhuri et al., 2003; Leng and Gurnis, 2011; Leng and Gurnis, 2015).

FIGURE 1. Model initial setting and boundary conditions. (A) The 2-D model is 4,000 $\times$ 1,000 km with an air layer 40 km thick at the top. At 3,300 km in the horizontal direction, we set up an MOR structure, and the two sides of the MOR are the foregoing oceanic plate and the following oceanic plate. The three white squares represent the sampling points placed in the model. Different colors in the model represent different lithologies, as shown at the bottom. $\theta$ is the initial dip angle of the foregoing plate. The initial slab length (L) and the imposed velocity ($V_x$) at the right boundary are the main control parameters in our model. (B) Initial temperature field. The white dashed line is the 1,300°C isotherm, which represents the bottom boundary of the lithosphere.

Our 2-D Cartesian model has a length of 4,000 km and a depth of 1,000 km (Figure 1). To simulate the top free-surface boundary, we set a sticky air layer 40 km thick at the top with zero density and minimum viscosity (Leng and Gurnis, 2011, 2015; Chen et al., 2020). The grid of our model is 513 $\times$ 129, which gives a resolution of 7.8 km in both the horizontal and vertical directions. To track the evolution and migration of the materials accurately, at the beginning of the model, we placed 32 particles uniformly distributed in each cell. Since our integration points for the finite element method are put on these tracking particles, our model can well resolve the compositional boundaries during model evolution (Moresi et al., 2003).

The upper, lower and left boundaries of the model are all set as free-slip boundary conditions, while at the right boundary, we impose a velocity boundary condition. A $tanh$ function of horizontal velocities is applied from the plate upper surface to the bottom of the model to ensure the conservation of mass during the evolutionary process of the model. The top and bottom boundary temperatures of the model are fixed at 0°C and the reference temperature (Table 1), respectively.

We divide the lithosphere into three parts in the model, i.e., the continental plate (x = 0–2,600 km), the foregoing oceanic plate (x = 2,600–3,300 km) and the following oceanic plate (x = 3,300–4,000 km) (Figure 1A). For the continental plate, the thickness of the lithosphere is 100 km with a 30-km thick crust. The initial temperature field of the continental lithosphere increases linearly from the surface to the bottom (Figure 1B). For the foregoing and following oceanic plates, the oceanic crust thickness is 7 km. The density of the continental and oceanic crust is 2,800 $\text{kg}/m^3$, and the density of the lithosphere is uniformly 3,300 $\text{kg}/m^3$. We set an initial slab with a length of L connecting to the foregoing oceanic plate, which has subducted into the mantle. The foregoing and the following oceanic plates are separated by the MOR at 3,300 km in the horizontal direction. We set the age of the oceanic plate at the MOR to 0 Myr, which linearly increases to 80 Myrs at both ends of the foregoing and following oceanic plates (i.e., at x = 2,600 and 4,000 km, respectively) (Figure 1B). The purple color in Figure 1A shows the oceanic plate thickness, which is determined by the thermal boundary layer thickness $y_L$ (Turcotte and Schubert, 2002):

where $\kappa$ is the thermal diffusivity and $t$ is the oceanic plate age. The initial temperature fields of the oceanic plate are computed with the half-space cooling model (Turcotte and Schubert, 2002). We take the 1,300°C isotherm as the boundary between the lithosphere and the mantle. The initial temperature field of the mantle increases linearly from 1,300°C at the bottom of the lithosphere to $T_0$ = 1,500°C at the bottom of the box (Figure 1B).

We use a simple model to estimate the mantle melting extent (Katz et al., 2003). The degree of mantle melting $F\left(P,T,X_{H_{2}O}^{bulk}\right)$ is calculated based on the temperature, pressure and water content:

where T is the temperature, $P$ is the pressure, $X_{H_{2}O}^{bulk}$ is the weight fraction of bulk water, $X_{H_{2}O}$ is the weight fraction of water in liquid phase, $T_{solidus}$ is the solidus temperature, $\Delta T\left(X_{H_{2}O}\left(X_{H_{2}O}^{bulk},P,T\right)\right)$ is the temperature decrease in the solidus at a water content and $T_{liquidus}^{lherz}$ is the lherzolite liquidus temperature.

For the water content part, we consider that the mantle area within 60 km above the slab surface is wet, and the water content of this layer is assumed to be 0.3% (Rüpke et al., 2004; Leng et al., 2012).

Results

The slab morphology is mainly affected by ridge push, slab pull, friction resistance, cohesion, among which the ridge push and slab pull force are the two main driving forces (Turcotte and Schubert, 2002). The slab pull force increases with the slab length. The imposed velocity is the main driving force of the following plate motion. So we run a total of 46 cases (Table 2) to study the effects of the initial slab length, imposed velocity, initial subduction angle and mantle viscosity structure on the evolution of trench-parallel MOR subduction. The initial length of the slab varies in the range of 200–800 km; the imposed velocity varies in the range of 1–10 cm/yr; the subduction angles are 30°, 45°, and 60°, respectively; and the depth-dependent viscosity structure includes the ViscT viscosity structure and ViscT&D viscosity structure, as indicated in the Methods section. We placed sampling points in the lithosphere of the foregoing and following oceanic plates to detect the relative motion of the two oceanic plates. To quantify the relative motion between the foregoing and following oceanic plates, we define the velocity difference $dv$ between the two oceanic plates as (the velocity direction is positive to the left):

where $v_1$ is the velocity of the foregoing oceanic plate and $v_2$ is the velocity of the following oceanic plate. $\Delta v=0$ indicates that the foregoing and the following oceanic plates have the same speed, whereas $\Delta v>0$ or $\Delta v<0$ indicates that the speed of the foregoing oceanic plate is greater or less than that of the following oceanic plate. Based on our model results, we find that $\Delta v$ is a key factor influencing the mode of trench-parallel subduction of the MOR.

TABLE 2. Case parameters and resulting modes of MOR subduction.

By analyzing and comparing all the model results, we divide the trench-parallel subduction of the MOR into three modes: 1) In the fast spreading mode, where the $\Delta v$ value is positive (case RS01, the red line in Figure 2A), representing a much faster velocity of the foregoing oceanic plate than that of the following oceanic plate. The MOR region experiences rapid expansion, and the melting area becomes wider (Figure 2B). 2) In the slow spreading mode, the $\Delta v$ value is close to 0 (case RS12, the blue dashed line in Figure 2A), indicating that the motion speed of the foregoing oceanic plate is almost the same as that of the following oceanic plate. The MOR area expands slowly, and the melting area is basically unchanged (Figure 2C). 3) In the extinction mode, the $\Delta v$ value is negative (case RS20, the black dot-dash line in Figure 2A). The following oceanic plate moves much faster than the foregoing oceanic plate, such that the following oceanic plate collides with the foregoing oceanic plate, causing the MOR structure to disappear. After ∼13.0 Myrs, the $\Delta v$ value increases to nearly 0, indicating that the foregoing and following oceanic plates merge and move at the same speed.

FIGURE 2. (A) Time evolution of $\Delta v$ for cases RS01, RS12 and RS20 with different MOR subduction modes. (B–D) show the zoomed-in area of the MOR for case RS01 at 5 Myrs, case RS12 at 18 Myrs and case RS20 at 15 Myrs. The white dashed lines represent the extent of melting with an interval of 5%. The red inverted triangle represents the location of the MOR.

According to the results of the model, the main factor controlling the subduction of the following plate is the imposed velocity $V_x$. Even if $V_x$ is a small value, such as $V_x$ = 1 cm/yr, the following plate can continue to subduct (e.g. cases RS10, RS11). However, when there are oceanic plateaus or microcontinents on the following plate, the imposed velocity $V_x$ can be reduced to 0 because of the negative density anomaly of the oceanic plateaus or microcontinent-arc collision. Therefore, when $V_x$ value is 0, the following plate lacks driving force and cannot subduct. For example, the Palawan subduction stopped at ∼20 Ma because of the microcontinent-arc collision (Keenan et al., 2016). This situation can be considered as a special case of the extinction mode. Later, we will discuss the other situation in detail where $V_x$ is greater than 0, that is, the following oceanic plate can subduct into the mantle.

The Fast Spreading Mode

Case RS01 represents a typical example of the fast spreading mode of MOR trench-parallel subduction. In this case, the length of slab L, the imposed velocity $V_x$, the initial dip angle of the slab $\theta$ and the viscosity structure of the mantle are set to 600 km, 3 cm/yr, 45° and the ViscT viscosity structure, respectively. Figure 3 shows the evolutionary process of case RS01. Due to the reference plastic strain and yield stress, the oceanic plate experiences strong weakening (Figures 4A–F). Within 0–6 Myrs (Figures 3A,B), the negative buoyancy of the slab is strong such that the foregoing oceanic plate is dragged by the subducted slab and quickly sinks into the mantle. Due to the release of water carried by the slab, the mantle wedge area experiences intense melting and weakening (Figures 3B, 4B). Therefore, the overriding continental plate experiences intense magmatic events. The topography of this area rapidly decreases due to slab subduction (the red line in Figure 5, 0–7 Myrs).

FIGURE 3. Model evolution of case RS01 (L = 600 km, $V_x$ = 3 cm/yr, $\theta$ = 45°) with the ViscT viscosity structure. (A–F) Snapshots at 0, 6, 7.5, 10, 14 and 36 Myrs. The white dashed line represents the extent of melting with an interval of 5%. The red inverted triangle represents the location of the MOR.

FIGURE 4. The viscosity structure for case RS01 (L = 600 km, $V_x$ = 3 cm/yr, $\theta$ = 45°) with the ViscT viscosity structure. (A–F) Snapshots at 0, 6, 7.5, 10, 14 and 36 Myrs.

FIGURE 5. Time evolution of the topography on the overriding continental plate for case RS01 (the fast spreading mode), case RS12 (the slow spreading mode) and case RS20 (the extinction mode). Three different background colors show the three stages for case RS01 of the fast spreading mode. Stage 1 is 0–7.5 Myrs, when the foregoing oceanic plate subducts into the mantle; stage 2 is 7.5–14 Myrs, when the MOR subducts; stage 3 is 14–65 Myrs, when the following oceanic plate subducts into the mantle.

In this case, the speed of the foregoing oceanic plate is much faster than that of the following oceanic plate; therefore, the MOR (the red inverted triangle in Figures 3B,C) has always been at the right boundary of the foregoing oceanic plate and moved forward with the foregoing oceanic plate. The rapidly expanding area between the MOR and the following oceanic plate is filled by the newly formed oceanic plate (Figures 3B, 4B, X = 2,800–3,100 km). The thickness of the newly formed oceanic plate is relatively thin because of the rapid extension of this area.

At 7.5 Myrs (Figures 3C, 4C), the foregoing oceanic plate has completely subducted into the mantle. The front segment of the following oceanic plate starts to contact the continental plate, marking the beginning of the subduction of the MOR. The buoyant MOR was strongly compressed with the overriding continental plate during the subduction process, resulting in a significant uplift of the continental plate near the trench (the red line in Figure 5, 7.5–14 Myrs). With the subduction of the MOR and the continuous forward movement of the following oceanic plate, the melting area of the mantle below the MOR gradually cools down, and the degree of melting decreases until melting disappears (Figures 3C–E).

From 10 to 14 Myrs (Figures 3D,E, 4D,E), the following oceanic plate starts to subduct, generating a new subduction zone. The thickness of the newly formed oceanic plate by the rapid expansion of the MOR is small, and the strength of the plate is weak; therefore, under the dual effect of the continuous compression of the following ocean plate and the obstruction of the front continental plate, the newly formed thin oceanic plate near the trench ruptures and gradually forms a new subduction (Figures 3E, 4E).

After 36 Myrs, a new subduction zone eventually forms. Mantle wedge melting resumes due to water transport from the following oceanic plate (Figures 3F, 4F). The continuous subduction of the following oceanic plate since 14 Myrs also results in the continuous reduction in the topography of the continental plate near the trench (the red line in Figure 5, 14–60 Myrs). During active MOR subduction, a gap of melting events and rapid topographic uplift occurs on the overriding continental plate. This melting gap divides the igneous events on the overriding plate into two separate periods.

The Slow Spreading Mode

Case RS12 is a typical case of the slow spreading mode, which has the same parameters as case RS01 except that the viscosity structure becomes a ViscT&D viscosity profile and the slab length increases to 700 km. Similar to case RS01, the foregoing oceanic plate in case RS12 continuously subducts due to slab pull. However, because of the change in the viscosity structure, the mantle resistance to the subducted plate increases, resulting in the decreased velocity of the foregoing oceanic plate. Therefore, the expansion of the MOR between the foregoing and following oceanic plates becomes slow, and the melting region of mantle beneath the MOR maintains a triangular shape (Figures 6A,B). Similar to case RS01, the topography on the continental plate near the trench also continuously decreases due to the subduction of the foregoing oceanic plate, except that the magnitude of the topography change is much smaller than that in case RS01 (the blue dashed line in Figure 5). At 18 Myrs (Figure 6B), the foregoing oceanic plate completely subducts into the mantle, and the MOR begins to contact the continental plate. Correspondingly, the topography of the overriding continental plate begins to rise due to the subduction of the buoyant MOR. Nevertheless, unlike case RS01, the newly formed oceanic plate is relatively thick and strong, and no fractures or distortions are observed during the formation of the new subduction zone (Figures 6C,D). With the continuous subduction of the following oceanic plate, the topography of the overriding continental plate changes from uplift to subsidence (the blue dashed line in Figure 5). The melting in case RS12 also experiences a gap during MOR subduction (Figures 6C,D), which divides the surface volcanism into two time periods.

FIGURE 6. Model evolution of case RS12 (L = 700 km, $V_x$ = 3 cm/yr, $\theta$ = 45°) with the ViscT&D mantle viscosity structure. (A–D) Snapshots at 5, 18, 33 and 50 Myrs. The white dashed line represents the extent of melting with an interval of 5%. The red inverted triangle represents the location of the MOR.

The Extinction Mode

Case RS20 represents a typical example of the extinction mode (Figure 7). In this case, the length of slab L, the imposed velocity $V_x$, the initial dip angle of the slab $\theta$ and the viscosity structure of the mantle are set to 200 km, 5 cm/yr, 45° and the ViscT&D viscosity structures, respectively. At the beginning, the length of the slab connecting to the foregoing oceanic plate is short (L = 200 km). The resulting negative buoyancy from the slab pull is less than the resistance force. Therefore, the foregoing plate cannot proceed to subduct into the mantle. On the other hand, the following oceanic plate is pushed to move under the imposed velocity boundary condition. Because the MOR has the characteristics of high temperature and low viscosity, the following oceanic plate starts to form new subduction at the MOR zone under strong compression, resulting in the closure and cooling of the MOR region (Figures 7A,B). The subducted part of the following oceanic plate is the youngest and thinnest plate; therefore, before it subducts into the deep mantle, it is coupled with the foregoing oceanic plate. The extinction of the original MOR forms a new plate with a thick coupled region (Figure 7C), which is driven by the velocity boundary condition to move forward as a whole and to subduct beneath the overriding continental plate (Figures 7C,D).

FIGURE 7. Model evolution of case RS20 (L = 200 km, $V_x$ = 5 cm/yr, $\theta$ = 45°) with the ViscT&D viscosity structure. (A–D) Snapshots at 0, 5, 15 and 33 Myrs. The white dashed line represents the extent of melting with an interval of 5%. The red inverted triangle represents the location of the MOR.

The subduction of the extinct MOR does not show a melting gap, as in case RS01 and case RS12 (Figure 7). The topography of the overriding plate has experienced continuous subsidence, and only a slight uplift occurs during the collision between the coupled region and the overriding plate (the black dotted-dashed line in Figure 5).

The Role of Slab Length (L) and the Imposed Velocity (${\mathit{V}}_{\mathit{x}}$)

The negative buoyancy of the slab connecting to the foregoing oceanic plate and the imposed velocity at the right boundary are the key driving factors for MOR subduction. We investigate the effects of different L and $V_x$ values on the MOR subduction (cases RS09, RS10, RS11, RS12, RS13, RS14, RS15, RS16, RS17, RS18, RS20, RS22, RS23, RS24, RS27, RS29, RS30 and RS32).

We run three cases (RS22, RS23 and RS24) with the same $V_x$ of 2 cm/yr but different initial slab lengths L of 200, 500 and 700 km. The model results show that when L becomes larger, the increased negative buoyancy of the slab leads to a larger subduction velocity of the foregoing plate. Consequently, the subduction mode of the MOR changes from the extinction mode to the slow spreading mode and eventually to the fast spreading mode (Figure 8). The oceanic plate age is also an important factor affecting the motion of plate. For the foregoing oceanic plate, the old plate can provide more negative buoyancy and drive the foregoing plate to subduct at a fast speed. The influence of old plate is similar to that of long initial slab length L. Therefore, we do not quantitatively explain the effect of oceanic plate age.

FIGURE 8. Influences of the slab length L and the imposed velocity $V_x$ on the modes of trench-parallel MOR subduction. All of these cases have a ViscT&D mantle viscosity structure and an initial slab dip angle of $\theta$ = 45°.

On the other hand, we fix the slab length L to 600 km and vary the imposed velocity $V_x$ from 1 to 3 cm/yr and 5 cm/yr for cases RS54, RS58 and RS60, respectively. The results show that the subduction mode of the MOR changes from the fast spreading mode to the slow spreading mode and then to the extinction mode (Figure 8). Therefore, the increased L and decreased $V_x$ effectively increase the magnitude of $\Delta v$, making the trench-parallel subduction of MOR move to the fast spreading mode. In contrast, the decreased L and increased $V_x$ decrease the magnitude of $\Delta v$, leading to the slow spreading mode of MOR subduction. When $\Delta v$ is small enough, we even observe shrinkage and extinction of the MOR. The following and foregoing oceanic plates strongly collide and merge into a coupled and thickened region (Figure 7C).

The Role of Other Parameters

Using case RS01 as a reference case, we refine the grid of our model to achieve a resolution of 2.0 km in both the horizontal and vertical directions. The results show that there is only a small difference from the previous model (Figure 2A). Then, we test the influence of the initial dip angle of the slab and the initial viscosity structure profile of the model on MOR subduction (cases RS02, RS03, RS04, RS05, RS06, RS07, RS08, RS19, RS21, RS25, RS26, RS28, RS31, RS33, RS34, RS35, RS36, RS37, RS38, RS40, RS41, RS42 and RS43). Cases RS40 and RS41 have the same parameters as case RS01 except that the initial dip angle $\theta$ changes to 30° and 60°, respectively. From the model results, the evolutionary processes of cases RS40, RS01, and RS41 are only slightly different from each other (Figure 9). These results indicate that under the same initial viscosity structure profile, even if the initial dip angles of the three models are different, the subduction angle of the slab becomes nearly vertical in a very short time. Therefore, the initial dip angle of the slab has little effect on the evolution of the subduction mode of the MOR.

FIGURE 9. The effects of different initial slab dip angles $\theta$ and mantle viscosity structures on the subduction of the MOR. We choose to show snapshots of the model when the foregoing oceanic plate is about to completely subduct. (A–C) show cases RS40, RS01 and RS41 at 6 Myrs when $\theta$ values are 30°, 45° and 60°, respectively, with the ViscT mantle viscosity structure. (D–F) show cases RS47, RS58 and RS19 at 21.5 Myrs when $\theta$ values are 30°, 45° and 60°, respectively, with the ViscT&D mantle viscosity structure.

We then vary the mantle viscosity structure between ViscT and ViscT&D and run case RS58. Cases RS01 and RS58 have the same parameters except for the initial viscosity structure. Since the change in the viscosity structure significantly affects the sinking velocity of the slab, the subduction mode of the MOR changes from the fast spreading mode (case RS01, Figure 9B) to the slow spreading mode (case RS58, Figure 9E). For the initial dip angle $\theta$ of 30°, when the mantle viscosity structure changes from ViscT to ViscT&D, the subduction mode of the MOR changes from the fast spreading mode to the extinction mode (Figures 9A,D). This similar variation also appears when the initial dip angle $\theta$ is 60° (Figures 9C,F). Therefore, the mantle viscosity structure may greatly affect the mode of trench-parallel subduction of the MOR.

Discussion

Implications for the Subduction of the Eastern and Western Ridges of the Pacific Plate

From our model results, for the fast-spreading subduction and slow-spreading modes, igneous rock event gaps and rapid topographic uplift on the overriding continental plate occur during active MOR subduction. Two typical subduction events of MORs are the subduction of the Western Pacific ridge (Izanagi-Pacific MOR) in the Late Cretaceous (Müller et al., 2016) and the subduction of the Eastern Pacific ridge (Pacific-Farallon MOR) from 30 Ma to the present (Ferrari et al., 2018).

The Izanagi-Pacific MOR that subducted toward the eastern Eurasian lithosphere underwent trench-parallel subduction at 55–43 Ma, as shown in Figure 10A (Müller et al., 2016). During the Late Cretaceous, the velocity of the Izanagi plate was 20 cm/yr, which was much larger than that of the Pacific plate moving at a speed of 7 cm/yr on the other side of the MOR (Wu and Wu, 2019). Similar to the model of the fast spreading mode, the velocity difference $\Delta v$ between the two sides of the MOR was much greater than 0. In this period, the young, hot and high-topographic MOR subducted under the Eurasian plate, which had strong positive buoyancy and enhanced coupling with the eastern edge of the Eurasian plate. The eastern Eurasian margin changed from regional extension to local compression (Liu et al., 2020; Zhang et al., 2020). During this period, the Northeast China regional topography also experienced obvious tectonic uplift (Wang et al., 2013; Song et al., 2014), similar to case RS01 at 7.5–14 Myrs when the overriding continental plate experienced rapid uplift (the red line in Figure 5). The distribution of igneous rocks since the Late Cretaceous is shown in Figures 10A,B. Based on the eruption time of the igneous rocks, the active volcanic activities in eastern Eurasia since the Late Cretaceous can be divided into two time periods: event 1 (100–55 Ma) and event 2 (45 Ma–present). Late Cretaceous igneous rocks are more widely distributed in NE Asia, but there are fewer igneous rocks after the Cenozoic. From the relationship between the $\epsilon Nd\left(t\right)$ values and ⁸⁷Sr/⁸⁶Sr values over time (Figures 10D,E), there is an obvious quiet igneous period between 55 and 45 Ma and a lack of igneous rock activity (Tang et al., 2018; Wu and Wu, 2019; Liu et al., 2020). Interestingly, the igneous rocks that formed before 55 Ma are significantly different from those formed after 45 Ma in terms of geochemical properties. The igneous rocks that formed before 55 Ma have low $\epsilon Nd\left(t\right)$ values and high ⁸⁷Sr/⁸⁶Sr values. The $\epsilon Nd\left(t\right)/\left({}^{87}Sr/{}^{86}Sr\right)$ values are mostly less than 0 (Figure 10E). The igneous samples that formed after 45 Ma have high $\epsilon Nd\left(t\right)$ values and low ⁸⁷Sr/⁸⁶Sr values, and the $\epsilon Nd\left(t\right)/\left({}^{87}Sr/{}^{86}Sr\right)$ ratios are mostly greater than 0 (Figure 10E). The difference in the geochemical properties of the igneous rocks in these two periods may be due to the change in the source of molten mantle material caused by the transformation of the subducting plate from the Izanagi plate to the Pacific plate. The quiet igneous period between 55 and 45 Ma in NE Asia can be compared to the disappearance of mantle melting above the subducting slab in our modeled case RS01 (Figures 3A–D).

FIGURE 10. Reconstructed model of the NE Eurasian plate and the West American plate and the $\epsilon Nd\left(t\right)$ and $\epsilon Nd\left(t\right)/\left({}^{87}Sr/{}^{86}Sr\right)$ versus ages for the NE Eurasian igneous rocks. (A,B) show the distribution maps of igneous rocks in northeastern Eurasia during the Late Cretaceous and Paleogene, respectively. The red dots show the locations of igneous rocks. The blue squares show the Hidaka belt where MORB-type basaltic rocks are dominant (Nanayama et al., 2019). The black lines with triangles are the subduction zones at that time. The red line is the MOR between the Izanagi plate and the Pacific plate. (C) Locations of the volcanic rocks after the Oligocene and the Sierra Madre Occidental silicic large igneous province (Ferrari et al., 2018). Plate reconstruction topologies are derived from Müller et al. (2016) and Vaes et al. (2019). (D,E) show $\epsilon Nd\left(t\right)$ and $\epsilon Nd\left(t\right)/\left({}^{87}Sr/{}^{86}Sr\right)$ values versus age for the northeast Eurasian igneous rocks from 100 to 20 Ma. The gray zone is a gap of continuous igneous events. Different color symbols [in (D,E)] represent igneous rock samples in different regions. Igneous rock data are from Wu and Wu (2019) and Tang et al. (2018).

On the other side of the Pacific plate, East Pacific ridge subduction began at 30 Ma and has continued to the present. The MOR has not yet fully subducted into the North American plate. We can observe that the East Pacific ridge is approximately parallel to the subduction zone (Figure 10C). In this period, the Sierra Madre Occidental silicic large igneous province with crustal melting and the Comondú Group of Baja California were formed, and resulting mafic magmas are widely distributed in western North America (Ferrari et al., 2018). However, the age of the subduction-related igneous rocks since 24 Ma can be divided into two stages (24–12 Ma and 12 Ma–present) (Figure 10C). There is also a short gap in the subduction-related igneous rocks at approximately 12 Ma (Martı́n et al., 2000; Pallares et al., 2007).

The magmatic activity gap in East Asia is longer than that in western North America (Figures 10A–C). We think it probably has two factors: the dip angle of the slab and the MOR subduction pattern, whether it is monolithic or segmented. The Western Pacific plate subducted under the East Asian plate at a high dip angle, caused strong backarc expansion, and formed a series of extensional basins on the edge of East Asia (e.g., Japan sea basin and Philippine basin). It maked the East Asian continent far away from the trench. The igneous rocks produced by oceanic plate subduction need a long time to migrate to the East Asian continent. As for the Farallon plate, it subducted under the North American plate with a low dip angle, and did not form extensional basins. In addition, the Izanagi-Pacific MOR is roughly parallel to the edge of the East Asian plate. The Izanagi-Pacific MOR subduction as a monolithic in a short period of time (within a few million years), while there was a small angle between the Farallon-Pacific MOR and the North American plate boundary. So the Farallon-Pacific MOR gradually subducted in segments. And the MOR subduction lasted for a long time (more than ten million years). This situation caused the lateral magmatic flow and resulted in a short igneous events gap.

Plate Coupling and Thickening in the Extinction Mode

From the results of the extinction mode, when the speed of the foregoing oceanic plate is much smaller than that of the following oceanic plate, the latter may be coupled with the foregoing oceanic plate. The extinction of the original MOR forms a new plate with a thickened coupling region, and the thickness of the coupling region can reach 150 km (Figures 7C,D). Therefore, the coupling region can be clearly distinguished from the normal thickness of the oceanic plate by seismic methods. Since the geological sequence of the coupling region is crust-lithosphere-crust-lithosphere (Figure 7C), the seismic waveform of this region can be relatively complicated. Moreover, the high temperature and high pressure applied to the coupled crustal region may generate metamorphic rocks and ophiolites, e.g., the Palawan ophiolite (Keenan et al., 2016). Due to the extinction of the active MOR, the following oceanic plate follows the foregoing oceanic plate to form continuous subduction. Therefore, a melting gap will not appear in the mantle wedge (Figures 7B–D). Continuous oceanic plate subduction causes continuous subsidence of the overriding continental plate topography (Figure 5).

The Shape and Mode of the MOR Subduction

Based on our results, the velocity difference $\Delta v$ between the foregoing and following oceanic plates determines the evolutionary modes and the shape of the MOR during its subduction. When $\Delta v$ is greater than 0, the velocity of the foregoing oceanic plate is much faster than that of the following oceanic plate, and the MOR region expands rapidly with a wide melting area (Figure 3B). When $\Delta v$ is close to 0, the foregoing and following oceanic plates move synchronously, the MOR area expands slowly, and the melting area appears to be normal and triangular. The center of the MOR is located in the center of the melting area (Figures 6A,B). However, when the MOR area is close to the trench, we find that the center of the MOR moves forward toward the left oceanic plate. This result occurs because when the foregoing oceanic plate is pulled by the slab and begins to subduct, due to the distance from the MOR, the age of the subducted foregoing oceanic plate gradually becomes younger and the density decreases as subduction continues. Because of continuous subduction, the length of the slab entering the mantle becomes longer, and the length of the unsubducted part of the foregoing oceanic plate becomes shorter. Therefore, the subduction velocity of the foregoing oceanic plate increases rapidly after a period of subduction (as shown in Figure 2A at 10–15 Myrs), resulting in the rapid forward movement of the center of the MOR, which is no longer at the center of the melting region. When $\Delta v$ is less than 0, the speed of the following oceanic plate is faster than that of the foregoing oceanic plate, and the MOR region gradually disappears (Figure 7B). Then, the foregoing and the following oceanic plates collide and merge into one plate with a thickened coupling region.

Nevertheless, previous studies have shown that many of the MORs are not completely parallel to the trench (Van Hunen et al., 2002; Gutiérrez et al., 2005; Rosenbaum et al., 2005; Espurt et al., 2008; Seton et al., 2015; Hu et al., 2016; Müller et al., 2016; Tang et al., 2018; Liu et al., 2020). The nonparallel subduction of MORs or asymmetric ridges may show different temporal and spatial distributions of igneous events and topographic evolution on the overriding plate. For these situations, our 2-D model has limitations, and 3-D models are required in future studies.

Conclusion

We construct a 2-D viscoelastoplastic model to study the modes and key parameters controlling the trench-parallel subduction of the mid-ocean ridge. Our model results show that the subduction modes of mid-ocean ridges can be primarily categorized into three types: the fast spreading mode, the slow spreading mode, and the extinction mode. The key factor controlling these subduction modes is the relative motion between the foregoing and following oceanic plates, which are separated by the mid-ocean ridge. The viscosity structure of the mantle significantly affects the sinking velocity of the preceding slab; therefore, the mantle viscosity structure may greatly affect the mode of trench-parallel subduction of the MOR. The volcanic gap and the topographic evolution observed in our model may successfully explain the geological observations for the subduction of the Izanagi-Pacific ridge at 55–45 Ma and the Pacific-Farallon ridge from 30 Ma to the present.

Data Availability Statement

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

Author Contributions

XS and WL contributed to conception and design of the study. XS and WL designed the model. XS ran the model and analyzed the results. XS 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 (41820104004, 41774105, 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.

Acknowledgments

We thank Jie Tang for providing igneous rock data in northeast Eurasia.

References