Plume–ridge interactions: ridgeward versus plate-drag plume flow

The analysis of mid-ocean ridges and hotspots that are sourced by deep-rooted mantle plumes allows us to get a glimpse of mantle structure and dynamics. Dynamical interaction between ridge and plume processes have been widely proposed and studied, particularly in terms of ridgeward plume flow. However, the effects of plate drag on plume–lithosphere and plume–ridge interaction remain poorly understood. In particular, the mechanisms that control plume flow towards vs. away from the ridge have not yet been systematically studied. Here, we use 2D thermomechanical numerical models of plume–ridge interaction to systematically explore the effects of (i) ridge-spreading rate, (ii) initial plume head radius and (iii) plume–ridge distance. Our numerical experiments suggest two different geodynamic regimes: (1) plume flow towards the ridge is favored by strong buoyant mantle plumes, slow spreading rates and small plume–ridge distances; (2) plume drag away from the ridge is in turn promoted by fast ridge spreading for small-to-intermediate plumes and large plume–ridge distances. We find that the pressure gradient between the buoyant plume and spreading ridge at first drives ridgeward flow, but eventually the competition between plate drag and the gravitational force of plume flow along the base of the sloping lithosphere controls the fate of plume (spreading towards vs. away from the ridge). Our results highlight that fast-spreading ridges exert strong plate-dragging force, which sheds new light on natural observations of largely absent plume–lithosphere interaction along fast-spreading ridges, such as the East Pacific Rise.

1 Introduction

Mid-ocean ridges (MORs) and hotspots are two main regions for deep-material recycling to the surface of the Earth. However, these two units are not always isolated; rather, they show strong interactions in some cases, termed as plume–ridge interaction (Morgan, 1978). Of up to 50 mantle plumes revealed by seismic tomography (French and Romanowicz, 2015; Montelli et al., 2004), more than 20 plumes are found to be associated with nearby ridges (Fig. 1a; Ito et al., 2003). Plume–ridge interaction is manifested by geophysical and geochemical anomalies along the ridge axis, e.g., high mantle potential temperature (Dalton et al., 2014), enriched radiogenic isotope anomalies (Cushman et al., 2004; Douglass et al., 1999; Yang et al., 2017) and adjacent lineations of seamounts (Ballmer et al., 2013b; Geissler et al., 2020; Lénat et al., 2009). Furthermore, plumes may promote migration of MOR spreading centers (Müller et al., 1998; Mittelstaedt et al., 2008, 2011; Whittaker et al., 2015), as evidenced by successive ridge jumps towards mantle plumes, e.g., at Iceland, Amsterdam–Saint Paul and Galapagos hotspots (Hardarson et al., 1997; Maia et al., 2011; Mittelstaedt et al., 2012). The interaction dynamics of a ridge with a ridge-centered and off-ridge plume has been widely studied and modeled in analogue and numerical experiments, revealing that the major controlling factors involve the ridge-spreading rate, plume buoyancy flux and their spatial distance (François et al., 2018; Ito et al., 1997; Kincaid et al., 1996; Ribe et al., 1995; Ribe, 1996; Sleep, 1997). Indeed, most plume–ridge interaction systems are associated with slow-spreading ridges, small mantle plumes and short plume–ridge distances (Fig. 1b). However, numerical studies systematically investigating the effects of these parameters on plume–ridge interaction and quantifying the controlling forces remain scarce.

As has been noted previously, buoyant plumes tend to spread ridgeward along the sloping base of the lithosphere (Morgan, 1978; Schilling, 1991; Small, 1995). Regions of divergent mantle flow beneath MORs represent the lowest dynamic-pressure regions in the oceanic asthenosphere and thus tend to suck ambient asthenospheric and plume materials towards the spreading center (Niu and Hekinian, 2004). On the other hand, the viscous drag at the base of the plate tends to convey the spreading plume material away from the MOR (Ribe and Christensen, 1994, 1999). Indeed, plume spreading at the base of the lithosphere is governed by the competition of trenchward viscous plate drag vs. ridgeward gravitational and pressure-driven forces (Kincaid et al., 1996). These gravitational and tectonic forces compete with each other to control the regime of plume–ridge interaction, but their balance remains to be quantified.

The different distribution of hotspots with classified plume–ridge interaction (ridgeward spreading) vs. no interaction (plate-drag spreading) also still remains enigmatic. Plume–ridge interaction is much more common near the Mid-Atlantic Ridge (MAR) than near the East Pacific Rise (EPR; Fig. 1a). Near the EPR, only the Pukapuka and Sojourn ridges display clear evidence of ridgeward flow of the magmatic source, but these volcanic ridges have been attributed to horizontally viscous differences or small-scale convection in the uppermost mantle and not to a mantle plume (Ballmer et al., 2013b; Clouard and Bonneville, 2005; Harmon et al., 2011). A previous study (Jellinek et al., 2003) proposed that fast-spreading ridges guide upwelling mantle flow towards the spreading center to convey the surrounding plumes from deep depths entirely into the MOR melting zone (Fig. 1c), resulting in the absence of hotspots adjacent to the EPR (see also Rowley et al., 2016; Rowley and Forte, 2022). However, fast plate spreading also tends to drag mantle plumes away from the MOR (Kincaid et al., 1995, 1996), leading to the typically parabolic shapes of hotspot swells, such as near Hawaii (Ribe and Christensen, 1994). Whether the increased spreading rates in the Pacific vs. Atlantic promote ridgeward vs. plate-drag plume flow remains an intriguing question.

The principal goal of this study is to investigate the process of plume–ridge interaction, with an emphasis on the effects of model parameters on the ridgeward vs. plate-drag plume spreading. We explore the effects of various model parameters, such as the size of the plume, ridge-spreading rate, and plume–ridge distance. Finally, we use our model results to interpret the differences between natural plume–ridge interaction systems in different oceans, particularly the striking difference between the East Pacific and Atlantic in this regard.

Figure 1 Global plume–ridge interaction systems. (a) Global distribution of mid-ocean ridges and mantle plumes. Residual bathymetry of the ocean basins comes from Straume et al. (2019). Mid-ocean ridges are painted in solid colored lines corresponding to half-spreading rate. Plumes not interacting with a ridge are shown by green circles, and hotspots linked to ridges are in red dots (Ito et al., 2003); size refers to the plume buoyancy flux from Hoggard et al. (2020). Black lines denote the regions of two LLSVPs (large low-shear velocity provinces) under South Africa and the Pacific Ocean (Torsvik et al., 2006). (b) Histograms of influential factors of plume–ridge interaction systems. Half the spreading rate and plume–ridge distance are taken from GPlates (Müller et al., 2016; Whittaker et al., 2015). Plume–ridge interaction systems link to slow-spreading ridges, small mantle plumes and short plume–ridge distances. (c) Sketches of ridgeward (top panel) and plate-drag plume flow (bottom panel) mode proposed, respectively.

Figure 2 Model setup. (a) Initial composition and boundary conditions. The oceanic plate consists of half-space cooling part and the thermal equilibrium part. A 50 Myr old mid-ocean ridge is set in the middle of the model based on half-space cooling-temperature structure. A thermal and chemical anormal mantle plume is located at 660 km. Different colors indicate the initial rock types and corresponding newly formed molten-rock types. Yellow arrows are the half-spreading rates imposed internally in the lithosphere (i.e., from 20 to 120 km in depth) to simulate ridge spreading. (b) Initial tested ridge and plume configurations. (c) Initial tested plume–ridge distances.

2 Numerical modeling

2.1 Modeling methods

To explore plume–lithosphere and plume–ridge interaction, we conduct numerical simulations utilizing the 2D thermo-mechanical code I2VIS, which is based on the staggered finite-difference method combined with marker-in-cell techniques (Gerya and Yuen, 2003, 2007). The equations of conservation of mass, momentum and energy (Eqs. 1–3, respectively) are solved in a fully staggered grid assuming an incompressible media:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{v}=\mathrm{0},\text{(2)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}{}^{\prime }}{\partial x_j}}-{\displaystyle \frac{\partial P}{\partial x_i}}+\mathit{\rho }{\mathit{g}}_i=\mathrm{0},\text{(3)}& {\displaystyle }\mathit{\rho }{C}_p\left({\displaystyle \frac{DT}{Dt}}\right)=-\mathrm{\nabla }\cdot \mathit{q}+H_{\mathrm{r}}+H_{\mathrm{a}}+H_{\mathrm{s}}+H_{\mathrm{l}},\end{array}$$

where v refers to the velocity, σ_ij^′ refers to the deviatoric stress tensor, P refers to the pressure, ρ refers to the density, g refers to the gravity acceleration, $\frac{D}{Dt}$ refers to the Lagrangian time derivative, C_p refers to the heat capacity, and q refers to the heat flux. Additionally, H_r, H_a, H_s and H_l are the radioactive, adiabatic, shear and latent heat productions, respectively.

We employ a non-Newtonian visco-plastic rheology (Gerya and Yuen, 2007) in the models. The viscous rheology depends on stress, temperature and pressure. The appropriate viscosity is expressed as that of a composite diffusion and dislocation-creep material (Eq. 4):

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{vis}}}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{diff}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{disl}}}},\end{array}$$

in which η_diff and η_disl are the diffusion- and dislocation-creep viscosity, respectively, and can be further computed using Eqs. (5) and (6):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{diff}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}A{\mathit{\sigma }}_{\mathrm{crit}}^{\mathrm{1}-n}\exp \left({\displaystyle \frac{P{V}_{\mathrm{a}}+E_{\mathrm{a}}}{RT}}\right),\text{(6)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{disl}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{\frac{\mathrm{1}}{n}}{\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{\frac{\mathrm{1}-n}{n}}\exp \left({\displaystyle \frac{P{V}_{\mathrm{a}}+E_{\mathrm{a}}}{nRT}}\right),\end{array}$$

where P is the pressure, T is the temperature; ${\dot{\mathit{\epsilon }}}_{\mathrm{II}}$ is the second invariant of the strain rate tensor; σ_crit is the diffusion- and/or dislocation-creep transition stress; and A, E_a, V_a and n are the strain rate pre-exponential factor, activation energy, activation volume and stress exponent, respectively. The plastic behavior η_pla is described by the Drucker–Prager yield criterion (Byerlee, 1978; Ranalli, 1995) according to Eqs. (7) and (8):

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\sigma }}_y=C+P\mathit{\phi },\text{(8)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{pla}}={\displaystyle \frac{{\mathit{\sigma }}_y}{\mathrm{2}{\dot{\mathit{\epsilon }}}_{\mathrm{II}}}},\end{array}$$

in which σ_y is the yield stress, C is the rock cohesion, and φ is the effective friction coefficient. The effective viscosity η_eff of rocks is thus constrained by both viscous and plastic deformation, where the rheological behavior depends on the minimum viscosity attained between ductile and brittle fields:

$$\begin{array}{}\text{(9)}& {\mathit{\eta }}_{\mathrm{eff}}=\min \left({\mathit{\eta }}_{\mathrm{vis}}{\mathit{\eta }}_{\mathrm{pla}}\right).\end{array}$$

Partial melting, melt extraction and percolation are also considered in the model in a simplified way (Gerya, 2013). The melt fraction (M₀) of the crust is assumed to increase with temperature and is calculated according to Eq. (10):

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{M}_{\mathrm{0}}& =\mathrm{0}\text{when}T\le T_{\mathrm{solidus}}\\ M_{\mathrm{0}}& ={\displaystyle \frac{(T-T_{\mathrm{solidus}})}{(T_{\mathrm{liquidus}}-T_{\mathrm{solidus}})}}\\ & \text{when}{T}_{\mathrm{solidus}}<T<T_{\mathrm{liquidus}}\\ M_{\mathrm{0}}& =\mathrm{1}\text{when}T\ge T_{\mathrm{liquidus}},\end{array}\end{array}$$

where T_solidus and T_liquidus are the solidus and liquidus temperature of different rock types, respectively, taken from Katz et al. (2003).

In our model, melt extraction is modeled indirectly and is considered to be an instantaneous process (Gerya et al., 2015). The extracted melt is assumed to move vertically from the molten source and then gets added to the bottom of the crust. Partial melt is extracted from the mantle and instantaneously displaced to the bottom of the crust and converted into hot mafic magma, obeying the conservation of material. The amount of extracted melt during the evolution of each experiment is traced by the Lagrangian markers (Gerya, 2013). The total amount of melt, M, for every marker excludes the amount of previously extracted melt according to Eq. (11):

$$\begin{array}{}\text{(11)}& M=M_{\mathrm{0}}-{\mathrm{\Sigma }}_n{M}_{\mathrm{ext}},\end{array}$$

where Σ_nM_ext refers to the total melt fraction extracted during the previous n melt extraction time steps.

The effective density of mafic magma and molten crust depends on its melt fraction and is calculated as follows (Gerya et al., 2015; Gülcher et al., 2020):

$$\begin{array}{}\text{(12)}& {\mathit{\rho }}_{\mathrm{eff}}={\mathit{\rho }}_{\mathrm{solid}}\left(\mathrm{1}-M+M{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0},\mathrm{molten}}}{{\mathit{\rho }}_{\mathrm{0},\mathrm{solid}}}}\right),\end{array}$$

where ρ_0,molten and ρ_0,solid are the reference densities of the molten and solid crust. ρ_solid is the crust density at a given pressure and temperature, which can be computed as follows:

$$\begin{array}{}\text{(13)}& {\mathit{\rho }}_{\mathrm{solid}}={\mathit{\rho }}_{\mathrm{0},\mathrm{solid}}\left[\mathrm{1}-\mathit{\alpha }\left(T-\mathrm{298}\right)\right]\left[\mathrm{1}+\mathit{\beta }\left(P-\mathrm{0.1}\right)\right],\end{array}$$

with thermal expansion $\mathit{\alpha }=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ K⁻¹ and compressibility $\mathit{\beta }={\mathrm{10}}^{-\mathrm{11}}$ Pa⁻¹.

Surface processes, such as erosion and sedimentation, are considered by solving the transport equation on the Eulerian nodes at each time step (Gerya and Yuen, 2003). Our erosion–sedimentation model uses gross-scale erosion–sedimentation rates which are independent of local elevation and topography (Burov and Cloetingh, 1997). We use constant and moderate rates of erosion (0.315 mm yr⁻¹) and sedimentation (0.0315 mm yr⁻¹), respectively, which fall within naturally observed ranges.

2.2 Model setup

The size of the model box is 6600×1200 km, with a nonuniform grid of 501×301 computational nodes in length and depth, respectively (Fig. 2). The densest grid is located in the center of the model domain (i.e., grid size decreases linearly from 20 km at the edges to 2 km at the ridge axis and increases linearly from 1 km at the top to 7 km at the bottom), where plume–ridge interaction would happen. The model consists of a 20 km thick sticky-air layer to accommodate crustal surface deformation. To reproduce the oceanic lithosphere, we choose a typical layered model, where the crust is composed of a water level (2 km), a sediment layer (1.5 km) and a basalt layer (7.5 km). The oceanic lithosphere and asthenosphere in the model are both modeled as dry olivine (the different colors for the mantle lithosphere and asthenosphere in the figures of this paper are only for better visualization). Besides, a 50 Myr old mid-ocean ridge is set on the central part of the lithosphere, splitting the model domain into two parts. At the depth of 660 km, a 200 km wide semicircular plume is located on the left of the model domain, corresponding to the onset of plume–ridge interaction from the mantle transition zone. Detailed rock parameters are listed in Table 1.

The thermal conditions at the top and bottom boundaries are fixed at 273 and 2513 K, respectively. The left and right boundaries are both insulating, with no external heat flow across them. The initial temperature structure of the mantle is adiabatic (0.5 K km⁻¹), which results in a temperature at 660 km depth of 1843 K. The initial temperature structure of the oceanic plate consists of a half-space cooling part and a thermal-equilibrium part (Fig. 2a). The half-space cooling model is used to describe the oceanic plate younger than 50 Myr, and the thermal-equilibrium structure is used to describe older oceanic parts. In other words, the thermal age of the lithosphere far away from the ridge is fixed at 50 Myr, with a constant plate thickness (i.e., ∼100 km). The hot plume is set a circular thermal and compositional (see Table 1) anomaly with an excess temperature of 250 K to trigger a rising thermochemical plume. All the velocity boundaries are free-slip boundaries. Additional internal boundary velocities are imposed at 500 km from each side boundary in the lithosphere to maintain the imposed half-spreading rate (Fig. 2a).

Table 1 Rock physical properties used in the numerical models.

^* 1 – Ranalli (1995); 2 – Turcotte and Schubert (2014). Other physical parameters used for all rocks include gas constant R=8.314 J K⁻¹ mol⁻¹, thermal expansion $\mathit{\alpha }=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ K⁻¹, compressibility $\mathit{\beta }=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{11}}$ Pa⁻¹ and heat capacity C_p=1000 J kg⁻¹ K⁻¹.

3 Model results

We conduct a series of numerical experiments to investigate ridge suction versus plate-drag acts on plumes. The effects of three major model parameters (i.e., the spreading rate of mid-ocean ridge, the initial plume head radius and the plume–ridge distance) are systematically studied. The typical dynamic evolution of models with ridgeward vs. plate-drag plume flow are demonstrated.

3.1 Model evolution with ridgeward plume flow

For models with dominant ridgeward flow, the typical model evolution is shown in Fig. 3 (the major model parameters used in this case are as follows: the half-spreading rate of 8 mm yr⁻¹, the initial plume head radius of 200 km and the off-axis distance of 800 km). In the early plume head stage, the buoyant mantle plume rises up rapidly in a mushroom-like shape and imposes dynamic stresses at the base of the overriding oceanic plate, leading to significant surface uplift (Fig. 3a, b). The ascending plume experiences extensive decompression melting at the base of the overriding plate and, due to the dynamic overpressure, spreads laterally, forming two branches that flow in opposite directions (Fig. 3c). A large amount of plume material is eventually entrained towards the spreading center, ponding underneath the ridge axis and significantly affecting the ridge dynamics. The entrainment of hot plume material promotes decompression melting (Fig. 3d, e) and increases the temperature beneath the ridge (Fig. S2 in the Supplement). Within the overlying lithosphere, the buoyant mantle plume leads to stress localization and strongly weakens the oceanic plate (Figs. S1, S3). As the plume eventually flows upward along the base of the gradually sloping plate near MOR, melting and crust production occurs (Fig. S1), forming an oceanic plateau of thickened crust. In addition to this gravitational force that guides plume material of the right branch ridgeward, plate spreading drags both branches in the opposite direction. Moreover, convective and tectonic stresses (“plume push” and “ridge suction”) affect both branches of the plume in a different way. As a consequence, the two branches evolve asymmetrically: the right branch that flows towards the ridge axis is more vigorous than the left branch, and the plume tail is also tilted towards the spreading center (Fig. 3c–e).

Figure 3 The evolution of the reference model M12 (see Table S1 in Supplement) with dominant ridgeward plume flow. The main model parameters employed in this case are as follows: half-spreading rate of 8 mm yr⁻¹, an initial plume head radius of 200 km and an off-axis distance of 800 km. (a) Surface topography over time along the flow path. (b–e) Snapshots of composition for the reference model (M12). (f) Profiles of the horizontal velocity component over time at the sections, as indicated (color-coded) in panel (c–e).

The mantle flow horizontal velocity profiles (Fig. 3f) further demonstrate the dominance of ridgeward plume flow, showing that plume flow is faster towards the spreading ridge than away from it. The velocity profiles elucidate dominant Poiseuille flow, with the maximum flow velocities in the middle of the asthenospheric channel. Such velocity profiles are well consistent with observations of seismic anisotropy at the Reunion plume (Barruol et al., 2019). The branches of the spreading plume head move significantly faster than the overriding plate. Because of the asymmetrical spreading of the plume head, the buoyancy flux carried by the right branch of the plume (density anomaly multiplied by horizontal velocity from Fig. 3f) is also much larger than that carried by the left branch.

3.2 Model evolution with plate-drag plume flow

For models with dominant plume flow away from the ridge (“plate-drag flow”), the typical model evolution is shown in Fig. 4. The controlling parameters of the representative model shown in Fig. 4 are the same as for the model shown in Fig. 3, except for a smaller radius (100 km) and faster spreading ridge (half-spreading rate: 45 mm yr⁻¹). At first, the ascending plume head spreads out similarly to the case described above and interacts with the overriding oceanic lithosphere. The largest surface uplift is sustained just above the plume head (Fig. 4a), slightly differently from the previous model, in which the highest surface elevation is observed on both sides of the plume conduit (Fig. 3a). Related to this spreading and uplift, the mantle plume changes the stress state of the overlying oceanic plate (Fig. S4), resulting in a slight weakening (Fig. S6). The plume head undergoes significant decompression melting near the deflection point (Fig. 4c). However, it is difficult for the small plume to extract molten plume due to the low melt flux (Fig. S4). As the plume cools, partially molten plume gets solidified speedily (Figs. 4d–e and S5). In contrast to the reference model from Sect. 3.1, this model displays most plume material flowing away from the ridge, likely due to dominant plate drag (Fig. 4c–e). Indeed, the left branch of the plume consistently displays larger buoyancy fluxes and maximum velocities than the right side over time (Fig. 4f).

Figure 4 Same as Fig. 3 for case M77 (i.e., the reference model for the plate-drag plume flow regime). The main model parameters employed in this case are as follows: half-spreading rate of 45 mm yr⁻¹, an initial plume head radius of 100 km and an off-axis distance of 800 km.

The underlying mechanism for dominant plate-drag plume flow is the frictional shear force of the moving plate, which is further demonstrated by the plume flow velocity profiles (Fig. 4f). In the early plume head stage (∼1.08 Myr), the plume spreads out faster than the plate velocity, which is primarily driven by the overpressure of the ponding plume head at this stage. After a short amount of time (∼2.32 Myr), however, plume spreading becomes significantly slower than plate velocity; hence, plate drag drives and controls the plume flow. Indeed, the flow mode in the asthenosphere rapidly shifts from Poiseuille flow (i.e., active plume flow) to Couette flow (i.e., passive plume flow; Fig. 4f), indicating the increasing role of plate drag in plume flow, soon after an initial plume head spreading.

Figure 5 Comparison between models with ridgeward vs. plate-drag plume flow. (a) Ridgeward flow with downwelling beneath the MOR (results from case M12, as in Fig. 3). Dashed white lines are streamlines; black arrows visualize the flow field. A schematic of flow is in the sub-panel on the right-hand side. (b) Plate-drag flow with upwelling mantle corner flow beneath the MOR (results from case M77, as in Fig. 4). (c) The dynamic-pressure and gravitational gradient of the plume marker (i.e., green circle in a) over time. (d) The dynamic-pressure and gravitational gradient of the plume marker (i.e., green circle in b) over time. The yellow box marks the location for the computation of average dynamic pressure at the ridge, needed for the calculation of the dynamic-pressure gradient (see text). The (c) and (d) represent the results of ridgeward and plate-drag flow, respectively.

Figure 6 Models varying initial plume head radii (model M53, M58 and M63; Table S1) shown by buoyancy flux and viscosity. (a–c) Buoyancy flux in spreading plume branches over time. Green and red triangles are markers used for buoyancy flux calculations for the left and right plume branches, respectively. U and R stand for half of the spreading rate and plume radius. (d–f) Viscosity snapshots of models with different plume head radii. Models with green circles represent plate-drag plume flow, and ridgeward plume flow is represented in red.

3.3 Two modes of plume–lithosphere interaction

The dominant ridgeward and dominant plate-drag plume flow regimes are two distinct modes of plume–plate interaction. The differences between these two regimes are further demonstrated in terms of mantle flow (Fig. 5a, b) and driving forces (Fig. 5c, d).

In the ridgeward-dominated models, a clockwise mantle develops from the plume to the spreading ridge (Fig. 5a). Molten plume material flows to the spreading ridge and occupies the space underneath the ridge axis, sustaining significant asymmetry of mid-ocean ridge melting (Conder et al., 2002). As a consequence of the continuous supply of the plume material, downward mantle flow forms beneath the ridge axis. This flow pattern dramatically differs from that shown in the plate-drag-dominated models, which show upward mantle flow underneath the ridge axis (Fig. 5b), as is typical for the flow beneath a MOR without the influence of a plume.

Figure 7 Models varying plume–ridge distances (model M57–M59; Table S1) shown by buoyancy flux and viscosity. (a–c) Buoyancy flux in spreading plume branches over time. Green and red triangles are markers used for buoyancy flux calculations for the left and right plume branches, respectively. U and D stand for half of the spreading rate and plume–ridge distance. (d–f) Viscosity snapshots of models with different plume–ridge distances. Models with green circles represent plate-drag plume flow, and the ridgeward plume flow is represented in red.

The distinct modes of plume–ridge interaction (ridgeward vs. plate-drag flow) are controlled by the competition of the tectonic (plate drag, ridge suction) and gravitational (plume buoyancy) driving forces. On the one hand, the moving plate drags sub-lithospheric plume material away from the ridge. On the other hand, the mechanism of ridgeward flow is twofold. First, the buoyant plume material flows along the sloping base of the lithosphere towards the shallow ridge along the gravitational gradient. Second, the plume is driven along the dynamic-pressure gradient from the pressure maximum (e.g., where the plume sustains dynamic topography) towards the pressure minimum beneath the diverging ridge. These gravitational (G_gv) and pressure-driven (G_dp) gradients are calculated by tracing plume markers (Fig. 5c, d) as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }{G}_{\mathrm{dp}}=(P_{\mathrm{mk}}-P_{\mathrm{r}})/L,\text{(15)}& {\displaystyle }{G}_{\mathrm{gv}}=({\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{mk}})\times \mathit{g}\times k,\end{array}$$

where P_mk is the dynamic pressure of the plume marker, and P_r is the averaged pressure in a 50 km box at the ridge center (Fig. 5b); L is the horizontal distance from the plume marker to the ridge axis; ρ_mk and ρ₀ are the plume marker density and initial density, respectively; g is the gravitational acceleration; k is the local slope of the base of the lithosphere.

In the early stage of model evolution, the plume head's dynamic overpressure is dominant, driving plume spreading in both directions (Fig. 5c, d), particularly in the direction of the low-pressure ridge. However, this pressure gradient systematically diminishes over time as the plume (head) spreads. Once the spreading plume approaches the ridge, the lithospheric slope increases. At some point, the gravitational gradient exceeds the dynamic-pressure gradient, taking over as the major driving force of guiding plume material towards the ridge. However, for low-flux plumes, the dynamic pressure and gravitation are not enough to support the flow to the ridge. Consequently, one of the essential conditions for plume–ridge interaction is that the plume must be able to reach the critical zone near the ridge, where the slope is sufficiently steep to take over for the ever-diminishing pressure gradient. This implies (1) that the plume buoyancy must overcome the shearing force of plate drag and (2) that the pressure gradient must be sustained long enough to reach the critical zone, in which the gravitational gradient can take over. It is noted that (1) shearing force scales with the rate of ridge spreading, and (2) the plume with high buoyancy flux is more likely to reach the critical region at a given plume–ridge distance.

3.4 Influence of model parameters

We have systematically investigated the effect of the three main model parameters (i.e., the spreading rate of the mid-ocean ridge, initial plume head radius and plume–ridge distance) on plume–ridge interaction. We explored half-spreading rates of the ridge of 8, 15, 30 and 45 mm yr⁻¹, corresponding to ultra-slow-, slow-, medium- and fast-spreading mid-ocean ridges, respectively (Gerya, 2012). We varied initial plume head radii in the range of 100 to 300 km. Further, we tested plume–ridge distance in the range of 600 to 1400 km.

3.4.1 Plume head radius

The size of the buoyant plume exerts an important control over plume–ridge interaction. Small plumes tend to be dragged away from the ridge, with typically larger lateral fluxes of the spreading plume in the left branch than in the right branch (Fig. 6a, b). The buoyancy flux, defined as the integral over a horizontal plane of the product of the vertical velocity and the density deficit within the plume, in each branch is calculated by multiplying the velocity of the markers in the plume pipe (Fig. 6d–f) by the density. The dynamic pressure decreases with decreasing plume size (Fig. S8a), and the pressures gradient is thus not strong enough for small plumes to reach the ridge. Plate shearing dominates plume flow soon after plume head spreading, and the moving plate then drags plume head material, leaving a tilted plume tail (Fig. 6d). In contrast, with larger initial plume head radius or buoyancy flux, the ponding plume spreads more vigorously (Fig. 6c) and sustains much higher overpressures at the base of the plate (Fig. S8a). This vigorous spreading can overcome plate drag to drive Poiseuille flow in both directions. Once the right plume branch approaches the spreading center, it is attracted and further accelerated by ridge suction. The plume tail is also markedly tilted towards the ridge axis due to asymmetric spreading (Fig. 6f). The larger the plume is, the more plume material gets entrained by the spreading center.

3.4.2 Plume–ridge distance

Plume–ridge distance also controls the regime of plume–ridge interaction. A plume at large distances spreads similarly to a plume at a small distance but is less likely to get affected by ridge suction (Fig. 7). The pressure gradient between the plume and ridge drives the ridgeward plume flow. However, the larger the plume–ridge distance, the smaller the pressure gradient would be (Fig. S8b), resulting in a lower buoyancy flux across the plume pipe (Fig. 7a–c). In the cases of distant plumes, the spreading of the plume head is strongly affected by plate drag (Fig. 7b, c).

3.4.3 Half-spreading rate of ridge

Another parameter that is worth investigating is the spreading rate of the ridge. The modeling results indicate that fast-spreading ridges promote plume flow away from the ridge due to the friction (Figs. 8 and 9a). With increasing spreading rate, the effect of plate shearing on plume–lithosphere interaction increases, as quantified by the spreading fraction. The spreading fraction γ (Eq. 14) is defined here as the ratio of ridgeward vs. plate-drag plume volume fluxes. We integrated the ridgeward plume volume flux (right branch), V_rw, and the plate-drag plume volume flux (left branch), V_tw. V_p is the total plume volume flux in the model. Ridgeward plume spreading is dominant for positive γ; plate-drag plume spreading is dominant for negative γ.

$$\begin{array}{}\text{(16)}& \mathit{\gamma }=(V_{\mathrm{rw}}-V_{\mathrm{tw}})/V_{\mathrm{p}}\end{array}$$

In the early stage (∼1 Myr), pressure-driven flow dominates in all models, and spreading fractions are positive, mainly driven by the expansion of the overpressured plume heads along the pressure gradient. After a certain time, the spreading fractions decrease dramatically with the decay of the mantle plume activity, representing the transition from the ridgeward to the plate-drag regime in some cases. The characteristic spreading fractions after an 8 Myr model time as a function of our model parameters are shown in Fig. 8. This compilation of our results reveals that the dominance of ridgeward flow decreases with increasing spreading rate and off-axis distance but significantly increases with plume size. For models with fast-spreading ridges, the parameter range of plate-drag flow-dominated models is expanded, indicating the critical role of plate drag in restricting ridgeward flow and plume–ridge interaction.

Figure 8 Parameter regime diagram of the contrasting modes of plume–ridge interaction. Spreading fractions γ (Eq. 14) at an ∼8 Myr model time. Each of the circles represents one of the numerical experiments, and sizes refer to γ. Circles in red and green represent models with dominant ridgeward plume flow and plate drag, respectively.

The transition from ridgeward (positive γ) to plate-drag (negative γ) flow in some of our cases is mainly determined by the competition between the effects of pressure-driven plume head spreading and plate shearing. The overpressure in the plume head drives plume materials towards the lower-pressure spreading center, while the moving plate shears plume away. Hence, we quantify the shear force of the overriding oceanic plate on the plume head using an integral approach:

$$\begin{array}{}\text{(17)}& {\mathbf{F}}_{\mathrm{s}}=\int {\mathit{\sigma }}_{xz}\mathrm{d}A.\end{array}$$

Equation (15) is employed to calculate the shear force, where F_s is the total shear force the spreading oceanic plate exerts on the uppermost part of the plume. σ_xz is the shear stress on each mantle plume grid cell, and A refers to the area of each grid cell. The pressure gradients, both gravitational and dynamic pressure, are calculated by tracing the plume markers according to Eqs. (12)–(13). As the plume material rises to the base of the lithosphere, the shear force exerted by the plate increases over time. We find that the integrated shear force between the spreading plate and the plume increases significantly as half-spreading rate increases (Fig. 9c).

Conversely, ridge-spreading rates control gravitational and pressure-driven plume driving forces (Fig. 9d). Increasing the spreading rate of the ridge implies a smaller dynamic-pressure gradient because the pressure gradient is related to the plate thickness difference at the ridge and plume, which is dependent on the spreading rate. A fast-spreading ridge also implies a smaller gravitational gradient because it leaves a more shallowly dipping lithospheric base. Thus, relatively strong plate shearing combined with relatively small pressure and gravitational gradients tend to advance plate-drag plume flow for high spreading rates.

Figure 9 Model results influenced by different half-spreading rates. (a) Effect of half-spreading rate and plume–ridge distance on ridgeward vs plate-drag. Viscosity snapshots are shown (models M7–M9, M82–M84; Table S1). The fast-spreading ridge promotes the dragging of plume material near the ridge. Models with green circles represent plate-drag plume flow, and ridgeward plume flow is represented in red. (b) Dynamic evolutions of ridgeward and plate-drag plume flow, revealed by defined ridge-spreading fraction (Eq. 14). (c) Shear force (F_s) between moving-plate and plume material under different spreading rates. (d) Pressure gradient between plume head and ridge center in different half-spreading rate models. The solid and dashed lines are the plume gravitation and dynamic-pressure gradient, respectively. R and D stand for half of the plume radius and plume–ridge distance.

4 Discussion

Natural observations show that there are only very few hotspots indicative of ridgeward plume flow close to the East Pacific Rise (EPR; Fig. 10a) in contrast to many such hotspots in the Atlantic and Indian oceans. A previous study (Jellinek et al., 2003) proposed that fast-spreading ridges such as the EPR efficiently convey any surrounding plumes into the spreading center from the deep mantle (Fig. 1c), which leads to fewer hotspots nearby fast-spreading ridges. However, based on our modeling results, fast-spreading ridges tend to promote plate-drag flow of the spreading plume material, providing an alternative explanation to the relative absence of hotspots along the EPR. We discuss the possibility of this potential explanation combined with geological and geophysical observations (Fig. 10).

Firstly, the plate-drag effect of fast-spreading ridges on plumes is evidenced by geophysical observations. We locate the positions of the mantle plumes at the core-mantle boundary (CMB) and the associated hotspots on the surface based on global seismic tomography (Jackson et al., 2021; Koppers et al., 2021). A lateral offset between the deep and surface positions of plumes is a common feature, indicating the deflection of plumes due to mantle flow. Specifically, a large portion of (7 of 14) plumes located in the Atlantic are tilted towards the mid-ocean ridge. However, only very few plumes (6 of 16) in the Pacific are tilted towards the mid-ocean ridge. Indeed, the majority of plumes are tilted away from the ridges, indicating the significant effect of plate drag on plumes beneath fast plates. Such observations are consistent with the predictions of our models with dominant plate-drag plume spreading.

Figure 10 A compilation of hotspots along with spreading ridges in the Atlantic and the Pacific. (a) Distribution of surface hotspots (circles) together with depth-projected source locations at CMB (blue dots) of the plumes based on Jackson et al. (2021). Plumes in magenta circles are mantle plumes interacting with ridges (Ito et al., 2003), and plumes not interacting with ridges are shown as green circles, whose size refers to the plume buoyancy flux (Hoggard et al., 2020). Blue arrows indicate the changes in the position of the plume at the top and bottom of the mantle. Yellow dots are MORB samples mapped in (b). (b) Plot of radioactive isotope ratios along ridge MORB samples. The data are downloaded from the PetDB database (http://portal.earthchem.org/, last access: 23 August 2021). The colored symbols refer to samples in different mid-ocean ridges. Main hotspots influencing MORBs are labeled with shaded bands. The dashed black lines are the mean MORB isotopes ratio from Gale (2013). Red and green lines are the mean ratios of the samples in the Mid-Atlantic Ridge and EPR, respectively.

Figure 11 Buoyancy flux, plate speed and plume–ridge distance of mantle plumes in different oceans. Mantle plumes in the Pacific, Atlantic and Indian oceans are shown in green, red and pink circles, respectively. Blue stars mark the ridge-interacted plumes according to Ito et al. (2003). (a) Plot of plume–ridge distance and plume buoyancy flux. Data are from Hoggard et al. (2020). (b) Plot of plume–ridge distance and plate speed at the location of plumes. Plume–ridge distance comes from GPlates (Müller et al., 2016; Whittaker et al., 2015), and plate speed data come from Becker et al. (2015).

Geochemical studies suggest that mantle plumes are enriched in light rare-earth elements (LREEs) and radiogenic isotopes of Sr and Pb but depleted in Nd isotopes. These geochemical anomalies are evident in MORB (mid-ocean-ridge basalt) at the sites of active plume–ridge interaction (Cushman et al., 2004; Douglass et al., 1999; Yang et al., 2017). We find that MORB sampled along both the Mid-Atlantic Ridge and the EPR indeed display geochemical anomalies (Fig. 10b), indicating ridgeward flow of plume material at specific locations. However, the Mid-Atlantic MORB dataset is slightly more heterogeneous than the East Pacific Rise in terms of geochemical isotopes. The EPR is basically characterized as normal oceanic basalt, along which only very few regions showing composition associated with nearby plumes. This contradicts the view (Jellinek et al., 2003) that mantle plumes are fully entrained into the central MOR melting zone at fast-spreading ridges.

Based on our modeling results, initial plume head radius and plume–ridge distance also control the mode of plume–ridge interaction. However, for different buoyancy fluxes B, there is only a small difference in the proportion of interacting plumes to non-interacting plumes. A small majority of major plumes (5 out of 8 with B>1.6 Mg s⁻¹) and a small minority of small-to-intermediate plumes (11 out of 25 for B<1.6 Mg s⁻¹) display interaction with the ridge (Fig. 11a). The underlying cause for this observation remains unclear but may be related to the distribution of large plumes globally, with many of them being located very far from MORs. Also note that our 2D models are limited in that plume material cannot spread in the out-of-plane direction, hence somewhat exaggerating the effects of buoyancy flux. In any case, the distribution of observed plume buoyancy fluxes (Hoggard et al., 2020) varies little across different oceans (Fig. 11a). Therefore, the effects of plume size are not a good candidate to explain the notable difference between the Atlantic and Pacific in terms of plume–ridge interaction mode.

On the other hand, compared with the Atlantic and Indian oceans, Pacific plumes are located significantly further from the mid-ocean ridge (Fig. 11b). Plume–ridge distances in the Pacific are mostly >2000 km, which exceeds the maximum plume–ridge interaction distance of 1400 km (Schilling, 1991). Most plumes in the Pacific exhibit the typical signatures of plume flow away from the ridge, such as parabolic swell shapes (e.g., Society, Marquesas and Hawaii plumes; Ballmer et al., 2013a; Ballmer et al., 2015; Cheng et al., 2015; Wolfe et al., 2009) and linear volcanic chains (Buff et al., 2021; Clouard and Bonneville, 2005; Jackson et al., 2010). Age-progressive hotspot trails indicate an absence of dominant ridgeward flow. By contrast, most plumes in the Atlantic have been close to the ridge since the opening of the ocean. These mantle plumes (e.g., Discovery, Iceland and Tristan–Gough plumes; O'Connor et al., 2012) have not moved much since the breakup of the Atlantic. One factor may be that the underlying plume-generation zone (i.e., the edge of the African LLSVP) round largely parallel to the Mid-Atlantic Ridge (Fig. 1) (Torsvik et al., 2006). In this case, plume–ridge distance may play a critical role in the plume–ridge interaction and could explain the striking difference between the Pacific and Atlantic in terms of the proportion of plume–ridge-interacting vs. non-interacting systems. In addition, the rapid movement of the Pacific plate tends to inhibit ridgeward plume flow at a given plume–ridge distance. The distribution of interacting (stars) vs non-interacting systems in Fig. 11b is almost exactly as predicted by our models for the coupled effects of plume–ridge distance and plate velocity. For example, we note that fast-spreading ridges can still interact with adjacent plumes under the appropriate conditions. In the case of very short plume–ridge distances, there is good evidence of plume–ridge interaction in the southern Pacific Ocean (e.g., Louisville plume; Conder et al., 2002; Toomey et al., 2002; Vlastélic and Dosso, 2005). Based on a series of numerical modeling, as well as geological and geophysical observations, we conclude that mantle plumes in the Pacific are more likely to spread away from the ridge and into the direction of plate motion than in the Atlantic and Indian oceans. However, the tendency of fast plate velocities to promote plume spreading away from the MOR through viscous drag may depend on the details of lithosphere–asthenosphere rheological coupling, such as the presence of a weak decoupling (e.g., melt) layer (Rychert et al., 2020). Further studies of plume spreading and plume–ridge interaction are needed to shed light on the coupling of the plate-mantle system.

5 Conclusions

In this study, we explore the evolution of plume–ridge interaction with 2D thermomechanical numerical models. Based on model results, we find the following:

1. Plume–ridge interaction is mainly governed by the competition between the effects of plume spreading (overpressure in the plume head stage), upward gravitationally driven flow of the plume along the base of the sloping lithosphere and plate shearing. These driving forces are controlled by plume size, plume–ridge distance and the spreading rate of the mid-ocean ridge.

2. MOR spreading does not only draw upwelling plumes into the spreading center but also tends to drag mantle plumes away from the ridge. Plume flow away from the ridge is favored by small and/or distant plumes, as well as by fast spreading rates, whereas plume flow towards the ridge is promoted by large and/or nearby plumes, as well as by slow spreading rates.

3. Considering the high plate velocity and typically large plume–ridge distances, mantle plumes in the Pacific are more likely to be dragged away from the EPR than they are to be drawn towards the ridge center.