Weak orogenic lithosphere guides the pattern of plume-triggered supercontinent break-up

The importance of nonrigid geological features (such as orogens) inside tectonic plates on Earth’s dynamic evolution lacks thorough investigation. In particular, the influence of continent-spanning orogens on (super)continental break-up remains unclear. Here we reconstruct global orogens and model their controlling effects on Pangea break-up. We show that while loci of Pangea break-up are linked to mantle plumes, development of continental rifts is guided by orogens. Rifting at Central Atlantic is driven by the modelled plume responsible for the Central Atlantic Magmatic Province (CAMP) within Pangea-forming orogens. South Atlantic rifting is controlled by necking between Pangea- and Gondwana-forming orogens with the assistance of plume-induced lithospheric weakening. Without CAMP-induced weakening, South Atlantic rifting fails between the West African and Amazonian cratons, but occurs between the West African and Saharan cratons instead. Our modeling on Pangea break-up is able to recreate present-day continental geometry through the combined effect of orogens and plume center-locations. The location of orogenic belts has a strong influence on the rifting of supercontinents after mantle plumes trigger the initial break-up and contribute to lithospheric weakening, according to numerical simulations of the break-up of Pangea.

T he break-up of supercontinents is a fundamental process of global tectonic evolution, and a key component of supercontinent cycles [1][2][3][4][5][6][7] . Recent studies focus primarily on the relative importance of the drag force due to subduction retreat [8][9][10] , the slab pull force due to far-field subduction 11,12 , and the push force of mantle plumes 7,13-21 . Other key factors such as non-rigid features that control the loci of continental rifting are not widely discussed, although it has been recognized that the loci of supercontinent break-up are influenced by the positions of both mantle plumes [18][19][20][21] and nonrigid features of the continental lithosphere, e.g. the global network of orogenic belts [22][23][24] . In this study, we model the process of supercontinent rifting implementing the position and geometry of pre-existing orogens, and compare the model results with geological records. The wellreconstructed geologic history of Pangea break-up makes it an ideal target for analyzing factors and processes controlling the loci of continental rifts during supercontinent breakup. Advances in 3D dynamic modeling allow the analysis of key factors influencing continental rifting during the break-up of Pangea with realistic geological boundary conditions and global orogenic belts.
The continental lithosphere consists of cratons surrounded by orogens. Cratons are essentially ancient continental lithosphere that has not been deformed for a significantly long period of time 25 . The physical properties of cratons such as strength, thickness, and viscosity are markedly different from those of orogens [26][27][28] . Conventionally, orogens are considered "weaker" than cratons due to anisotropy anomalies and/or major lithospheric shear zones, high crustal fluids, or thinner lithosphere. The latter is either intrinsic to the location or due to convective or gravitational removal of mantle lithospheric roots [29][30][31][32][33][34] . Orogens are therefore likely locations for subsequent continental rifting. For example, 45% of the~25,000 km of Gondwana rift margins developed along orogens 29,30 . During the opening of the Central and South Atlantic Ocean, most of the rifting occurred along orogens; e.g., the formation of the Central Atlantic occurred along the Alleganian-Variscan-Hercynian orogen, separating eastern North America and northwest Africa 28 , whereas the opening of the South Atlantic occurred roughly along with the Pan-African-Brasiliano belts such as the Araçuaí-West Congo orogen between the São Francisco and Congo cratons 32,35 (Fig. 1). Nonetheless, the opening of the northern part of the South Atlantic Ocean failed to follow the orogens. Instead, it splits the West African and Amazonian cratons, which were considered a coherent craton since at least the Mesoproterozoic 36 .
Another important factor that influences the rift loci is the location of mantle plumes 16,18,37 . Thermal erosion caused by plume impingement beneath the continental lithosphere can weaken the lithosphere 18,20 . The push force of uprising subcontinental mantle plumes may also provide a dynamic driving force for continental break-up 38 . In general, there is a clear correlation between downward projected large Igneous Provinces (LIPs) eruption sites of the past 300 Ma and the margins of the Large Low Shear Velocity Provinces (LLSVPs) at the base of the mantle [39][40][41][42] . The majority of LIPs are thought to have been generated by plumes that rose from the edges of the LLSVPs. As a result, many rifted continental margins formed during Pangea break-up have records of plume-induced LIPs, such as the Parana, North Atlantic Igneous Province (NAIP), Central Atlantic Magmatic Province (CAMP) and Karoo LIPs 19,20 . However, it is unresolved whether plumes alone can explain the loci of continental rifting as multiple mantle plumes (and associated LIPs) of similar ages exist along failed continental rifts (e.g., the Benue Trough, the Oslo Rift, and the Baikal Rift).
Here, we evaluate the roles of different factors, especially the geometry of global orogenic belts, in the break-up of Pangea through a 3D free-slip spherical global dynamic modeling. The model result demonstrates the combination of orogens of different ages, plume push force, and plume-induced lithospheric weakening control the initial rifting of Pangea. Although subduction retreat may also affect supercontinent break-up [8][9][10] , its effect in the opening the South Atlantic is relatively minor as no subduction retreat has been documented for the region during the modeled time period. This study represents a step forward in understanding how a continent's lithospheric heritage influences plume-induced continental break-up.

Results
Modeling the orogen-plume interaction and lithosphere rifting. The spherical model is constrained and driven by reconstructed subduction history. The model's simulation time starts at 250 Ma (see "Methods") when supercontinent Pangea is surrounded by a subduction girdle 2 . We focus on the model outcomes after 200 Ma when plumes break out and continental rifting occurs widely. Continental lithosphere older than 0.7 billion years ago (Ga) is treated with the same physical property as old cratons, whereas continental lithosphere younger than 0.7 Ga is treat as weak zones (orogens) (see Methods for further information). Our model also considers both partial melting 43-47 due to mantle plume and rheological weakening of the lithosphere subjected to melt percolation [48][49][50] (see "Methods"). Our standard model (Case 1) accounting for orogenic belts (Fig. 2) shows that the subducting slabs drive hot instabilities to rise at the boundary of African LLSVP. The generated mantle plumes first arrive under the central Pangean lithosphere between Laurasia and Gondwana ( Fig. 2a), causing rifting in the Central Atlantic region along Alleganian-Variscan-Hercynian orogens, as shown by the diverging plate motion velocity vectors (Fig. 2a) and high values of the second invariant of strain-rate tensor (_ ε II ) (Fig. 2b). This process continues until the plume between Laurasia and Gondwana disappears (Fig. 2c). Magma produced by the plume results in the thermal and melt weakening of the cratonic lithosphere between the West African and Amazonian cratons, shown as an increase in _ ε II (Fig. 2c, d). Later, the southern thermal instability beneath South America (Fig. 2a) grows to form a new plume beneath South Atlantic (Fig. 2e). Due to this plume event and its induced melt migration into the northern tip of the Pan-African-Brasiliano orogens, cratons in Africa start to move to different directions relative to South America as represented by the diverging motion vectors. South Atlantic opening initializes at the southern part as _ ε II in this area is evidently higher than the other areas (Fig. 2e, f). With the continuing influence of this mantle plume from the south, plume-induced melts along the western boundary of the orogens, and the retreating slabs at the west, the motion direction of Africa become distinctly different from that of South America (Fig. 2g). The cratonic lithosphere between the Amazonian and West African cratons is destroyed and _ ε II in this area increases sharply, leading to the break-up between the Amazonian and West African cratons as the final stage of South Atlantic opening (Fig. 2g, h). Therefore, our model produces a two-stage opening for the Central and South Atlantic Ocean: (1) the opening of Central Atlantic by the separation of Laurasia from Gondwana ( Fig. 2a-d) and (2) the south to the north opening of the South Atlantic between the South America and Africa cratons (Fig. 2e-h).
We next compare the effects of orogens of different ages on the Atlantic rifting (Case 2). In this case, we assign orogens of <1.6 Ga ages into two groups (1.6-0.7 Ga and <0.7 Ga, Supplementary  Fig. 1b) with various physical properties (see "Methods") and repeated the calculation. The result shows a similar opening process from Central to South Atlantic. This illustrates that orogens older than 0.7 Ga did not have much influence on Pangea's break-up, especially the rifting along Amazon and West Africa craton, because orogens of 1.6-0.7 Ga are not weak enough to alter the plume and rifting behaviors.
Our standard case shows that the South Atlantic rifting follows the Pan-African-Brasiliano orogens but interestingly migrates across the joint between the West African and Amazonian cratons instead of the orogen between West African and Sahara cratons. We hence add an exploration model (Case 3) to investigate the influence of plume-induced thermal and melt weakening on the rifting process. Case 3 is the same as Case 1 except that we removed the effect of CAMP plume-induced lithospheric weakening effect through heat and melts. Case 3 shows a different evolution path after 140 Ma. Similar to Case 1, South Atlantic starts to open at the middle part as _ ε II in this area is evidently higher than the other areas (Fig. 3b). However, the rifting migrates to the northeast along the weak zone (orogen) between West African and Sahara cratons (Fig. 3c), and eventually cuts through between the two cratons and extends to the Tethys Ocean (Fig. 3d).

Discussion
Triggering effect of the plume in Pangea's initial break-up. Plumes with the associated partial melting appear to be necessary to initiate the break-up of Pangea and determine the timing and locations of continental rifting. First, our model suggests that the locations where supercontinent break-up happens are determined by plume locations (Fig. 2). Our model clearly shows three plume pulses, corresponding to the CAMP (Fig. 4a), Karoo (Fig. 4b) and Parana ( Fig. 4c) LIPs, respectively. At the first stage of our model, the curved subduction boundary spans from North America to South America 51-53 (Fig. 4a), inducing intensive downwellings and a strong plume at Central Atlantic. During the second stage, the southern thermal instability beneath South America (Fig. 2a) grew to form a new plume at South Atlantic (Fig. 2e). This plume also benefits from the subduction zone curvature along the west margin of South America 51-53 but it forms after the first plume eruption at Central Atlantic. The second plume pulse initializes the rifting of South Atlantic at first (Fig. 2e). All thermal instabilities show a linear distribution influenced by the subduction zone and the LLSVP boundary. We note that the modeled location of the Karoo LIP is further southeast (Figs. 4b, 5b) than its actual location, and the modeled time for the Karoo and Parana LIPs are~10 Myr later than in the geological record. Furthermore, despite that the lithosphere above the Karoo plume is weakened and that its viscosity has a sharp drop after Karoo plume's eruption (Figs. 4b, 6 and Supplementary Fig. 2a), the model results show a slower break-up than the geological record. This inconsistency is due likely to the less well-constrained subduction location and angles at southern Pangea. Additional data parameterizing this subduction system are required to produce a more realistic break-up of southern Pangea.
The plume associated with the CAMP erupts right beneath the orogens between Laurasia and Gondwana at around 195 Ma in Case 1 (Figs. 4a, 5a). As it approaches the surface, the plume causes melt production which further weakens the orogens, leading to the eventual opening of Central Atlantic 19,54 . Similarly, the initiation of South Atlantic rifting occurs where the plume corresponding to the Parana LIP erupts 55,56 . This is due to the proximity of the plume and orogens there at~135 Ma 19,57 (Figs. 4c, 5c), causing the ascending plume head to be funneled toward the thinned orogenic lithosphere and the eruption of the flood basalts, and triggering South Atlantic rifting. Our model results are reasonably consistent with the positions of the actual LIPs (Figs. 4a-c). The velocity distribution after the opening of Central Atlantic (Supplementary Fig. 2a) and South Atlantic in our model ( Supplementary Fig. 2b) is similar to the reconstruction result of Muller et al. 52 , with the mean angle differences of velocity distribution between our model and the reconstruction being smaller than 50° (Supplementary Table 3). Furthermore, the Case 1 modeled root mean square velocities for 200-100 Ma (Supplementary Fig. 2c) are also resemble that of the reconstructions 58 .
Guiding role of young orogens in Pangea rifting. Our results also suggest that orogens younger than the previous supercontinent play the main guiding role to the rifting and break-up of a supercontinent, and largely determine the shapes of the new continental margins. Previous work 23,24 demonstrated that orogens, as weak zones in continental lithosphere compared to cratons, are preferred locations of continental rifting and break-up. In this study, the rifting of the Central and South Atlantic Ocean occurs along the orogens after plume eruptions (Fig. 2). Especially, for the opening of South Atlantic, the plume generating the Parana LIP is located beneath the Sao Francisco craton rather than the orogens between Africa and South America, yet the rifting still occurs along the orogens (Figs. 2g, 5d). Furthermore, a comparison of Case 1 and Case 2 demonstrates that orogens younger than the previous supercontinent (generally younger than 600 Ma) are the preferred locations for rifting and break-up because they have thinner lithosphere, lower viscosity, smaller yield stress (Supplementary Table 2), and are prone for plume-induced thermal and melt weakening.
An interesting geological observation is that the South Atlantic rifting cut through the West Africa-Amazon craton but not along the pre-existing orogen between the West African and Sahara cratons (the Dahomeyide-Pharuside belt). In our standard model (Case 1), the necking between CAMPand Parana-induced lithospheric weakening reproduced such  Necessity of melt-induced weakening for the lithosphere rifting. By calculating melts induced by plumes, the effect of rheological weakening on continental lithosphere as observed in many cratons 25 is introduced into our model. We find that partial melting induced by uprising plumes at the bottom of the lithospheric is necessary for continental break-up. For example, the northern South Atlantic opening occurred between Amazonian and West African cratons primarily caused by the weakening efforts of CAMP plume-induced melts (Fig. 6a). The CAMP partial melts accumulated along the curvature of Alleganian-Variscan-Hercynian belts first weakens the northern cratonic lithosphere of the West African and Amazonian cratons during the opening of Central Atlantic (Fig. 6b).
During the eruption of Parana plume at South Atlantic, a significant amount of melts appears in the South Atlantic. Both the melts and the rift migrate from south to north, and eventually connect to the earlier weakened lithosphere between the West African and Amazonian cratons by CAMP-induced melts, leading to the final separation of the Amazonian craton from the West African craton (Fig. 5d). We demonstrate through Case 3 (Fig. 3a-d) that if we remove the effect of lithospheric weakening by melts generated by the CAMP plume, then the model is unable to regenerate the opening of northern South Atlantic between West African and Amazon cratons. In model Case 4 (Fig. 3e-h) we remove the effect of melts generated by both the CAMP and the Parana plumes. In this case, the South Atlantic fails to open all together. This further demonstrates the importance of plume-induced melting in the opening of supercontinents. The continental margins produced by our model is in good agreement with the actual continental coastlines at Central and South Atlantic, indicating that orogens not only guide the positions of supercontinent break-up but also determine the geometry of continental margins. Overall, our work demonstrates that orogens and the loci of plume eruption control the loci of Atlantic rifting 23,33 . The break-up of Pangea is triggered by mantle plumes and guided by old lithospheric scars (<0.7 Ga orogens), which together provide the dominant influence on the shapes of a large portion of present-day continents. With more sophisticated global orogens, LIP reconstructions 60,61 , and plate motion models, geodynamical models should be able to more realistically reproduce the shapes of all continents through supercontinent cycles.

Methods
Model set-up. Our three-dimensional spherical mantle convection model with continental blocks uses CitcomS 62 , assuming an infinite Prandtl number. We use the Boussinesq approximation to solve the conservation equations of mass, momentum, energy and composition. The non-dimensional governing equations for mantle with different compositions 62-64 are: where u is the velocity vector, P the dynamic pressure, η the viscosity, _ ε the strain-rate tensor, δT the temperature perturbation, g the gravitational acceleration vector, T the temperature, t the time, H the internal heat generation rate and C i the composition used to define continental blocks and LLSVPs (0 ≤ C i ≤ 1; 0 for the mantle, 1 for the chemically distinct materials-cratons if C 1 = 1, orogens if C 2 = 1, LLSVPs if C 3 = 1). Ra and Ra c_i , the Rayleigh number and the chemical Rayleigh number, respectively, are defined as: where α is the thermal expansivity, ρ the density, ΔT the temperature difference between the top and bottom boundaries, R the radius of the Earth, η ref the reference viscosity, κ the thermal diffusivity, and Δρ c_i the density contrast between chemically distinct materials (cratons, orogens and LLSVPs) and oceanic mantle (Supplementary Table 1). Notice that the definition of Rayleigh number is based on Earth's radius rather than the thickness of the mantle, with that based on the former being~10 times of that based on the latter. The supercontinent, which includes cratons and orogens, is modeled as a chemically distinct material with different density and viscosity. The viscosity formula used in our model relies on depth, temperature, composition and melt fraction 65 .
where η r (r) is the depth-dependent viscosity prefactor, η c (C i ) the compositional prefactor for chemically distinct materials (100 for cratons, 10 (or 30) for orogens and 1 for LLSVPs), f i the fraction of the ith composition in an element, E the activation energy, α F the melt fraction factor, F the melt fraction which defines the degree of melting expected by equilibrium batch melting prior to the exhaustion of clinopyroxene (cpx). The other rheological parameter used is plasticity 66 , which is controlled by the yield stress. If the convective stress exceeds the yield stress, the viscosity is reduced to the yielding viscosity. It follows the equation of: where η y ; σ y and _ ε II are the nondimensional yielding viscosity, yield stress, and the second invariant of strain-rate tensor, respectively. The final effective viscosity η e is a combination of Eqs. (7) and (8), given by: This plasticity formulation has been commonly used in previous studies [67][68][69] , and the approach approximates the effect of strain-related softening 70,71 .
The buoyancy number is defined as the ratio between compositional density anomaly and density anomaly caused by expansion: where Δρ is compositional density anomaly compared with background mantle. Non-dimensional activation energy E is set to 9.21 to produce temperatureinduced viscosity variation of 10 4 between the temperature at 1 and 0 and mobilelid convection with relatively strong lithosphere 72 . We choose η r (r) = 1/30 for the upper mantle except for the lithosphere (depth < 100 km) where η r (r) = 1 and the viscosity prefactor in lower mantle is fixed at 1. Such a viscosity profile prefers a long-wavelength planform of mantle convection 72 .

Key variables and implications
Internal heat generation rate H. We give a non-dimensional internal heat generation rate (H = 100, corresponding to 9.22 × 10 −12 W/kg) which lead to a~60 internal heating 43,73,74 . We use Cases 5 and 6 to investigate the effect of varying H on melt generation and continental break-up. Case 5 has H = 50 (or 4.61 × 10 −12 W/kg) and Case 6 has H = 150 (or 13.83 × 10 −12 W/kg) (Supplementary Table 2). The results demonstrate that the internal heating rate changes the mean temperature of the mantle ( Supplementary Fig. 3a), which leads to varying melt amounts and depths for a mantle plume ( Supplementary Fig. 3b). In Case 5, no melt is produced since the mantle temperature remains lower than the solidus temperature. On the contrary, Case 6 generates excessive melts which leads to a disordered break-up process of Pangea ( Supplementary Fig. 3c).
Large low shear velocity provinces (LLSVPs). We defined the two dense chemical piles above the core-mantle boundary (CMB) to follow the same shape of the LLSVPs defined by the SMEAN share wave velocity anomaly (−1% contour) model at 2800 km depth 75,76 . Our simplified flat-top LLSVPs have a thickness of 200 km for the standard case (Case 1) as well as for cases 2-6 (Supplementary Table 2), as by previous work 73 . This thickness is the lower bound of estimated LLSVP thickness [77][78][79] . We also test the effect of changing LLSVPs thickness (e.g., 400 km in Case 7, and 600 km in Case 8; Supplementary Table 2) on the modeling results. Buoyancy number (parameter B) is set at 0.5 which corresponds to a 125 kg/m³ density contrast between the LLSVPs and the surrounding mantle. The compositional pre-factor of η c (C 3 ) = 1 gives the LLSVPs the same viscosity as the surrounding mantle. Since it has been widely recognized that the LLSVPs are relatively stable since 200 Ma 80 , we set our model with stable LLSVPs by setting the horizontal velocity of LLSVPs to zero.
Our calculations show that thicker LLSVPs accelerate the growth of the thermal instabilities originating from the edges of the LLSVPs. In Case 7, both the CAMP and Karoo plumes erupt at ca. 199 Ma ( Supplementary Fig. 4a), but the eruption of the Parana plume forms~20 Myr too early (Supplementary Table 2). In Case 8, these plumes all erupt at ca. 199 Ma (Supplementary Fig. 4b). We interpret that a thickness range of 200-400 km for LLSVPs is appropriate for our free-slip boundary model to reproduce the plume events in Central and South Atlantic. When the LLSVP thickness increases to beyond 600 km, more parameterization for the slope of the LLSVPs are required to reduce the heat transport to the deep plume conduits and obtain the desired timing of plume eruptions.
Rayleigh number. Rayleigh number also influences the convection vigor 81,82 and hence the plume eruption times. The thermal Rayleigh number Ra in our model Cases 1-8 is fixed to 1.5 × 10 8 (i.e., Ra = 1.37 × 10 7 if the mantle thickness is used to define Ra) with corresponding reference viscosity of 5 × 10 21 . In Cases 9 and 10 the Rayleigh numbers are three times and 1/3 of that for Cases 1-8, respectively. These two cases produce the eruption of the Parana plume over 30 Myrs earlier and 60 Myrs later than that in Case 1 (Supplementary Table 2), respectively.
Subduction zones. Subduction zones play the most important role in controlling the plume locations in our free-slip boundary model. A 200-km-wide weak zone at the supercontinent's boundary is prescribed to make it possible for the surrounding oceanic lithosphere to subduct 17,23 . The thickness of the weak zone is 150 km, the same as the thickness of the oceanic lithosphere. The viscosity of the weak zone is reduced by two orders of magnitude relative to the oceanic lithosphere 13 . The shapes of the reconstructed subduction zones vary with time following the plate reconstruction models used [51][52][53] .
It appears that subduction zones and LLSVPs determine the timing and location of plume eruptions 83 . The modeled CAMP and Parana plumes are both generated near the curvature of subduction zones in western Pangea (Fig. 2a-b, e-f). Though ca. 2000 km away from the actual location, the modeled Karoo plume is also generated by a curved subduction zone along the tip of South America and Antarctic [51][52][53] . A better parameterization on the subduction zone's shape and history is required to generate a more realistic generation of the Karoo plume. Free-slip model boundary conditions. The non-dimensional radii for the top and bottom boundaries are 1 and 0.55, respectively. Free-slip and isothermal boundary conditions are applied at the surface and the CMB in all calculations. Therefore, the supercontinent is free to move and deform, and to interact with mantle flows. The mantle is divided into 12 caps and each cap is further divided into 96 elements in three directions. The top and CMB thermal boundaries are refined to achieve a resolution of 20 km.
Due to the free-slip boundary condition and a lack of refined parameterization for the Farallon slab 84 , the plume beneath North Atlantic (corresponding to the Iceland plume) is not generated in an expected timing for the opening of the North Atlantic. We hypothesize that the opening of the North Atlantic should also be related to the trigger of plumes (the NAIP) and the guidance of orogens (the Caledonian and Alleganian-Variscan-Hercynian orogens 19,[85][86][87][88][89] . Classification and parameterization of cratons, orogens, and oceanic lithosphere. In our model, we use the IGCP 648 database for cratons and orogens >0.7 Ga 4 (Fig. 1). Orogens of 2.2-1.6 Ga are those related to the assembly of the supercontinent Nuna. Orogens of 1.6-0.7 Ga are generally those related to the assembly of the Rodinia supercontinent. Orogens of 0.7-0 Ga ages are those related to the assembly of Gondwana and Pangea, and those formed after Pangea breakup. Supplementary Fig. 1 represents simplified continental crustal ages for our modeling purpose. For our standard model Case 1, we give Archean cratons and Proterozoic orogens of >1.6 Ga stronger physical properties featuring larger thickness (200 km 90 ), and orogens of <0.7 Ga ages are defined as weak zones (young orogens) which have a thinner continental lithosphere (100 km thick) than the cratons (Supplementary Fig. 1a). In Case 2, we alternatively assigned orogens <1.6 Ga into two groups (1.6-0.7 Ga and <0.7 Ga; Supplementary Fig. 1b) with varying physical properties (140 and 100 km thick, respectively).
Compositional prefactor η c (C i ), buoyancy ratio B and yield stress σ y are used to describe the supercontinent. Continental cratons are much more viscous than the oceanic lithosphere 91 . However, too high a pre-factor may introduce numerical convergence difficulties 92 . We set the η c (C 1 ) (i.e., cratons) to 100, which is consistent with previous works 17,23,73,93 . The compositional prefactor of orogens (i.e., η c (C 2 )) is 10 whereas the compositional prefactor for orogens of 1.6-0.7 Ga (i.e., η c ðC 0 2 Þ) is 30 (Case 4). Parameter B in our model is all set at −0.2, which corresponds to a −50 kg/m 3 density contrast between the continental lithosphere and oceanic lithosphere. This number is the lower bound for the chemical buoyancy 26 . In our model, we neglect the localized yielding in the oceanic lithosphere in order to focus on the stress distribution within the continental lithosphere and get a faster computational convergence. The yield stress of cratons is set at 100 MPa which is consistent with laboratory experiments 8 and previous numerical works 17,23,94 . According to previous work 95 , the yield stress of orogens is at least a factor of 10 less than that of cratons, so we set yield stress for orogens of <0.7 Ga at 10 MPa. Yield stress for orogens of 1.6-0.7 Ga in Case 4 is set at 30 MPa.
Partial melting and melt extraction modeling. Partial melting of mantle plumes and the asthenosphere is implemented in a simplified manner 48-50 and we calculate the melt fraction for each time step. Melting in the geodynamic model depends on pressure and temperature, expressed as F = F (T, P) 96,97 . It is calculated based on the parameterization for batch melting of anhydrous peridotite 96,97 . This work uses a simple equations which represents the parameterized relationship between F, T and P 94 , allowing a fast and relatively accurate estimation of F: Here the pressure dependence is hidden in T solidus and T liquidus . T solidus and T liquidus are solidus and liquidus temperatures (in Kelvin), respectively. When T < T solidus , F = 0; when T > T liquidus , F = 1. Our model also considers the effect of latent heating when calculating melt fraction, and it is treated as a heat sink in the energy Eq. (3). Temperature T is calculate by adding a pre-factor (1 þ L C p ΔT 0 ) 98,99 before the first term of the left-hand side in Eq. (3), where L is the latent heat (assume to be 640 kJ/kg 43,100 ) and ΔT′ is the difference between solidus and liquidus temperatures. This implementation is benchmarked to a previous melt model 43 . Since adiabatic heating is ignored in our model under Boussinesq approximation, we need to consider the effect of adiabatic temperature when calculating melt fraction. Here we use adiabatic temperature gradient of 0.4 K km −1 (Supplementary Table 1).
Initial conditions. In order to achieve a reasonable initial mantle structure beneath the supercontinent, we perform a three-step calculation based on previous models 17, 23 as well as geological constraints. For the first step, we first set our model using an initial temperature profile obtained from the statistical steady-state thermal field of a pre-calculation ( Supplementary Fig. 5a-b). Secondly, we use the steady-state temperature field profile from the purely thermal convection model as the initial temperature condition, and introduce reconstructed subductions as weak zones. We then use the reconstructed Pangea supercontinent at 200 Ma with orogens and a pre-existing Pacific LLSVP, and run the model for up to 150 Myr to generate a global 3D mantle structure and an appropriate upper-mantle thermal status under the supercontinent. This will avoid the influence of the initial 3D return flow structures. In the third step, the Africa LLSVP is now introduced to promote the plumes under Pangea. When the first plume reaches the continental lithosphere ( Supplementary Fig. 5c), our model clock is set at 200 Ma, and the continental yield stress is applied to the supercontinent. In the second step, we examine the effects of the subduction history on the plume locations using three alternatively models with subduction starting from 300, 250, and 200 Ma, respectively. For the initial condition tests with subduction starting from 300 Ma, we run the subduction from 300 Ma for one transit time (a time period for a particle to migrate from the surface to the CMB with the average surface horizontal velocity of 5 cm/yr in this study) and update the subduction history from 300 to 200 Ma for every 1 Myr. The subduction migration is achieved by manually relocating the low viscosity weak zone for every 1 Myr. Similar runs are applied to test model with subduction starting from 250 Ma and 200 Ma. The produced plume locations for CAMP and Parana show little difference between the three test models, although the locations of other plumes vary between models due to their close relationships with subductions in the Tethys oceans and along the boundary of Siberia.

Data availability
All the data that are necessary (including orogens, pre-process Matlab script, and our source code) for evaluating the findings of this study are available at https://data. mendeley.com/datasets/ky4j7h7z8z/3/.

Code availability
The version of the CitcomS code used to simulate mantle convection can be accessed at https://data.mendeley.com/datasets/ky4j7h7z8z/3/ . The code used to make the figures can be accessed at www.generic-mapping-tools.org/ and www.paraview.org.