# Intraplate Volcanism Driven by Slab-Plume Interaction: Numerical Modeling and Its Application to the Eifel, Massif Central and Hainan Volcanic Areas

Yingying Li, Bernhard Steinberger, Sascha Brune, Eline Le Breton, Anne

Glerum, Michaël Pons

<figure><img src="image_1.png" /></figure>

## To cite this version:

Yingying Li, Bernhard Steinberger, Sascha Brune, Eline Le Breton, Anne Glerum, et al. Intraplate Volcanism Driven by Slab-Plume Interaction: Numerical Modeling and Its Application to the Eifel, Massif Central and Hainan Volcanic Areas. Journal of Geophysical Research: Solid Earth, 2026, 131 (1), pp.e2025JB032799. 10.1029/2025JB032799. insu-05461349

HAL Id: insu-05461349

https://insu.hal.science/insu-05461349v1

Submitted on 16 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2025JB032799

**Key Points:**

* Slabs can trigger secondary plumes from the base of the transition zone, where lower-mantle plumes stall

* Slab geometry, and the temperature and heating duration of the 670-km thermal boundary layer control plume morphology and strength

* Secondary plumes triggered by slabs might explain the formation of intraplate volcanism in the Eifel, Massif Central and Hainan

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

Y. Li,
yingying.li@gfz.de

**Citation:**

Li, Y., Steinberger, B., Brune, S., Le Breton, E., Glerum, A., & Pons, M. (2026). Intraplate volcanism driven by slab-plume interaction: Numerical modeling and its application to the Eifel, Massif Central and Hainan volcanic areas. *Journal of Geophysical Research: Solid Earth*, 131, e2025JB032799. https://doi.org/10.1029/2025JB032799

Received 22 AUG 2025
Accepted 25 DEC 2025

**Author Contributions:**

**Conceptualization:** Yingying Li,
Bernhard Steinberger, Sascha Brune

**Formal analysis:** Yingying Li

**Investigation:** Yingying Li

**Methodology:** Yingying Li,

Bernhard Steinberger, Sascha Brune,

Michaël Pons

**Resources:** Sascha Brune

**Software:** Yingying Li,

Bernhard Steinberger, Sascha Brune,

Anne Glerum

**Supervision:** Bernhard Steinberger,
Sascha Brune, Anne Glerum

**Validation:** Yingying Li,

Bernhard Steinberger, Sascha Brune,

Eline Le Breton

© 2026. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_3.png" /></figure>

## Intraplate Volcanism Driven by Slab-Plume Interaction: Numerical Modeling and Its Application to the Eifel, Massif Central and Hainan Volcanic Areas

Yingying Li¹,² ID, Bernhard Steinberger¹ ID, Sascha Brune¹,³ ID, Eline Le Breton³,⁴ ID, Anne Glerum¹ ID, and Michaël Pons¹ ID

¹GFZ Helmholtz Centre for Geosciences, Potsdam, Germany, ²Institute of Geological Sciences, Freie Universität Berlin, Berlin, Germany, ³Institute of Geosciences, University of Potsdam, Potsdam-Golm, Germany, ⁴University Rennes, CNRS, Geosciences Rennes – UMR 6118, Rennes, France

**Abstract** Intraplate volcanism has long been linked to deep mantle plumes. However, recent studies showed that intraplate magmatism can originate from transition zone dynamics, where lower-mantle plumes might be ponding, creating a Thermal Boundary Layer (TBL). Inspired by intraplate volcanoes in Eifel, Massif Central and Hainan that are distributed near tips of stagnant slabs imaged at transition zone depth, we hypothesize that subducting slabs might destabilize the TBL and trigger upper mantle plumes (secondary plumes), leading to intraplate volcanism. So far, the generation of such secondary plumes and the influence of slabs on plumes remain poorly understood. In this study, we perform 2D upper-mantle geodynamic models with a TBL imposed at 670 km depth interacting with a slab of an intra-oceanic subduction zone. The effects of various slab geometries (rollback, rollover and intermediate), TBL temperature and heating time are tested. Our models show that slabs of all geometries can destabilize the TBL, initiating secondary plumes ahead of and behind the slab. All plumes are deflected by the slab-induced mantle flow and a sinking slab may even suppress plumes beneath it. However, a higher TBL temperature and a longer pre-subduction heating duration may increase buoyancy flux of secondary plumes, making them more resistant to slab-driven flow. Under all conditions explored in this study, the strength of secondary plumes produced in our models is comparable to that of the Eifel plume. This paper elucidates slab-plume interaction and their impact on intraplate volcanism with applications to the Eifel, Massif Central and Hainan volcanic areas.

**Plain Language Summary** Intraplate volcanoes are located far away from plate boundaries. Earth's hot core heats the mantle, causing it to rise toward the surface in the form of narrow, hot columns known as mantle plumes. These mantle plumes have long been thought to explain intraplate volcanism. However, recent studies showed that these plumes might stall at 670 km depth, forming a hot layer. Shallower plumes, namely secondary plumes, are suspected to be generated at this depth, yet the process is not fully understood. We found in this study that subducting tectonic plates, sinking into the mantle, can induce secondary plumes from 670 km depth. Plume shapes and their upwelling intensity are affected by the shape of the subducting plate and the temperature and heating time of the bottom layer. For example, a subducting plate may weaken the underlying plumes. However, a hotter layer and a longer heating time leads to more rapidly rising plumes, which may overcome the suppressing effect of a slab. Overall, the secondary plumes modeled here ascend less vigorously than plumes generated at great depth. Secondary plumes triggered by subducting slabs might explain intraplate volcanism in the Eifel, Massif Central and Hainan volcanic areas.

## 1. Introduction

Intraplate volcanism has often been attributed to “primary” plumes, which are rooted in the deep mantle (Morgan, 1971). However, Courtillot et al. (2003) suggested that less than 10 out of 49 intraplate volcanic areas (hotspots) originate from the deep mantle plumes, while the other hotspots may have shallower origins: (a) They might be upper-mantle features controlled by the lithosphere (Anderson, 2000, 2001; Foulger, 2002) and do not need the excess temperature brought by plumes from the lower mantle. Other upper-mantle processes and features include edge driven convection (King, 2007; King & Anderson, 1995, 1998), small-scale sublithospheric convection (Ballmer et al., 2007, 2009), shear-driven upwelling (Conrad et al., 2010, 2011), downwelling drips (Elkins-Tanton, 2007), wet plumes (Kim et al., 2022; Long et al., 2019; Richard & Bercovici, 2009; Richard & Iwamori, 2010), and slab-driven upwellings, either sourced from the hydrous mantle transition zone
**Visualization:** Yingying Li, Michaël Pons

**Writing – original draft:** Yingying Li

**Writing – review & editing:**

Bernhard Steinberger, Sascha Brune,

Eline Le Breton, Anne Glerum,

Michaël Pons

(Lee et al., 2025; Yang & Faccenda, 2020) or simply a slab return flow (Faccenna et al., 2010) (b) They might also be associated with secondary plumes (Liu & Leng, 2020), referring to plumes rooted at the bottom of the transition zone between the upper and the lower mantle (~670 km in depth). The formation and dynamics of secondary plumes are the focus of this study.

At the base of the transition zone, a plume originating from the core-mantle boundary might be stagnant, forming a thermal boundary layer (TBL) that may generate secondary plumes above it (Courtillot et al., 2003). Seismic tomography has shown that mantle plumes can stagnate at specific depths: French and Romanowicz (2015), for instance, observed that plumes get laterally deflected at 1,000 km depth. Plume stagnation might be related to multiple processes. First, the phase transition with a negative Clapeyron slope at 670 km tends to impede the upwelling of plumes (Faccenda & Dal Zilio, 2017; Schubert et al., 2001). Second, entrainment of the recycled subducted oceanic crust in the form of eclogite (Ballmer et al., 2013; Dannberg & Sobolev, 2015) might cause plume stagnation at ~300–410 km depth. Third, the phase transition toward the mineral assemblage of wadsleyite, majorite, and ferropericlase with a negative thermal expansivity (Dannberg et al., 2022; Kerr & Stegman, 2024; R. Li et al., 2025) might reduce plume buoyancy at ~440–580 km depth, leading to plume stagnation. Finally, slab stagnation in the transition zone, as revealed by seismic tomography (Fukao et al., 2009; Fukao & Obayashi, 2013), can also occur due to the negative Clapeyron slope and the viscosity jump at this depth (e.g., Fukao et al., 1992; Torii & Yoshioka, 2007). These stagnant slabs may act as an additional barrier to plume ascent. In this work, we selected 670 km as a representative stagnation depth for modeling secondary plume generation. The dynamics of the secondary plumes generated at other depths are likely to be broadly similar to the plumes generated at 670 km.

While secondary plumes have been proposed to explain intraplate volcanism, their formation mechanisms remain poorly understood. Liu and Leng (2020) demonstrate that clusters of such plumes can be generated due to the combined effect of weak layers in the asthenosphere and mantle transition zone. Koptev et al. (2021) investigate secondary plume longevities and the role they play in lithospheric break-up. Cloetingh et al. (2022) examine the conditions under which they can penetrate the overlying continental lithosphere. Negredo et al. (2022) evaluate the possibility that the Canary Islands originated from a secondary plume combined with edge-driven convection. Dasgupta and Lee (2024) explore the interaction between a secondary plume and an overlying continental plate, taking into account plume-driven lithospheric weakening. Among these studies, secondary plumes are generated self-consistently only in Liu and Leng (2020). In other numerical models, the secondary plumes are imposed, leaving their initiation mechanism unresolved.

In several regions worldwide, intraplate volcanism appears to be spatially and temporally connected to subduction zones. For example, the intraplate volcanoes in western Europe are located above the NW edge of a fast seismic velocity anomaly underneath Europe (Figure 1a), interpreted as cold subducted lithosphere conveyed into the Alpine-Mediterranean area along different subduction zones (Piromallo & Morelli, 2003). A seismic tomography study also suggested the existence of a lower mantle plume under central Europe feeding the secondary plumes in the upper mantle (Goes et al., 1999). Other intraplate volcanoes, such as those near Hainan Island in southeast Asia, are also horizontally located close to the tips of slab imaged at depth (Zhao et al., 2021), indicating a possible link between subduction and plume generation. These observations motivate our hypothesis that secondary plumes might be triggered by or interact with the subducting plate.

Earlier analog and numerical models have explored slab-plume interaction, demonstrating a two-way interaction between them. Mantle plumes can modify the geometry of slabs, leading to the formation of a slab window within the subducted slab (Betts et al., 2012) or slab tear at the trench (Leonard & Liu, 2016). They can also influence the subduction rate and trench motions (Betts et al., 2012; Mériaux et al., 2015, 2016). Heilman and Becker (2022, 2024) have shown that a plume can even shut off subduction in both 2D and 3D numerical models. In turn, plumes can also be deflected by toroidal and poloidal flow induced by the subducting slab (Chang et al., 2016; Druken et al., 2014; Kincaid et al., 2013; Mériaux et al., 2015, 2016; Strak & Schellart, 2018; K. Wang et al., 2025). These studies have advanced our understanding of the mutual interaction between slabs and plumes. However, unlike models simulating plume generation from the lower mantle (e.g., Burke et al., 2008; Heron et al., 2015, 2024; Roy et al., 2025; Tan et al., 2002), the possibility of plume generation at the bottom of the transition zone and the subsequent interaction between the slab and the secondary plume remains poorly explored.

Seismic tomography (Lallemand et al., 2001; Van der Voo et al., 1999; Wortel & Spakman, 2000) and previous forward modeling studies (Bellahsen et al., 2005; Di Giuseppe et al., 2009; Funiciello et al., 2008; Xue
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) The distribution of intraplate volcanoes in Western Europe and the Central-Western Mediterranean region over the last 30 Ma, modified from Faccenna et al. (2010). The background map is a seismic tomographic model of P-wave velocity anomaly at 550 km depth from Piromallo and Morelli (2003). Note that the volcanoes in western Europe are distributed above the NW edge of the high-velocity anomaly in the transition zone. (b–d) Sketches of slabs with various geometries. Velocities applied at the plate boundaries to generate each slab geometry in our models are indicated. On each lateral boundary, inflow is compensated by prescribed outflow as shown in panel (e). (e) Boundary velocity, lithospheric composition and temperature distribution for Model R1 at 2 My. The prescribed boundary velocity is taken from R1, other velocities are used for alternative models. (f) The zoomed-in plate interface shows various compositional fields of the presented models. (g) Temperature profile and (h) viscosity profile through the subducting plate taken at *x* = 3,000 km.</figcaption></figure>

et al., 2020) have distinguished three slab geometries: (a) slab rollback (lazy “S” shaped slab; Figure 1b), (b) slab rollover (“U” shaped slab with a ~90° backward rotation; Figure 1c), (c) intermediate slab geometries with the slab outlining a lazy “S” shape while the tip of the slab is backward folded (Figure 1d). We include all three slab geometries in our model and examine their impact on plume dynamics.

This study aims to explore the following two questions: (a) Can a subducting slab trigger secondary plumes from a heated 670 km discontinuity? (b) How do variations in slab geometries affect morphology and strength of plumes that are either triggered by the slab or result from TBL heating? To this end, we perform 2D upper-mantle numerical models with a TBL imposed at the bottom and a kinematically driven subduction occurring on top. The various geometries of the slab are achieved by varying the velocity prescribed on the subducting and overriding plates. We also systematically test the effect of the temperature and TBL heating time on the morphology and strength of the plume. The detailed methodology and model set-up are shown in Section 2. The evolution of the
models is presented in Section 3. Lastly, in Section 4, we compare our model results with previous studies and discuss possible applications to the Eifel, Massif Central and Hainan volcanic areas.

## 2. Methods

### 2.1. Governing Equations

We perform 2D slab-plume interaction models utilizing ASPECT (Advanced Solver for Planetary Evolution, Convection and Tectonics, version 2.6.0-pre; Bangerth et al., 2023; Heister et al., 2017; Kronbichler et al., 2012) built on deal.II 9.4.0 (Arndt et al., 2022). The code solves the equations for the conservation of momentum (Equation 1), mass (Equation 2) and energy (Equation 3) as well as the advection equation for each compositional field (Equation 4):

$$ -\nabla \cdot (2\eta\dot{\epsilon}) + \nabla P = \rho g \quad (1) $$

$$ \nabla \cdot u = 0 \quad (2) $$

$$ \rho C_p \left( \frac{\partial T}{\partial t} + u \cdot \nabla T \right) - \nabla \cdot k \nabla T = 0 \quad (3) $$

$$ \frac{\partial c_i}{\partial t} + u \cdot \nabla c_i = q_i \quad (4) $$

with $\eta$ the effective viscosity, $\dot{\epsilon}$ the deviatoric strain rate tensor (defined as $\frac{1}{2}(\nabla u + \nabla u^T)$), $u$ the velocity, $P$ the pressure, $\rho$ the density, $g$ the gravity acceleration, $C_p$ the specific heat, $T$ the temperature, $k$ the thermal conductivity, $c_i$ the compositional fields (i.e., upper and lower oceanic crust, lithospheric mantle of overriding plate and subducting plate as well as total strain), $i$ the index indicating different fields, and $q_i$ the reaction rate (which is zero for all compositional fields, except for total strain, which we also treat as a “compositional field,” and where it consists of the current strain rate in each point multiplied by the time step size). Our simulations assume incompressibility under the Boussinesq approximation, where the density varies as:

$$ \rho = \rho_0(1 - \alpha(T - T_0)) \quad (5) $$

with $\rho_0$ representing the reference density at $T = T_0$, and $\alpha$ the thermal expansivity. The values for $C_p$, $k$, $\rho_0$, $\alpha$, $T_0$ of each compositional field can be found in Table S1 in Supporting Information S1.

### 2.2. Rheology

Our models apply a visco-plastic rheology (Glerum et al., 2018). In the viscous regime, we use a composite diffusion and dislocation creep flow law. The viscosity for diffusion creep and dislocation creep is defined as:

$$ \eta_{\text{eff}}^{\text{diff|dis}} = \frac{1}{2} A_{\text{diff|dis}}^{\frac{1}{n}} d_{n}^{m} \dot{\epsilon}_{II}^{\frac{1-n}{n}} \exp \left( \frac{E_{\text{diff|dis}} + \text{PV}_{\text{diff|dis}}}{\text{nRT}} \right) \quad (6) $$

where $A$ is the prefactor, $n$ the stress exponent, $d$ the grain size, $m$ the grain size exponent, $\dot{\epsilon}_{II}$ the square root of the deviatoric strain rate tensor second invariant, $E$ the activation energy, $P$ the pressure, $V$ the activation volume, $R$ the gas exponent, and $T$ the temperature. Diffusion creep, where $n = 1$ and $m \neq 0$, is strain rate independent but notably dependent on grain size. Dislocation creep, on the contrary, is grain size independent with $n > 1$ and $m = 0$. The contributions from dislocation creep and diffusion creep are averaged into a composite effective viscosity (van den Berg et al., 1993; Equation 7).

$$ \eta_{\text{eff}}^{\text{com}} = \left( \frac{1}{\eta_{\text{eff}}^{\text{diff}}} + \frac{1}{\eta_{\text{eff}}^{\text{dis}}} \right)^{-1} \quad (7) $$

We apply the flow law parameters of wet quartzite (Ranalli, 1995) to the altered upper oceanic crust with sediments, plagioclase (An75) (Ranalli, 1995) to the basaltic lower oceanic crust, wet olivine (Hirth &
Kohlstedt, 2003) to the weak subduction zone interface (hereafter referred to as the “weak zone”), dry olivine (Hirth & Kohlstedt, 2003) to the lithospheric mantle and the rest of the upper mantle. A representative viscosity profile is shown in Figure 1h. Details of the model layering are provided in Section 2.3. The values for A<sub>diffldis</sub>, n, d, m, E<sub>diffldis</sub>, V<sub>diffldis</sub> of each compositional field can be found in Table S1 in Supporting Information S1.

When viscous stresses exceed the yield stress, plastic yielding occurs. We apply the Drucker–Prager yield criterion (Davis & Selvadurai, 2002), where the effective plastic viscosity is defined as:

$$ \eta_{\text{eff}}^{\text{pl}} = \frac{C \cdot \cos(\Phi) + P \cdot \sin(\Phi)}{2\dot{\epsilon}_{\text{II}}} \qquad (8) $$

with C the cohesion and Φ the angle of internal friction. A cohesion value of 5 MPa is prescribed to the weak zone and the oceanic crust of the subducting plate. Internal friction angles of 5° and 2° are prescribed to the weak zone and oceanic crust, respectively. Strain weakening based on the total strain accumulation over time is also included in our models. The internal friction angles of the upper and lower crust of the subducting plate are decreased linearly from 5° to 2° over the accumulated total strain interval of 0–1. These low values are chosen to approximate a weak, lubricated plate interface which is required to maintain subduction (Duarte et al., 2015). Refer to Table S1 in Supporting Information S1 for C and Φ for other compositional fields.

## 2.3. Model Setup

This study focuses on the influence of subduction on the formation and evolution of secondary plumes in the upper mantle. Our model box has an extent of 5,000 km × 670 km (Figure 1). For simplicity, we model subduction of oceanic plates to avoid lithospheric complexities associated with continental collision. The subducting plate is located on the right side of the model domain. With a total thickness of 95 km, it comprises a 3 km thick upper oceanic crust, a 5 km thick lower oceanic crust and an 87 km thick lithospheric mantle. Since older, thicker, and denser oceanic plates generally sink beneath younger and thinner ones (Earle, 2015), we adopt a slightly thinner overriding plate in our model. The overriding plate is 82 km thick and consists of a 2 km thick upper oceanic crust, a 5 km thick lower crust and a 75 km thick lithospheric mantle. The overriding and subducting plates are initially separated by a 14 km thick weak zone, dipping 35° toward the overriding plate and extending to 82 km depth. Such weak zones allow the initiation of subduction and are commonly used in subduction models (e.g., Erdős et al., 2021; Quinquis & Buiter, 2014). Figures 1e and 1f show plate and weak zone geometries after 2 Myr of model evolution.

The thickness for the subducting and overriding plate corresponds to ages of 40 and 53 Myr, respectively, according to the half-space cooling model (Turcotte & Oxburgh, 1967). For simplicity, the plates are assumed to have laterally uniform ages, resulting in uniform initial thermal structures defined by the plate model (Carslaw & Jaeger, 1959), formulated as:

$$ T = T_0 + (T_1 - T_0) \left[ \frac{y}{y_{L_0}} + \frac{2}{\pi} \sum_{n=1}^{2} \exp \left( \frac{-\kappa n^2 \pi^2 t_a}{y_{L_0}^2} \right) \sin \left( \frac{n\pi y}{y_{L_0}} \right) \right] \qquad (9) $$

with T<sub>0</sub> the surface temperature (273.15 K), T<sub>1</sub> the temperature at the base of the plates (1573.15 K for both plates), y the depth, y<sub>L<sub>0</sub></sub> the depth of the base of the plates, κ the thermal diffusivity (10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup>), t<sub>a</sub> the age of the plates. The upper mantle beneath the lithosphere has a thermal gradient of 0.25 K km<sup>-1</sup>. Using this thermal gradient under the incompressibility assumption might slightly enhance mantle convection and increase the vigor of modeled plumes. As the temperature increase is gradual, the effect is expected to be less than if it was concentrated in the thermal boundary layer, which will be introduced in the following paragraph. Given that the two plates with 13 km thickness difference feature the same imposed basal temperature, the linear thermal gradient leads to a 3.25 K temperature difference at x = 2,500 km, at any given depth below the plates (Figure 1e, the kink of the 1400°C isotherm). The temperature difference helps to destabilize the TBL (described below) and generate plumes within a reasonable geological timescale, allowing us to test the scenario in which the slab interacts with a pre-existing plume. Such thermal variations and other heterogeneities, for example, lithospheric drips, slab fragments, commonly exist in the upper mantle and provide perturbations that can trigger the formation of small-scale convection cells.
<table><caption>**Table 1**
Summary of Key Model Parameters in This Study</caption><thead><tr><th>Model name</th><th>V<sub>o</sub> (cm/yr)</th><th>V<sub>s</sub> (cm/yr)</th><th>Slab geometry</th><th>t<sub>h</sub> (Myr)</th><th>T<sub>e</sub> (K)</th></tr></thead><tbody><tr><td>R1 (Reference)</td><td>2</td><td>1</td><td>Rollback</td><td>0</td><td>200</td></tr><tr><td>R2 (Reference)</td><td>0</td><td>3</td><td>Rollover</td><td>0</td><td>200</td></tr><tr><td>R3 (Reference)</td><td>\</td><td>\</td><td>\</td><td>\</td><td>200</td></tr><tr><td>RB20-T200</td><td>2</td><td>1</td><td>Rollback</td><td><strong>20</strong></td><td>200</td></tr><tr><td>RB80-T200</td><td>2</td><td>1</td><td>Rollback</td><td><strong>80</strong></td><td>200</td></tr><tr><td>RO20-T200</td><td>0</td><td>3</td><td>Rollover</td><td><strong>20</strong></td><td>200</td></tr><tr><td>RO80-T200</td><td>0</td><td>3</td><td>Rollover</td><td><strong>80</strong></td><td>200</td></tr><tr><td>RB00-T300</td><td>2</td><td>1</td><td>Rollback</td><td>0</td><td><strong>300</strong></td></tr><tr><td>RO00-T300</td><td>0</td><td>3</td><td>Rollover</td><td>0</td><td><strong>300</strong></td></tr><tr><td>IN00-T200</td><td><strong>1</strong></td><td><strong>2</strong></td><td><strong>Intermediate</strong></td><td>0</td><td>200</td></tr><tr><td>IN80-T200</td><td><strong>1</strong></td><td><strong>2</strong></td><td><strong>Intermediate</strong></td><td><strong>80</strong></td><td>200</td></tr></tbody></table>

Note. Values in bold font indicate the variations with respect to the reference models. $V_o$, $V_s$—velocity imposed for the overriding and subducting plate, $t_h$—the TBL heating time before onset of convergence, $T_e$—the excess temperature of the thermal boundary layer. Naming convention: We identify models with their key parameters following the format “[Slab geometry] [t_h]-T[T_e].” RB-rollback, RO-rollover, IN-intermediate slab geometry.

The temperatures of the top, left and right boundaries are fixed at the initial temperature determined by Equation 9 and the upper mantle thermal gradient. For the bottom boundary, in order to simulate the heating from a mantle plume stagnant at depth 670 km, we impose a thermal boundary layer (TBL). For simplicity and to avoid introducing additional free parameters, such as the width and position of the TBL, we apply the TBL along the entire bottom boundary. However, it does not imply a stagnant plume head with a lateral extent of 5,000 km. We merely use a relatively broad TBL to maintain the generality of the model and ensure that our model set up is applicable to a range of natural scenarios, in which the TBL width is poorly constrained. Also, we only interpret the central 3,000 km of the model. Judging from the flow fields, it is only the thermal boundary layer in the vicinity of the plume that destabilizes. Therefore, it does not matter there is a thermal boundary layer elsewhere, as long as it exists where the plume develops. The excess temperature of the TBL is 200 K in the reference models. It is varied from 200 to 400 K when testing the influence of the excess temperature (see Section 2.4). Figure 1g shows the temperature profile of the model taken at a distance of 3,000 km from the left boundary.

The mesh of our model is dependent on compositional field, temperature and depth. The weak zone, upper and lower crust of the subducting plate are prescribed with the highest resolution of 2.5 km × 2.5 km at a shallow depth of <200 km to allow subduction, while the resolution of the subducting plate can reach up to only 5 km × 5 km when it descends below 200 km. The mesh

cells of the remaining compositional fields (the whole upper plate and the lithospheric mantle of the subducting plate) are always kept 5 km × 5 km. The upper mantle below the lithosphere is mostly resolved at 10 km × 10 km, while it can also be finer resolved at 5 km × 5 km where temperature variations occur, that is, where the plumes are generated. This makes sure that the modeled mantle plumes are well resolved.

## 2.4. References Models and Key Model Parameters

The following key parameters are tested in order to study their influence on secondary plume formation and dynamics: (a) slab geometry, (b) TBL heating time before convergence ($t_h$), and (c) temperature of the TBL ($T_e$). The reference models we present in this work are based on the two end-members slab geometries: a slab rollback geometry in Model R1 (Figure 1b) and a slab rollover geometry in Model R2 (Figure 1c). A third reference case, Model R3 without subduction, is also included.

In Model R1, the slab rollback geometry is obtained by imposing a 2 cm/yr inflow at the left side of the upper plate and a 1 cm/yr inflow at the right side of the subducting plate (Figures 1b and 1e; Table 1). Model R2 has a 3 cm/yr inflow only from the right side of the subducting plate (Figure 1c and Table 1). The total convergence rate of 3 cm/yr is close to the present-day median subduction rate (Young et al., 2019) and lies between the rates of the natural cases discussed in Section 4: (a) 5.2–7.3 cm/yr for the Manila subduction zone (Schellart et al., 2011) near the Hainan plume, (b) less than 1 cm/yr for the Alpine subduction (Le Breton et al., 2021) and around 2.5 cm/yr for Calabrian subduction (Faccenna et al., 2001). To balance the model volume, outflow is prescribed at the lower part of the boundary where there is an inflow at the top (in Model R2, the outflow is only imposed on the right boundary). To avoid strong shearing on the lateral boundaries, a linear velocity gradient from in-to outflow is prescribed over a depth of 20 km (Figure 1e). Model R3 has no subduction (Table 1) with its left and right boundaries featuring free-slip conditions. The bottom boundary of all our models has non-slip conditions to simulate the effect of the high-viscosity lower mantle (Steinberger & Calderwood, 2006). The top boundary is set to free slip. Additionally, intermediate slab geometries between slab rollback (Model R1) and slab rollover (Model R2) were tested by varying the velocity of the overriding plate to 1 cm/yr and the subducting plate velocity to 2 cm/yr (models “IN” in Table 1, Figure 1d).

The impact of TBL heating time before convergence ($t_h$) was studied by varying the time before applying boundary velocities to the plates. Indeed, in models R1 and R2, the convergence is initiated at the beginning of the models ($t_h = 0$ Myr; Table 1). However, when testing $t_h$, we vary the time before we start moving the plates to 20
and 80 Myr (RB20/80-T200 and RO20/80-T200 in Table 1). A broader range of $t_h$ values ranging from 20 to 80 Myr was also explored but not shown in Section 3; these results are instead included in the regime diagram of Section 4.1. Before the initiation of convergence, while the TBL keeps heating the model domain, the lateral two boundaries are set to free slip. Finally, we tested the influence of the excessive temperature of the TBL, with a temperature ranging from 200 K in the reference models (models “T200” in Table 1), 300 K (models “T300” in Table 1) to 400 K (Section 4).

## 2.5. Calculation of the Buoyancy Flux

We quantify the strength of the plume by computing its buoyancy flux in a postprocessing step for which the script can be found in Y. Li (2025). The buoyancy flux is calculated at a depth of 300 km. Here the target plume still exists as a narrow conduit and has not yet spread beneath the lithosphere. We sample upwelling velocity and density at 10 km horizontal intervals. We then integrate plume buoyancy flux over its diameter.

$$B = 160,000 * \int_{x_0 - \frac{D}{2}}^{x_0 + \frac{D}{2}} (\rho_m - \rho_p) * V_y dx \quad (10)$$

with $x_0$ the plume center, $D$ the plume diameter, $\rho_m - \rho_p$ the density contrast of ambient mantle versus plume at 300 km depth, $V_y$ the upwelling velocity, at each point inside the plume, $dx$ the interval between the sampled points (10 km). The plume center $x_0$ represents the location where the plume reaches the highest temperature. To avoid including buoyancy flux from other plumes, we first manually check the x-axis range of motion of the target plume over the entire model, then extract $x_0$ at each timestep within this range. The diameter $D$ of the plume is initially assumed to be 160 km in our 2D simulation plane. This value is chosen to be large enough to capture the full extent of the plume as closely as possible. It can also be reduced automatically to exclude downwelling (negative $V_y$), or densities higher than that of the ambient mantle. We also assume a plume diameter of 160 km in the third dimension, which is incorporated as a multiplier in front of the integral in Equation 10. In this way, the unit of the buoyancy flux is consistent with that used in observations.

## 2.6. Model Limitations

Our models are designed for studying the generation and dynamics of secondary plumes, assuming that plume activity might potentially lead to intraplate volcanism. We note that there are 5 main simplifications of our models: (a) Partial melting is not included. Melt generation depends not only on plume strength but also on other factors such as lithospheric thickness (Ball et al., 2021; Davies et al., 2015), which are not explicitly tested in our models. Therefore, we do not directly simulate whether the melt reaches the surface and induces volcanism. Instead, we focus on the conditions under which secondary plumes can form, which may subsequently lead to melt generation and volcanism. (b) Our model domain comprises a depth range from 0 to 670 km, representing the upper mantle. Thus, the eventual penetration of slabs into the lower mantle is not simulated. However, seismic observations and numerical modeling suggest that slabs may reside in the transition zone for tens to hundreds of millions of years before descending into the lower mantle (Goes et al., 2017). Therefore, our model results remain applicable to natural systems during this prolonged stagnation stage. (c) To mimic heating by stalled plumes, we impose a TBL at 670 km depth, while the possible material influx from the lower mantle is not considered. A realistic plume excess temperature at the surface ranging from 100 to 300 K corresponds to higher excess temperatures up to 400 K at the base of transition zone (Albers & Christensen, 1996; Steinberger et al., 2025). Notably, even a plume with a relatively high excess temperature at the base of transition zone, for example, 400 K, may still be stalled. For instance, Dannberg and Sobolev (2015) demonstrate that dense thermochemical plumes, compared to thermal plumes, need higher temperatures to advance to the base of lithosphere, which implies that some relatively hot thermochemical plumes might stall. Although we do not simulate thermochemical plumes in the upper mantle, it might be possible that the thermal secondary plumes in the upper mantle have a thermochemical source in the lower mantle supplying heat for them. Moreover, stagnant slabs in the transition zone can also hinder the plume upwelling, where the heat from the plume might be trapped, leading to a higher temperature in the transition zone, similar to how plumes may be trapped by slabs near the core-mantle boundary (Tan et al., 2002). (d) We adopt a 2D setup where the toroidal flow occurring around the edge of the slab is not included. The possible entrainment of the plume material from the slab back to the mantle wedge is therefore not accounted
for. However, such entrainment portion might be rare (only 10% as shown in Mériaux et al. (2016)). Additionally, the 2D set up requires us to assume a fixed plume diameter in the third dimension, which would develop self-consistently in a 3D model. (e) The role of the strength of both plates is not explored. As we focus on plume dynamics, we use a strong overriding plate with high cohesion and friction angle (Figure S1 and Table S1 in Supporting Information S1) to minimize plate deformation complexities and reduce computation time. Variations in the subducting plate strength, density contrast with the ambient mantle, and the resulting sinking velocity were also not considered in our models. These factors could modulate the efficiency of slabs disturbing the TBL and triggering plumes.

### 3. Results

In this section, we investigate how plumes are generated and how they interact with slabs of various geometries. We first present the results of the three reference models (a) with slab rollback geometry, (b) with slab rollover geometry, and (c) without subduction. We then present models varying (a) TBL heating time, (b) TBL excess temperature, and (c) slab geometry and compare them with the reference models. In what follows, the terms “front” and “rear” plume refer to the location of the plumes with respect to the subducting plate (“front plume”: in front of the slab, beneath the overriding plate, and when two plumes are present in front of the slab, it denotes the one closest to the slab; “rear plume”: behind the slab, beneath the subducting plate). The buoyancy fluxes of the front and rear plume are both shown. Plumes first reach their maximum buoyancy flux and then gradually weakens. In the following, we refer to the plume buoyancy flux in this post-peak weakening stage as “tail buoyancy flux.”

#### 3.1. Reference Models

##### 3.1.1. Model With Slab Rollback

The reference Model R1 has a slab rollback geometry with a relatively fast trench retreat of ~2 cm/yr due to the imposed boundary velocity (See Movie S1 for its full dynamic evolution). By 18 Myr, the subducting plate penetrates halfway into the upper mantle, distorting the isotherm with the input of cold materials (Figure 2a). In this case, the poloidal flow ahead of the slab is way stronger than the poloidal flow behind it. The hot material right under the slab is displaced and preferably accumulated ahead of the slab (Figures 2a and 2b). The mantle plume ahead of the slab is subsequently triggered, reaching its highest buoyancy flux (~380 kg/s) at ~34 Myr and spreading below the bottom of the lithosphere (Figures 2c and 2e). With the further retreat of the trench, more and more space is created toward the mantle wedge, allowing the head of the front plume being advected backwards to the trench (Figure 2d). At this stage the plume becomes weaker with its buoyancy flux decreasing to 68 kg/s (Figure 2e). A weaker plume together with the poloidal return flow induced by the slab might both contribute to the plume head motion toward the trench, while the bottom of the plume is pushed progressively away from the trench by the slab (Figure 2d). Unlike the front plume, which ascends to the base of the lithosphere, there is no well-developed rear plume in this model. The thermal structure below the slab, including the 1400°C isotherm (Figure 2d), only uplifts slightly and does not form a distinct plume-like thermal anomaly. Therefore, the buoyancy flux measured at the depth of 300 km in the rear plume region remains zero throughout the simulation (Figure 2f).

##### 3.1.2. Model With Slab Rollover

The reference Model R2 features a slab rollover geometry and an approximately fixed trench (Movie S2). At ~16 Myr when the slab has sunk more than halfway through the upper mantle, the slab tip is backward reclined, cooling the ambient mantle and bending the 1400°C isotherm in a symmetric pattern around the slab (Figure 3a). When the slab reaches the bottom of the upper mantle, the hot material is almost equally pushed forward and backward by the slab (Figure 3b). Following that, the front and rear plumes are generated and reach their highest buoyancy flux at almost the same time—41 Myr for the front plume (Figures 3c and 3e) and 44 Myr for the rear plume (Figures 3d and 3f). The highest buoyancy flux for the front plume (~390 kg/s; Figure 3e) is larger than for the rear plume (~214 kg/s; Figure 3f) though. There is no trench retreat or advance in this case and both plume heads are titled away from the trench (Figures 3c and 3d).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Evolution of Model R1 with a slab rollback geometry. (a–d) Evolution of Model R1 showing the various compositional fields, temperature distribution (transparent to red), isotherms (black lines) and the velocity field (gray arrows). *t* and *t*<sub>c</sub> refer to model time and time since the onset of convergence, respectively. The isotherm initially at the depth of ~490 km corresponds to the 1400°C isotherm. Velocity arrows indicating the total velocity are scaled with the velocity magnitude. Note that the front plume is triggered by the slab while the rear plume is absent. (e, f) The buoyancy flux evolution for the front plume triggered by the slab (e) and the rear plume which is not generated (f). The four snapshot times for (a–d) are marked on (e, f), illustrating key evolutionary stages: (a) slab reaches the mid-upper mantle, (b) slab approaches the model base, (c) front plume reaches its maximum buoyancy flux and (d) trench retreats further after front plume reaches maximum buoyancy flux.</figcaption></figure>

3.1.3. Model Without Subduction

Model R3 has no subduction (Movie S3). The heat from the TBL is conducted upward into the ambient mantle, making the thermal boundary layer thicker over time (Figure 4a). The front plume in this case, is generated close to the middle of the model due to the initial horizontal temperature change in the upper mantle (see Section 2.3; Figure 4b). It reaches its highest buoyancy flux at ~ 64 Myr and its strength gradually decreases to 193, 100, 67 kg/s at 10, 20 and 30 Myr after the maximum buoyancy flux (Figures 4c and 4e). The downwelling next to the plume, forming a convection cell with the upwelling within the plume, brings cool material to deeper levels and bends the isotherm downward (Figures 4b and 4c). As a result, the hot material is accumulated adjacently at 2,900–3,300 km. Following that, another plume, termed the rear plume in this model, is generated there, achieving its highest buoyancy flux of 670 kg/s at 81 Myr and dropping to 174, 90, and 49 kg/s over the next 30 million years. Note that, due to a sufficient time of heating, the plumes generated by the TBL in this case are way stronger than the plumes triggered by the slabs in model R1 and R2. At 81 Myr, another mantle plume is about to form under the overriding plate in the left part of the domain (x ~ 1,400 km). This plume is not analyzed, as it lies further from the plate boundary than the front plume.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Evolution of model R2 with slab rollover geometry. Both plumes are triggered by the slab. The rear plume is weaker than the front one. Panels (a–d) illustrate key evolutionary stages: (a) slab reaches the mid-upper mantle and (b) slab reaches the upper mantle base, (c) front and rear plume attain their maximum buoyancy flux, respectively. (e, f) The buoyancy flux evolution of the front and rear plumes. Refer to Figure 2 for legend and notation explanations.</figcaption></figure>

## 3.2. Effect of the TBL Heating Time ($t_h$)

### 3.2.1. Slab Rollback Models With Varying $t_h$

In model RB20-T200 (see the footnote of table 1 for the naming convention), the convergence is initiated at 20 Myr model time, and the thermal boundary layer is already quite thick upon the onset of the convergence. At 38 Myr (18 Myr after the initiation of convergence), the poloidal flow ahead of the slab becomes more vigorous (Figure 5a) than that in Model R1 (Figure 2a). In this case, it only takes 27 Myr for the front plume to reach its peak strength of 543 kg/s (Figures 5b and 5d). The strength of the plume then decreases to 213, 112, 68 kg/s over the following 10, 20, 30 Myr. The highest and the tail buoyancy flux of the front plume are, thus, both increased upon the increase of the $t_h$. Although the conduit of the rear plume is tilted by the retreating slab, it still manages to reach the bottom of the lithosphere. The strength of the rear plume is smaller than that of the front plume (Figures 5d and 5e).

In model RB80-T200 where the convergence is initiated at 80 Myr model time, the convergence occurs later than the generation of the mantle plumes. The peak buoyancy flux and its occurrence time for the front and rear plume are, as a result, not influenced by the subduction and are the same as in model R3 (Figures 5i and 5j). However, the subduction does influence the tail buoyancy flux of the plumes. The front plume collides with the slab instantly upon the initiation of the convergence (Figure 5f), which reduces the strength of the plume (Figure 5i) with respect to R3. When the slab goes deeper, the plume head, specifically the side close to the subducting slab, is detached from the main upwelling structure (Figure 5g). At this stage, the subducting plate generates a poloidal flow ahead of it (Figure 5g), advecting the rest of the head of the front plume backward toward the trench, and enhancing the strength of the remaining front plume (Figures 5g and 5i). After 40 Myr of convergence, the base of the front plume is driven further from the trench, while the poloidal flow continues to slightly enhance the strength of the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolution of model R3 without subduction. Plumes are generated due to heating from the TBL. With sufficient heating, the plumes are more buoyant than plumes triggered by slabs in Models R1 and R2. Key stages of plume evolution shown in panels (a–d): (a) ambient mantle above the TBL warms before plume generation, (b) front plume ascends, and (c, d) front and rear plumes reach their maximum buoyancy flux, respectively. (e, f) Evolution of buoyancy flux for the front and rear plumes. Refer to Figure 2 for legend and notation explanations.</figcaption></figure>

front plume (Figure 5h). This results in the front plume in model RB80-T200 remaining stronger after the initial collision than in model R3, which has no subduction (Figure 5i). The rear plume keeps being tilted by the retreating trench, but its strength remains unaffected (Figures 5g and 5j). The strength of the plume is only affected after 105 Myr of model time when the trench is located above the plume and the retreating slab comes into contact with the rear plume, followed by a tilt of the plume conduit and a reduction in the strength of the plume (Figures 5h and 5j).

Overall, the morphology of the plumes in these two models is quite similar, with the geometry of the front plumes also similar to model R1.

### 3.2.2. Slab Rollover Models With Varying $t_h$

In the case of a slab rollover model, when the convergence is initiated at 20 Myr, the hot material pushed aside by the subducting plate is thicker and hotter when the slab reaches the bottom of the model than in model R2. It is especially obvious for the hot material ahead of the slab (Figure 6a). The front plume then reaches a peak strength of 555 kg/s (Figure 6b) higher than Model R2 (390 kg/s; Figure 6d), This occurs at 55 Myr, only 35 Myr after the initiation of the convergence, which is shorter than in Model R2 (41 Myr). The rear plume reaches its highest strength 7 Myr later (Figure 6c). Although this time difference is a bit larger than for Model R2, it is still a short time interval compared to models with slab rollback (e.g., the difference is 29 Myr for Model RB20-T200). The rear plume also attains its peak strength more rapidly after the onset of convergence (42 Myr) than model R2 (44 Myr), although the peak strength (~189 kg/s) is slightly smaller than for Model R2 (~214 kg/s; Figure 6e).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Evolution of slab-rollback models with varying convergence initiation time. (a–e) Evolution of model RB20-T200, including (d, e) the comparison of buoyancy flux with model R1. A 20 Myr pre-subduction heating from the TBL makes the slab-triggered front plume stronger than that of Model R1 and allows the rear plume to reach the base of lithosphere. Key stages of model evolution shown in panels (a–c): (a) slab reaches the mid-upper mantle, and (b, c) front and rear plume achieve their maximum buoyancy flux, respectively. (f–j) Evolution of model RB80-T200, including (i, j) the comparison of buoyancy flux with model R3. An 80 Myr pre-subduction heating leads to plume formation before the onset of subduction. The strength of the front plume is affected by the slab, whereas the rear plume, although deflected by trench retreat, remains unaffected in strength until 105 Myr. Key stages of model evolution shown in panels (f–h): (f, g) the tail flux of the front plume reaches its lowest point and highest point, respectively, (h) slab comes into contact with rear plume and the buoyancy flux of the rear plume gets reduced. Refer to Figure 2 for legend and notation explanations.</figcaption></figure>

The rollover-type slab initiated at 80 Myr has ample contact with the front plume. The tail buoyancy flux of the front plume is thus greatly affected by the subducting plate. At 87 Myr that is, 7 Myr after the initiation of the convergence, the subducting plate bumps into the front plume, splitting it in two. A part of the plume head remains in front of the slab, while the other part of the plume stays at the rear side of the slab (Figure 6f). Such collision between plume head and slab tip leads to a sudden drop of the plume strength evaluated at depth 300 km (Figure 6i). When the slab reaches the bottom of the box at 108 Myr, the rolled-over slab then pushes the hot material backwards, inducing a second pulse of the rearward portion of the front plume (Figures 6g and 6i). With more and more subducting plate reclined backward at the depth of 670 km, the rearward portion of the front plume is finally shoveled and shut down by the subducting plate 40 Myr after onset of convergence (Figures 6h and 6i). The rear plume is gradually pushed further back by the poloidal flow induced by the subducting plate since the onset of subduction. However, its strength is not influenced by the subducting plate, that is, it remains almost the same as in Model R3 over the evolution of the whole model (Figure 6j).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Evolution of slab-rollover models with varying convergence initiation time. (a–e) Evolution of Model RO20-T200, including (d, e) plume strength comparison with Model R2. A 20 Myr pre-subduction heating leads to a higher buoyancy flux for the front plume. The rear plume strength remains suppressed. Representative model stages shown in panels (a–c): (a) slab reaches the model base, (b, c) front and rear plume achieve their maximum buoyancy flux, respectively. (f–j) Evolution of Model RO80-T200, including (i, j) the plume strength comparison with Model R3. An 80 Myr pre-subduction heating leads to plume formation before the onset of subduction. The front plume, split in two by the slab, is strongly affected by the slab, whereas the rear plume, although deflected, remains unaffected in strength throughout the model run. Representative model stages shown in panels (f–h): (f, g) the tail flux of the front plume drops to its lowest point and reaches a maximum, respectively, (h) front plume is shoveled by the slab, its buoyancy flux drops again. Refer to Figure 2 for legend and notation explanations.</figcaption></figure>

### 3.3. Effect of the TBL Excess Temperature ($T_e$)

#### 3.3.1. Slab Rollback Models With Elevated $T_e$

The temperature of the thermal boundary layer in Model RB00-T300 is set to 300 K, 100 K higher than in the reference Model R1. A broader and more intense convection cycle is formed ahead of the slab owing to the hotter thermal boundary layer (Figure 7a). The front plume is generated at approximately the same position as in Model R1, but the peak strength of the plume (1,050 kg/s) is nearly 2.8 times stronger (379 kg/s for Model R1; Figures 7b and 7d).

The rear plume is generated (Figure 7c), in contrast with Model R1 where the rear plume is totally suppressed. It reaches its highest buoyancy flux at 47 Myr since convergence, when the trench just moves over the top of the plume. Thus, the conduit of the plume is only slightly tilted (Figure 7c), with the slab exerting a minor influence on its strength. The hot thermal boundary layer governs the strength of the plume which peaks at 1373 kg/s and surpasses the strength of the front plume (Figure 7e).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Evolution of the slab rollback (a–e) and rollover (f–j) models with elevated thermal boundary layer temperature. Plume buoyancy flux evolution and comparison with reference models are shown in panels (d, e) and (i, j). Note that the maximum value for the y axis in panels (d, e) and (i, j) is increased to 1500 kg/s. Plume strength in both models strongly increases due to the TBL temperature increase. Slab suppression of the rear plume is removed for rollback models (e), not for rollover models (j). Key model stages shown for models RB00-T300 (a–c) and RO00-T300 (f–h): (a) slab has sunk halfway through the upper mantle, (b, c) front and rear plume reach their maximum buoyancy flux, respectively, (f) slab reaches the bottom of the model, (g, h) front and rear plume reach their maximum buoyancy flux, respectively. Refer to Figure 2 for legend and notation explanations.</figcaption></figure>

### 3.3.2. Slab Rollover Models With Elevated $T_e$

For Model RO00-T300, the slab-rollover model with a warmer thermal boundary layer, both plumes' strengths are strongly enhanced, and their peak buoyancy flux occurs earlier (Figures 7f–7j), as also seen in Model RB00-T300. The difference between the peaking time of the two plumes remains small (Figures 7g and 7h). The rolling-back slab still has a suppressing effect on the rear plume despite the increased temperature of the thermal boundary layer, causing that the rear plume is still weaker than the front plume (Figures 7i and 7j).

### 3.4. Effect of the Slab Geometry

Model IN00-T200 has a moderate trench retreat (Table 1) and a slab geometry intermediate between slab rollover and slab rollback. At 18 Myr, the subducting plate penetrates midway into the upper mantle with a nearly vertical dip angle (Figure 8a). When it reaches the depth of 670 km, the tip of the slab is folded, outlining a slab rollover geometry (Figure 8b). The front plume is advected away from the trench, similar to Model R2 (Figure 3). Afterward, the trench retreats and the slab outlines a rollback geometry overall (Figures 8c and 8d). The morphology of the plumes thus becomes more similar to that shown in Model R1 (Figures 2 and 8d). The front plume reaches its highest buoyancy flux of 393 kg/s at 37 Myr (Figure 8e), of which the magnitude and timing are both similar to Model R1 (379 kg/s, 34 Myr) and R2 (391 kg/s, 41 Myr) with pure slab rollback and rollover geometry. In Model
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Evolution of the models with intermediate slab geometry. (a–f) Model IN00-T200 with convergence initiated at 0 Myr. The slab first folds, then lies flat. Plume morphology varies with slab geometry. The rear plume is not suppressed by the slab. Important model stages shown in panels (a–d): (a) slab has sunk halfway through the upper mantle, (b, c) front and rear plume reach their maximum buoyancy flux, respectively (d) trench retreats further. (g–l) Model IN80-T200 with convergence initiated at 80 Myr. The front plume prevents slab folding and plume morphology is closer to Model R1. The slab affects the tail flux of front plume, while the slightly tilted rear plume remains unaffected in strength. Important model stages shown in panel (g–j): (g) the head of the front plume is separated (h) the slab-induced poloidal flow enhances front plume (i) slab reaches the bottom of the model (j) trench retreats further. Refer to Figure 2 for legend and notation explanations.</figcaption></figure>

R1, the rear plume is prevented altogether by the fast trench retreat. The rollover-style slab in Model R2 also tends to reduce the rear plume strength. In contrast, in the intermediate case, the rear plume isn't affected by the slab and reaches its peak of 480 kg/s at 51 Myr (Figure 8f), exceeding the strength of the front plume.

In Model IN80-T200, the front plume formed prior to the onset of subduction prevents the folding of the plate, leading to a rollback inclined slab geometry. As a result, the geometries of the plumes are also similar to that in Model R1 and RB80-T200 with the front plume advected toward the trench and the rear plume away from the trench. The influence of the slab on the buoyancy flux of the front plume is similar to that in RB80-T200. The slab first separates the plume head (the side close to the subducting slab) from the front plume and reduces the strength of the plume (Figures 8g and 8k). The poloidal flow ahead of the slab then enhances the strength of the front plume (Figures 8h and 8k). The intensity of the poloidal flow ahead of the slab is closely associated with the
velocity of the subducting plate. At 114 Myr when the slab encounters resistance at the 670 km discontinuity, the plume ahead also weakens considerably (Figures 8i and 8k). As the trench in this case retreats less than in RB80-T200, the rear plume is less titled and its strength is less reduced compared to RB80-T200 (Figures 8j and 8l).

# 4. Discussion

## 4.1. Comparison With Previous Work

### 4.1.1. Plume Generation

Although our models employ an impenetrable, no-slip bottom boundary, the resulting plume behavior is broadly consistent with the model of Filippo and Faccenda (2025), who use a deeper domain (1,000 km in depth) incorporating a ten-fold viscosity jump and a phase transition with a negative Clapeyron slope of -2.5 MPa/K at 660 km depth but otherwise similar model setup. Their results show that due to the viscosity jump and phase transition at 660 km, there is only a limited amount of lower mantle material entrained into the transition zone, which supports the validity of our upper-mantle-only setup.

Our model results show that in both slab rollback and slab rollover cases (e.g., Figures 2 and 3), the descending slab depresses the thickness of the underlying TBL and pushes the hot mantle sideways. As a result, neighboring areas with a thickened TBL become unstable, eventually giving rise to the rear and front plumes. When the slab reaches the base of the transition zone, it slides horizontally and sweeps the root of the secondary plume. Therefore, the plume (front/rear plume in the slab rollback/rollover case, respectively) is always situated at the slab tip. A similar phenomenon has also been observed when the slab interacts with the TBL at the core-mantle boundary (Tan et al., 2002), despite the great difference in temperature, viscosity and other physical conditions between the upper and lower mantle.

### 4.1.2. Plume Geometry

Secondary plumes appear in various shapes in slab rollback and rollover models. In the slab rollback models, driven by poloidal flow and enabled by the additional space created by trench retreat, the head of the front plume and the upper part of the conduit are increasingly advected toward the slab and the trench, forming an elongated plume head (Figures 2 and 5), similar to that shown in earlier models (Filippo & Faccenda, 2025; Kincaid et al., 2013; Mériaux et al., 2016). For the rear plume, our models corroborate previous results (Druken et al., 2014; Mériaux et al., 2015, 2016) that trench retreat-induced mantle flow governs the geometry of the plume reaching the base of the lithosphere. The conduit of the rear plume is swept away from the slab while the portion of the plume head adjacent to the slab is carried downward by the slab (Figure 5h, *x* = 2,900–3,200 km). Beyond this, we find that the slab might even prevent the generation of the rear plume from the transition zone under the condition of fast trench retreat and slab rollback (Figure 2).

The interaction between plumes and slabs with a rollover geometry has rarely been studied before. K. Wang et al. (2025) show that the rollover-type slab tilts the rear plume away from the trench before the slab sinks into the lower mantle, which is also corroborated by our models (Figures 3 and 6). For the front plume, our simulations additionally demonstrate that, due to a lack of space between the plume and slab, the front plume is also tilted away from the trench (Figures 3 and 6, rollover models; Figure 8b, the rollover stage of the intermediate models). This is distinct from the front plume in the rollback models. Although, in reality, the slab will eventually penetrate the 670 km discontinuity, plume head tilting and plume migration may still occur while the slab remains stagnant at this depth. Depending on the position, the plume located right below the slab in this case can also be split (Figure 6f), in agreement with Leonard and Liu (2016).

### 4.1.3. Effect of Boundary Conditions

An important parameter we tested is the temperature of the TBL. The temperature increases the buoyancy of the plume, making the plumes stronger across all slab geometries (Figures 9a, 9c and 9e). It also leads to an earlier occurrence of upwelling (Figures 9b, 9d and 9f) consistent with K. Wang et al. (2025). Moreover, the stronger and faster plume resulting from the elevated TBL temperature exhibits greater resistance to the slab-driven flow, in agreement with Druken et al. (2014). This is evident in models where plumes are triggered by the slab (models left of the dashed line in Figure 9). The greater resistance to slab-driven flow manifests as follows: (a) In the slab rollback models, a higher TBL temperature allows the formation of the rear plume (Figure 7), thus its buoyancy
<figure><img src="image_12.png" /></figure>
flux and timing of occurrence can be detected (Figures 9a and 9b). Whereas in Model R1, the rear plume is not generated (Figure 2). (b) Rear plumes are not weakened by slab rollback at a high TBL temperature compared to the front plumes. The rear plumes in the 200-K models are consistently suppressed by the slab, with lower buoyancy fluxes than the front plumes in the slab rollback and rollover models (Figures 9a and 9c), and they are only slightly stronger in the intermediate models (Figure 9e). The time at which the rear plume reaches its highest buoyancy flux is also delayed in Model RB20-T200, almost as much as in RB40-T200 (Figure 9b). On the contrary, in the 300-K and 400-K models, the rear plumes are mostly stronger than the front plumes and they are not delayed.

In addition, we test the effect of the TBL heating time at the transition zone. This has not yet been explored by previous studies. When the plumes are triggered by the slab (left to the dashed line in Figure 9), TBL heating time plays a similar role as temperature in strengthening plumes (Figures 9a, 9c and 9e) and promoting earlier plume ascent (Figures 9b, 9d and 9f). This also makes plumes more resistant to slab-driven flow. For example, the rear plume remains weak in RO00-T300, while it is stronger and no longer suppressed by the slab in RO20-T300. When the plumes are generated prior to subduction (those right of the dashed line in Figure 9), the plots are not significantly different from the corresponding no-subduction models in most cases, indicating that the slab in general does not have a significant impact on the highest strength of the plume and its occurrence time. One exception is Model RO40-T300, where the subducting slab happens to collide with the rising front plume, splitting it in two and largely reducing the strength of the plume (Figure 9c). Subduction does not have any impact on the plumes' highest buoyancy flux for the models in the shaded areas of Figure 9, while it does influence the tail buoyancy flux of the plume as described in Section 3. However, across all tested TBL excess temperatures, the tail buoyancy flux of the plumes remains relatively low overall (Figure S2 in Supporting Information S1). These values are comparable with upper mantle plumes in nature, which will be discussed in Section 4.2.

### 4.1.4. Comparison With Other Intraplate Volcanism Models

Our models account for both slabs and the TBL that may be caused by a stagnant plume, highlighting the role of the subducting plate in destabilizing the TBL and facilitating secondary plume generation. The stagnant plume itself can branch into clusters of secondary plumes, which move, merge then branch again due to the combined effect of the low-viscosity layer in the asthenosphere and the transition zone (Liu & Leng, 2020). Our models show that without the weak layers, clusters of secondary plumes can still be observed (e.g., Model R3) while they remain stationary. Slabs can promote the secondary plume generation and drive lateral plume migration.

Slab-related hypotheses for the formation of intraplate volcanism have been proposed by earlier studies. Faccenna et al. (2010) found that the subducted slab itself can trigger a focused upwelling in front of it, after which decompression melting might occur, leading to intraplate volcanism. Their models do not however include a stagnant plume from the lower mantle. The maximum vigor of the poloidal flow in front of the slab is generated before the slab reaches the base of the transition zone. Our models show that in the presence of a TBL at 660 km depth, the highest buoyancy flux of the front plume is obtained when the slab is close enough to the TBL, pushing the hot material ahead. While melt generation is not calculated in our models nor in Faccenna et al. (2010), it is clear that the presence of a TBL would enhance melt generation because of the excess temperatures it introduces.

Our models do not account for water, neither from the slab's uppermost layer nor from the mantle transition zone. Nevertheless, small plumes rise from the low-viscosity uppermost layer of the slab (e.g., Movie S1, 90 Myr onwards). These structures reproduce the dynamics of wet plumes (Kim et al., 2022; Richard & Bercovici, 2009; Richard & Iwamori, 2010) that are at the top surface of the slab, unlike the slab-triggered secondary plumes that are located at the slab tip. The expelling of water from the transition zone by the slab can also lead to intraplate volcanism (Yang & Faccenda, 2020). The location of the melt generation in this case is similar to the plume location in our model, that is, in front of the slab tip and behind the slab.

**Figure 9.** The regime diagram showing the maximum buoyancy flux (a, c, e) and the time at which plumes reach their maximum buoyancy fluxes (b, d, f) versus thermal boundary layer (TBL) heating time before subduction of the slab rollback models (a, b), rollover models (c, d) and intermediate models (e, f). The models without subduction are plotted in the rightmost column of each subfigure ("No subd."). The dashed line separates the plumes triggered by the subducting slab to the left and the plumes generated by the TBL to the right. The symbols enclosed by the gray shadow—indicating where the peak flux occurs earlier than subduction initiation— represent plumes whose maximum buoyancy flux is not influenced by subduction.
Geophysical observations might help to distinguish our proposed mechanism from other slab-related upper mantle processes. For instance, in seismic tomographic images, our proposed mechanism is expected to produce a secondary plume located at the tip of a stagnant slab, above a stagnant lower-mantle plume. Electrical conductivity can provide constraints on the water content in the transition zone (Karato, 2011; Yoshino & Katsura, 2013), which helps to assess whether water plays a role in melt generation. However, due to the uncertainty related to the resolution of the seismic tomography and the different relationships between water content and electrical conductivity proposed by different authors (e.g., Dai & Karato, 2009; Romano et al., 2009; Yoshino & Katsura, 2012), it is overall difficult to distinguish between different mechanisms. All slab-related upper mantle plumes are expected to share similar geochemical signatures, for example, low $ ^{3}He/$^{4}He ratio characteristic of upper-mantle sources, high $ \delta^{13}C $ and high $ ^{238}U/$$ ^{204}Pb $ ratio indicative of recycled oceanic crust, which makes it difficult to distinguish between different mechanism by solely using geochemical constraints.

## 4.2. Application to Eifel and Massif Central Volcanic Areas

The generation of the intraplate volcanoes in the Eifel volcanic field (EVF) and the Massif Central (MC), part of the central European volcanic province (CEVP), has long been associated with plume-slab interactions (Goes et al., 1999). Low P-velocity anomalies can be imaged down to at least 410 km depth below the EVF (Figures 10b and 10c; Ritter, 2007; Ritter et al., 2001; Zhu et al., 2012) and MC (Figure 10d; Granet et al., 1995; Sobolev et al., 1997). A root in the lower mantle which might support the heat or material for the Eifel plume in the upper mantle can also be visualized (Figures 10b and 10c; Goes et al., 1999; Zhao, 2007). Nobel gas isotopes studies support that the EVF, along with some other central European volcanoes, has a deep mantle plume origin (Bräuer et al., 2013; Buikin et al., 2005; Caracausi et al., 2016), while Bekaert et al. (2019) show that the mantle source beneath EVF and CEVP is akin to the upper mantle reservoir.

The plume beneath the EVF in the upper mantle seems detached from its lower mantle root (Figures 10b and 10c), which is also reflected in the depth of the discontinuity: Only the 410-km discontinuity is depressed while the 660-km discontinuity remains unchanged (e.g., Budweg et al., 2006). It has been suggested that the detachment of the lower mantle root and the upper mantle plume is related to the stagnation of subducted lithosphere of the Alpine and Mediterranean subduction zones (Figures 10b–10d; Budweg et al., 2006; Goes et al., 1999). The evolution of these subduction zones is quite complex. Indeed, the subduction changes polarity between the Alps and the Apennines (west of Adriatic Plate)-Dinarides (east of Adriatic Plate, Figure 10a). The European Plate is the subducting plate beneath the Alps but it is the overriding plate of the Western Mediterranean subduction zones (Apennines-Calabria-Maghrebides-Betic, Figure 10a). The cross-sections A-A' and B-B' (Figure 10) show the location of the EVF and MC with respect to these subduction zones. The NW-SE profile A-A' passes through the EVF, the frontal thrusts of the Alps, ~400 km southeast of the EVF, and then cuts obliquely through the Apennines and reaches the Calabrian Arc. Seismic tomography models along this profile (Figures 10b and 10c) show that the Eurasian slab beneath the Alpine subduction appears to be connected at depth to subducted Alpine Tethys oceanic lithosphere (“Tethys relicts” as interpreted in e.g., Handy et al., 2010, 2021), forming a slab rollover geometry. Note that the slabs might also be juxtaposed rather than forming a continuous slab. However, as the subducted Eurasian slab goes down, it will inevitably push Tethys relicts situated at the Eurasian slab tip within the transition zone. Thus, the slabs can also be “dynamically connected” in a shape mimicking a slab rollover geometry. Although a continuous plate with a slab rollover geometry is rarely observed at present, and is commonly attributed to trench advance (Di Giuseppe et al., 2009; Xue et al., 2020), the “dynamically connected” slabs in the shape of a slab rollover geometry could occur even if a slab retreats. Also, the profile A-A' cuts obliquely through the Adriatic slab beneath the Apennines and through the long Calabrian slab. Due to the obliquity of the profile, the large positive anomaly at depth most likely represents a mix of Adriatic slab, Calabrian slab and relics of subducted oceanic lithosphere (Figures 10b and 10c). The cross-section B-B' (Figure 10) passes through the MC on the Eurasian plate and orthogonally through the Calabrian slab, which depicts a clear slab rollback geometry associated with the rollback of the Western Mediterranean subduction zones (e.g., Faccenna et al., 2001).

Our models show that the slab of the rollover and rollback geometries can both push the hot material in the direction of slab tip motion, leading to plume formation (Figures 2 and 3). The subducted Eurasian plate together with the putative Alpine-Tethys slab within the transition zone outlines a slab rollover geometry (Figures 10b and 10c). The total Adriatic-European convergence ranges from 560 km (Van Hinsbergen et al., 2020) to 680 km (Le Breton et al., 2021) to 930 km (Handy et al., 2015), which might allow the subducted Eurasia slab and Alpine-
<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a) Geological overview of the region surrounding the Eifel volcanic field and the Massif Central. Base map from GeoMapApp (www.geomapapp.org)/CC BY/CC BY (Ryan et al., 2009) (b, c) Cross sections of P-velocity tomography from (b) UUP07 (Amaru, 2007) and (c) MITP08 (C. Li et al., 2008) through the Eifel volcanic field and multiple thrusts to the south, with their location shown in panel (a). Zoom-ins of the upper 1000 km are shown on the right with the name and the contours of the slabs, interpreted by for example, Handy et al. (2010, 2021). (d) Cross section of P-velocity anomaly through the Massif Central and Calabrian subduction zone from Piromallo and Morelli (2003) and Faccenna et al. (2010), location shown in panel (a). Abbreviations: AFT, Alpine Frontal Thrust; CFT, Calabrian Frontal Thrust; EVF, Eifel volcanic field; MC, Massif Central; LG, Limagne Graben; BG, Bresse Graben; URG, Upper Rhine Graben; LRG, Lower Rhine Graben; HG, Hessian grabens; EG, Eger (Ohre) Graben.</figcaption></figure>

Tethys slab to reach the vicinity of the EVF and MC and generate the upper mantle plume to its northwest as shown in Model R2 (Figure 3, rear plume). Conversely, the Calabrian slab depicts a slab rollback geometry (Figures 10b–10d). It has a length of ~1,200 km, that is, related to the slow Africa–Europe convergence and the fast trench retreat (Faccenna et al., 2001), and it could push the material accumulated in the transition zone from other subduction zones forward, contributing to the generation of the MC (Figure 2, front plume). This interpretation aligns with Faccenna et al. (2010), who attributed the generation of MC to the emplacement of cold
material in the transition zone. Note that a burst of magmatism is expected to occur in front of the tip of the slab when the highest buoyancy flux of the slab-driven plume is achieved, not at the present-day slab tip location. This means it is difficult to find unambiguous evidence in nature in the present day. In terms of plume geometry, in our slab rollback models, the head of the front plume is advected toward the trench due to the space created by trench retreat. In this case, the upper mantle seismic anisotropy also suggests the southeastward deflection of the plume under MC due to slab rollback (Barruol & Granet, 2002). Since the Eurasian slab lies between the EVF and the Calabrian subduction zone, the Calabrian subduction might not contribute substantially to the EVF.

The Alpine subduction initiated at ~84 Myr ago (Ma) (Le Breton et al., 2021), predating the Quaternary intraplate volcanism in the Eifel and the major magmatic activity in the Massif Central (~15 Ma) for more than 69 Myr. In our models, the time required for a plume to reach the base of the lithosphere after the onset of subduction is 17–56 Myr, depending on factors such as temperature and slab geometry (Figures 9b, 9d and 9f). Given that the Alpine subduction (0.71–1.26 cm/yr) and collision rate (maximum 0.6–0.89 cm/yr) are relatively low (data from Handy et al., 2015; Le Breton et al., 2021; Van Hinsbergen et al., 2020), the actual time needed in the study area is likely longer than in our models. Nevertheless, the more than 69 Myr interval remains sufficient for a plume to arrive at the lithosphere base. The older phases of the volcanism in Eifel (46 Ma) and Massif Central (65 Ma) are temporally too close to the initiation of the Alpine subduction, and considering the slow subduction rate they might not be associated with secondary plumes triggered by slabs. This supports previous geological, geochronological and isotopic studies, which suggest that only the Quaternary volcanism in Eifel and the post-15 Ma volcanism in Massif Central are plume-related, whereas the older phases are pre-rift or rift-related magmatic event (Fekiacova, Mertz, & Hofmann, 2007; Fekiacova, Mertz, & Renne, 2007; Michon & Merle, 2001). The rollback subduction and fast trench retreat of the Calabrian subduction zone started at ~35 Ma with a subducting velocity of ~2.5 cm/yr (Faccenna et al., 2001, 2003). Prior to this, because the slow and oblique Africa–Europe convergence had likely already driven ~300 km of Calabrian slab subduction at a shallow dip since ~80 Ma, the slab was able to reach the base of the transition zone at 15 Ma (Faccenna et al., 2001, 2003). When the slab reached the base of the transition zone, a front plume triggered by the slab could reach the surface, leading to intraplate volcanism in the MC, as suggested by our slab rollback models (Figures 2c and 5b). The older phases of volcanism in MC (65 Ma) occurred before the fast rollback subduction of the Calabrian slab and the slab did not reach a sufficient depth due to its shallow subduction dip by then. Therefore, the older phase of volcanism in MC might not be related to the secondary plume triggered by the Calabrian subduction either.

The Eifel plume is classified as a secondary plume in the upper mantle (Courtillot et al., 2003). Its buoyancy flux has been estimated using different methods, yielding values of 10, 40 and 120 kg/s (King & Adam, 2014). It is a lot weaker than the primary plumes originating from the core mantle boundary with buoyancy fluxes exceeding 1,000 kg/s (Courtillot et al., 2003). This is supported by our result that all the upper-mantle plumes are quite weak regardless of the slab geometry. Twenty million years after the peak flux, the buoyancy flux of the plumes in all the models drops to less than ~250 kg/s, including the models with the highest excess temperature of 400 K (Figure S2 in Supporting Information S1). Previous studies indicate an excess temperature of ~200 to ~300 K for the Eifel plume (Budweg et al., 2006) and the plume under Massif Central (Sobolev et al., 1997). In our models, plumes of these excess temperatures are even weaker than the 400-K models. The buoyancy fluxes drop down to less than ~150 kg/s at 20 Myr and ~100 kg/s at 30 Myr after the peak flux, in agreement with the estimate by King and Adam (2014). For the plume beneath Massif Central (also the Hainan plume which we discuss in the next section), no observational estimates are available, and therefore a direct comparison with our modeled buoyancy flux is not possible. However, our model results establish the basis for future comparison once observational buoyancy-flux estimates for these plumes are obtained.

By accounting for subducting slabs and a stagnant plume, which have both been visualized under Eifel and Massif Central, our models provide a plausible mechanism for the Eifel and Massif Central intraplate volcanism. However, we cannot rule out other possibilities suggested by earlier studies, such as wet plumes formed due to the detachment of the low-viscosity portion of the uppermost subducting slab (Kim et al., 2022; Long et al., 2019; Richard & Bercovici, 2009; Richard & Iwamori, 2010; this can also be seen from our models, e.g., Movie S1), edge-driven convection (El-Sharkawy et al., 2024; both Eifel and Massif Central are located at the transition from thick to thin lithosphere), slab-driven upwelling (Faccenna et al., 2010; Yang & Faccenda, 2020) and TBL heating without slab (Liu & Leng, 2020). Note that these mechanisms are not mutually exclusive and may operate simultaneously. Further geological, geophysical and geochemical constraints are needed to distinguish their relative importance.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** (a) Major subduction zones and Holocene volcanoes in the vicinity of the Hainan plume (located under Hainan Island). Base map from GeoMapApp (www.geomapapp.org)/CC BY/CC BY (Ryan et al., 2009). (b) E-W cross sections through P-wave seismic tomography models UU-P07 (Amaru, 2007) and MITP08 (C. Li et al., 2008) which cut across the Hainan plume and the Manila trench. (c) P-wave tomographic image from Lei et al. (2009), showing the tilted Hainan plume. Abbreviations: HNI, Hainan Island; LP, Leizhou Peninsula; PP, Pacific plate; MNT, Manila trench; PLT, Philippine trench; ELT, East Luzon trench; IBT, Izu-Bonin trench; PSCS, proto-South China Sea; PSP, Philippine Sea plate; PP, Pacific plate; SCS, South China Sea.</figcaption></figure>

4.3. Application to the Hainan Plume – South China Sea Slab

The Hainan plume is evidenced by the Holocene intraplate volcanoes on Hainan Island, the Leizhou Peninsula and the Indochina block (Figure 11a). Cenozoic basalts sharing common intraplate oceanic basalt-like features, of which the ages range from 28.5 to <0.1 Ma, are also widely distributed in these areas and the South China Sea basin (Gu et al., 2019; Ho et al., 2000; Hoang et al., 1996; Yan et al., 2018). Multiple subduction zones, created by the subduction of the South China Sea plate, Philippine plate and Pacific plate to the east and the Indo-Australian plate to the southwest, surround the Hainan plume (Figure 11a). Seismic tomography models imaged that (a) the Hainan plume is rooted in the lower mantle (Figure 11b), as also shown in other tomographic models (J. Huang & Zhao, 2006; Montelli et al., 2006; Xia et al., 2016; Zhao, 2007); and (b) it is blocked by the relicts of the slabs in the transition zone, thus ponding below the depth of 660 km (Figure 11b; also shown in Zhao, 2007; Zhao et al., 2021). Xia et al. (2016) show a shallower ponding depth of 350–500 km, above which the plume branches out into multiple smaller conduits. The slab relicts have been identified by many tomographic studies (Hall & Spakman, 2015; Z. Huang et al., 2015; Lallemand et al., 2001; Obayashi et al., 2013; Rangin et al., 1999; Zahirovic et al., 2016), however their origins remain uncertain. Possible origins include the proto-South China Sea “north” slab (Wu & Suppe, 2018), the (Paleo-) Philippine Sea plate and the Pacific remnants (Hall, 2002; Hall
& Spakman, 2015; Yin, 2010). Overall, the proximity of the plume to subduction zones and its ponding suggest interaction between the slabs and the plume, which is also supported by the coexistence of the FOZO (lower mantle component) and recycled oceanic crust component in basaltic rocks of the region (X.-C. Wang et al., 2013; Yu et al., 2022).

Among the four subducting plates mentioned above, the deep subduction of the Indian and Pacific plates down to the core-mantle boundary (He & Wen, 2011) was proposed to be associated with the birth of the Hainan plume (Wei et al., 2012; Zhang & Li, 2018; Zhao et al., 2011). The interaction between the deep slab and the plume is out of the scope of our models. However, the plume did not reach the surface until ~28.5 Ma, when the earliest intraplate basalt was generated in the Leizhou peninsula, possibly due to the blocking by the stagnant slabs (Yu et al., 2022). The stagnant Indo-Australia, Pacific and proto-South China slabs may all contribute to plume ponding. Thus, one can expect the formation of secondary plumes on top of the ponding plume in the upper mantle. Zhao et al. (2021) proposed that the opening of the South China Sea, starting at ~33 Ma and ending at ~15 Ma (C.-F. Li et al., 2014), strengthened the convection in the upper mantle which might induce the upwelling of the secondary plume from the transition zone. Here we show that subduction in the vicinity might also contribute to the upwelling of the Hainan plume.

No buoyancy flux estimates for the Hainan plume are available. Nevertheless, the spatial relationship between the Hainan plume and subducted slabs in the vicinity allows a qualitative comparison. The subduction of the South China Sea (SCS) plate along the Manila Trench and the Philippine Sea plate along the East Luzon Trough creates a double subduction zone to the east of the Hainan plume. The subducted SCS slab dips eastward with a steep angle of about 80° (Figure 11b). Driven by a westward mantle flow (Steinberger & Calderwood, 2006; Steinberger & Holme, 2008), it might be mechanically coupled with the stagnant slab to its west, pushing it further west as subduction continues. This process may facilitate the generation of the secondary plume further west, regardless of the slab origin. The subducted SCS slab together with the stagnant slab thus outlines a slab rollover geometry as a whole, with the Hainan Plume located on the rear side of the slab. As shown in our model with this slab geometry and a fixed trench, the rear plume head can be deflected away from the trench due to a lack of space between the slab and plume (e.g., Figures 3c and 3d). The Manila trench retreats westward with a trench-normal speed of 3–5 cm/yr from south to north (Schellart et al., 2011). Such a retreat might further reduce the space between the slab and the Hainan plume at its rear, causing the NW tilting of the plume observed in the seismic image (Figure 11c). This provides an alternative explanation for the Hainan plume tilting in addition to the slab rollback-induced toroidal flow (Mériaux et al., 2015), and the rising of a tilted plume conduit from the lower mantle (Lei et al., 2009).

The subduction of the SCS lithosphere along the Manila trench was initiated in the Mid-Miocene (~16 Ma to 11.6 Ma; Taylor & Hayes, 1983). Therefore, it could have triggered a secondary plume that resulted in the younger group of OIB at the surface, particularly the most recent ones that formed long after subduction initiation. Wu and Suppe (2018) proposed a northward subduction of the Proto-South China Sea (PSCS) “north” slab that initiated at ~34 Ma. This happened before subduction of the SCS, which generated the slab relicts observed in the transition zone in Figure 11b. With its trench retreating southward, the slab might follow a slab rollback geometry. We speculate that the rollback-style PSCS slab might have advected the head of Hainan plume to the SSE toward its trench before the plume was governed by mantle flow induced by SCS subduction and tilted toward WNW. PSCS subduction might also have generated earlier OIB-type basalts, although the oldest basalts (~28.5 Ma) in this area, which were generated only ~5.5 Myr after subduction initiation might be unrelated. The East Luzon Trough initiated only at 3–5 Ma (Cardwell et al., 1980; Hamburger et al., 1983; Karig, 1975; Ozawa et al., 2004), and the Philippine Sea plate that subducted there only extends to quite shallow depths, characterized by a Wadati-Benioff zone of up to 80 km (Defant et al., 1989). As a result, it might not have had a huge impact on the Hainan plume at a distance of ~1,300 km away from it. Altogether, our models suggest that the subduction of the SCS and PSCS slabs might help to destabilize the TBL and promote the upwelling of the Hainan plume. However, other explanations, such as the opening of the South China Sea (Zhao et al., 2021), might complement the mechanism we propose.
5. Conclusion

Secondary plumes originating from the transition zone can result in intraplate volcanism. This study explores the
effect of subducting slabs of various geometries on the generation and morphology of such plumes. Our main
findings are:

1. Subducting slabs can trigger secondary plumes from the base of the mantle transition zone. These plumes can migrate with the motion of the slab along the bottom of the transition zone. Their morphologies are modified in different ways depending on the slab geometry and trench motion. Driven by slab rollback and the associated mantle flow, the rear plume head (beneath the subducting plate) in the slab rollback scenario is tilted away from the trench, while the space created in the mantle wedge allows the head of the front plume (beneath the overriding plate) to be advected toward the trench. In the slab rollover scenario, the heads of rear and front plumes are both advected away from the trench, as the trench location is fixed and no space is created toward the trench.

2. The secondary plumes generated in our models are overall quite weak (with buoyancy fluxes up to 100–250 kg/s), consistent with natural observations. Slab geometries, TBL excess temperature ($T_e$) and TBL heating time ($t_h$) can all influence the plumes' strength. The slab tends to dampen the rear plume in both slab rollback and rollover cases, while in the intermediate model, the slab does not appear to significantly influence the rear plume. $T_e$ plays a key role in enhancing plume strength, thereby increasing its ability to resist slab-driven flow. $t_h$ can only increase the strength of secondary plumes that are triggered by the slab.

3. We propose a new mechanism for the generation of intraplate volcanism which can be applied to at least 3 natural cases. The subducted Eurasian plate beneath the Alps and Calabrian slab beneath the Western Mediterranean, dynamically connected with stagnant slabs of subducted oceanic lithosphere (Tethys) in the transition zone, might trigger secondary plumes including the Eifel plume, leading to the Quaternary volcanic activities in Eifel and the major magmatic phase of Massif Central that began at ~15 Ma. The subducted South China Sea plate, coupled with the stagnant slab in the transition zone, and the possible subduction of the proto-South China Sea “north” slab might have contributed to the further upwelling of the Hainan plume in the upper mantle. The rollover geometry outlined by the subducted South China Sea slab and the stagnant slab, combined with the Manila trench retreat, might explain the northwest tilt of the Hainan plume.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The input files, additional plugins, postprocessing script, and the modified version of ASPECT used in this paper
can be found in the Zenodo repository at: https://doi.org/10.5281/zenodo.16884447 (Y. Li, 2025). The modified
version of ASPECT (version 2.6.0-pre, commit ca24c3a, dealii 9.4.0) is based on ASPECT 2.6.0-pre, commit
025144e25. Figures were generated using ParaView 5.11.2 and Python 3.11.4.

References

Albers, M., & Christensen, U. R. (1996). The excess temperature of plumes rising from the core-mantle boundary. *Geophysical Research Letters*, *23*(24), 3567–3570. https://doi.org/10.1029/96GL03311

Amaru, M. L. (2007). *Global travel time tomography with 3-D reference models* (Vol. 274). Utrecht University.

Anderson, D. L. (2000). The thermal state of the upper mantle; No role for mantle plumes. *Geophysical Research Letters*, 27(22), 3623–3626.
https://doi.org/10.1029/2000GL011533

Anderson, D. L. (2001). Top-down tectonics? *Science*, 293(5537), 2016–2018. https://doi.org/10.1126/science.1065448

Arndt, D., Bangerth, W., Feder, M., Fehling, M., Gassmöller, R., Heister, T., et al. (2022). The deal.II library, Version 9.4. *Journal of Numerical Mathematics*, 30(3), 231–246. https://doi.org/10.1515/jnma-2022-0054

Ball, P. W., White, N. J., Maclennan, J., & Stephenson, S. N. (2021). Global influence of mantle temperature and plate thickness on intraplate volcanism. *Nature Communications*, 12(1). https://doi.org/10.1038/s41467-021-22323-9

Ballmer, M. D., Ito, G., Wolfe, C. J., & Solomon, S. C. (2013). Double layering of a thermochemical plume in the upper mantle beneath Hawaii. *Earth and Planetary Science Letters*, 376, 155–164. https://doi.org/10.1016/j.epsl.2013.06.022

Ballmer, M. D., van Hunen, J., Ito, G., Bianco, T. A., & Tackley, P. J. (2009). Intraplate volcanism with complex age-distance patterns: A case for small-scale sublithospheric convection. *Geochemistry, Geophysics, Geosystems*, 10(6). https://doi.org/10.1029/2009GC002386

Ballmer, M. D., van Hunen, J., Ito, G., Tackley, P. J., & Bianco, T. A. (2007). Non-hotspot volcano chains originating from small-scale sublithospheric convection. *Geophysical Research Letters*, 34(23). https://doi.org/10.1029/2007GL031636

Acknowledgments

We thank the Computational Infrastructure
for Geodynamics (geodynamics.org),
which is funded by the National Science
Foundation under award EAR-0949446
and EAR-1550901, for supporting the
development of ASPECT. Y. Li.
acknowledges financial support from the
China Scholarship Council (CSC). B.
Steinberger is funded by the innovation
pool of the Helmholtz Association through
the Advanced Earth System Modelling
Capacity (ESM) activity. S. Brune is
funded by the European Union (ERC,
EMERGE, 101087245). A. Glerum is
funded by the Helmholtz Association
(Recruitment Initiative 0316). M. Pons is
funded by the ERC MEET Grant (856555).
We gratefully acknowledge the computing
time granted by the Resource Allocation
Board and provided on the supercomputer
Lise at NHR@ZIB and Emmy at NHR-
Nord@Göttingen as part of the NHR
infrastructure. The calculations for this
research were conducted with computing
resources under the project bbp00039 and
bbp00064. The authors would like to thank
the reviewers Changyeol Lee, Ana M.
Negredo, the associate editor, and editor
Mark Behn for their constructive feedback
on this manuscript. Open access funding
enabled and organized by Project DEAL.
Bangerth, W., Dannberg, J., Fraters, M., Gassmoeller, R., Glerum, A., Heister, T., et al. (2023). Geodynamics/aspect: ASPECT 2.5.0 [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.8200213

Barruol, G., & Granet, M. (2002). A Tertiary asthenospheric flow beneath the southern French Massif Central indicated by upper mantle seismic anisotropy and related to the west Mediterranean extension. *Earth and Planetary Science Letters*, 202(1), 31–47. https://doi.org/10.1016/S0012-821X(02)00752-5

Bekaert, D. V., Broadley, M. W., Caracausi, A., & Marty, B. (2019). Novel insights into the degassing history of Earth's mantle from high precision noble gas analysis of magmatic gas. *Earth and Planetary Science Letters*, 525, 115766. https://doi.org/10.1016/j.epsl.2019.115766

Bellahsen, N., Faccenna, C., & Funiciello, F. (2005). Dynamics of subduction and plate motion in laboratory experiments: Insights into the “plate tectonics” behavior of the Earth. *Journal of Geophysical Research*, 110(B1). https://doi.org/10.1029/2004JB002999

Betts, P. G., Mason, W. G., & Moresi, L. (2012). The influence of a mantle plume head on the dynamics of a retreating subduction zone. *Geology*, 40(8), 739–742. https://doi.org/10.1130/G32909.1

Bräuer, K., Kämpf, H., Niedermann, S., & Strauch, G. (2013). Indications for the existence of different magmatic reservoirs beneath the Eifel area (Germany): A multi-isotope (C, N, He, Ne, Ar) approach. *Chemical Geology*, 356, 193–208. https://doi.org/10.1016/j.chemgeo.2013.08.013

Budweg, M., Bock, G., & Weber, M. (2006). The Eifel Plume—Imaged with converted seismic waves. *Geophysical Journal International*, 166(2), 579–589. https://doi.org/10.1111/j.1365-246X.2005.02778.x

Buikin, A., Trieloff, M., Hopp, J., Althaus, T., Korochantseva, E., Schwarz, W. H., & Altherr, R. (2005). Noble gas isotopes suggest deep mantle plume source of late Cenozoic mafic alkaline volcanism in Europe. *Earth and Planetary Science Letters*, 230(1–2), 143–162. https://doi.org/10.1016/j.epsl.2004.11.001

Burke, K., Steinberger, B., Torsvik, T. H., & Smethurst, M. A. (2008). Plume Generation Zones at the margins of Large Low Shear Velocity Provinces on the core–mantle boundary. *Earth and Planetary Science Letters*, 265(1–2), 49–60. https://doi.org/10.1016/j.epsl.2007.09.042

Caracausi, A., Avice, G., Burnard, P. G., Füri, E., & Marty, B. (2016). Chondritic xenon in the Earth's mantle. *Nature*, 533(7601), 82–85. https://doi.org/10.1038/nature17434

Cardwell, R. K., Isaacks, B. L., & Karig, D. E. (1980). The spatial distribution of earthquakes, focal mechanism solutions, and subducted lithosphere in the Philippine and northeastern Indonesian islands. *Geophysical Monograph Series*, 23, 1–35. https://doi.org/10.1029/gm023p0001

Carslaw, H. S., & Jaeger, J. C. (1959). *Conduction of heat in solids* (2nd ed.). Oxford University Press.

Chang, S.-J., Ferreira, A. M. G., & Faccenda, M. (2016). Upper- and mid-mantle interaction between the Samoan plume and the Tonga–Kermadec slabs. *Nature Communications*, 7(1), 10799. https://doi.org/10.1038/ncomms10799

Cloetingh, S., Koptev, A., Lavecchia, A., Kovács, I. J., & Beekman, F. (2022). Fingerprinting secondary mantle plumes. *Earth and Planetary Science Letters*, 597, 117819. https://doi.org/10.1016/j.epsl.2022.117819

Conrad, C. P., Bianco, T. A., Smith, E. I., & Wessel, P. (2011). Patterns of intraplate volcanism controlled by asthenospheric shear. *Nature Geoscience*, 4(5), 317–321. https://doi.org/10.1038/ngeo1111

Conrad, C. P., Wu, B., Smith, E. I., Bianco, T. A., & Tibbetts, A. (2010). Shear-driven upwelling induced by lateral viscosity variations and asthenospheric shear: A mechanism for intraplate volcanism. *Physics of the Earth and Planetary Interiors*, 178(3–4), 162–175. https://doi.org/10.1016/j.pepi.2009.10.001

Courtillot, V., Davaille, A., Besse, J., & Stock, J. (2003). Three distinct types of hotspots in the Earth's mantle. *Earth and Planetary Science Letters*, 205(3–4), 295–308. https://doi.org/10.1016/S0012-821X(02)01048-8

Dai, L., & Karato, S. (2009). Electrical conductivity of wadsleyite at high temperatures and high pressures. *Earth and Planetary Science Letters*, 287(1), 277–283. https://doi.org/10.1016/j.epsl.2009.08.012

Dannberg, J., Gassmöller, R., Li, R., Lithgow-Bertelloni, C., & Stixrude, L. (2022). An entropy method for geodynamic modelling of phase transitions: Capturing sharp and broad transitions in a multiphase assemblage. *Geophysical Journal International*, 231(3), 1833–1849. https://doi.org/10.1093/gji/ggac293

Dannberg, J., & Sobolev, S. V. (2015). Low-buoyancy thermochemical plumes resolve controversy of classical mantle plume concept. *Nature Communications*, 6(1), 6960. https://doi.org/10.1038/ncomms7960

Dasgupta, R., & Lee, C. (2024). Lithospheric weakening by a small-scale plume and its geodynamic implications. *Earth and Planetary Science Letters*, 626, 118514. https://doi.org/10.1016/j.epsl.2023.118514

Davies, D. R., Rawlinson, N., Iaffaldano, G., & Campbell, I. H. (2015). Lithospheric controls on magma composition along Earth's longest continental hotspot track. *Nature*, 525(7570), 511–514. https://doi.org/10.1038/nature14903

Davis, R. O., & Selvadurai, A. P. S. (2002). *Plasticity and geomechanics*. Cambridge University Press. https://doi.org/10.1017/CBO9780511614958

Defant, M. J., Jacques, D., Maury, R. C., de Boer, J., & Joron, J.-L. (1989). Geochemistry and tectonic setting of the Luzon arc, Philippines. *Geological Society of America Bulletin*, 101(5), 663–672. https://doi.org/10.1130/0016-7606(1989)101<0663:gatsot>2.3.co;2

Di Giuseppe, E., Faccenna, C., Funiciello, F., Van Hunen, J., & Giardini, D. (2009). On the relation between trench migration, seafloor age, and the strength of the subducting lithosphere. *Lithosphere*, 1(2), 121–128. https://doi.org/10.1130/L26.1

Druken, K. A., Kincaid, C., Griffiths, R. W., Stegman, D. R., & Hart, S. R. (2014). Plume–slab interaction: The Samoa–Tonga system. *Physics of the Earth and Planetary Interiors*, 232, 1–14. https://doi.org/10.1016/j.pepi.2014.03.003

Duarte, J. C., Schellart, W. P., & Cruden, A. R. (2015). How weak is the subduction zone interface? *Geophysical Research Letters*, 42(8), 2664–2673. https://doi.org/10.1002/2014GL062876

Earle, S. (2015). *Physical geology*. BCcampus. Retrieved from https://openlibrary-repo.ecampusontario.ca/jspui/handle/123456789/402

Elkins-Tanton, L. T. (2007). Continental magmatism, volatile recycling, and a heterogeneous mantle caused by lithospheric gravitational instabilities. *Journal of Geophysical Research*, 112(B3). https://doi.org/10.1029/2005JB004072

El-Sharkawy, A., Hansteen, T. H., Clemente-Gomez, C., Fullea, J., Lebedev, S., & Meier, T. (2024). Cenozoic volcanic provinces and shallow asthenospheric volumes in the Circum-Mediterranean: Evidence from magmatic geochemistry, seismic tomography, and integrated geophysical-petrological thermochemical inversion. *Geochemistry, Geophysics, Geosystems*, 25(9), e2023GC011017. https://doi.org/10.1029/2023GC011017

Erdös, Z., Huismans, R. S., Faccenna, C., & Wolf, S. G. (2021). The role of subduction interface and upper plate strength on back-arc extension: Application to mediterranean back-arc basins. *Tectonics*, 40(8), e2021TC006795. https://doi.org/10.1029/2021TC006795

Faccenna, M., & Dal Zilio, L. (2017). The role of solid–solid phase transitions in mantle convection. *Lithos*, 268–271, 198–224. https://doi.org/10.1016/j.lithos.2016.11.007

Faccenna, C., Becker, T. W., Lallemand, S., Lagabrielle, Y., Funiciello, F., & Piromallo, C. (2010). Subduction-triggered magmatic pulses: A new class of plumes? *Earth and Planetary Science Letters*, 299(1–2), 54–68. https://doi.org/10.1016/j.epsl.2010.08.012

Faccenna, C., Becker, T. W., Lucente, F. P., Jolivet, L., & Rossetti, F. (2001). History of subduction and back arc extension in the Central Mediterranean. *Geophysical Journal International*, 145(3), 809–820. https://doi.org/10.1046/j.0956-540x.2001.01435.x
Faccenna, C., Jolivet, L., Piromallo, C., & Morelli, A. (2003). Subduction and the depth of convection in the Mediterranean mantle. *Journal of Geophysical Research, 108*(B2). https://doi.org/10.1029/2001JB001690

Fekiacova, Z., Mertz, D. F., & Hofmann, A. W. (2007). Geodynamic setting of the Tertiary Hocheifel volcanism (Germany), Part II: Geochemistry and Sr, Nd and Pb isotopic compositions. In J. R. R. Ritter & U. R. Christensen (Eds.), *Mantle plumes: A multidisciplinary approach* (pp. 207–239). Springer. https://doi.org/10.1007/978-3-540-68046-8_7

Fekiacova, Z., Mertz, D. F., & Renne, P. R. (2007). Geodynamic setting of the Tertiary Hocheifel volcanism (Germany), Part I: 40Ar/39Ar geochronology. In J. R. R. Ritter & U. R. Christensen (Eds.), *Mantle plumes: A multidisciplinary approach* (pp. 185–206). Springer. https://doi.org/10.1007/978-3-540-68046-8_6

Filippo, Z., & Faccenda, M. (2025). *Interaction between subducting slabs and an anomalous lower mantle as source of intraplate volcanism* (Master thesis). University of Padova. Retrieved from https://hdl.handle.net/20.500.12608/92416

Foulger, G. R. (2002). Plumes, or plate tectonic processes? *Astronomy and Geophysics, 43*(6), 6.19–6.24. https://doi.org/10.1046/j.1468-4004.2002.43619.x

French, S. W., & Romanowicz, B. (2015). Broad plumes rooted at the base of the Earth's mantle beneath major hotspots. *Nature, 525*(5767), 95–99. https://doi.org/10.1038/nature14876

Fukao, Y., & Obayashi, M. (2013). Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity. *Journal of Geophysical Research: Solid Earth, 118*(11), 5920–5938. https://doi.org/10.1002/2013JB010466

Fukao, Y., Obayashi, M., Inoue, H., & Nenbai, M. (1992). Subducting slabs stagnant in the mantle transition zone. *Journal of Geophysical Research, 97*(B4), 4809–4822. https://doi.org/10.1029/91JB02749

Fukao, Y., Obayashi, M., Nakakuki, T., & Deep Slab Project Group. (2009). Stagnant slab: A review. *Annual Review of Earth and Planetary Sciences, 37*(1), 19–46. https://doi.org/10.1146/annurev-earth.36.031207.124224

Funiciello, F., Faccenna, C., Heuret, A., Lallemand, S., Di Giuseppe, E., & Becker, T. W. (2008). Trench migration, net rotation and slab–mantle coupling. *Earth and Planetary Science Letters, 271*(1–4), 233–240. https://doi.org/10.1016/j.epsl.2008.04.006

Glerum, A., Thieulot, C., Fraters, M., Blom, C., & Spakman, W. (2018). Nonlinear viscoplasticity in ASPECT: Benchmarking and applications to subduction. *Solid Earth, 9*(2), 267–294. https://doi.org/10.5194/se-9-267-2018

Goes, S., Agrusta, R., van Hunen, J., & Garel, F. (2017). Subduction-transition zone interaction: A review. *Geosphere, 13*(3), 644–664. https://doi.org/10.1130/GES01476.1

Goes, S., Spakman, W., & Bijwaard, H. (1999). A lower mantle source for Central European volcanism. *Science, 286*(5446), 1928–1931. https://doi.org/10.1126/science.286.5446.1928

Granet, M., Wilson, M., & Achauer, U. (1995). Imaging a mantle plume beneath the French Massif Central. *Earth and Planetary Science Letters, 136*(3–4), 281–296. https://doi.org/10.1016/0012-821X(95)00174-B

Gu, X.-Y., Wang, P.-Y., Kuritani, T., Hanski, E., Xia, Q.-K., & Wang, Q.-Y. (2019). Low water content in the mantle source of the Hainan plume as a factor inhibiting the formation of a large igneous province. *Earth and Planetary Science Letters, 515*, 221–230. https://doi.org/10.1016/j.epsl.2019.03.034

Hall, R. (2002). Cenozoic geological and plate tectonic evolution of SE Asia and the SW Paci®c: Computer-based reconstructions, model and animations. *Journal of Asian Earth Sciences*.

Hall, R., & Spakman, W. (2015). Mantle structure and tectonic history of SE Asia. *Tectonophysics, 658*, 14–45. https://doi.org/10.1016/j.tecto.2015.07.003

Hamburger, M. W., Cardwell, R. K., & Isacks, B. L. (1983). Seismotectonics of the northern Philippine island arc. *Geophysical Monograph Series, 27*, 1–22. https://doi.org/10.1029/gm027p0001

Handy, M. R., M. Schmid, S., Bousquet, R., Kissling, E., & Bernoulli, D. (2010). Reconciling plate-tectonic reconstructions of Alpine Tethys with the geological–geophysical record of spreading and subduction in the Alps. *Earth-Science Reviews, 102*(3–4), 121–158. https://doi.org/10.1016/j.earscirev.2010.06.002

Handy, M. R., Schmid, S. M., Paffrath, M., & Friederich, W., & the AlpArray Working Group. (2021). Orogenic lithosphere and slabs in the greater Alpine area – Interpretations based on teleseismic P-wave tomography. *Solid Earth, 12*(11), 2633–2669. https://doi.org/10.5194/se-12-2633-2021

Handy, M. R., Ustaszewski, K., & Kissling, E. (2015). Reconstructing the Alps–Carpathians–Dinarides as a key to understanding switches in subduction polarity, slab gaps and surface motion. *International Journal of Earth Sciences, 104*(1), 1–26. https://doi.org/10.1007/s00531-014-1060-3

He, Y., & Wen, L. (2011). Seismic velocity structures and detailed features of the D" discontinuity near the core-mantle boundary beneath eastern Eurasia. *Physics of the Earth and Planetary Interiors, 189*(3–4), 176–184. https://doi.org/10.1016/j.pepi.2011.09.002

Heilman, E., & Becker, T. W. (2022). Plume-Slab interactions can shut off subduction. *Geophysical Research Letters, 49*(13), e2022GL099286. https://doi.org/10.1029/2022GL099286

Heilman, E., & Becker, T. W. (2024). Plume-Driven subduction termination in 3-D mantle convection models. *Geochemistry, Geophysics, Geosystems, 25*(7), e2024GC011523. https://doi.org/10.1029/2024GC011523

Heister, T., Dannberg, J., Gassmöller, R., & Bangerth, W. (2017). High accuracy mantle convection simulation through modern numerical methods – II: Realistic models and problems. *Geophysical Journal International, 210*(2), 833–851. https://doi.org/10.1093/gji/ggx195

Heron, P. J., Gün, E., Shephard, G. E., Dannberg, J., Gassmöller, R., Martin, E., et al. (2024). The role of subduction in the formation of Pangaean oceanic large igneous provinces. *Geological Society, London, Special Publications, 542*(1), 105–128. https://doi.org/10.1144/SP542-2023-12

Heron, P. J., Lowman, J. P., & Stein, C. (2015). Influences on the positioning of mantle plumes following supercontinent formation. *Journal of Geophysical Research: Solid Earth, 120*(5), 3628–3648. https://doi.org/10.1002/2014JB011727

Hirth, G., & Kohlstedt, D. (2003). Rheology of the upper mantle and the mantle wedge: A view from the experimentalists. In J. Eiler (Ed.), *Geophysical monograph series* (Vol. 138, pp. 83–105). American Geophysical Union. https://doi.org/10.1029/138GM06

Ho, K., Chen, J., & Juang, W. (2000). Geochronology and geochemistry of late Cenozoic basalts from the Leiqiong area, southern China. *Journal of Asian Earth Sciences, 18*(3), 307–324. https://doi.org/10.1016/S1367-9120(99)00059-0

Hoang, N., Flower, M. F. J., & Carlson, R. W. (1996). Major, trace element, and isotopic compositions of Vietnamese basalts: Interaction of hydrous EM1-rich asthenosphere with thinned Eurasian lithosphere. *Geochimica et Cosmochimica Acta, 60*(22), 4329–4351. https://doi.org/10.1016/S0016-7037(96)00247-5

Huang, J., & Zhao, D. (2006). High-resolution mantle tomography of China and surrounding regions. *Journal of Geophysical Research, 111*(B9). https://doi.org/10.1029/2005JB004066

Huang, Z., Zhao, D., & Wang, L. (2015). P wave tomography and anisotropy beneath Southeast Asia: Insight into mantle dynamics. *Journal of Geophysical Research: Solid Earth, 120*(7), 5154–5174. https://doi.org/10.1002/2015JB012098
Karato, S. (2011). Water distribution across the mantle transition zone and its implications for global material circulation. *Earth and Planetary Science Letters, 301*(3–4), 413–423. https://doi.org/10.1016/j.epsl.2010.11.038

Karig, D. E. (1975). Basin genesis in the Philippine Sea. Initial Reports of the Deep Sea Drilling Project (Vol. 31, pp. 857–879).

Kerr, M. C., & Stegman, D. R. (2024). Mantle avalanches in a Venus-like stagnant lid planet. *Earth and Planetary Physics, 8*(5), 686–702. https://doi.org/10.26464/epp2024062

Kim, H., Lee, Y., Kim, D., & Lee, C. (2022). Roles of partial hydration of mantle transition zone and wadsleyite-olivine phase transition in wet plume development: Origin of Quaternary intraplate volcanoes in Korea. *Tectonophysics, 837*, 229467. https://doi.org/10.1016/j.tecto.2022.229467

Kincaid, C., Druken, K. A., Griffiths, R. W., & Stegman, D. R. (2013). Bifurcation of the Yellowstone plume driven by subduction-induced mantle flow. *Nature Geoscience, 6*(5), 395–399. https://doi.org/10.1038/ngeo1774

King, S. D. (2007). Hotspots and edge-driven convection. *Geology, 35*(3), 223–226. https://doi.org/10.1130/G23291A.1

King, S. D., & Adam, C. (2014). Hotspot swells revisited. *Physics of the Earth and Planetary Interiors, 235*, 66–83. https://doi.org/10.1016/j.pepi.2014.07.006

King, S. D., & Anderson, D. L. (1995). An alternative mechanism of flood basalt formation. *Earth and Planetary Science Letters, 136*(3–4), 269–279. https://doi.org/10.1016/0012-821X(95)00205-Q

King, S. D., & Anderson, D. L. (1998). Edge-driven convection. *Earth and Planetary Science Letters, 160*(3–4), 289–296. https://doi.org/10.1016/S0012-821X(98)00089-2

Koptev, A., Cloetingh, S., & Ehlers, T. A. (2021). Longevity of small-scale ('baby') plumes and their role in lithospheric break-up. *Geophysical Journal International, 227*(1), 439–471. https://doi.org/10.1093/gji/ggab223

Kronbichler, M., Heister, T., & Bangerth, W. (2012). High accuracy mantle convection simulation through modern numerical methods. *Geophysical Journal International, 191*(1), 12–29. https://doi.org/10.1111/j.1365-246X.2012.05609.x

Lallemand, S., Font, Y., Bijwaard, H., & Kao, H. (2001). New insights on 3-D plates interaction near Taiwan from tomography and tectonic implications. *Tectonophysics, 335*(3–4), 229–253. https://doi.org/10.1016/S0040-1951(01)00071-3

Le Breton, E., Brune, S., Ustaszewski, K., Zahirovic, S., Seton, M., & Müller, R. D. (2021). Kinematics and extent of the Piemont–Liguria Basin – Implications for subduction processes in the Alps. *Solid Earth, 12*(4), 885–913. https://doi.org/10.5194/se-12-885-2021

Lee, C., Cerpa, N. G., Lee, J., & Huh, M. (2025). Intraplate volcanism induced by magmatic solitary waves ascending from the mantle transition zone. *Geology, 40*(12), 538001.

Lei, J., Zhao, D., Steinberger, B., Wu, B., Shen, F., & Li, Z. (2009). New seismic constraints on the upper mantle structure of the Hainan plume. *Physics of the Earth and Planetary Interiors, 173*(1–2), 33–50. https://doi.org/10.1016/j.pepi.2008.10.013

Leonard, T., & Liu, L. (2016). The role of a mantle plume in the formation of Yellowstone volcanism. *Geophysical Research Letters, 43*(3), 1132–1139. https://doi.org/10.1002/2015GL067131

Li, C., van der Hilst, R. D., Engdahl, E. R., & Burdick, S. (2008). A new global model for P wave speed variations in Earth's mantle. *Geochemistry, Geophysics, Geosystems, 9*(5). https://doi.org/10.1029/2007GC001806

Li, C.-F., Xu, X., Lin, J., Sun, Z., Zhu, J., Yao, Y., et al. (2014). Ages and magnetic structures of the South China Sea constrained by deep tow magnetic surveys and IODP Expedition 349. *Geochemistry, Geophysics, Geosystems, 15*(12), 4958–4983. https://doi.org/10.1002/2014GC005567

Li, R., Dannberg, J., Gassmöller, R., Lithgow-Bertelloni, C., & Stixrude, L. (2025). How phase transitions impact changes in mantle convection style throughout Earth's history: From stalled plumes to surface dynamics. *Geochemistry, Geophysics, Geosystems, 26*(2), e2024GC011600. https://doi.org/10.1029/2024GC011600

Li, Y. (2025). yingyil06/Paper-Li-Intraplate-volcanism-driven-by-slab-plume-interaction: V1.1.0 [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.16884447

Liu, H., & Leng, W. (2020). Plume-Tree structure induced by low-viscosity layers in the upper mantle. *Geophysical Research Letters, 47*(1), e2019GL086508. https://doi.org/10.1029/2019GL086508

Long, X., Ballmer, M. D., Córdoba, A. M.-C., & Li, C.-F. (2019). Mantle melting and intraplate volcanism due to self-buoyant hydrous upwellings from the stagnant slab that are conveyed by small-scale convection. *Geochemistry, Geophysics, Geosystems, 20*(11), 4972–4997. https://doi.org/10.1029/2019GC008591

Mériaux, C. A., Duarte, J. C., Schellart, W. P., & Mériaux, A.-S. (2015). A two-way interaction between the Hainan plume and the Manila subduction zone. *Geophysical Research Letters, 42*(14), 5796–5802. https://doi.org/10.1002/2015GL064313

Mériaux, C. A., Mériaux, A.-S., Schellart, W. P., Duarte, J. C., Duarte, S. S., & Chen, Z. (2016). Mantle plumes in the vicinity of subduction zones. *Earth and Planetary Science Letters, 454*, 166–177. https://doi.org/10.1016/j.epsl.2016.09.001

Michon, L., & Merle, O. (2001). The evolution of the Massif Central Rift; spatio-temporal distribution of the volcanism. *Bulletin de la Societe Geologique de France, 172*(2), 201–211. https://doi.org/10.2113/172.2.201

Montelli, R., Nolet, G., Dahlen, F. A., & Masters, G. (2006). A catalogue of deep mantle plumes: New results from finite-frequency tomography. *Geochemistry, Geophysics, Geosystems, 7*(11). https://doi.org/10.1029/2006GC001248

Morgan, W. J. (1971). Convection plumes in the lower mantle. *Nature, 230*(5288), 42–43. https://doi.org/10.1038/230042a0

Negredo, A. M., van Hunen, J., Rodríguez-González, J., & Fullea, J. (2022). On the origin of the Canary Islands: Insights from mantle convection modelling. *Earth and Planetary Science Letters, 584*, 117506. https://doi.org/10.1016/j.epsl.2022.117506

Obayashi, M., Yoshimitsu, J., Nolet, G., Fukao, Y., Shiobara, H., Sugioka, H., et al. (2013). Finite frequency whole mantle P wave tomography: Improvement of subducted slab images. *Geophysical Research Letters, 40*(21), 5652–5657. https://doi.org/10.1002/2013GL057401

Ozawa, A., Tagami, T., Listanco, E. L., Arpa, C. B., & Sudo, M. (2004). Initiation and propagation of subduction along the Philippine Trench: Evidence from the temporal and spatial distribution of volcanoes. *Journal of Asian Earth Sciences, 23*(1), 105–111. https://doi.org/10.1016/S1367-9120(03)00112-3

Piromallo, C., & Morelli, A. (2003). P wave tomography of the mantle under the Alpine-Mediterranean area. *Journal of Geophysical Research, 108*(B2). https://doi.org/10.1029/2002JB001757

Quinquis, M. E. T., & Buiter, S. J. H. (2014). Testing the effects of basic numerical implementations of water migration on models of subduction dynamics. *Solid Earth, 5*(1), 537–555. https://doi.org/10.5194/se-5-537-2014

Ranalli, G. (1995). *Rheology of the Earth*. Springer Science & Business Media.

Rangin, C., Spakman, W., Pubellier, M., & Bijwaard, H. (1999). Tomographic and geological constraints on subduction along the eastern Sundaland continental margin (South-East Asia). *Bulletin de la Societe Geologique de France, 170*, 775–788.

Richard, G. C., & Bercovici, D. (2009). Water-induced convection in the Earth's mantle transition zone. *Journal of Geophysical Research, 114*(B1). https://doi.org/10.1029/2008JB005734
Richard, G. C., & Iwamori, H. (2010). Stagnant slab, wet plumes and Cenozoic volcanism in East Asia. *Physics of the Earth and Planetary Interiors, 183*(1-2), 280-287. https://doi.org/10.1016/j.pepi.2010.02.009

Ritter, J. R. R. (2007). The seismic signature of the Eifel plume. In J. R. R. Ritter & U. R. Christensen (Eds.), *Mantle plumes: A multidisciplinary approach* (pp. 379-404). Springer. https://doi.org/10.1007/978-3-540-68046-8_12

Ritter, J. R. R., Jordan, M., Christensen, U. R., & Achauer, U. (2001). A mantle plume below the Eifel volcanic fields, Germany. *Earth and Planetary Science Letters, 186*(1), 7-14. https://doi.org/10.1016/S0012-821X(01)00226-6

Romano, C., Poe, B. T., Tyburczy, J., & Nestola, F. (2009). Electrical conductivity of hydrous wadsleyite. *European Journal of Mineralogy, 21*(3), 615-622. https://doi.org/10.1127/0935-1221/2009/0021-1933

Roy, P., Steinberger, B., Faccenda, M., & Pons, M. (2025). Modeling anisotropic signature of slab-induced mantle plumes from thermochemical piles in the lowermost mantle. *Geophysical Research Letters, 52*(10), e2024GL113299. https://doi.org/10.1029/2024GL113299

Ryan, W. B. F., Carbotte, S. M., Coplan, J. O., O'Hara, S., Melkonian, A., Arko, R., et al. (2009). Global Multi-Resolution Topography synthesis. *Geochemistry, Geophysics, Geosystems, 10*(3). https://doi.org/10.1029/2008GC002332

Schellart, W. P., Stegman, D. R., Farrington, R. J., & Moresi, L. (2011). Influence of lateral slab edge distance on plate velocity, trench velocity, and subduction partitioning. *Journal of Geophysical Research, 116*(B10), B10408. https://doi.org/10.1029/2011JB008535

Schubert, G., Turcotte, D. L., & Olson, P. (2001). *Mantle convection in the Earth and planets*. Cambridge University Press.

Sobolev, S. V., Zeyen, H., Granet, M., Achauer, U., Bauer, C., Werling, F., et al. (1997). Upper mantle temperatures and lithosphere-asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations. *Tectonophysics, 275*(1-3), 143-164. https://doi.org/10.1016/S0040-1951(97)00019-X

Steinberger, B., & Calderwood, A. R. (2006). Models of large-scale viscous flow in the Earth's mantle with constraints from mineral physics and surface observations. *Geophysical Journal International, 167*(3), 1461-1481. https://doi.org/10.1111/j.1365-246X.2006.03131.x

Steinberger, B., & Holme, R. (2008). Mantle flow models with core-mantle boundary constraints and chemical heterogeneities in the lowermost mantle. *Journal of Geophysical Research, 113*(B5). https://doi.org/10.1029/2007JB005080

Steinberger, B., Roy, P., Pons, M., & Jopke, M. P. (2025). Why are plume excess temperatures much less than the temperature drop across the lowermost-mantle thermal boundary layer? *Journal of Geophysical Research: Solid Earth, 130*(4), e2024JB030111. https://doi.org/10.1029/2024JB030111

Strak, V., & Schellart, W. P. (2018). A subduction and mantle plume origin for Samoan volcanism. *Scientific Reports, 8*(1), 10424. Article 1. https://doi.org/10.1038/s41598-018-28267-3

Tan, E., Gurnis, M., & Han, L. (2002). Slabs in the lower mantle and their modulation of plume formation. *Geochemistry, Geophysics, Geosystems, 3*(11), 1-24. https://doi.org/10.1029/2001GC000238

Taylor, B., & Hayes, D. E. (1983). Origin and history of the South China Sea basin. In D. E. Hayes (Ed.), *Geophysical monograph series* (Vol. 27, pp. 23-56). American Geophysical Union. https://doi.org/10.1029/GM027p0023

Torii, Y., & Yoshioka, S. (2007). Physical conditions producing slab stagnation: Constraints of the Clapeyron slope, mantle viscosity, trench retreat, and dip angles. *Tectonophysics, 445*(3), 200-209. https://doi.org/10.1016/j.tecto.2007.08.003

Turcotte, D. L., & Oxburgh, E. R. (1967). Finite amplitude convective cells and continental drift. *Journal of Fluid Mechanics, 28*(1), 29-42. https://doi.org/10.1017/S0022112067001880

van den Berg, A. P., van Keken, P. E., & Yuen, D. A. (1993). The effects of a composite Non-Newtonian and Newtonian rheology on mantle convection. *Geophysical Journal International, 115*(1), 62-78. https://doi.org/10.1111/j.1365-246x.1993.tb05588.x

Van der Voo, R., Spakman, W., & Bijwaard, H. (1999). Tethyan subducted slabs under India. *Earth and Planetary Science Letters, 171*(1), 7-20. https://doi.org/10.1016/S0012-821X(99)00131-4

Van Hinsbergen, D. J. J., Torsvik, T. H., Schmid, S. M., Mañenco, L. C., Maffione, M., Vissers, R. L. M., et al. (2020). Orogenic architecture of the Mediterranean region and kinematic reconstruction of its tectonic evolution since the Triassic. *Gondwana Research, 81*, 79-229. https://doi.org/10.1016/j.gr.2019.07.009

Wang, K., Li, K., Cai, K., Wang, X., & Hu, J. (2025). Subduction-stalled plume tail triggers Tarim large igneous province. *Earth and Planetary Science Letters, 651*, 119150. https://doi.org/10.1016/j.epsl.2024.119150

Wang, X.-C., Li, Z.-X., Li, X.-H., Li, J., Xu, Y.-G., & Li, X.-H. (2013). Identification of an ancient mantle reservoir and young recycled materials in the source region of a young mantle plume: Implications for potential linkages between plume and plate tectonics. *Earth and Planetary Science Letters, 377-378*, 248-259. https://doi.org/10.1016/j.epsl.2013.07.003

Wei, W., Xu, J., Zhao, D., & Shi, Y. (2012). East Asia mantle tomography: New insight into plate subduction and intraplate volcanism. *Journal of Asian Earth Sciences, 60*, 88-103. https://doi.org/10.1016/j.jseas.2012.08.001

Wortel, M. J. R., & Spakman, W. (2000). Subduction and slab detachment in the Mediterranean-Carpathian Region. *Science, 290*(5498), 1910-1917. https://doi.org/10.1126/science.290.5498.1910

Wu, J., & Suppe, J. (2018). Proto-South China Sea plate tectonics using subducted slab constraints from tomography. *Journal of Earth Sciences, 29*(6), 1304-1318. https://doi.org/10.1007/s12583-017-0813-x

Xia, S., Zhao, D., Sun, J., & Huang, H. (2016). Tele seismic imaging of the mantle beneath southernmost China: New insights into the Hainan plume. *Gondwana Research, 36*, 46-56. https://doi.org/10.1016/j.gr.2016.05.003

Xue, K., Schellart, W. P., & Strak, V. (2020). Effect of plate length on subduction kinematics and slab geometry: Insights from buoyancy-driven analog subduction models. *Journal of Geophysical Research: Solid Earth, 125*(11), e2020JB020514. https://doi.org/10.1029/2020JB020514

Yan, Q., Shi, X., Metcalfe, I., Liu, S., Xu, T., Kornkanitnan, N., et al. (2018). Hainan mantle plume produced late Cenozoic basaltic rocks in Thailand, Southeast Asia. *Scientific Reports, 8*(1), 2640. https://doi.org/10.1038/s41598-018-20712-7

Yang, J., & Faccenda, M. (2020). Intraplate volcanism originating from upwelling hydrous mantle transition zone. *Nature, 579*(7797), 88-91. https://doi.org/10.1038/s41586-020-2045-y

Yin, A. (2010). Cenozoic tectonic evolution of Asia: A preliminary synthesis. *Tectonophysics, 488*(1), 293-325. https://doi.org/10.1016/j.tecto.2009.06.002

Yoshino, T., & Katsura, T. (2012). Re-evaluation of electrical conductivity of anhydrous and hydrous wadsleyite. *Earth and Planetary Science Letters, 337-338*, 56-67. https://doi.org/10.1016/j.epsl.2012.05.023

Yoshino, T., & Katsura, T. (2013). Electrical conductivity of mantle minerals: Role of water in conductivity anomalies. *Annual Review of Earth and Planetary Sciences, 41*(1), 605-628. https://doi.org/10.1146/annurev-earth-050212-124022

Young, A., Flament, N., Maloney, K., Williams, S., Matthews, K., Zahirovic, S., & Müller, R. D. (2019). Global kinematics of tectonic plates and subduction zones since the late Paleozoic Era. *Geoscience Frontiers, 10*(3), 989-1013. https://doi.org/10.1016/j.gsf.2018.05.011

Yu, X., Liu, Z., Zeng, G., Cao, W., Meas, R., Hoang, L. V., & Sang, P. N. (2022). Mantle plume-stagnant slab interaction controls the generation of a mixed mantle source for continental intraplate basalts. *Lithos, 426-427*, 106795. https://doi.org/10.1016/j.lithos.2022.106795
Zahirovic, S., Matthews, K. J., Flament, N., Müller, R. D., Hill, K. C., Seton, M., & Gurnis, M. (2016). Tectonic evolution and deep mantle structure of the eastern Tethys since the latest Jurassic. *Earth-Science Reviews*, **162**, 293–337. https://doi.org/10.1016/j.earscirev.2016.09.005

Zhang, N., & Li, Z.-X. (2018). Formation of mantle “lone plumes” in the global downwelling zone—A multiscale modelling of subduction-controlled plume generation beneath the South China Sea. *Tectonophysics*, **723**, 1–13. https://doi.org/10.1016/j.tecto.2017.11.038

Zhao, D. (2007). Seismic images under 60 hotspots: Search for mantle plumes. *Gondwana Research*, **12**(4), 335–355. https://doi.org/10.1016/j.gr.2007.03.001

Zhao, D., Toyokuni, G., & Kurata, K. (2021). Deep mantle structure and origin of Cenozoic intraplate volcanoes in Indochina, Hainan and South China Sea. *Geophysical Journal International*, **225**(1), 572–588. https://doi.org/10.1093/gji/ggaa605

Zhao, D., Yu, S., & Ohtani, E. (2011). East Asia: Seismotectonics, magmatism and mantle dynamics. *Journal of Asian Earth Sciences*, **40**(3), 689–709. https://doi.org/10.1016/j.jseaes.2010.11.013

Zhu, H., Bozdağ, E., Peter, D., & Tromp, J. (2012). Structure of the European upper mantle revealed by adjoint tomography. *Nature Geoscience*, **5**(7), 493–498. https://doi.org/10.1038/ngeo1501