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

Abstract

Trench-parallel subduction of mid-ocean ridges occurs frequently in plate motion history, such as along the western boundary of the Pacific plate in the early Cenozoic and along the eastern boundary of the Pacific plate at present. Such subduction may strongly alter the surface topography, volcanic activity and slab morphology in the mantle, whereas few studies have been conducted to investigate its evolutionary process. Here, we construct a 2-D viscoelastoplastic numerical model to study the modes and key parameters controlling trench-parallel subduction of mid-ocean ridges. 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 the following oceanic plates, which are separated by the mid-ocean ridge. Different subduction modes exert different surface geological expressions, which may explain specific evolutionary processes related to mid-ocean ridge subduction, such as topographic deformation and the eruption gap of volcanic rocks in East Asia within 55–45 Ma and in the western North American plate during the late Cenozoic.

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 is the strain rate, is the shear modulus, is the time variation rate of deviatoric stress, is the viscosity, and is the deviatoric stress. and represent spatial indices.

The viscosity equation we use is non-Newtonian and is temperature- and depth-dependent (Karato and Wu, 1993):where is the reference viscosity, is the second invariant of the deviatoric strain rate tensor, is the reference strain rate, is the strain exponent, is the activation energy, is the gas constant, is the absolute temperature (in Kelvin), is the reference temperature, is the pre-exponential factor of depth-dependent viscosity, is the model box thickness, and 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, . For the ViscT&D viscosity structure, .

The plasticity of the material is applied with the Drucker-Prager criterion through the yield stress:where is the yield stress, is the coefficient of friction, is the pressure and 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 and linearly decrease when the accumulated plastic strain increases:where is the initial coefficient of friction, is the accumulated plastic strain, is the reference plastic strain, is the initial cohesion and is the minimum cohesion (Table 1). The material no longer weakens when 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, , , and , 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

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 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 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 , and the density of the lithosphere is uniformly 3,300 . 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 (Turcotte and Schubert, 2002):where is the thermal diffusivity and 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 = 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 is calculated based on the temperature, pressure and water content:where T is the temperature, is the pressure, is the weight fraction of bulk water, is the weight fraction of water in liquid phase, is the solidus temperature, is the temperature decrease in the solidus at a water content and 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 between the two oceanic plates as (the velocity direction is positive to the left):where is the velocity of the foregoing oceanic plate and is the velocity of the following oceanic plate. indicates that the foregoing and the following oceanic plates have the same speed, whereas or 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 is a key factor influencing the mode of trench-parallel subduction of the MOR.

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 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 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 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 value increases to nearly 0, indicating that the foregoing and following oceanic plates merge and move at the same speed.

FIGURE 2

According to the results of the model, the main factor controlling the subduction of the following plate is the imposed velocity . Even if is a small value, such as = 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 can be reduced to 0 because of the negative density anomaly of the oceanic plateaus or microcontinent-arc collision. Therefore, when 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 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 , the initial dip angle of the slab 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

FIGURE 4

FIGURE 5

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

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 , the initial dip angle of the slab 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

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

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

On the other hand, we fix the slab length L to 600 km and vary the imposed velocity 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 effectively increase the magnitude of , making the trench-parallel subduction of MOR move to the fast spreading mode. In contrast, the decreased L and increased decrease the magnitude of , leading to the slow spreading mode of MOR subduction. When 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 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

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 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 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 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 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 values and high ⁸⁷Sr/⁸⁶Sr values. The values are mostly less than 0 (Figure 10E). The igneous samples that formed after 45 Ma have high values and low ⁸⁷Sr/⁸⁶Sr values, and the 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

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 between the foregoing and following oceanic plates determines the evolutionary modes and the shape of the MOR during its subduction. When 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 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 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.

References

• 1

  BuckW. R.PoliakovA. N. B. (1998). Abyssal hills Formed by Stretching Oceanic Lithosphere. Nature392 (6673), 272–275. 10.1038/32636

  • CrossRef
  • Google Scholar

• 2

  BurkettE. R.BillenM. I. (2009). Dynamics and Implications of Slab Detachment Due to ridge‐trench Collision. J. Geophys. Res. Solid Earth114, B12402. 10.1029/2009jb006402

  • CrossRef
  • Google Scholar

• 3

  ChenM.ShenX.LengW.ChenL. (2020). Destruction of Cratonic Lithosphere Induced by Oceanic Subduction Initiation. Geophys. Res. Lett.47 (15), e2020GL089140. 10.1029/2020GL089140

  • CrossRef
  • Google Scholar

• 4

  ChristensenU. R. (1996). The Influence of Trench Migration on Slab Penetration into the Lower Mantle. Earth Planet. Sci. Lett.140 (1-4), 27–39. 10.1016/0012-821x(96)00023-4

  • CrossRef
  • Google Scholar

• 5

  ČížkováH.BinaC. R. (2013). Effects of Mantle and Subduction-Interface Rheologies on Slab Stagnation and Trench Rollback. Earth Planet. Sci. Lett.379, 95–103. 10.1016/j.epsl.2013.08.011

  • CrossRef
  • Google Scholar

• 6

  ČížkováH.van den BergA. P.SpakmanW.MatyskaC. (2012). The Viscosity of Earth’s Lower Mantle Inferred from Sinking Speed of Subducted Lithosphere. Phys. earth Planet. Interiors200, 56–62. 10.1016/j.pepi.2012.02.010

  • CrossRef
  • Google Scholar

• 7

  EspurtN.FunicielloF.MartinodJ.GuillaumeB.RegardV.FaccennaC.et al (2008). Flat Subduction Dynamics and Deformation of the South American Plate: Insights from Analog Modeling. Tectonics27 (3), TC3011. 10.1029/2007tc002175

  • CrossRef
  • Google Scholar

• 8

  FerrariL.Orozco-EsquivelT.BryanS. E.López-MartínezM.Silva-FragosoA. (2018). Cenozoic Magmatism and Extension in Western Mexico: Linking the Sierra Madre Occidental Silicic Large Igneous Province and the Comondú Group with the Gulf of California Rift. Earth-Science Rev.183, 115–152. 10.1016/j.earscirev.2017.04.006

  • CrossRef
  • Google Scholar

• 9

  GarelF.GoesS.DaviesD. R.DaviesJ. H.KramerS. C.WilsonC. 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. 10.1002/2014gc005257

  • CrossRef
  • Google Scholar

• 10

  GeorgievaV.GallagherK.SobczykA.SobelE. R.SchildgenT. F.EhlersT. A.et al (2019). Effects of Slab-Window, Alkaline Volcanism, and Glaciation on Thermochronometer Cooling Histories, Patagonian Andes. Earth Planet. Sci. Lett.511, 164–176. 10.1016/j.epsl.2019.01.030

  • CrossRef
  • Google Scholar

• 11

  GroomeW. G.ThorkelsonD. J. (2009). The Three-Dimensional Thermo-Mechanical Signature of ridge Subduction and Slab Window Migration. Tectonophysics464 (1-4), 70–83. 10.1016/j.tecto.2008.07.003

  • CrossRef
  • Google Scholar

• 12

  GurnisM.HagerB. H. (1988). Controls of the Structure of Subducted Slabs. Nature335 (6188), 317–321. 10.1038/335317a0

  • CrossRef
  • Google Scholar

• 13

  GutiérrezF.GioncadaA.FerranO. G.LahsenA.MazzuoliR. (2005). The Hudson Volcano and Surrounding Monogenetic Centres (Chilean Patagonia): an Example of Volcanism Associated with ridge–trench Collision Environment. J. Volcanology Geothermal Res.145 (3-4), 207–233. 10.1016/j.jvolgeores.2005.01.014

  • CrossRef
  • Google Scholar

• 14

  HuJ.LiuL.HermosilloA.ZhouQ. (2016). Simulation of Late Cenozoic South American Flat-Slab Subduction Using Geodynamic Models with Data Assimilation. Earth Planet. Sci. Lett.438, 1–13. 10.1016/j.epsl.2016.01.011

  • CrossRef
  • Google Scholar

• 15

  IsozakiY.AokiK.NakamaT.YanaiS. (2010). New Insight into a Subduction-Related Orogen: A Reappraisal of the Geotectonic Framework and Evolution of the Japanese Islands. Gondwana Res.18 (1), 82–105. 10.1016/j.gr.2010.02.015

  • CrossRef
  • Google Scholar

• 16

  KaratoS.-i.WuP. (1993). Rheology of the Upper Mantle: A Synthesis. Science260 (5109), 771–778. 10.1126/science.260.5109.771

  • CrossRef
  • Google Scholar

• 17

  KatzR. F.SpiegelmanM.LangmuirC. H. (2003). A New Parameterization of Hydrous Mantle Melting. Geochem. Geophys. Geosystems4 (9), 1–9. 10.1029/2002gc000433

  • CrossRef
  • Google Scholar

• 18

  KeenanT. E.EncarnaciónJ.BuchwaldtR.FernandezD.MattinsonJ.RasoazanamparanyC.et al (2016). Rapid Conversion of an Oceanic Spreading center to a Subduction Zone Inferred from High-Precision Geochronology. Proc. Natl. Acad. Sci. USA113 (47), E7359–E7366. 10.1073/pnas.1609999113

  • CrossRef
  • Google Scholar

• 19

  KingS. D.FrostD. J.RubieD. C. (2015). Why Cold Slabs Stagnate in the Transition Zone. Geology43 (3), 231–234. 10.1130/g36320.1

  • CrossRef
  • Google Scholar

• 20

  LengW.GurnisM.AsimowP. (2012). From Basalts to Boninites: The Geodynamics of Volcanic Expression during Induced Subduction Initiation. Lithosphere4 (6), 511–523. 10.1130/l215.1

  • CrossRef
  • Google Scholar

• 21

  LengW.GurnisM. (2011). Dynamics of Subduction Initiation with Different Evolutionary Pathways. Geochem. Geophys. Geosystems12 (12), Q12018. 10.1029/2011gc003877

  • CrossRef
  • Google Scholar

• 22

  LengW.GurnisM. (2015). Subduction Initiation at Relic Arcs. Geophys. Res. Lett.42 (17), 7014–7021. 10.1002/2015gl064985

  • CrossRef
  • Google Scholar

• 23

  LiS. M.WangQ.ZhuD. C.CawoodP. A.SternR. J.WeinbergR.et al (2020). Reconciling Orogenic Drivers for the Evolution of the Bangong‐Nujiang Tethys during Middle‐Late Jurassic. Tectonics39 (2), e2019TC005951. 10.1029/2019TC005951

  • CrossRef
  • Google Scholar

• 24

  LingM.-X.WangF.-Y.DingX.HuY.-H.ZhouJ.-B.ZartmanR. E.et al (2009). Cretaceous ridge Subduction along the Lower Yangtze River belt, Eastern China. Econ. Geology.104 (2), 303–321. 10.2113/gsecongeo.104.2.303

  • CrossRef
  • Google Scholar

• 25

  LiuK.ZhangJ.XiaoW.WildeS. A.AlexandrovI. (2020). A Review of Magmatism and Deformation History along the NE Asian Margin from Ca. 95 to 30 Ma: Transition from the Izanagi to Pacific Plate Subduction in the Early Cenozoic. Earth-Science Rev.209, 103317. 10.1016/j.earscirev.2020.103317

  • CrossRef
  • Google Scholar

• 26

  MaedaJ. I.KagamiH. (1996). Interaction of a Spreading ridge and an Accretionary Prism: Implications from MORB Magmatism in the Hidaka Magmatic Zone, Hokkaido, Japan. Geol24 (1), 31–34. 10.1130/0091-7613(1996)024<0031:ioasra>2.3.co;2

  • CrossRef
  • Google Scholar

• 27

  MartinodJ.GuillaumeB.EspurtN.FaccennaC.FunicielloF.RegardV. (2013). Effect of Aseismic ridge Subduction on Slab Geometry and Overriding Plate Deformation: Insights from Analogue Modeling. Tectonophysics588, 39–55. 10.1016/j.tecto.2012.12.010

  • CrossRef
  • Google Scholar

• 28

  Martı́nA.FletcherJ. M.López-Martı́nezM.Mendoza-BorundaR. (2000). Waning Miocene Subduction and Arc Volcanism in Baja California: the San Luis Gonzaga Volcanic Field. Tectonophysics318 (1-4), 27–51.

  • Google Scholar

• 29

  MoresiL.DufourF.MühlhausH.-B. (2003). A Lagrangian Integration point Finite Element Method for Large Deformation Modeling of Viscoelastic Geomaterials. J. Comput. Phys.184 (2), 476–497. 10.1016/s0021-9991(02)00031-1

  • CrossRef
  • Google Scholar

• 30

  MuhuriS. K.DewersT. A.ScottT. E.JrRechesZ. E. (2003). Interseismic Fault Strengthening and Earthquake-Slip Instability: Friction or Cohesion?Geol31 (10), 881–884. 10.1130/g19601.1

  • CrossRef
  • Google Scholar

• 31

  MüllerR. D.SetonM.ZahirovicS.WilliamsS. E.MatthewsK. J.WrightN. M.et al (2016). Ocean basin Evolution and Global-Scale Plate Reorganization Events since Pangea Breakup. Annu. Rev. Earth Planet. Sci.44, 107–138. 10.1146/annurev-earth-060115-012211

  • CrossRef
  • Google Scholar

• 32

  NanayamaF.YamasakiT.IwanoH.DanharaT.HirataT. (2019). Zircon U-Pb Ages of Sedimentary Complexes in the Hidaka Belt. Jour. Geol. Soc. Jpn.125 (6), 421–438. 10.5575/geosoc.2019.0011

  • CrossRef
  • Google Scholar

• 33

  PallaresC.MauryR. C.BellonH.RoyerJ. Y.CalmusT.Aguillón-RoblesA.et al (2007). Slab-tearing Following ridge-trench Collision: Evidence from Miocene Volcanism in Baja California, México. J. Volcanology Geothermal Res.161 (1-2), 95–117. 10.1016/j.jvolgeores.2006.11.002

  • CrossRef
  • Google Scholar

• 34

  PoliakovA. N.BuckW. R. (1998). Mechanics of Stretching Elastic-Plastic-Viscous Layers: Applications to Slow-Spreading Mid-ocean Ridges. Geophys. Monogr. Ser.106, 305–324. 10.1029/GM106p0305

  • CrossRef
  • Google Scholar

• 35

  PütheC.GeryaT. (2014). Dependence of Mid-ocean ridge Morphology on Spreading Rate in Numerical 3-D Models. Gondwana Res.25 (1), 270–283. 10.1016/j.gr.2013.04.005

  • CrossRef
  • Google Scholar

• 36

  QingJ.LiaoJ.LiL.GaoR. (2021). Dynamic Evolution of Induced Subduction through the Inversion of Spreading Ridges. J. Geophys. Res. Solid Earth126 (3), e2020JB020965. 10.1029/2020JB020965

  • CrossRef
  • Google Scholar

• 37

  RosenbaumG.GilesD.SaxonM.BettsP. G.WeinbergR. F.DubozC. (2005). Subduction of the Nazca Ridge and the Inca Plateau: Insights into the Formation of Ore Deposits in Peru. Earth Planet. Sci. Lett.239 (1-2), 18–32. 10.1016/j.epsl.2005.08.003

  • CrossRef
  • Google Scholar

• 38

  RüpkeL. H.MorganJ. P.HortM.ConnollyJ. A. (2004). Serpentine and the Subduction Zone Water Cycle. Earth Planet. Sci. Lett.223 (1-2), 17–34. 10.1016/j.epsl.2004.04.018

  • CrossRef
  • Google Scholar

• 39

  SetonM.FlamentN.WhittakerJ.MüllerR. D.GurnisM.BowerD. J. (2015). Ridge Subduction Sparked Reorganization of the Pacific Plate‐mantle System 60-50 Million Years Ago. Geophys. Res. Lett.42 (6), 1732–1740. 10.1002/2015gl063057

  • CrossRef
  • Google Scholar

• 40

  SongY.RenJ.StepashkoA. A.LiJ. (2014). Post-rift Geodynamics of the Songliao Basin, NE China: Origin and Significance of T11 (Coniacian) Unconformity. Tectonophysics634, 1–18. 10.1016/j.tecto.2014.07.023

  • CrossRef
  • Google Scholar

• 41

  SteinbergerB.CalderwoodA. R. (2006). Models of Large-Scale Viscous Flow in the Earth's Mantle with Constraints from mineral Physics and Surface Observations. Geophys. J. Int.167 (3), 1461–1481. 10.1111/j.1365-246x.2006.03131.x

  • CrossRef
  • Google Scholar

• 42

  TangJ.XuW.WangF.GeW. (2018). Subduction History of the Paleo-Pacific Slab beneath Eurasian Continent: Mesozoic-Paleogene Magmatic Records in Northeast Asia. Sci. China Earth Sci.61 (5), 527–559. 10.1007/s11430-017-9174-1

  • CrossRef
  • Google Scholar

• 43

  ToriiY.YoshiokaS. (2007). Physical Conditions Producing Slab Stagnation: Constraints of the Clapeyron Slope, Mantle Viscosity, Trench Retreat, and Dip Angles. Tectonophysics445 (3-4), 200–209. 10.1016/j.tecto.2007.08.003

  • CrossRef
  • Google Scholar

• 44

  TurcotteD. L.SchubertG. (2002). Geodynamics. Cambridge: Cambridge University Press.

  • Google Scholar

• 45

  VaesB.HinsbergenD. J. J.BoschmanL. 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. Tectonics38 (4), 1367–1413. 10.1029/2018tc005164

  • CrossRef
  • Google Scholar

• 46

  van HunenJ.van den BergA. P.VlaarN. J. (2002). The Impact of the South-American Plate Motion and the Nazca Ridge Subduction on the Flat Subduction below South Peru. Geophys. Res. Lett.29 (14), 35–41. 10.1029/2001gl014004

  • CrossRef
  • Google Scholar

• 47

  WangC.DingW.LiJ.DongC.FangY.TangL.et al (2019). Effects of Trench‐perpendicular ridge Subduction on Accretionary Wedge Deformation: Clues from Analogue Modelling. Geol. J.54 (4), 2665–2678. 10.1002/gj.3317

  • CrossRef
  • Google Scholar

• 48

  WangC.FengZ.ZhangL.HuangY.CaoK.WangP.et al (2013). Cretaceous Paleogeography and Paleoclimate and the Setting of SKI Borehole Sites in Songliao Basin, Northeast China. Palaeogeogr. Palaeoclimatol. Palaeoecol.385, 17–30. 10.1016/j.palaeo.2012.01.030

  • CrossRef
  • Google Scholar

• 49

  WuJ. T.-J.WuJ. (2019). Izanagi-Pacific ridge Subduction Revealed by a 56 to 46 Ma Magmatic gap along the Northeast Asian Margin. Geology47 (10), 953–957. 10.1130/g46778.1

  • CrossRef
  • Google Scholar

• 50

  YangT.MoresiL.ZhaoD.SandifordD.WhittakerJ. (2018). Cenozoic Lithospheric Deformation in Northeast Asia and the Rapidly-Aging Pacific Plate. Earth Planet. Sci. Lett.492, 1–11. 10.1016/j.epsl.2018.03.057

  • CrossRef
  • Google Scholar

• 51

  ZhangB.LiuS.LinC.ShenW.LiX. (2020). Reconstruction of the Stress Regime in the Jiaolai Basin, East Asian Margin, as Decoded from Fault-Slip Analysis. J. Struct. Geology.141, 104190. 10.1016/j.jsg.2020.104190

  • CrossRef
  • Google Scholar