The role of edge-driven convection in the generation ofvolcanism – Part 2: Interaction with mantle plumes, applied to the Canary Islands

In the eastern Atlantic Ocean, several volcanic archipelagos are located close to the margin of the African continent. This configuration has inspired previous studies to suggest an important role of edge-driven convection (EDC) in the generation of intraplate magmatism. In a companion paper (Manjón-Cabeza Córdoba and Ballmer, 2021), we showed that EDC alone is insufficient to sustain magmatism of the magnitude required to match the volume of these islands. However, we also found that EDC readily develops near a step of lithospheric thickness, such as the oceanic–continental transition (“edge”) along the western African cratonic margin. In this work, we carry out 3D numerical models of mantle flow and melting to explore the possible interactions between EDC and mantle plumes. We find that the stem of a plume that rises close to a lithospheric edge is significantly deflected ocean-ward (i.e., away from the edge). The pancake of ponding hot material at the base of the lithosphere is also deflected by the EDC convection cell (either away or towards the edge). The amount of magmatism and plume deflection depends on the initial geometric configuration, i.e., the distance of the plume from the edge. Plume buoyancy flux and temperature also control the amount of magmatism, and influence the style and extent of plume–EDC interaction. Finally, comparison of model predictions with observations reveals that the Canary plume may be significantly affected and deflected by EDC, accounting for widespread and coeval volcanic activity. Our work shows that many of the peculiar characteristics of eastern Atlantic volcanism are compatible with mantle plume theory once the effects of EDC on plume flow are considered.

1 Introduction

Volcanism exerts a major control for material flux between the interior of the Earth and the surface–atmosphere system. Volcanic activity along mid-ocean ridges and subduction zones is readily explained by plate tectonics. However, in the absence of nearby plate boundaries, plate tectonics cannot account for intraplate volcanism.

Several models have been proposed to explain the origin of such magmatism. The leading hypothesis is mantle plume theory, in which a deep columnar thermal anomaly rises from the core–mantle boundary to the base of the lithosphere in order to support localized hotspot volcanism (Wilson, 1963; Morgan, 1971). Still, several predictions of plume theory are not fulfilled at many locations worldwide (e.g., linear age progressions consistent with plate velocity; Courtillot et al., 2003) and other models have been put forward: small-scale convection (SSC; Richter, 1973; Parsons and McKenzie, 1978; Huang et al., 2003; Dumoulin et al., 2005; Ballmer et al., 2007), shear-driven upwelling (SDU; Conrad et al., 2010), or edge-driven convection (EDC; King and Anderson, 1995, 1998).

In the eastern Atlantic, several volcanic archipelagos are located on the ocean floor near continental lithosphere. At these locations, many of the predictions of plume theory are not met. In the Canary Islands (where volcanism is as recent as the 2021 eruption of La Palma) volcano ages do not follow a consistent linear age–distance relationship, with coeval volcanism occurring across several hundreds of kilometers (Abdel-Monem et al., 1971, 1972; Thirlwall et al., 2000; Geldmacher et al., 2005), the plume swell is nearly absent (Sleep, 1990; King and Adam, 2014) (although see Huppert et al., 2020), and the duration of volcanism at a single island is longer than expected in comparison with other chains (e.g., Carracedo, 1999). Besides, a cogenetic relation of these volcanoes with Alboran Domain volcanism has been suggested due to tectonism (Doblas et al., 2007) and with the northwest Africa Cenozoic volcanism as part of the same upwelling (Duggen et al., 2009). These inconsistencies have led several authors to reject the plume model for these islands (e.g., Doblas et al., 2007; Martínez-Arevalo et al., 2013). Similar arguments against the plume model have been made for Cabo Verde (King and Ritsema, 2000; Helffrich et al., 2010) or the Cameroon Volcanic Line (Fitton, 1980; Déruelle et al., 2007; Milelli et al., 2012), both of which also formed near the African continental margin.

Of the alternative models put forward to substitute mantle plume theory, EDC is the only one that has been proposed for the three aforementioned volcanic regions (King and Anderson, 1998; King and Ritsema, 2000; Milelli et al., 2012). The EDC model postulates that a convection cell is generated due to the juxtaposition of two lithospheric sections of different age or structure: the related density difference is sufficient to generate a downwelling and an associated upwelling. In theory, the return upwelling flow would be enough to generate magma to sustain ocean island volcanism, provided that the overlying lid was sufficiently thin to facilitate decompression melting. EDC has been proposed for other regions of the globe as well. For example Vogt (1991) suggested that EDC could be related to the Bermuda Rise, although several studies (Shahnas and Pysklywec, 2004; Ramsay and Pysklywec, 2011) showed that EDC would produce a maximum upwelling (and a related increase in topography) closer to the margin than the Bermuda Rise, Davies and Rawlinson (2014) have proposed EDC as the mechanism for the southeastern Australian volcanic province, and Afonso et al. (2016) have proposed EDC as an important mechanism for the central–western US volcanism.

Nonetheless, in a previous paper (Manjón-Cabeza Córdoba and Ballmer, 2021), we quantitatively tested the hypothesis of edge-driven convection as an origin of oceanic intraplate volcanism near continental margins, and our results showed that, by itself, EDC can only support minor magmatism even under the most favorable conditions and is clearly insufficient to generate long-lived island-building volcanism. While other studies have shown that a very steep oceanic–continental transition (Kim and So, 2020; Negredo et al., 2022) or additional geometrical complexities (Duvernay et al., 2021) could increase the amount of EDC-related melting calculated by the companion study, all of them agree that magmatism is very restricted to account for volcanism in the Canary Islands. Furthermore, we speculated that due to the prevalence of EDC with Earth-like mantle properties, most of EDC-related flow and melting should occur near mid-ocean ridges in young lithospheres, which was previously observed by geological (Ligi et al., 2011) and geodynamic (Buck, 1986; Boutilier and Keen, 1999; Sleep, 2007) studies alike, and not in old lithospheres (as is the case for the Canary Islands).

On the other hand, recent seismic-tomography studies provide evidence of deeply rooted mantle plumes in the eastern Atlantic by imaging continuous near-vertical low-velocity anomalies in the mantle (French and Romanowicz, 2015) or broad upwellings just below these archipelagos (Civiero et al., 2021). In addition, additional geophysical evidence points to the presence of thermal upwellings (plumes) at least from the base of the transition zone (Liu and Zhao, 2014; Saki et al., 2015).

In the light of the evidence gathered along these lines, we here explore the dynamics of mantle flow and melting related to plumes that rise near a continental margin (Fig. 1). We hypothesize that the interaction between plumes and EDC can explain (at least some of) the discrepancies between the predictions of plume theory and observations, as already suggested by Geldmacher et al. (2005). To study the interaction between plumes and EDC, we carry out three-dimensional (3D) numerical models of flow and melting near the transition between the oceanic and the continental lithosphere. Negredo et al. (2022) carried out a preliminary study of plumes and their potential interactions with EDC in 2D, finding that EDC could be responsible for plume migration in the Canaries. However, plumes are inherently 3D, and including this third dimension allows us to include the effects of plate velocity. We explore the parameters that control plume flow (e.g., plume buoyancy flux, plume excess temperature) and EDC (e.g., mantle viscosity, distance of the plume from the continental margin).

Figure 1 Schematic of a plume interacting with edge-driven convection. In this work, we build on the models of Manjón-Cabeza Córdoba and Ballmer (2021) and add a plume in the form of a temperature anomaly at the bottom and a plate velocity v_plate=2 cm yr⁻¹ consistent with the age–distance track to the north of the Canary Islands. An idealized lithosphere–asthenosphere boundary (LAB) is labeled in the figure.

2 Methods

We run 3D Cartesian numerical models using the same version of the finite-element code CITCOM (Moresi and Solomatov, 1995; Moresi and Gurnis, 1996; Zhong et al., 2000) as in our previous paper (Manjón-Cabeza Córdoba and Ballmer, 2021). The conservation equations of mass, momentum, and energy are solved on the finite-element mesh according to the “extended Boussinesq approximation” (Christensen and Yuen, 1985) (although we impose the adiabat as a constant gradient with depth of 0.3 K km⁻¹); composition is tracked by passive Lagrangian particles (i.e., tracers). Three-dimensional geometry of the model box is chosen due to the intrinsic 3D nature of the problem (see Fig. 1) and the related complex flow patterns. To make our models comparable with the 2D cases in the companion paper (Manjón-Cabeza Córdoba and Ballmer, 2021), we use the same model-box depth z_box and width x_box. The total extent of our computational domain is $\mathrm{2640}\times \mathrm{1980}\times \mathrm{660}$ km (x_box, y_box, and z_box, respectively). This domain is resolved by a grid of $\mathrm{384}\times \mathrm{288}\times \mathrm{96}$ elements with a uniform spacing. Resolution tests in the companion paper (Manjón-Cabeza Córdoba and Ballmer, 2021) confirmed that this is enough to accurately model EDC. Our plumes do not feature viscosities that are too low due to our low activation energies (Appendix, Fig. A2) and are similarly resolved as regional or global models (e.g., Ballmer et al., 2013; Davies and Davies, 2009) for hotter (less viscous) plumes.

Free slip is imposed at the side boundaries (x=2640 and x=0 km); no slip is imposed at the bottom. To model Atlantic plate motion and achieve a steady state for plume inflow and outflow, we impose a plate velocity at the top boundary layer parallel to the y direction (the new dimension added, parallel to the edge) of v_plate=2 cm yr⁻¹ and a related Couette flow at the inflow boundary (y=0 km) that is consistent with the viscosity profile. We acknowledge that the real absolute African plate motion could be oblique to the African margin near the Canaries today, but the volcanic track reflects a history of motion nearly parallel to the African margin (Geldmacher et al., 2005); in any case, most frames of references depict a plate movement parallel to the margin (Schellart et al., 2008; Martín et al., 2014). The corresponding outflow velocity boundary (at y=1980 km) remains unconstrained to allow free exit of material, but we impose all flow perpendicular to the boundary (that is, no slip in the directions parallel to the boundary). We also open an unconstrained circular “hole” at the bottom of the box and y=660 km to allow free inflow at the plume location (Ballmer et al., 2011).

The top boundary is fixed at T_surf=0 ^∘C, while the bottom boundary is fixed at T_ref = 1350 ^∘C (+198 ^∘C is added corresponding to the adiabatic gradient increase of 0.3 K km${}^{-\mathrm{1}}\times \mathrm{660}$ km); the x-normal boundaries are reflective. The models are internally heated as well ($H=\mathrm{7.75}\times {\mathrm{10}}^{-\mathrm{12}}$ W kg⁻¹). At the inflow boundary, the thermal distribution corresponds to the initial condition, which is identical to that of the 2D profile of the previous paper (Fig. 2 in Manjón-Cabeza Córdoba and Ballmer, 2021), including a continental “edge” at x=1320 km (Fig. 2). In nearly all cases of this study, the initial thermal age of the juxtaposed continental and oceanic lithospheres is τ_c=100 Ma and τ_o=40 Ma, respectively, except when otherwise specified. This choice of τ_o results in an age of ${\mathit{\tau }}_{\text{o},y=\mathrm{660}}=\mathrm{73}$ Ma for the oceanic lithosphere right above the plume anomaly.

Figure 2 Two-dimensional sections depicting the initial thermal (potential temperature) and compositional (depletion for the “depleted component”) profiles for the models in this work. Depletion (F) is defined as the amount of melt extracted from the mantle. Plate motion is imposed at the top boundary from left to right. The left side shows the initial thermal and inflow boundary conditions. The right side shows the initial compositional and inflow boundary conditions. The front and back sides show the thermal side boundary conditions. Solid and discontinuous arrows represent constrained (imposed) and unconstrained (open boundary) velocity boundaries, respectively. For further details on the calculations of the initial profiles, see text.

The transition between the two lithospheric thicknesses is linearly interpolated for both temperature and composition (Fig. 2). The width of the transition is w=264 km in all cases. In addition, we impose a circular plume thermal anomaly of radius r_plume centered at y=660 km (i.e., sufficiently far away to avoid artifacts due to the proximity to the inflow boundary) and variable distances from the edge $D_{\text{plume}}=\mathrm{1320}\text{km}-x_{\text{plume}}$. The plume thermal anomaly at the bottom is described by the following condition:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }T=\mathrm{\Delta }{T}_{\text{plume}}\cdot e^{\frac{-r^{\mathrm{2}}}{r_{\text{plume}}^{\mathrm{2}}}},\end{array}$$

where ΔT is the difference between the plume temperature and the background temperature and r is the distance from x_plume. Plume buoyancy flux is kept nearly constant during the simulation by automatically adjusting r_plume every 50 time steps. For example, if the plume buoyancy flux B(t) (measured at the bottom boundary) is different from the target value B_plume for a given model, r_plume is adjusted by a factor of ${\frac{B_{\text{plume}}}{B\left(t\right)}}^{\mathrm{0.5}}$. This approach keeps B(t) practically constant through much of the simulation but renders the ratio between the radius of the opening at the bottom of the model and r_plume variable between cases. Nonetheless, we make sure that this ratio remains between 3.5 and 4 for all models in the statistical steady state. In this statistical steady state, r_plume does not vary (see Supplement).

Statistical steady state is evaluated by analyzing the vVrms (root mean squares of the vertical velocity component) and the melt flux (see below). Models are only evaluated when these properties do not change systematically over time (small changes are expected due to random thermal noise). The exception to this evaluation is the models with a D_plume=400 km (see below, Supplement), which featured periodic behavior. In those cases, evaluation occurs as soon as the cycles are statistically symmetric. For examples of vVrms and melt-flux plots used to evaluate the statistical steady state, see Appendix, Fig. A1.

The mantle source consists of a mechanical mixture of three different lithological components (depleted/dry peridotite, enriched/hydrous peridotite, pyroxenite), which make up 82 %, 15 %, and 3 % of the volume of the mantle, respectively. We assume that these lithologies are in thermal equilibrium but in chemical disequilibrium due to their fine-scale nature (i.e., smaller than the finite-element mesh). Each of these lithologies has a different density and is subject to a different melting law (Katz et al., 2003; Pertermann and Hirschmann, 2003; Manjón-Cabeza Córdoba and Ballmer, 2021; Ballmer et al., 2009). Initially, the lithosphere is depleted in all of the lithologies and hence is buoyant and does not melt immediately. This lithological depletion is pre-calculated from 2D models of flow and melting at a mid-ocean ridge (Manjón-Cabeza Córdoba and Ballmer, 2021), using the same method and parameters (e.g., potential temperature and melting laws) as the current study. Such an approach restricts excess melting that may otherwise occurs due to arbitrary model choices (e.g., an increased potential temperature of the models would be met with an increased melt depletion at mid-ocean ridges, therefore limiting further melting).

Progressive melting during the simulation affects the relevant densities due to melt retention and depletion of the residue. The driving forces (density anomalies) further depend on temperature and composition. As in the companion paper, the linearized density formulation is

$$\begin{array}{}\text{(2)}& \mathit{\rho }={\mathit{\rho }}_{\text{ref}}-\mathit{\alpha }\cdot {\mathit{\rho }}_{\text{rhef}}\cdot (T-T_{\text{ref}})+F\cdot \mathrm{\Delta }{\mathit{\rho }}_F+\mathit{\varphi }\cdot \mathrm{\Delta }{\mathit{\rho }}_{\mathit{\varphi }}.\end{array}$$

This equation takes into account thermal expansivity (α) as well as the density effects of melt fraction (F) and melt retention (ϕ). For the values of the characteristic constants of these density effects (Δρ_F and Δρ_ϕ), see Table 1.

In turn, the resisting forces (viscosities) do not depend on melt retention or depletion. Our Newtonian viscosity formulation is temperature- and depth-dependent:

$$\begin{array}{}\text{(3)}& \log \mathit{\eta }=log{\mathit{\eta }}_{\mathrm{0}}+{\displaystyle \frac{E_{\text{a}}+P\cdot V_{\text{a}}}{R\cdot T}}-{\displaystyle \frac{E_{\text{a}}}{R\cdot T_{\text{ref}}}},\end{array}$$

with E_a and V_a the activation energy and volume, respectively. The chosen value for activation energy (Table 1) is smaller than values obtained for fitting experimental data (Kohlstedt and Hansen, 2015) but is useful to obtain lithospheric thicknesses close to those obtained with non-Newtonian rheologies (Christensen, 1984; van Hunen et al., 2005).

Table 1 Relevant parameters for the models described in this paper. Values outside and within parentheses provide the reference value and the explored parameter space, respectively.

3 Results

In the 2D models of the companion paper (Manjón-Cabeza Córdoba and Ballmer, 2021), we find that EDC starts right at the onset of the model evolution with a dominant downwelling below the continental side of the edge (or ocean–continent transition) and a return-flow upwelling below the oceanic side. The upwelling sustains erosion of the lithosphere, creating a “bump” or “dent” at its base. This dent is characteristic of every model in this work and also present when no plume is imposed (Appendix, Fig. A3). Ultimately, SSC also occurs at the base of the oceanic lithosphere far from the edge. We refer to SSC as a thermal boundary layer instability that (in contrast to EDC) is not immediately triggered by the presence of a nearby edge but rather typically occurs as soon as the boundary layer (nearly) reaches its critical thickness (Richter, 1973; Parsons and McKenzie, 1978).

In this study, the test cases without a mantle plume confirm that the results of Manjón-Cabeza Córdoba and Ballmer (2021) are robust and hold in our 3D geometry: EDC begins right after the material enters the model box, promoting a convection cell and related sublithospheric erosion above the upwelling on the oceanic side. SSC develops in our 3D models, appearing sooner (i.e., closer to the inflow boundary) near the edge than far away from it, with convection cells typically aligned parallel to plate motion as “Richter rolls” (Richter, 1973; Richter and Parsons, 1975; Marquart, 2001; Huang et al., 2003). The development of Richter rolls is stable even for our low v_plate=2 cm yr⁻¹ due to our high Rayleigh number (Korenaga and Jordan, 2003). Note that no EDC melting is found in the case without a plume (Fig. A3) nor in 2D cases with the same parameters as in this work (Manjón-Cabeza Córdoba and Ballmer, 2021).

Figure 3 shows the results of the reference case, which includes a plume with ΔT_plume=150 ^∘C, B_plume=100 kg s⁻¹, and D_plume=0 km. Compared to other geodynamic studies of mantle plumes (Ribe and Christensen, 1994; Ballmer et al., 2011), the most evident characteristic of this model is the lateral deflection of the plume conduit. Instead of ascending vertically, the plume conduit is displaced towards the oceanic side with thinner lithosphere. This displacement suggests some interaction of plume flow with EDC-related flow.

Figure 3 Steady-state temperature field and melting for the reference case with ΔT_plume=150 ^∘C, B_plume=100 kg s⁻¹, D_plume=0 km, and ${\mathit{\eta }}_{\mathrm{0}}=\mathrm{8.29}\times {\mathrm{10}}^{\mathrm{18}}$ Pa s. (a) Front perspective of the reference case. Cross sections of potential temperature (at the margins of the model) are shown in red-to-blue colors. The light blue and white surfaces mark isotherms (as labeled), showing the base of the lithosphere and the plume, respectively. The purple isosurface outlines the region of active melting, while the orange isosurface outlines the region of finite melt presence, including where active melt re-freezing occurs. Dashed and solid lines depict the dynamic and thermal PET (see text), respectively. (b) Side view of the model with stream lines. (c) Top view with cross section of the temperature field at z=106 km. Black arrows depict plate motion. Note that the color scale is different than in (a) and (b).

The plume ponds at the base of the lithosphere as a pancake of hot material. Before melting depths, the temperature excess of the plume with respect to the initial potential temperature is nearly the same (within 2 ^∘ depending on the depth) as the thermal anomaly at the bottom, due to our approximation of the adiabat (note that at some depths, shear heating and internal heating may increase the temperature excess of the plume above ΔT_plume, but this was found to be the exception compared to a slightly lower value). The hottest central part of the plume pancake is located at the minimum thickness of the oceanic lithosphere (i.e., at the aforementioned dent or bump). Without further analysis, however, it remains unclear whether the plume is conveyed to this minimum thickness created by EDC or if the plume actively creates a dent and EDC reorganizes accordingly. The plume pancake and melting zones are slightly asymmetric, but, again, it remains unclear whether this asymmetry is due to the spreading of the pancake at the base of a lithosphere with variable thickness at the ocean–continent transition or whether it is caused by EDC-related flow.

The plume acts to efficiently erode the imposed edge at the base of the lithosphere, displacing the thermal boundary layer. This erosion also creates a plume erosion track (PET) that is observed in all our models (Ribe and Christensen, 1994). The PET can be defined thermally, as the region were heat flow is increased due to plume erosion (discontinuous line in Fig. 3), or dynamically, as the erosion of the lithosphere that would contribute to dynamic topography (continuous line in Fig. 3). In the reference case (Fig. 3), the PET is mostly parallel to the direction of plate motion (independently of how it is defined). The plume also displaces the main EDC downwelling continent-wards: this effect starts around the plume pancake, but it continues downstream. In the reference case, it is difficult to clearly distinguish between the downwelling associated with the PET (i.e., plume curtain) and the main EDC downwelling (compare Figs. A3 and 3).

To better quantify the lateral displacement of the plume stem and the pancake, we calculate a plume deflection index (PDI) defined simply as the inverse of the slope ($\frac{\mathrm{\Delta }x}{\mathrm{\Delta }z}$) between two temperature maxima at two different depths. PDI_stem is a proxy for the plume-stem displacement, calculated as the lateral distance between the plume stem (and related thermal maxima) at z=220 km and at z=660 km (divided by the difference between both these depths, i.e., 440 km); in addition, we also report a PDI_pancake, which is calculated from the thermal maxima at depths of z=110 km and z=220 km and otherwise analogously to PDI_stem. In this work, we arbitrarily define positive values of PDI as distortions of upwelling flow “away from the edge” and negative values as “towards the edge”. In the reference case, PDI_stem=0.143 and PDI_pancake=0.109. Both the stem and the pancake are deflected towards the oceanic side. These values correspond to absolute displacements of the plume towards the oceanic domain of 63 km from 660 to 220 km depth and another 12 km from 220 to 110 km depth. In particular, the lateral displacement of the plume stem is significant. We will discuss the relevance of these values in comparison to other cases below.

We also investigate the compositional origin of mantle melts as a proxy for their geochemical signature. To do this, we evaluate the total melt volume flux M (i.e., melt produced in the mantle) and total volcanic volume flux V (i.e., melt extracted from the mantle), along with the melt flux and volcanic flux that is related to pyroxenite melting only: M_PX and V_PX, respectively. These metrics provide a compositional index for mantle melting, $\frac{M_{\text{PX}}}{M}$, and melt extraction, $\frac{V_{\text{PX}}}{V}$. The latter is the compositional origin of volcanism explicitly predicted by our models. Note, however, that this specific prediction of lithological origin of volcanism depends on the critical porosity explicitly assumed here (1 %) and on the style of melt extraction. For example, if pyroxenite-derived and peridotite-derived melts were already pooled in the mantle (instead of in a shallow magma chamber) and then were extracted together or if all melts were efficiently extracted (i.e., for fully fractional melting), $\frac{M_{\text{PX}}}{M}$ would provide a more appropriate geochemical proxy than $\frac{V_{\text{PX}}}{V}$. In other words, both $\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$ provide reasonable bounds for the compositional origin of predicted lavas.

In the reference case, $\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$ are 0.774 and 0.994, respectively. Such a dominance of pyroxenite-derived melting and volcanism is mostly explained by the relatively low plume excess temperatures and large relative seafloor ages modeled here (and relevant for eastern Atlantic volcanism; Müller et al., 2008). The related large lithospheric thicknesses restrict extensive peridotite melting, even though peridotite is the most abundant component in the plume source. Also note that pyroxenite melting starts at greater depths than peridotite melting and efficiently extracts latent heat, such that the ascent of peridotitic material is sub-adiabatic (less melting) and the ascent of pyroxenite material is super-adiabatic (more melting; Hirschmann and Stolper, 1996).

3.1 Effects of plume temperature

We conduct a series of cases with variable plume excess temperature ΔT_plume and constant buoyancy flux B_plume. We find that ΔT_plume only has minor effects on the overall flow patterns at a given B_plume. As B_plume is kept constant, the radius of a hotter plume is implicitly smaller than that of a cooler plume. As a consequence of this implicit effect of ΔT_plume on plume radii, the plume pancake and the related PET tend to be wider for smaller ΔT_plume. While the base of the lithosphere is eroded more efficiently for large ΔT_plume because a hotter plume sustains a lower-viscosity pancake, the differences are extremely small as far as B_plume is kept constant (Appendix, Fig. A4).

There is no indication that changing ΔT_plume while keeping B_plume constant systematically changes the effect of EDC-related flow on plume ascent (or the effect of plumes on EDC). The lateral displacement of the plume by EDC is similar across all our models with different ΔT_plume, as evidenced by the nearly flat trends of PDI_stem and PDI_pancake (Fig. 4a). The only noticeable difference between the models is that the plume pancake is more asymmetric for the case with ΔT_plume=100 ^∘C than for greater ΔT_plume. Note also that all PDIs in Fig. 4a are positive, implying that the plume is consistently deflected away from the edge at all depths.

Figure 4 Diagrams showing the sensitivity of several output parameters as a function of ΔT_plume. (a) Plume distortion index (PDI) for the different models (see text for explanation). Note that all values are positive (deflection away from the edge). (b) Melt volumes fluxes (M and V) and melt compositional index ($\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$) for the different models. Predictions in terms of melt production (M, $\frac{M_{\text{PX}}}{M}$) are given as dashed lines; predictions in terms of melt extraction (V, $\frac{V_{\text{PX}}}{V}$) as dotted lines.

Melt fluxes (i.e., volume fluxes of melts produced in the mantle) and volcanic fluxes (i.e., volume fluxes of melts extracted from the mantle) systematically increase with ΔT_plume (Fig. 4b). This result is intuitive, and consistent with previous work (Ribe and Christensen, 1994; Ballmer et al., 2011). In terms of the compositional origin of magmas, $\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$ decreases with the amount of melt produced and therefore decreases with increasing ΔT_plume (Fig. 4b). This result is expected, as PX (pyroxenite)-derived melts are diluted by peridotite-derived melts for increasing degrees of melting.

3.2 Effects of plume buoyancy flux

We also explore the influence of B_plume on model results. Figure 5 shows steady-state model predictions for cases with different B_plume but otherwise the same parameters as in the reference case. Increasing B_plume implicitly increases the radius of the plume. Thereby, the width and volume of the melting zone and of the plume pancake also increase, as does the area of the PET. The PET remains mostly parallel to the plate velocity vector, as for the reference case.

Figure 5 Steady-state snapshots of representative cases with different B_plume but otherwise the same parameters as in the reference case (Fig. 3). (a) B_plume=50 kg s⁻¹; (b) B_plume=200 kg s⁻¹. For reference to colors of surfaces and cross sections, see Fig. 3 caption and legend.

The lateral deflection of the plume stem is less evident for cases with higher B_plume than for cases with lower B_plume. Indeed, the high buoyancy flux plume is straighter while rising through the model box than the plume in the reference case. In fact, PDI_stem tends to 0 as B_plume increases (Fig. 6a, c), providing evidence of a limitation of the ability of EDC (or of SSC in general) to affect the rise of plumes: efficient displacement is restricted to plumes with moderate-to-low buoyancy fluxes. Nevertheless, the melting zone and the plume pancake display subtle asymmetry also in the case with the highest B_plume modeled here. As for PDI_stem, PDI_pancake also tends to decrease for increasing B_plume but remains positive.

Figure 6 Diagrams showing the variation in key output parameters to changing B_plume. (a) PDI_stem and PDI_pancake for models with variable B_plume and ΔT_plume=150 ^∘C. (b) M, V, $\frac{M_{\text{PX}}}{M}$, and $\frac{V_{\text{PX}}}{V}$ for models with variable B_plume and ΔT_plume=150 ^∘C. (c) Same as (a) for models with ΔT_plume = 200 ^∘C. (d) Same as (b) for models with ΔT_plume=200 ^∘C.

As far as the position of the main downwelling of EDC is concerned, increasing B_plume increases the displacement of the main EDC downwelling towards the continent around the plume pancake. Once the plume pancake erodes the continent–ocean transitional lithosphere, the main EDC downwelling is not shifted towards the continent side further downstream, even if PET is wider for cases with higher B_plume. In other words, while the plume controls the position of the main EDC downwelling close to the hotspot, the PET does not have an active effect on the position of the main EDC downwelling.

Due to the aforementioned radius increase as a function of B_plume, M and V both systematically increase with increasing B_plume. Regarding melt compositions, $\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$ display a shape that mirrors melt volumes (Fig. 6b, d), decreasing with increasing B_plume. Similar to the effects of plume excess temperature (see Fig. 4), the trends of melt volumes and compositions as a function of B_plume mirror each other because $\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$ decrease with increasing degrees of melting of the dominant lithology, peridotite. The influence of B_plume on magma compositions decreases at higher buoyancy fluxes, probably because the extent of vertical sublithospheric erosion becomes nearly independent of B_plume at some point. Note that the convex upward shape of the dotted gray line in Fig. 6b is due to the saturation of PX contributions at ∼100 %. In Fig. 6, the difference in composition between the produced melts and the extracted melts is greater for lower ΔT_plume, which is not evident in Fig. 4. This is explained by the much higher productivity of pyroxenite melting (and hence smaller sensitivity to an extraction threshold) than for peridotite melting at high B_plume. Among all the parameters explored in this work, B_plume shows the strongest effect on plume vigor and related melting.

3.3 Effects of distance of the plume from the edge

Next, we analyze the effects of the distance of the base of the plume from the edge, D_plume, on model results. The effects of this parameter are a good indicator of plume–EDC interaction because D_plume changes the spatial relationship between the plume and the edge while leaving intrinsic plume parameters unchanged. Figure 7 shows 3D snapshots of mantle temperature and melting as a function of D_plume for two sets of ΔT_plume and B_plume (in panels a and b, for a relatively weak plume with parameters such as in the reference case: ΔT_plume=150 K and B_plume=100 kg s⁻¹; in panels c and d, for a moderately strong plume with ΔT_plume=200 K and B_plume=200 kg s⁻¹). Interaction of the plume with the EDC convection cell and topography at the base of the lithosphere causes systematic changes in the flow patterns and related melting characteristics. For D_plume=200 km (Fig. 7a, c) the plume stem is deflected in a similar way as in the reference case (Fig. 3), for which D_plume=0 km. For D_plume=400 km, the plume stem is instead generally less affected by the presence of the edge and related EDC.

Regarding the plume pancake, we find two significant changes in the behavior predicted by our cases with variable D_plume with respect to the reference case (Fig. 3). First, as the plume is shifted away from the edge (i.e., for increasing D_plume), the plume pancake is deflected towards the edge. This transition happens at a different D_plume depending on plume properties (∼25 km in Fig. 8a, ∼125 km in Fig. 8c), but it happens nonetheless. After this rather sudden transition, the edge-ward deflection of the pancake decreases progressively with increasing D_plume. At a distance of D_plume=400 km, another notable phenomenon occurs: vigorous SSC appears in the plume pancake with dominant transverse rolls (i.e., perpendicular to the edge, Fig. 7b, d). This peculiar geometry of SSC separates the plume-fed melting zone into two distinct melting zones (Fig. 7b, d). This separation is transient, however: as the SSC downwellings move with the plate, the two melting zones are separated and merged periodically (see Supplement). Dominant transverse rolls are a specific prediction for cases with D_plume=400 km. They occur neither for cases with any other D_plume nor for test cases without a plume (see Appendix, Fig. A3). This prediction highlights the subtle effects that EDC may have on plume–lithosphere interaction as a function of D_plume and EDC flux. In any case, SSC transverse rolls have also been found in studies of plumes without a nearby edge (Ballmer et al., 2011) or studies of EDC without a plume present (Kaislaniemi and Van Hunen, 2014).

Figure 7 Steady-state snapshots of representative cases with variable D_plume. (a, c) Cases with D_plume=200 km. (b, d) Cases with D_plume=400 km. Panels (a) and (b): models with a relatively weak plume with ΔT_plume=150 ^∘C and B_plume=100 kg s⁻¹ are shown. Panels (c) and (d): models with a relatively strong plume with ΔT_plume=200 ^∘C and B_plume=200 kg s⁻¹ are shown. For reference to colors of surfaces and cross sections, see Fig. 3 caption and legend. The 3D perspective is the same as in Fig. 3.

Figure 8 shows the effect of D_plume on quantitative characteristics of plume ascent. Note that the significant changes shown in Fig. 8, both in terms of PDI and melt fluxes as a function of D_plume, are exclusively due to plume–EDC interaction (intrinsic plume parameters remain unchanged in each row of Fig. 8). PDI_stem is generally positive but highly variable. It peaks at D_plume=50 km and D_plume=150 km for the relatively weak and strong plumes shown in the top and bottom rows of Fig. 8, respectively. For higher D_plume, PDI_stem systematically decreases with D_plume. In turn, PDI_pancake becomes strongly negative for the D_plume at which PDI_stem peaks and progressively less negative for any higher D_plume. These results emphasize the strong effects of plume–EDC interaction and its diversity as a function of D_plume (and for plumes with different ΔT_plume and/or B_plume). The switch to dominantly transverse rolls in the pancake for plumes far from the edge (i.e., at D_plume=400 km) does not seem to strongly affect the deflection of the plume stem or shallow pancake.

As the plume pancake, the PET changes with changing D_plume, becoming more asymmetric (with respect to the plate velocity vector) and being deflected continent-wards whenever PDI_pancake is negative. As a result, the continent-side limit of the PET (right side of the panels in Figs. 3 and 7) always remains close to the main EDC downwelling. The PET is also affected by SSC in the plume pancake perpendicular to plate motion. Very likely, these predictions have implications for dynamic topography and swell geometry. While the PET and the plume pancake are greatly affected by D_plume, this is not the case for the geometry of the main downwelling of EDC. This downwelling is displaced by the plume in a similar way in all cases with variable D_plume.

Figure 8 Diagrams showing the sensitivity of key output parameters to D_plume. (a, c) PDI_stem (light brown) and PDI_pancake (blue) in the steady state for models with variable D_plume. (b, d) M, V, $\frac{M_{\text{PX}}}{M}$, and $\frac{V_{\text{PX}}}{V}$, for models with variable D_plume. Panels (a) and (b): results for a relatively weak plume with ΔT_plume=150 ^∘C and B_plume=100 kg s⁻¹ are given. Panels (c) and (d): results for with a relatively strong plume with ΔT_plume=200 ^∘C and B_plume=200 kg s⁻¹ are given. Legend as in Figs. 4, 6.

Compared to those on PDI, the effects of D_plume on melt fluxes and compositions are less severe. Figure 8b and d show the trends of melting-related parameters as a function of D_plume at the same scale as other figures (e.g., Fig. 6). The effects on melt fluxes appear small, which is mostly due to the logarithmic scale of the figure; nonetheless several “regimes” or different behaviors can be distinguished on the basis of distance of the plume to the edge. Similar to the PDI figure, there is an initial regime (D_plume=0 km for the cases with B_plume=100 kg s⁻¹, and D_plume=0–100 km for the cases with B_plume=200 kg s⁻¹) with lower PDI, and the melting volumes remain mainly flat (at least in Fig. 8d). Then, at greater distances, EDC interacts strongly with the plume, resulting in slightly lower melting volumes with a smooth peak around D_plume=200 km. Finally, at D_plume=400 km, melting volumes increase substantially due to SSC, but volcanism remains practically the same (suggesting that the main volume of melting still happens at the top of the plume conduit).

In general, plume deflection, as caused by the effects of EDC, tends to systematically decrease the amounts of hotspot magmatism for a given plume vigor and/or temperature. The least negative (or most constructive) plume–EDC interaction occurs near D_plume=0 km and D_plume=200 km. These locations roughly reflect the intrinsic pattern of the EDC-related and neighboring “triggered” SSC-related upwellings, as predicted by the companion paper, but not exactly so. The differences are likely due to the effect of the presence of a plume on EDC and SSC patterns. That the distance of the change in regime (from weak to strong) of influence of EDC on the plume depends on plume vigor is also related to the effects of the plume (and plume pancake) on the wavelength of EDC. Thus, EDC appears to affect plume ascent and vice versa.

3.4 Effects of mantle viscosity

Finally, we explore models with different reference viscosities. Figure 9 shows data for cases with variable viscosity (D_plume=200 and 300 km) and otherwise the same parameters as in the reference case. Similar to the effects of plume temperature, the width of the plume stem is implicitly smaller with decreasing reference viscosity. At D_plume=200 km, one of the most striking characteristics of these cases is that the deflection of the plume stem is less severe for the high-viscosity and the low-viscosity case than for the intermediate-viscosity case, shown in Fig. 7c (same D_plume and same B_plume but intermediate η₀). Indeed, PDI indexes (Fig. 9a) display a maximum in terms of plume deflection for the intermediate viscosity value of the reference case ${\mathit{\eta }}_{\mathrm{0}}=\mathrm{8.29}\times {\mathrm{10}}^{\mathrm{18}}$ Pa s (Fig. 3).

However, at D_plume=300 km the trends in Fig. 9 depict a more systematic behavior with cases with lower viscosity featuring higher deflection. The differences between the cases at D_plume=200 and 300 km may be analogous to the differences between cases with different B_plume, in which the critical D_plume, where PDI_pancake switches from positive to negative, is increased with increasing buoyancy flux. However, note that in this case, the change is not from negative to positive PDI_pancake. Instead, we find that the greater spread of the plume pancake in cases with lower reference viscosity (Fig. 10) leads to the plume pancake being affected by the topography of the edge, decreasing PDI_pancake at D_plume to smaller than 300 km.

Since the vigor of EDC decreases with increasing η₀ (Sleep, 2007; Till et al., 2010; Davies and Rawlinson, 2014; Manjón-Cabeza Córdoba and Ballmer, 2021), EDC–plume interaction also becomes less important. This is evident at higher viscosities, where both PDI_stem and PDI_pancake decrease (Fig. 9).

Figure 9 Diagrams showing the sensitivity of selected output parameters to η₀ in the steady state. Note that the PDI change is different for different D_plume.

Figure 10 Steady-state model snapshot for the low-viscosity cases (${\mathit{\eta }}_{\mathrm{0}}=\mathrm{5.53}\times {\mathrm{10}}^{\mathrm{18}}$ Pa s). (a) Case with D_plume=200 km and B_plume=200 kg s⁻¹. (b) Case with D_plume=300 km and B_plume = 200 kg s⁻¹. Note that the deflection of the pancake at 300 km is greater than at 200 km (Fig. 9) as opposed to the case for intermediate and high viscosity.

The position of the main EDC downwelling is weakly affected by changing viscosity. However, the relationship between the main EDC downwelling and the PET changes substantially. In cases with decreased viscosity (Fig. 10), the PET expands, and the continent-side limit of the dynamic PET crosses the main downwelling of EDC, which is not observed in any of our other models (i.e., changing any other property). This implies that in models with low viscosity the main EDC downwelling is situated inside the PET. We also find that the symmetry of the PET with respect to the vector of plate velocity is higher for the cases with low viscosity than for the cases with intermediate and high viscosity. In other words, the PET is more asymmetric (continent-ward) even if the plume pancake features small PDIs. These model predictions are explained by the very low viscosity of the pancake (in cases with low η₀), which promotes spreading independent of nearby features, such as the edge.

In turn, model predictions in terms of melting as a function of η₀ are as expected. Both M and V increase with decreasing η₀ and, along with this decrease, $\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$ decrease (Fig. 9b). In addition to plume-related hotspot melting, melting away from the hotspot (i.e., directly due to EDC or SSC) appears for some cases with low η₀. This melting is minor and consistent with the low-η₀ cases of the companion paper (Manjón-Cabeza Córdoba and Ballmer, 2021).

4 Discussion

We ran a wide range of 3D numerical models to systematically explore the interaction between EDC and mantle plumes. One of our main results is that the plume geometry, its interaction with the lithosphere, and the extent of related melting depends on the distance of the plume from the edge, being altered by EDC. Despite these important effects, the buoyancy flux of the plume still remains the main influence on the characteristics of plume–lithosphere interaction and hotspot magmatism.

We quantify the deflection of plumes by two characteristic parameters: the deflection of the plume stem and the deflection of the shallow plume conduit and the plume pancake. The plume stem is systematically deflected away from the edge. This may provide an explanation as to why hotspot tracks in the Atlantic preferentially occur near and sub-parallel to the continental margin but rarely across it (an exception to this is the Cameroon Volcanic Line;  Fitton, 1980; Déruelle et al., 2007). On the other hand, plume pancake deflection commonly (but not always) occurs towards the edge. This prediction may explain why some hotspot tracks (such as the Canaries) do not strictly align with plate velocity, and volcanism is widespread with more activity far from the continental margin than near it (e.g., La Palma vs. Gran Canaria).

The nature of the deflection of the plume pancake (and also of the plume stem) systematically changes in our models with the distance of the plume to the edge. However, this deflection is generally predicted to decrease with increasing B_plume (Fig. 6a, c) relative to the EDC flux. Test cases with a greater step of lithospheric thickness at the continental margin confirm this prediction (Appendix, Fig. A5). Such a configuration implies greater downwelling fluxes of EDC (see Manjón-Cabeza Córdoba and Ballmer, 2021), and leads to greater deflections of the plume and more asymmetric plume–lithosphere interaction than the reference case. In terms of PDI, the absolute value of PDI_stem is 0.12, as opposed to 0.10 in the case with regular edge and similar parameters (Fig. 7c); PDI_pancake is −0.07, as opposed to −0.05. A higher value for both indices suggests that the plume–EDC interaction depends on the ratio between the downward flux of materials due to EDC and the upward flux of materials due to plume activity.

The dependence of plume–EDC interaction to these two fluxes is consistent with previous work such as Ballmer et al. (2011), who found that, even for high plume fluxes, sublithospheric convective instabilities can have an effect on the surface expression of mantle plumes. In turn, our results challenge the opinion that strong external fluxes will overprint or even ignore EDC (King and Anderson, 1998; Till et al., 2010; Kaislaniemi and Van Hunen, 2014), a statement that may only hold true for high buoyancy flux plumes.

We find that weak-to-intermediate plumes can be strongly deflected by EDC with significant effects on plume-related magmatism (Figs. 7, 8). Plumes can be laterally deflected by >100 km (PDI > 0.2) as a function of plume–EDC geometrical configuration. While the specific D_plumes for which maximum plume deflection occurs depends on plume flux and asthenospheric viscosity, these extents of deflection are of the same order than the apparent displacements of the location of the main EDC downwelling.

There is a clear difference in the position of the main EDC downwelling between cases with and without plumes (see Fig. A1 in this paper and  Manjón-Cabeza Córdoba and Ballmer, 2021). The downwelling is systematically pushed continent-wards near the hotspot and also further downstream. Nonetheless, this displacement is only influenced little by changing plume properties or the rheology of the mantle. Of the parameters studied, only B_plume efficiently affects the position of the main EDC downwelling (see also  Burov et al., 2007). More surprisingly, the main EDC cell is not strongly bound to the plume flow downstream of the hotspot (and the PET), as evidenced by our low-η₀ models. On the contrary, we find a high coupling between the plume pancake and stem and the EDC cell near the hotspot. Overall, this evidence points towards a greater effect of EDC on plumes than vice versa. However, this affirmation is likely dependent on parameter choices (such as plume buoyancy flux), as well as initial conditions (such as edge geometry).

The presence of the plume can also affect the vigor of EDC. In particular, the vertical component of the main EDC downwelling is faster in the cases with negative PDI_pancake than in cases with positive PDI_pancake. For example, the reference case shows a vertical velocity component of up to ∼1.3 cm yr⁻¹ near the plume pancake, while an equivalent case of D_plume=200 km features a vertical velocity component of ∼1.5 cm yr⁻¹. This again highlights the effects of constructive vs. destructive plume–EDC interaction as a function of geometry (or plume–edge distance).

Our models stand in contrast to other studies that focused on very high buoyancy flux plumes and plume arrival (plume heads). For example, Burov et al. (2007) and François et al. (2018) observed a deflection of the plume head ocean-wards, not unlike our models with D_plume=0 km. In such a scenario, the plume completely disrupts the EDC cell (favoring craton removal), and plume flow is conditioned by the topography of the lithosphere. For very strong plumes and plume heads, it is likely that plume flow dominates over EDC. Concerning the influence of the plume head, our models are not suited to evaluate its influence because they are focused on a steady-state plume setting (Ribe and Christensen, 1999). Moreover, various studies have shown that different lithospheric strength models will interact with different kinds of plume heads in radically different ways, everything else being equal (Gerya et al., 2015; Koptev et al., 2021). It is unlikely, however, that the plume head of the Canaries still has an influence on the archipelago since the hotspot track is >60 Myr old.

In contrast to the EDC-only cases in the first part of our work (Manjón-Cabeza Córdoba and Ballmer, 2021), volcanic (or melt) volume fluxes are significantly higher, displaying strong variations as a function of plume parameters and moderate variations as a function of D_plume. Therefore, a subset of our models can account for the volumes of eastern Atlantic hotspots. However, both B_plume and ΔT_plume also affect the geochemistry of the melts. In our models, the composition of melting due to plume + EDC is generally less enriched than for melting due to EDC-only (Manjón-Cabeza Córdoba and Ballmer, 2021). Any increase in volcanism due to higher B_plume or ΔT_plume is associated with a decrease in enrichment (i.e., the fraction of melting products from enriched lithologies such as pyroxenite; Figs. 4b, 6b, d), with only minor effects on geochemical proxies as a function of D_plume (Fig. 8).

In the companion paper (Manjón-Cabeza Córdoba and Ballmer, 2021), we clearly showed that EDC alone is insufficient to generate the Canary Islands magmatism and that the contribution from a mantle plume (or equivalent source) is required. The Canary Islands feature two islands in shield-building stages: for El Hierro, Carracedo et al. (1998) estimates a minimum volcanic flux of $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{1}}$ km³ kyr⁻¹; for La Palma, the inferred volcanic rates are around 1 km³ kyr⁻¹ (Day et al., 1999). For the whole archipelago, independent estimates place extrusion volumes at 1–10 km³ kyr⁻¹, depending on whether the 18th century Timanfaya eruption(s) are considered an anomalous event (Longpré and Felpeto, 2021). These numbers can easily be doubled when considering underplating and plutonism beneath the islands (Klügel et al., 2005). Overall, these values are higher than any published estimate of EDC-generated volume, including work that considered additional geometric complexity for the ocean–continent transition (Duvernay et al., 2021; Negredo et al., 2022).

Our models predict that relatively weak plumes with parameters similar to that of the reference case are sufficient to generate these amounts of magmatism. However, it remains difficult to pinpoint plume parameters, e.g., as plume temperature and buoyancy flux are traded off against each other (Figs. 4b, 6b, d). This result implies that the Canary plume must be of either medium/low flux (i.e., B_plume<200 kg s⁻¹) or medium/low temperature (i.e., ΔT_plume<200 ^∘C) but not both, as the melt volumes would be too low then (i.e. < 0.1 km³ yr⁻¹). Note that M and V in the relevant model cases are of the same order of magnitude, suggesting that our results are robust despite model limitations and simplifications in terms of modeling mantle melting and extraction. Our predictions for $\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$ further constrain the properties of the plume: the occurrence of shield-like magmas in the Canary Islands (Abdel-Monem et al., 1971, 1972; Carracedo et al., 1998) strengthens the suggestion that B_plume and ΔT_plume cannot (both) be lower than for our reference case, as this would generate melts that are too enriched, i.e., with PX contributions much higher that inferred from petrological work (Day et al., 2009; Day and Hilton, 2011). More likely, one of these parameters (or both) must be slightly higher than in the reference case.

The distance of the Canary hotspot from the African passive margin is ∼250–300 km. Considering the PDI values predicted for relatively weak plumes (Fig. 8b), we estimate that the Canary plume at 660 km depth is centered ∼50–100 km closer to the African margin than the hotspot (which is located near El Hierro). Likewise, we estimate that the plume at 410 km depth is centered ∼30–70 km closer to the margin. This prediction is consistent with the receiver functions study of Saki et al. (2015): the location of the shallowest 410 km discontinuity is shifted from the hotspot at El Hierro eastward towards Lanzarote. In addition, the recent tomographic study by Civiero et al. (2021) shows an arcuate upper-mantle plume beneath the Canaries not unlike the plumes from our study. Attending to this evidence, we estimate that the cases with $\mathrm{150}\le D_{\text{plume}}\le \mathrm{250}$ km best match the configuration of the Canary plume. This finding implies that plume–EDC interaction (e.g., as quantified by PDI values) is significant for the Canary hotspot. The plume stem is pushed to the west (away from the edge) by about 80–110 km and the plume pancake is pulled back to the east (towards the edge) by about 25–35 km. If the volcanic flux at the Canaries is significantly higher than estimated by Carracedo et al. (1998), e.g., due to unaccounted magmatic crustal underplating, we reach the same conclusions, predicting very similar PDI values (Fig. 8c).

In addition, we find that several key characteristics of the Canary Islands are matched by our models. The Canaries present active volcanism far from the inferred deflection point of the plume stem near El Hierro (in fact, all islands are currently active with the exception of La Gomera; Abdel-Monem et al., 1971, 1972; Carracedo, 1999; Geldmacher et al., 2005). Several of our cases predict deflection of the plume pancake and the melting zones toward the continental margin (Fig. 5), including the cases with $\mathrm{150}\le D_{\text{plume}}\le \mathrm{250}$, which would explain the shape of the whole archipelago and the geographic distribution of volcanism. Even with a plate velocity that produces a volcanic track parallel to the ocean–continent transition, given the right distances to the edge, the plume pancake may not necessarily be parallel to the plate movement.

While our models are able to explain coeval volcanism across the islands, they do not reproduce the general westward progression of the main shield stage over time (Geldmacher et al., 2005). Indeed, the Canary hotspot may have moved westward with respect to the African margin in the last few million years (Wang et al., 2018). Accordingly, the distance of the hotspot relative to the African plate may have changed, rendering plume–EDC and plume–lithosphere interaction a transient phenomenon, which cannot be explicitly addressed by our steady-state model setting. However, Figs. 7 and 8 provide an indication of how the geometry of the plume and plume–lithosphere interaction may have changed during such a movement. As the plume moved away from the margin, the effects of changes in plume–EDC interaction may have extended and deflected the pancake and, therefore, extended the area of volcanism from a single track to a wide zone, consistent with Geldmacher et al. (2005).

Recently, Negredo et al. (2022) provided a mechanism by which plume migration due to EDC can also occur. While the deflection of the plume pancake continent-ward is not predicted by their models, it is possible that a combination of (their and our) phenomena can explain the history of the Canaries. Our model predicts that, eventually, further movement of the plume away from the edge may decrease the extent of plume deflection. Alternatively, if the vigor of the plume (or of EDC) has recently changed, plume–EDC interaction and plume displacement would also have changed and the recent movement of the hotspot relative to the African margin would be potentially unrelated to any movement of the deep plume stem or even the plate movement.

The application of our models to other hotspots in the eastern Atlantic is less obvious. Lodhia et al. (2018) and King and Ritsema (2000) have suggested a link between the Cabo Verde plume and the downwelling at the African margin near Cabo Verde. From our models, however, a significant effect of EDC on plume ascent over such long distances (over 1000 km) is not justified. It is true that higher mantle viscosities may result in larger EDC cells, but they will also result in lower EDC-related fluxes as already shown in Fig. 5a of Manjón-Cabeza Córdoba and Ballmer (2021). Regardless, the models presented here include a plate velocity that is not fully consistent with the Cabo Verde “near-zero” plate velocity. In fact, Patriat and Labails (2006) detected a “bulge” or “bump” along the continental–oceanic transition between the Canary Islands and Cabo Verde. Whether this bulge is related to an EDC upwelling is difficult to determine, but a topography high is expected in the area of maximum sublithospheric erosion above an EDC upwelling (Manjón-Cabeza Córdoba and Ballmer, 2021). Such a relationship is also consistent with lithospheric models that detect EDC-related erosion at the bottom of the lithosphere beneath the Canaries (see for instance model c2 in Fig. 7 in Fullea et al., 2015), which would imply that the main EDC upwelling happens very close to the edge. Such a proximity is also suggested by our 3D models here. The aforementioned bulge is consistent with the location of the eastern islands of the Canary archipelago but does not seem to be related directly to the Cabo Verde hotspot.

The Cameroon volcanic line is very close to the continental margin and even crosses it, such that an influence of EDC is expected. From our models, it can be inferred that at least part of the volcanic “track” with widespread volcanism perpendicular to the ocean–continent transition is consistent with plume–EDC interaction. However, the geometry of the edge near the Cameroon volcanic line is considerably more complicated than in our models, and Duvernay et al. (2021) have shown that the patterns of volcanism due to EDC can change considerably with complex 3D configurations. To better understand this volcanic field, more specific work taking into consideration the shape of the ocean–continent transition in this region is required.

Regarding future work, our models can be expanded to address limitations. While we use a decreased activation energy to mimic dislocation creep, we expect the explicit effects of a composite rheology including non-Newtonian viscosity laws will quantitatively modify our results. In particular, shear-thinning rheologies are expected to localize EDC flow and, potentially, plume ascent. In addition, the continental lithosphere in our models is not rheologically stabilized. Differences between models with and without rheological stabilization (e.g.,  Kaislaniemi and Van Hunen, 2014; Manjón-Cabeza Córdoba and Ballmer, 2021) suggest that stabilization may decrease EDC downwelling flux and protect the continental lithosphere against erosion. As a consequence, the migration of the EDC main downwelling may decrease for a given lithospheric age; PDI absolute values also decrease for a given τ_c.

Our current models also rely on a simplified compositional approximation. In our work, we considered purely thermal plumes with a composition equal to the “background” mantle. In turn, this is inconsistent with CO₂-rich volcanism in the Canary Islands (Allegre et al., 1971; Taracsák et al., 2019). However, explicitly modeling CO₂-related melting in geodynamic models is not possible to date. Moreover, the fraction of pyroxenite in the source is unlikely to be constant along the plume path (Day et al., 2010). Future models are needed to explore the transient changes along the Canary hotspot track.

Moreover, as said above, in this work we did not explore in-depth geometrical considerations of the craton and the continental–oceanic transition. In previous works, several authors noted the importance of these geometrical aspects on conditioning the mantle flow around the edge (Till et al., 2010; Kim and So, 2020; Manjón-Cabeza Córdoba and Ballmer, 2021; Duvernay et al., 2021). In this article, we focused on plume characteristics instead, but we acknowledge that the quantitative inferences made could change when changing craton characteristics, although preliminary tests confirm that results will be qualitatively consistent (Appendix, Fig. A5).

One of the main characteristics of the Canary Islands hotspot is the near absence of a hotspot swell. While the PET provides a proxy for the potential location of the plume swell, analysis of dynamic topography using geodynamic models may answer the question of whether deflection of the plume and the pancake by EDC can blur the dynamic-topography signal of the plume (Huppert et al., 2020). Another area of potential further work is to better constrain the geochemical fingerprint of our model magmas. In our models, ΔT_plume and B_plume have very similar effects on the compositional proxies of volcanism used here ($\frac{M_{\text{PX}}}{M}$ and $\frac{V_{\text{PX}}}{V}$). However, as soon as several geochemical systems are considered (e.g., major and trace elements, isotopes), the effects of ΔT_plume should have a distinct effect on the geochemistry of magmas from increasing B_plume. Unfortunately, additional assumptions in terms of starting composition of PX and peridotite (or their potential interactions; Ballmer et al., 2013; Jones et al., 2017) are required to explicitly predict trace-element and isotopic signatures (Bianco et al., 2008). To date, no practical melting parameterization is available to realistically predict major element compositions from geodynamic models. Future work will focus on a new melting parameterization that can help to discriminate between parameters in this setting and other geodynamic models.

5 Conclusions

We studied the effects of edge-driven convection (EDC) on low-to-intermediate buoyancy flux plumes. The following points summarize the main findings of this study:

• Low and intermediate buoyancy flux plumes interact with shallow mantle flow related to sublithospheric convective instability, which causes the plume to be deflected with important effects on the volume flux (and composition) of hotspot melting.

• The interaction of the plume with edge-driven convection highly depends on the distance of the plume to the ocean–continent transition, but the distance for which EDC has the strongest influence varies with physical properties of the mantle and plumes. For example, weaker plumes (lower buoyancy flux, lower temperature) are most affected closer to the edge (i.e., continental margin) than more vigorous plumes (higher buoyancy flux, higher temperature).

• The ratio of the buoyancy flux of the plume with respect to the flux of material from EDC is one of the most important factors to control plume–EDC interaction at a given plume–edge distance, including deflection of the plume stem away from the continental edge and of the pancake towards the edge.

• In the Canary Islands, a plume of low buoyancy flux and high temperature or, alternatively, a plume with moderate buoyancy flux and low temperature may be rising at 200 km from the continental margin, being deflected and creating the complex age progression and widespread volcanism.

Appendix A: Additional figures

Figure A1 Example of trends used to evaluate steady state in our models. Shown are the reference case (a and b; see also Fig. 3) and the case of D_plume=400 km, ΔT_plume=150 ^∘C and B_plume=100 kg s⁻¹ (c and d; see also Fig. 7b and Supplement). Panels (a) and (c) show the total melting rate (integrated across the whole model) and (b) and (d) show the root mean square of the vertical velocity (vVrms). All axes are in non-dimensional scale.

Figure A2 Viscosity profiles representative of the reference case (Fig. 3). As dashed lines, viscosity profiles for the inflow boundary conditions: in black, oceanic side; in gray, continental side (see also Manjón-Cabeza Córdoba and Ballmer, 2021). As solid lines, two representative viscosity profiles of the plume at x=1254 km (green line) and x=1287. Note that due to the plume deflection, a single vertical profile representing the whole plume is not possible (see also Fig. A4).

Figure A3 Vertical cross sections of a case without a plume perpendicular to plate velocities. (a) y=220 km; (b) y=660 km; (c) y=1320 km. Compared with Fig. 3 of Manjón-Cabeza Córdoba and Ballmer (2021), changes in time of EDC are well represented at depth, while plate velocities do not disrupt the general planforms of EDC.

Figure A4 Vertical cross sections parallel to the edge (perpendicular to the x direction) and x=1254 km. The isotherm of T_pot=1215 ^∘C is drawn as a black line. (a) Reference case. (b) Case with ΔT = 200 ^∘C and other values as in the reference case. (c) Closeup of the reference case with the isotherm of the case depicted in (d) (dashed yellow line). Note that due to plume deflection, the x coordinate for maximum lithospheric erosion does not correspond to the location of the thermal anomaly at the bottom.

Figure A5 Snapshot of test case with an increased continental lithospheric thickness (τ_c=200 Ma), D_plume=200 km, and other properties as in the reference case. See text for details.