The role of subduction in the formation of Pangaean oceanic large igneous provinces Open Access

Supplementary material: Numerical code input parameters and animations are available at https://github.com/heronphi/Heron_Super22 [last accessed 28 July 2023]

Pangaea is the most recent supercontinent that has assembled and dispersed over the course of Earth's history (Worsley et al. 1984; Nance et al. 1986; Dalziel 1991; Hoffman 1991, 1997; Moores 1991; Murphy and Nance 1991; Torsvik 2003; Li et al. 2004, 2008; Murphy et al. 2006, 2009; Nance and Murphy 2013; Mitchell et al. 2021). In particular, Pangaea was formed during the Carboniferous (c. 320 Ma; Pastor-Galán 2022) and broke up from the Jurassic–Cretaceous (c. 175 Ma) onwards. Although there is no agreement on a formal definition for what constitutes a supercontinent (Pastor-Galán et al. 2019; Heron et al. 2021), the process by which continents amalgamate is well established as part of the life cycle of an ocean (e.g. the Wilson Cycle; Wilson 1966; Burke and Dewey 1975). The cycle of supercontinent formation and dispersal, in the context of whole-mantle flow, can be viewed as involving the following steps:

1. Closure of oceans by subduction causes large-scale continental collisions (Fig. 1a). Any interior subduction that yields the assembly of the newly amalgamated landmass terminates (Fig. 1b).

2. Termination of interior subduction related to supercontinent formation changes the dynamics both within and below the supercontinent. The initiation and development of circum-supercontinent exterior subduction zones, and the related descent of oceanic lithosphere into the mantle, disrupts the regional mantle dynamics and produces a return flow that can potentially trigger deep mantle plumes (Tan et al. 2002; McNamara and Zhong 2005; Heron and Lowman 2010; Davies et al. 2012, 2015; Heron et al. 2015; Dannberg and Gassmöller 2018; Crameri et al. 2019; Heyn et al. 2020).

3. Coupled lithosphere and mantle dynamics ultimately break the landmass up after a period of 100–400 myr (Condie 2022), leading to continental dispersal. These dynamics may be related to plate driving forces, including a mantle push force from newly formed plumes beneath the supercontinent and/or slab-pull from circum-supercontinent subduction zones, as well as mantle traction and insulating temperature effects (e.g. Burke and Dewey 1973; McKenzie 1977; White and McKenzie 1989; Hill 1991; Courtillot et al. 1999; Tan et al. 2002; Gurnis et al. 2004; Lenardic et al. 2005, 2011; Coltice et al. 2007, 2009; Brandl et al. 2013).

Ancient deep mantle plumes are thought to have been manifested at the Earth's surface by voluminous magmatism erupted over relatively short geological timescales (typically 1–5 myr) to form large igneous provinces (LIPs) (Burke and Torsvik 2004), and these LIPs can be found in both oceanic and continental lithosphere. A deep mantle source for such plumes can be inferred from different geological and geophysical lines of evidence including regional domal uplift, which suggests the sub-crustal arrival of buoyant plume material (e.g. Sleep 1990; Davies 1999; Rainbird and Ernst 2001; Sengör 2001; Saunders et al. 2007). Further, a deep mantle source hypothesis is supported by tomographic images beneath modern hotspots (French and Romanowicz 2015), deep mantle geochemical signatures including recycled oceanic lithosphere (Sobolev et al. 2007; Williams et al. 2019) and a potential link to low shear-wave velocity seismic features above the core–mantle boundary (Burke et al. 2008). Alternative interpretations of the origin and source of LIPs have been suggested but remain contested (White and McKenzie 1989; Coffin and Eldholm 1994; Ernst et al. 2005; Campbell and Kerr 2007; Foulger 2007; Ernst 2014). Although other mechanisms have the ability to produce LIPs, in this study we concentrate on the role of a deep mantle plume in the formation of a LIP and explore the role of plate tectonics in this process.

Several studies have investigated the temporal, spatial and geodynamic relationships between LIP activity and plate tectonics (e.g. Yale and Carpenter 1998; Ernst et al. 2005; Ernst and Bleeker 2010; Sobolev et al. 2011; Buiter and Torsvik 2014; Hastie et al. 2014; Heron et al. 2015; Glišović and Forte 2017; Marzoli et al. 2018; Beccaluva et al. 2020; Heyn et al. 2020; Condie et al. 2021). Studies suggest that relatively little LIP activity is associated with the amalgamation stage of the supercontinent cycle (e.g. Yale and Carpenter 1998; Ernst et al. 2005; Ernst and Bleeker 2010). However, after a supercontinent has existed for a period of time (100–200 myr; Sobolev et al. 2011), the number of LIPs increases on a global scale (e.g. Yale and Carpenter 1998; Ernst et al. 2005; Ernst and Bleeker 2010; Sobolev et al. 2011). For example, c. 26 LIPs have been emplaced globally since the formation of Pangaea (Table 1).

A number of papers have discussed the role of subducting slabs of oceanic lithosphere controlling mantle flow and therefore the position and timing of plumes (Tan et al. 2002; McNamara and Zhong 2005; Heron and Lowman 2010; Davies et al. 2012, 2015; Heron et al. 2015; Dannberg and Gassmöller 2018; Heyn et al. 2020). These studies show that sinking slabs of oceanic lithosphere have the power to stir mantle flow and control the thermal evolution of the mantle, which is of particular interest during amalgamation and advanced break-up phases of the supercontinent cycle, when subduction patterns change on a global scale.

Most of the discussion on the supercontinent cycle and the formation of plumes has been directed at the formation of LIPs and/or the break-up of a (super)continent (see Heron 2019 for a review). The discussion has focused in particular on the formation of the Central Atlantic Magmatic Province (CAMP) (Marzoli et al. 2018), Karoo (Hastie et al. 2014) and Parana-Etendeka (Beccaluva et al. 2020) LIPs, which may have impacted the opening of the Atlantic, and the Deccan Traps LIP, which may have played a role in India's separation from Madagascar (Glišović and Forte 2017). Although mostly located within oceans today, at the time of their eruption all of these LIPs formed within the interior of Pangaea and so were first manifested on Earth's surface within Pangaea's circum-supercontinent subduction girdle (Fig. 2). Indeed, focused long-term subduction may have been capable of forcing a return flow in the mantle beneath Pangaea, thereby producing sub-supercontinental plumes as discussed in a number of previous studies (Zhong et al. 2007; Yoshida 2010; Li and Zhong 2009; Heron and Lowman 2010, 2011). In addition, anomalous material near the core–mantle boundary beneath Africa and the Pacific Ocean, characterized by low shear velocities in seismic tomography (see below), has been indicated as a location for plume generation (e.g. Torsvik et al. 2006, 2008, 2010; Burke et al. 2008). In particular, the margins of the African Large Low Shear Velocity Provinces (LLSVPs) have been linked to LIP formation within the interior of Pangaea (e.g. Torsvik et al. 2010).

However, a number of LIPs that formed after supercontinent amalgamation do not lie within this interior framework. Those in the Pacific Ocean (oceanic LIPs or plateaux), for example, formed outside the supercontinent subduction girdle (i.e. in areas that were exterior to the supercontinent) on plates that were subducting beneath it. LIPs such as the Ontong Java Plateau (OJP) (Fig. 1c), Caribbean LIP (CLIP) (Fig. 1d) and the Shatsky Rise all formed in provinces exterior to Pangaea (Fig. 1), yet investigations into their formation in the context of circum-supercontinent subduction are limited.

Here, we test the hypothesis that the formation of mantle plumes and resulting LIPs that are interior (i.e. on the overriding plate) and exterior (i.e. on the subducting plate) to a supercontinent could be generated by the same subduction setting. This dual-mechanism for plume formation can be conceptualized in the following way: the subducting slab pushes oceanic material in the direction of its dip beneath a continent (also referred to as Andean-style subduction) to generate an interior mantle plume through a return flow, while the same subduction back flow can also trigger one or several plumes on the oceanward side (i.e. exterior scenario) (Fig. 3). Both interior and exterior plumes may be linked to the margin of an LLSVP, and in this study we prescribe a thermochemical layer that is passive to the mantle flow (rather than fixed in place) in order to isolate the impact of subduction on plume formation.

Several mechanisms have been discussed regarding the formation of exterior (oceanic) LIPs, in particular for the OJP (Tarduno et al. 1991; Zhang et al. 2020; Isse et al. 2021) and Shatsky Rise (Sager et al. 2016). As a contribution to this discussion, we test the theory for exterior plume formation caused by circum-supercontinent subduction by using numerical experiments of 2D and 3D thermochemical mantle convection that investigate the effects of subduction processes related to the supercontinent cycle. We analyse the thermal response of the mantle driven by the evolution of surface tectonics to identify the timing and location of exterior plumes, and link this distribution to the formation of oceanic LIPs following Pangaea amalgamation.

In this section we outline the base numerical model setup, followed by the specific supercontinent process for 2D and 3D modelling. A full description of the numerical model setup can be found in the ASPECT input files provided in the Supplementary material.

The equations are solved for time t, velocity u, pressure p, temperature T and chemical composition C (e.g. harzburgite or recycled ocean crust). η is the viscosity, $ε˙$ is the deviatoric strain rate $(12(∇u+∇uT)−13(∇⋅u)δij)$⁠, g is the gravitational acceleration, ρ is the density, C_p is the specific heat capacity (at constant pressure), k is the thermal conductivity, H is the radiogenic heat production and α is the thermal expansivity.

Our models begin from an initial condition that has an average mantle temperature gradient and composition that is appropriate for the present-day mantle (Fig. 4) as determined by previous studies (e.g. Dannberg and Gassmöller 2018; Zhang and Li 2018). The initial mantle temperature is laterally homogeneous, following an adiabatic profile with thermal boundary layers at the top and bottom (Dannberg and Gassmöller 2018). We fix the surface temperature to 273 K and core–mantle boundary temperature to 3700 K (Fig. 4). Although 3700 K is a higher temperature than used in some previous studies (e.g. Zhang and Li 2018), it is well within the range of uncertainty for the exact present-day value (Nomura et al. 2014). There is greater uncertainty in historical core–mantle boundary temperature values, with previous studies indicating a higher basal temperature in the Hadean (e.g. 4000 K; O'Neill and Zhang 2018). However, the change in temperature during the geological timeframe we studied is not significant given the uncertainty in its exact value (Nomura et al. 2014). Consequently, we can prescribe a fixed a value of 3700 K for all times during the simulation (e.g. Lay et al. 2008).

As this viscosity law produces very strong lateral and vertical viscosity variations over the temperature range appropriate for Earth (ΔT = 3427 K) (Fig. 4), we limit the viscosity range to four orders of magnitude to reduce the computation time and make model runs feasible (1 × 10²⁰–5 × 10²³ Pa s). In Earth, lateral viscosity contrasts are likely to be larger than in our models, and we test the robustness of our base model in a series of 2D simulations with an altered viscosity regime. Furthermore, our model setup does not include plastic yielding.

As in previous studies (e.g. Dannberg and Gassmöller 2018), we assume an average mantle composition of 82% harzburgite and 18% recycled oceanic crust (Xu et al. 2008; Nakagawa et al. 2010), and compute the material properties ρ, α and C_p using PerpleX (Connolly 2009) and a mineral physics database (Stixrude and Lithgow-Bertelloni 2011). Accordingly, all material properties include the effects of phase transitions, with thermal expansivity and specific heat taking into account the corresponding latent heat release and consumption (Nakagawa et al. 2009).

Here Δρ is the density contrast between the anomalously dense material of the thermochemical pile and the average mantle composition directly above the layer. In the buoyancy equation, ΔT is the difference between the actual core–mantle boundary temperature (fixed at 3700 K) and the adiabatic profile (at a reference depth of 2500 km). At the core–mantle boundary, the buoyancy number of the thermochemical layer is c. 1.9, where ΔT is 1284 K and the density difference is an increase of 2% (ρ = 5551 kg m⁻³).

In 2D models, we test different parameters (Table 2) against the generation of plumes following supercontinent formation. Table 3 provides a list of the different models tested, which include running simulations that deviate from the base model in viscosity profile, thermochemical pile thickness and resolution. In order to simulate global-scale tectonic processes, we prescribe plate motions at the surface of the model to re-create a simplified mechanism for the formation of a supercontinent (e.g. a degree-1 mantle flow): (1) we first apply 4 cm a⁻¹ of total surface motion convergence (2 cm a⁻¹ per side) with an initial large-scale downwelling (subduction) simulating the assembly of plates (Fig. 5a), and subsequently (2) the cessation of this initial downwelling is followed by the initiation of subduction zones on the margins of a stationary portion of the surface (i.e. circum-supercontinent subduction, Fig. 5b–d). Following the cessation of initial subduction, the surface velocity within the model supercontinent is set to zero, while the convergent velocity on its margins remains at 2 cm a⁻¹. This two-step process is shown in Figure 5, where three different sizes of the supercontinent are provided (i.e. portions of stationary material covering 20%, 30% and 40% of the surface). The first step lasts for 400 myr in order to allow time for the formation of plumes and slabs from the laterally homogeneous initial condition (e.g. Fig. 4).

Following the formation of the model supercontinent, the simulation runs for another 300 myr (i.e. 700 myr in total) and the location of a plume arriving at the base of the lithosphere in regions exterior to the supercontinent (i.e. in the oceanic domain) is analysed. Given the timeframe of the formation of the supercontinent Pangaea to the present-day is c. 300 myr, we have used this temporal range as a cut-off for our simulations. Figure 6 shows the process for analysing the full 2D simulation, with Figure 6a showing the subduction zone that would bring together a model supercontinent. The overall northward-dipping subduction of the Palaeotethys that contributed to Pangaea formation is considered as the interior subduction (Fig. 1a; Murphy et al. 2009). After 400 myr, this subduction zone ceases and is replaced by circum-supercontinent subduction around a stationary region covering 30% of the Earth's surface (e.g. Fig. 6b). Exterior plume positions forming after supercontinent formation are then calculated using the great circle distance from the edge of the circum-supercontinent subduction zone to the centre of the plume tail (Fig. 6c, d). As an example, the base model with a 30% stationary surface generates exterior plumes 3781 km away from the edge of the supercontinent 136 myr after supercontinent formation (Fig. 6c) and 1334 km from the supercontinent edge 249 myr after supercontinent formation (Fig. 6d). The exterior plume timing is manually calculated based on the visual inspection of when the plume head starts to spread out beneath the surface boundary layer (our model lithosphere).

For this study, the simulations and analysis of plumes forming after supercontinent formation (i.e. after the cessation of interior, Palaeotethyan analogous subduction) are further split into two categories based on their location. In the simplified 2D simulations, we base this classification on the initial subduction polarity, i.e. the polarity of the angled downwelling (e.g. Fig. 1), which also relates to the geographical longitude (i.e. eastern and western hemispheres). In particular, the initial subduction that forms the supercontinent can be split into an overriding and a subducting plate: the western hemisphere of Pangaea (Laurentia, Baltica, Siberia, Amazonia) is described as the overriding plate and the eastern hemisphere is the subducting plate (Fig. 1a). In other words, we identify a western and an eastern margin of Pangaea, where the western side corresponds to the North and South American subduction boundary and the eastern refers to the Australasian subduction zones (Fig. 1b). This definition is used to geographically label the overriding and subducting plates in our 2D models (Fig. 6) to allow for comparison with Pacific LIPs forming closer to the Australasian (eastern) or American (western) subduction zones.

The parameters of the low-resolution 2D setup (Model 6) are also applied to a 3D global mantle convection model. The changes between Model 6 and the base model include a reduction in resolution in the bottom 250 km and top 80 km of the model, and the removal of particles to track the compositions within the simulation. The computational time for a 3D spherical model spanning supercontinent formation and dispersal is significant, and therefore the use of particles to track composition was removed to reduce computational time. However, tests in the 2D setup using this reduced resolution profile did not provide significant changes to the model results (see 2D results section). As a result of this 2D testing and the computational expense of the global models, we only provide one 3D model here.

The 3D simulation setup is similar to the 2D setup in that an initial formation period is applied with one interior subduction zone at 4 cm a⁻¹ total convergence (Fig. 4a). However, in the 3D models this time period is 300 myr, which is 100 myr shorter than the 2D simulations in order to save computational time. The goal of this initial formation period in 3D is to provide sufficient time for the mantle to mix from its initial condition. After 300 myr, the boundary condition at the surface is then provided by global plate reconstructions of Pangaean formation and dispersal, rather than the simple 2D setup of a stationary supercontinent. In this study, we use the Matthews et al. (2016) plate reconstructions, which extend back in time 410 myr (i.e. to prior to supercontinent formation) to the present day and thereby encompass the formation and break-up of Pangaea.

We apply the Matthews et al. (2016) surface velocities at 1 myr intervals, with a resolution of 1° × 1°. As the boundary condition at the surface is provided by plate reconstructions and is therefore kinematic (i.e. it does not change in accordance with the force balance of the underlying mantle), we can show the impact of subduction-driven mantle flow on the positioning of mantle plumes. The 3D simulation is then analysed for plumes in and around the potential sites of the Shatsky Rise (170 Ma), Ontong Java (122 Ma) and Caribbean (95 Ma) LIPs.

As our study is focused on the impact of subduction on plume generation and applies time-dependent boundary conditions related to plate reconstruction histories to determine mantle flow, we do not need to prescribe continental and oceanic material at the surface as in previous studies (e.g. Coltice et al. 2007, 2009; Phillips et al. 2009; Phillips and Coltice 2010; Yoshida 2010; Rolf et al. 2012; Arnould et al. 2020). Although more realistic, the impact of chemically distinct continental material would not be a primary driver of deep mantle plumes (e.g. Heron and Lowman 2010, 2011, 2014; Yoshida 2013; Heron et al. 2015) and therefore has not been applied. However, the model uses temperature-dependent viscosity, which forms strength differences within the top boundary layer. These are a proxy for the outlines of plates bounded by the ridges and subduction zones produced from the applied plate reconstruction velocities. In that respect, we can generate areas of reduced heat flux due to cold material that thickens in response to plate tectonic processes. As a result, our simulations yield areas of the surface that could be described as a model continent (or interior of a supercontinent) and a model ocean with mid-ocean ridges (exterior of a supercontinent) (see 3D results section for an example).

In this section we outline the results from the models shown in Table 3, starting with the 2D simulations.

Figure 7 shows the position of the exterior (i.e. sub-superoceanic) plumes for the base model with supercontinent coverage of 20% (Fig. 7a, d), 30% (Fig. 7b, e) and 40% (Fig. 7c, f). In these models, the change in mantle flow due to the formation of the circum-supercontinent downwelling (subduction) culminates with the production of the first exterior plume near the eastern boundary for all supercontinent coverage percentages c. 140 myr after the final supercontinent formation (i.e. shutdown of interior subduction) (Fig. 7a–c). These eastern plumes form at c. 4400 km (20% coverage, Fig. 7a), 3700 km (30%, Fig. 7b) and 3800 km (40%, Fig. 7c) from the peripheral subduction zone. The western plumes form later than their eastern counterparts (Fig. 7d–f), with their distance from the downwelling being more variable (20%: 3600 km; 30%: 1300 km; 40%: 3600 km). The delay in the initiation of the plumes may relate to the western side of the supercontinent being on the overriding plate during supercontinent formation (e.g. Fig. 1). Therefore, plume formation beneath the western margin is suppressed by the cold subducted interior ocean lithosphere covering the core–mantle boundary in that area. In addition, eventual overturn of the subducted slab on the eastern margin allows for a stronger interaction between the eastern slab and the sub-superocean thermochemical structure.

The other models performed in this study (Table 3) yield a similar behaviour to the 2D base model (Fig. 7). Figure 8 shows an overview of the timing of the exterior plume arrival plotted against distance from the peripheral subduction zone for both the eastern and western regions. Plumes that formed due to the initial condition, and plumes that would have occurred if no supercontinent had formed (e.g. plumes forming due to the initial subduction flow), are not included in our summary. Plumes included in this study are only those related to the new subduction initiation. Eastern plumes typically form before western plumes, and show a smaller variance in distance from the downwelling. The average timing of the exterior plumes post-supercontinent formation is 157 ± 6 myr for eastern plumes (n = 21) and 229 ± 9 myr for western plumes (n = 12). The standard deviation for the eastern and western ages is 28 and 32 myr, respectively. The average distance from the subduction for the 2D models was 4132 ± 221 km for eastern plumes and 3252 ± 388 km for western plumes. The standard deviation for the eastern and western locations is 766 and 1347 km, respectively. The timing cut-off for measuring plumes is set at 300 myr after supercontinent formation. There are consequently fewer western plumes during this timeframe due to the delay in plume formation (with plumes forming after the cut-off time of 300 myr; Fig. 7d).

Figure 9 shows the western plume positions for the different models described in Table 3, for 30% supercontinent coverage (Model 7 is not shown since no plume(s) formed near the western margin within 300 myr of supercontinent formation, only forming after this cut-off period). As illustrated in Figure 9, there is variation in the plume positions, with Model 3 (decrease in dense layer thickness; Fig. 9c) and Model 5 (increase in maximum viscosity; Fig. 9e) yielding plumes with the greatest distance from the peripheral subduction zones. Decreasing the resolution in Model 6 (Fig. 9f) produces a similar result to that of Model 1 (Fig. 9a).

Figure 10 shows the eastern plume position for the different models with 30% coverage, and highlights the summary shown in Figure 8 (earlier plume arrivals and less variance between the distances than western plumes). The only notable difference between these models is that the decrease in dense material initial thickness in Model 3 (Fig. 10c) produces more vigorous mantle convection. The initially thin dense layer becomes more readily mixed into the mantle, producing more plumes emanating from the core–mantle boundary. However, despite these different dynamics, the plume position for Model 3 forms closest to the eastern boundary but is similar in timing to the other models (Fig. 11, diamonds).

When changing the continental coverage, no plumes form for the 20% coverage on the western side during the 700 myr model run (Fig. 12). Furthermore, the 30% continental coverage produces western plumes at shorter distances from the subduction zone than the 40% coverage (Fig. 12). The interaction of the newly formed circum-supercontinent subduction on the western side either captures a plume beneath the supercontinent (20%), or pushes the thermal instabilities out beneath the exterior on the oceanside (30%, 40%). Conversely, the consistency of the eastern side results in producing a narrow band of plume timings and distances relates to the absence of subducted lithosphere from the interior ocean, allowing the new circum-supercontinent subduction to disrupt the mantle flow and to form deep thermal instabilities that lead to mantle plumes (regardless of initial model setup; Figs 11 & 12).

Figure 13 shows the evolution of the 3D numerical simulation that prescribes 410 myr of plate reconstruction history as a time-dependent velocity boundary condition (Matthews et al. 2016). We present snapshots of the thermal structure at a depth of 70 km from the surface, which is sufficiently shallow to resolve the colder continental lithosphere and sufficiently deep to capture any thermal plume arriving beneath the oceanic lithosphere from the core–mantle boundary. The snapshots show the southern Pacific in the left column and the western boundary of the Americas in the right column, both in the present-day reference frame. Figure 13a, b shows the formation of our model Pangaea at 300 Ma, bounded by peripheral subduction zones (Fig. 1b). We then analyse the models at four key times relating to the formation of oceanic LIPs: 170 Ma, 145 Ma (Shatsky Rise), 122 Ma (Ontong Java) and 95 Ma (Caribbean LIP). Figure 14 shows temperature cross-sections that correspond to these four key times.

Figure 13d shows a warm thermal anomaly in the cold lithosphere off the coast of South America at 170 Ma (as highlighted by the white dashed circle). This thermal anomaly can be seen to be related to the arrival of a plume head (Fig. 14a). The location of this plume is in the approximate position of the present-day Galapagos hotspot. The plume is also in this location during the timing of the CLIP (Figs 13j & 14d). Our models indicate that the exterior plume activity off the coast of South America is related to circum-supercontinent subduction and occurs at a time period that links the Galapagos hotspot to the CLIP, as previously discussed in Courtillot and Renne (2003) and Hoernle et al. (2004).

Figure 13e shows the formation of a thermal anomaly in the Pacific at the approximate formation timing and location of the Shatsky Rise (145 Ma). Analysing a thermal cross-section at this time period shows the formation of a plume interacting with a subduction zone situated c. 7200 km away at the eastern boundary of the Pacific (Fig. 14b). Furthermore, the dense anomaly on the core–mantle boundary below the Pacific is clearly visible and shows a plume rising from its margin (Torsvik et al. 2010) despite our thermochemical piles being mobile and not fixed in location (Zhang et al. 2010; Zhong and Rudolph 2015; Heron et al. 2021; Flament et al. 2022). As indicated by the white lines in Figure 14b, we interpret this plume to be an exterior upwelling related to the nearby eastern subduction.

Figure 13g, h shows the southern Pacific and western America boundary at 122 Ma and highlights a thermal anomaly arriving near the Antarctica subduction zone. In the cross-section of the mantle's thermal field (Fig. 14c), a plume related to this lithospheric anomaly is shown forming along the southern margin of the dense layer. We interpret this interaction to be due to the exterior return flow related to the Antarctic subduction zone, which occurs at the time of the OJP LIP. Owing to the angle of the cross-section and the 3D nature of the models, it is difficult to show the closeness of the plume and subduction zone. However, the great circle distance between the surface expression of our OJP plume in Figures 13g and 14c and the Antarctica subduction zone is c. 4300 km.

Results from our modelling show the impact of modifying global subduction patterns on the location and timing of plumes exterior to a supercontinent (Fig. 3). Although a number of studies have previously modelled or discussed plume development after supercontinent formation (Sobolev et al. 2011; Biggin et al. 2012; Heron et al. 2015; Flament et al. 2017; Dal Zilio et al. 2018; Dannberg and Gassmöller 2018; Zhang and Li 2018), the mechanisms that affect the timing and location of exterior plumes are not well understood. Our study offers a formal analysis of the generation of exterior plumes following supercontinent formation and their relationship to circum-supercontinent subduction (Fig. 3). We provide an estimate of the timing and location of such plumes and their variability (Fig. 8). Through the 3D simulations, we also visualize this subduction–plume interaction on a global scale. Below we discuss the potential relevance of these findings to the exterior LIPs that currently remain in the Pacific Ocean (Fig. 2): the western LIPs of Ontong Java and the Caribbean, and the eastern Shatsky Rise LIP.

The OJP is the world's largest oceanic plateau and is linked to a massive volcanic event that significantly altered global climate and ocean oxygen levels, which was responsible for a mass extinction event (Zhang et al. 2020). The position and geometry of the OJP has changed significantly over time associated with the drifting Pacific plate. The three currently separated plateaux of Ontong Java, Manihiki and Hikurangi were originally formed in a single complex known as the greater OJP (Taylor 2006; Chandler et al. 2012). While the greater OJP is thought to have been emplaced by two major magmatic events at c. 122 and 90 Ma (Mahoney et al. 1993; Tejada et al. 1996), the cause of the magmatism is still debated (Isse et al. 2021). Two main hypotheses relate OJP formation to a mantle plume (i.e. melting anomalies arising from hot diapirs ascending from the deep mantle: Tarduno et al. 1991; Isse et al. 2021) or to a plate boundary interaction (i.e. melting anomalies fuelled by shallow mantle processes associated with plate tectonics: Nakanishi et al. 1999; Frey et al. 2000; Sager et al. 2019).

Although a number of studies have analysed the viability of a mantle plume origin, including studies using seismic evidence (Isse et al. 2021), there is little discussion of the mechanism that would have generated the upwelling. The Louisville hotspot, for example, is often discussed as a potential source of the OJP (Tarduno et al. 1991; Chandler et al. 2012), but the mechanism that caused the Louisville plume to form at its specific time and location is not known.

The reconstructed initial emplacement position of the OJP complex is at c. −38° latitude and 219° longitude (Torsvik et al. 2006, 2008). The closest subduction zone to this region would have been that of Antarctica's western boundary near the South Pole (−75°, 219°), as seen in Figure 1c (Seton et al. 2012). The distance from the OJP centre to this subduction margin is calculated to be c. 4000 km (using a great circle path). In our 3D model that uses plate reconstruction velocities to move the surface, we form a plume at 122 Ma and at a distance of c. 4300 km from the Antarctic subduction zone.

To put our 2D models into geological context, we need to modify the OJP emplacement age to a timescale relevant to Pangaea formation. The formation of Pangaea occurred during the Late Carboniferous period and it is difficult to apply one single age to the formation. However, we assume the cessation of the interior subduction zone that helped form Pangaea occurred at c. 320 Ma (Pastor-Galán 2022). We use this single time as a reference in order to create an approximate timing of our plumes relative to LIP magmatism. The approximate time lag to OJP magmatism following Pangaea formation is 198 myr (320–122 Ma). However, given the uncertainty related to Pangaea formation, the timing has an uncertainty in the range of ± 10 myr. In our 2D models presented above, the western plumes (which would be most dynamically relatable to the OJP) have an average timing and position (with standard deviation) of 229 ± 32 myr and 3252 ± 1347 km (Fig. 8), which fits well with the observed spatial and temporal range for the OJP.

The CLIP has a similar chemical composition to the OJP (Hauff et al. 2000), with the main volcanic emplacement occurring between 95 and 88 Ma (Fig. 1d). Compared to other LIPs that formed after Pangaea formation, the CLIP is among the most recent (Torsvik et al. 2006, 2008). Although suggested formation mechanisms for the CLIP vary, the Galapagos hotspot has been proposed as a potential deep mantle plume source (Courtillot and Renne 2003). Our 3D models show the formation of a plume in the approximate location of the CLIP from 170 Ma (Figs 13c, d & 14a).

The reconstructed initial emplacement position of the CLIP was at c. 0° latitude and 250° longitude (Scaife et al. 2017). At that time, the closest subduction zone to this was along the western boundary of North and South America (0°, 250°), as seen in Figure 1d (Scaife et al. 2017). The great circle distance from the CLIP to this subduction margin is c. 2200 km and, using 320 Ma for the cessation of interior subduction, the age of the CLIP with respect to Pangaea formation is 225 ± 10 myr (320–95 Ma). The CLIP therefore appears relatively close to the subduction zone, but within the margins of plumes observed in our 2D models, and with a timing that is consistent with our model results (Fig. 8). Our 3D model generates a plume c. 1000 km from the subduction zone of South America and further south of the equator than the CLIP emplacement.

The Shatsky Rise is the Earth's third largest oceanic plateau (Geldmacher et al. 2014), forming in the exterior to Pangaea at c. 145 Ma (Mahoney et al. 2005). Like the OJP, two main hypotheses exist for the formation of the plateau: a plume origin or plate boundary interaction (seafloor spreading) (Sager et al. 2013, 2019; Heydolph et al. 2014; Li et al. 2016).

The reconstructed initial emplacement position of the Shatsky Rise was at c. −4° latitude and 219° longitude (Torsvik et al. 2006, 2008). The closest subduction zone to this region would have been in the eastern Pacific region on the margins of Eurasia (−4°, 150^o). The great circle distance from the Shatsky Rise to the subduction margin is therefore nearly 8000 km. The timing post-Pangaea for this LIP is 175 ± 10 myr and therefore consistent with our 2D results showing earlier eastern plume formation (Fig. 8). While almost all of our 2D models also feature plumes at similarly substantial distances from the subduction margin (e.g. Figs 6c, d & 7), they also highlight that this distance is likely too great for the plume to be caused by the direct impact of the circum-supercontinent subduction and it instead may have formed as a natural instability of the thermal boundary layer (which we did not analyse here). However, our 3D models still show a potential interaction with an eastern subduction zone, producing a mantle return flow and an exterior plume at a distance of c. 7200 km (Fig. 14b). For the formation of the Shatsky Rise, Sager et al. (2016) suggested a similar return flow and exterior plume mechanism related to the subducted slab and its interaction with the core–mantle boundary, which is demonstrated in our model results.

The 2D and 3D models presented here are not intended to be a definitive study of the range in the potential timing and location of different oceanic plumes. There are uncertainties in the material properties of the mantle that may affect the calculation of the plume to subduction zone distance, as well as in the timing of the plume head arrival. Moreover, the timings and locations of the real-world examples presented here (OJP, CLIP and Shatsky Rise) also have uncertainty associated with them. In addition, our 2D models feature stationary subduction zones, which is not realistic when compared to the evolution of a typical subduction zone (e.g. Matthews et al. 2016). Although the western subduction boundary is relatively stationary in a plate tectonic reference frame during the Pangaea life span, it is not perfectly stationary as in our 2D models. This is also evidenced by the different timing and location of plumes in our 3D models, where ancient subducted oceanic lithosphere is distributed more widely in the lowermost mantle. However, by using a simplified set of surface boundary conditions we are able to focus on the first-order mechanisms linking supercontinents and plumes.

Despite these uncertainties, our 2D and 3D models are consistent with the timing and location of the three exterior plumes presented here. The findings therefore provide a viable link between the timings and locations of our modelled exterior plumes and observed oceanic plateaux. As such, our study successfully highlights subduction as a potential mechanism for the formation of oceanic LIPs exterior to a supercontinent.

The modelling presented here proposes a new mechanism explaining the timing and location of exterior plumes in response to supercontinent formation. However, there are still various aspects related to the numerical experiments and geological analyses that would be intriguing to explore further. For example, we present a limited number of 2D numerical experiments and there are a number of key parameters that could be tested to expand the breadth of these interpretations. In particular, our subduction zones are static in these 2D models. Trench advance or retreat rate, in addition to the plate convergence rate, would impact the slab morphology and therefore the timing and positioning of these plumes. This may be a factor in the 3D numerical modelling results for the Shatsky Rise, which fall outside the range provided in the 2D simulations (i.e. at a distance >7000 km). Furthermore, our 2D setup presents a simple model of Pangaean tectonics without any post-formation interior subduction. Plate reconstruction history does indicate subduction within the interior of Pangaea related to the Palaeotethys and Neotethys oceans. An avenue for further exploration would be to determine the potential impact of such events on both interior and exterior plume formation. In addition, including chemically distinct continental lithosphere within the supercontinent area would be beneficial in creating more realistic subduction dynamics, which could impact the interaction of the subduction zone and the mantle flow. Despite these caveats, however, we maintain that our simplified 2D models provide a primary step in investigating the origin of plumes exterior to the supercontinent.

Although we present models with and without dense piles of material near the base of the mantle, future work could explore the evolution and impact of these piles on long-wavelength mantle convection. In the 2D models presented here, dense piles form under the exterior ocean realm rather than below the supercontinent, which may impact plume formation and timing. In addition, we do not explore the potential impact on our findings of the density and geometry of the deep LLSVP material (and only present the end-member cases). Denser and less mobile piles of material could have a control on plume location and plume generation, given their potential impact on the core–mantle heat flow. Furthermore, the impact of changing the initial condition of our numerical models (e.g. temperature profile) would also assist in placing our findings within the context of Earth's secular evolution.

Future work could additionally explore and compare the forces exerted by both exterior and interior plumes on the plates that they impinge on (e.g. Stern 2007). This might address the role of exterior plumes in supercontinent break-up through far-field forces.

Once a further understanding of the mechanisms involved in the formation of plumes from subduction-related mantle flow is achieved (Arnould et al. 2020), it would be possible to characterize all available LIP and hotspot data.

The numerical models presented here show, for the first time, that subduction related to the supercontinent cycle can reproduce the location and timing of exterior (intra-ocean plate) mantle plumes, including those related to the LIPs (or oceanic plateaux) of the Ontong Java, Caribbean, and Shatsky Rise formations (Figs 1, 8, 13 & 14). The results offer insights into the formation of other oceanic LIPs, and suggest that the impact of circum-Pangaean subduction zones can extend away from the supercontinent into the mantle beneath the surrounding ocean (Fig. 3). This study adds to a body of research that highlights the importance of considering how a supercontinent impacts mantle convection at every stage of the supercontinent cycle.

We would like to thank Daniel Pastor-Galán and an anonymous reviewer for their thoughtful comments on this paper. A number of figures in this manuscript were generated using Scientific Colour Maps (Crameri 2021), the online figure repository s-ink.org (Crameri et al. 2022) and with accessibility in mind (Crameri et al. 2020).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

PJH: conceptualization (lead), data curation (lead), formal analysis (lead), funding acquisition (lead), investigation (lead), methodology (lead), project administration (lead), resources (lead), software (lead), supervision (lead), validation (lead), visualization (lead), writing – original draft (lead); EG: conceptualization (supporting), investigation (supporting), methodology (equal), visualization (lead), writing – original draft (supporting); GES: conceptualization (equal), data curation (equal), methodology (equal), software (equal), writing – original draft (supporting); JD: methodology (equal), writing – original draft (supporting); RG: data curation (equal), methodology (equal), software (equal), writing – original draft (supporting); EM: writing – original draft (supporting); AS: investigation (equal), writing – original draft (supporting); RNP: writing – original draft (supporting); RDN: conceptualization (supporting), investigation (supporting), writing – original draft (equal); JBM: conceptualization (equal), investigation (supporting), writing – original draft (equal).

This research was enabled in part by support provided by Compute Ontario (computeontario.ca) and the Digital Research Alliance of Canada (alliancecan.ca) for allocations to PH and RP. We also acknowledge support from a Natural Sciences and Engineering Research Council of Canada Discovery Grants for PH, RP and JBM. We thank the Computational Infrastructure for Geodynamics (geodynamics.org), which is funded by the National Science Foundation under the awards EAR-2149126, EAR-0949446 and EAR-1550901, for supporting the development of ASPECT. GES acknowledges support from the Research Council of Norway through its Centres of Excellence funding scheme, Project Number 223272, and through its Young Research Talent scheme for ‘POLARIS – Evolution of the Arctic in Deep Time’, Project Number 326238. JD and RG were supported by the NSF grants EAR-2054605 and EAR-1925677.

Data and input files related to the numerical modelling of the current study are available in the GitHub repository, https://github.com/heronphi/Heron_Super22.