LIP formation and protracted lower mantle upwelling induced by rifting and delamination

Large Igneous Provinces (LIPs) are commonly attributed to mantle plumes, hot upwellings from the deep lower mantle, apparently unrelated to plate motions. However, LIPs often form in association with rifting and breakup. Using numerical modelling, we introduce a novel idea that explains plume-like mantle upwelling by plate tectonic processes. Our model indicates that rifting-induced delamination of orogenic lithosphere can perturb the thermochemical mantle stratification and induce lower mantle upwelling which causes syn-rift LIP formation followed by protracted and enhanced mid ocean ridge basalt (MORB) generation. Our model provides an explanation for the geographical correlation between the Caledonian suture, the North Atlantic Igneous Province (NAIP) and present-day Icelandic magmatism.

Continental breakup between North America and Europe in the late Palaeocene was accompanied by the formation of the NAIP, one of the classic LIPs. More than 1.3 × 10 6 km 2 of the present day North Atlantic region were covered by up to 1 × 10 7 km 3 of flood basalts 1,2 . Magmatism occurred in pulses between 62 and 54 Ma with a major outburst within the continental rift system that eventually became the North Atlantic Ocean 2,3 . Following continental breakup, the North Atlantic remained magmatically hyperactive, and thick (possibly >30 km) igneous crust was produced along the Greenland-Iceland-Faeroe Ridge (GIFR) 1 . Present-day Iceland similarly represents a mid-ocean ridge with anomalous topography, melt productivity and composition 2,4-6 .
The concept of mantle plumes, active upwellings with distinct chemical characteristics sired in the deep lower mantle 7 explains many first-order properties of the NAIP as well as other LIPs 3,8 . High melt productivity, trace element and isotopic compositions of melts reflect elevated mantle potential temperatures (T p ) and fertile melt sources, possibly with components of recycled crust 5 and primordial mantle 6,8 . Finally, low velocity anomalies in seismic tomography reaching into the lower mantle have been interpreted as plume conduits 9 . These observations are widely seen as evidence for plumes originating from the deep lower mantle independently of plate motion. However, most LIPs formed in the context of continental rifting and breakup 10 and are thus causally related to plate tectonic processes. The NAIP developed just where a Caledonian suture was dissected by Palaeogene rifting and breakup 11,12 , and the unlikeliness of a plume coinciding with exactly this intersection has been invoked as an argument against a deep origin 13 . Alternative 'top-down' explanations for LIP generation include edge-driven convection 14 , recycling of subducted oceanic crust 15 and delamination of lithospheric mantle 16 . However, these non-plume models have difficulties explaining the particular chemical signature of melts and the long-lived high melt productivity that prevails today in Iceland 55 Myr after initial LIP formation 17 .
We hypothesise that orogenic lithosphere with a dense crustal root in eclogite facies 18 may delaminate upon rifting 19 , rapidly sink through the upper mantle, penetrate the lower-upper mantle boundary (LUMB), and perturb the metastable stratification across this thermochemical boundary layer [20][21][22] . This would induce plume-like upwelling from the hot lower mantle and thus produce voluminous magmatism.
To explore this hypothesis, we use 2D-thermo-mechanical modelling. Our model accounts for pressure-and temperature dependent viscous, elastic and plastic flow in the lithosphere-mantle system 23,24 , and for decompression mantle melting. Relevant physical properties associated with mineralogical changes are tracked using tabulated thermodynamic, free energy minimised models of density, heat capacity and entropy (methods). The modelling domain is 2000 × 2000 km (Fig. 1, Supplementary Figure S1) and encompasses the upper and lower mantle. Initially, a 600 km wide orogen comprising thickened continental crust with 20 km of mafic lower crust is placed in the centre of the model. A 200 °C T p step is assumed across the LUMB, and a 40 km thick basaltic layer Colour striping is used to discriminate between crust (red/tan), mantle lithosphere (light/dark green), upper mantle (light blue/tea green) and lower mantle (cream/pink), and to show strain since model onset. Purple indicates mafic material, orange indicates regions were melts are currently forming, and grey tones with white text show total degree of melt depletion of mantle material that melted at earlier stages of the model run. Thin black lines are isotherms, and bold black lines are minimum and maximum temperature as functions of depth. For comparison, grey curves show theoretical adiabats for T p 's of 1300 and 1500 °C, respectively. (g-l) close ups of central crustal region (dashed boxes in a-f). Within 2 Myr, the delaminated lithospheric bodies sink through the upper mantle, reach the MTZ and penetrate the LUMB. This perturbs the metastable thermomechanical boundary layer and hot, buoyant lower mantle material rises to the rift zone causing voluminous melting (Figs 1c,d,i,j and 2f,g). The time between delamination and the onset of decompression melting of lower mantle upwellings is 6 Myr, and within 0.3 Myr, melt productivity (here defined as the volume of melt that is generated in a given time step normalised by extension rate; see also methods) reaches 280 km, 40 times average mid ocean ridge productivity 28 (Fig. 3). The high productivity is a consequence of both a large upwards flux of lower mantle material and its high T p .
Lower mantle upwelling is caused by the delaminating material that sinks into the lower mantle and shears the MORB graveyard-layer so that it locally thickens towards regions of downwelling and thins further away. Where the MORB layer is thinned, Rayleigh-Taylor instabilities of buoyant lower mantle material form and start penetrating into the upper mantle. The MORB graveyard-layer, that hitherto had suppressed convective mixing between upper and lower mantle 29 for 300 Myr, is thereafter breached as the lower mantle convective instabilities continue to grow, and as the delaminated material continues to sink into the lower mantle while dragging a wake of cold upper mantle material along ( Fig. 1c-e). Plume-like lower mantle upwelling continues to increase melt productivity for the 20 Myr that follow the arrival of the first upwelling and melting of lower mantle material, and average igneous crustal thickness reaches ~35 km, with peaks 5-10 times higher (Fig. 3). The high productivity is, apart from a higher lower mantle temperature and a fast upwelling rate, also a consequence of the thinned lithosphere (due to rifting and delamination), allowing for a relatively shallow upper limit of the melt of the melt-forming region of ~35 km, close to the base of the continental crust (Fig. 1j,k). The transient variability of melt productivity is a consequence of asthenospheric small-scale convention that modulates the upwelling flow from the lower mantle and pulsation of lower mantle upwelling. At 330-340 Myr, melt productivity is relatively low (5-10 km), because lower mantle upwelling slows down at this time and is not entrained into the rift zone. From ca. 342 Myr, a new pulse of lower mantle upwelling arrives into the rift system and increases melt productivity to 20-40 km for the next 10 Myr (Figs 1e,k and 3). In the following, pulsing and ascending hot lower mantle material continues to flow towards the thinned lithosphere, but gradually the sublithospheric, originally upper mantle, material becomes replaced with hot lower mantle material to depths greater than the maximum depth of melting at ~160 km ( Fig. 1f,l). Consequently, decompression melting of a hot lower mantle source continues, but flow within the melt zone is, in turn, predominantly driven by the passive flow associated with plate separation, rather than buoyantly enhanced active upwelling (Fig. 1d,j). The elevated and stable T p in the melt-generating zone leads to an average crust productivity of 16 km at times later than 50 Myr after the onset of extension.
Our model implies that plume-like upwelling of lower mantle material can be a direct cause of plate tectonic processes. It provides an explanation of why LIPs develop just where rifting intersects orogenic lithosphere and sutures. We propose that the pre-LIP magmatism in the NAIP at ca. 61-55 Ma is a consequence of rifting and delamination, while the intense activity during and following LIP formation at ~55 Ma was caused by actively rising lower mantle following a breach of the LUMB. This is consistent with models that point towards elevated temperatures and enriched or primordial melt source 5,6,8 . Also, low seismic velocities 9 , a high geoid and elevated topography 30 would be consistent with our model. Our model predicts magmatic pulsation due to chaotic interaction between active upwelling and upper mantle convection. We suggest that such upper mantle modulation could explain observations that have been attributed to plume pulsation 4,31 .
The initially elevated lower mantle T p , assumed in our model, is a possible consequence of extended periods of stagnation where the lower mantle would gradually heat up due internal radiogenic heating and/or conductive heat flow from the core 32 . Different estimates of radiogenic heat production would correspond to heating rates up to 160 to K per Gyr, ignoring heat loss and heating from the core (see supplementary note on radiogenic heating rate).
Large low shear velocity provinces at the core-mantle boundary have been suggested to correlate with historical positions of LIP formation, and it has been hypothesised that they control the generation of plumes and LIPs 33 . Another explanation involves the effect of thermal insulation beneath supercontinents 34 . Our model of a top-down control on lower mantle upwelling does not preclude the presence of either these mechanisms. Instead, it provides a triggering mechanism for rapid tapping and destabilisation of a regionally heated lower mantle reservoir. Such destabilisation would likely manifest as plume-like upwelling of various magnitude or as a regional mantle overturn, depending on the amount of trapped heat and its distribution. The latter would also determine the duration of upwelling, as the process of delamination-induced flow from the lower mantle would be limited by the amount of hot lower mantle available. The generic model presented here assumes that the whole lower mantle is uniformly hot and upwelling therefore continues until the end of the simulation at 100 Myr after the onset of extension (Fig. 3). The development of an improved model with a melt productivity evolution more consistent to that of the NAIP, or other LIPs, would require knowledge on the specific mechanisms that heat the lower mantle and its thermal state prior to the onset of extension. The model presented here requires a process that maintains thermal stratification with an elevated T p of the lower mantle relative to the upper mantle for extended periods of time. It has been shown in earlier works 21,22 that, depending on the somewhat controversial 35 magnitude of the negative Clapeyron slope of the Ringwoodite-out reaction near the LUMB, mantle convection is possibly layered due to the convection-inhibiting effect of this phase transition. For the thermodynamic database employed in the present study 36 , the Clapeyron slope is half that of several recent estimates, and the associated convection-inhibiting effect is therefore relatively small (supplementary discussion). However, the possible presence of a MORB layer at the LUMB 25,26 , as assumed here, additionally to sustain mantle stratification 29 (compare with the model presented in Supplementary Figure S7 where no such MORB is assumed).
An additional vital component of the model presented here is the presence of a dense orogenic root that facilitates delamination during rifting. To illustrate this point, a similar simulation without an orogenic root is presented in Supplementary Figure 3. In that case, rifting does not trigger delamination, and no upwelling from the lower mantle is induced. Instead, mantle melting is predominantly driven by passive upwelling of the upper mantle, and melt production only reaches 4-5 km.
While numerous studies consider lower mantle upwelling as the cause for rifting and continental breakup, the results presented here suggest the reverse relationship, in which rifting and continental breakup trigger lower mantle upwelling. This means that plume activity does not need to be complementary to plate tectonic processes to explain anomalous magmatism. If a certain lower mantle region is not cooled by descending slabs and has not been tapped by mantle upwelling for sufficiently long time, a regional thermal anomaly that is preserved/protected by the convection-suppressing nature of the LUMB can develop. A regional mantle overturn can then be triggered by plate tectonic processes like the above-described mechanism. Our generalised model can explain the large rates of melt production associated with the formation of the NAIP and why excess magmatism still occurs in present-day Iceland. An improved model would require additional constraints on structure of the mantle beneath the North Atlantic region (e.g. ref. 37 ), should be tested against geophysical and geochemical observations and would possibly also need to account for 3D effects involved in the process of delamination-induced mantle upwelling. Several other LIPs than the NAIP formed concurrently with continental breakup at locations of earlier suturing 10 . The model presented here provides a potential cause for this apparently common correlation. The implication would then be that mantle upwelling is not necessarily a process that operates outside of the context of plate tectonics, but is instead intimately related to the latter. Future work will address this question.

Methods
Overview. The employed numerical method for modelling visco-elastic-plastic flow is generally similar to that described by Gerya and Yuen 38 , but differs by employing a multigrid-based approach 24,39,40 which allows for high resolution simulation in both space and time. The presented method also differs from the abovementioned references in that we employ a quantitative petrological model for the calculation of density and latent heat effects 41,42 .
The simulations of the present paper have a grid resolution of ~2 km for the entire 2000 km × 2000 km modelling domain and are run for 400 Myr.
The method employs a combined particle-in-cell approach 38 for the coupled equations for the conservation of energy, momentum and mass in 2 dimensions, k is thermal conductivity, T is temperature, C P is isobaric heat capacity, ρ is density and DT Dt is material time derivative of temperature. The source terms, H r + H a + H S + H L represent the total rate of enthalpy change due to radiogenic, adiabatic, shear, and solid state phase change heating, respectively. DF Dt is the material time derivative of the total melt depletion of the mantle source, and L is enthalpy change due to melting. σ′ ij is the deviatoric stress tensor defined as σ σ δ ′ = + P ij ij ij , where σ ij is the stress tensor, and 22 is pressure. v j is the velocity field vector that is constitutively related to stress through a viscous-elastic-plastic strain rate tensor (see the section on the rheological model in the below). Petrologic and thermodynamic model. To account for solid state phase changes at given pressure and temperature conditions, phase equilibria are calculated using the Gibb's free energy minimisation software, Perple_X 42 together with a thermodynamic database 36 for common mantle mineral assemblages. We employ 3 general lithologies with different compositions in the Na 2 O-CaO-FeO-MgO-Al 2 O 3 -SiO 2 (NCFMAS) system. These include (1) pyrolite 43 for the ambient upper and lower mantle (with the exception of the mantle lithosphere and pyrolite that is subject to melt depletion during model evolution -see below equations 9 and 10, (2) mid-ocean ridge basalt, MORB 44,45 for orogenic lower crust and mafic material in the transition zone 25 and (3) Continental crust 44,46 . Compositions are given in Supplementary Table S2.
All endmembers and solution models are employed from the Stixrude and Lithgow-Bertelloni-database 36 are employed, assuming the 'Gt' solution model for garnet.
To reduce computation time, the phase-equilibria are pre-calculated in Perple_X, and lookup tables for density (ρ), isobaric heat capacity (C P ), and entropy of formation (S) are calculated on a grid in P-T-space with a resolution of 75 MPa and 10 K, respectively. In our numerical model, the tables are pre-loaded, and bilinear interpolation in P-T-space allows for efficient look-up of the relevant physical properties. To account for solid state phase changes and melting, we employ a method that can be regarded as an extension of any given thermo-mechanical method that otherwise does not account for such phase changes.
Density and isobaric heat capacity are calculated at each Lagrangian particle that has one of the three abovementioned lithologies and pressure, P 0 , and temperature, T 0 . Each time step is broken into two sub-steps where effects of solid state phase changes are accounted for. A third step accounts for effects of mantle melting.
Isentropic advection. The coupled momentum and continuity equations (2 and 3) are solved, and particles are then moved according to the calculated velocity field. This step involves a pressure change at each marker: dP = P 1 − P 0 , where P 1 is the new pressure. The associated temperature change due to adiabatic heating is found by assuming that entropy is conserved during the pressure change: S(P 0 , T 0 ) = S(P 1 , T 1 ). This equation is solved for T 1 by Newton-Raphson iterations that are evaluated by calculating S(P 1 , T 1 ) from repeated table interpolations with varying temperature. If Newton-Raphson iterations do not converge (e.g. in the case of steep entropy gradients with respect to temperature), a solution is found by using bisection starting at an interval that corresponds to the entire temperature range of the entropy table.
Isobaric thermal diffusion. This step involves the solution of the energy equation (1) using a procedure that has been employed in earlier works 24,39 , but differs by assuming that initial temperature is T 1 .
First, the energy equation (1) is solved for a time step dt, by assuming an initial temperature field T 1 and that source terms for adiabatic and latent heating (due to both melting and solid state phase changes) are 0. The resulting time integrated temperature, T 2 , accounts for thermal diffusion, radiogenic heating and shear heating. The updated temperature T 2 also implicitly accounts for adiabatic heating as this was calculated in the above step (1). To also account for latent heat effects due to solid state phase changes, a correction is applied by calculating the entropy change dS corresponding to the temperature difference between T 2 and T 1 : This entropy change does not necessarily correspond to a change of sensible heat, only, but could also be associated with latent heating. Therefore the entropy table is employed to solve S(P 1 , T 1 ) + dS = S(P 1 , T 3 ) for T 3 . Like in the advection step, Newton-Raphson iterations or bisection are employed together with repeated entropy table interpolations.
Isobaric mantle melting. The thermodynamic database employed for the calculations in the above does not include a melt phase. Instead melt productivity is calculated following a procedure described elsewhere 39,47 . It is summarized here: We employ a parameterization of mantle melting, where, the temperature T m of a mantle peridotite source in equilibrium with melt is given as a function of pressure P 1 and degree of melt depletion F. Here, we specifically employ the T m -parameterisation for batch melting given by Katz et al. 48 , and assume a water content of 0 %, and an initial modal clinopyroxene concentration of 17 % for unmolten solid mantle. The change in the degree of melting is calculated by assuming an isobaric process where an excess temperature, ΔT = T 3 − T m , due to the above steps 1-2, provides an entropy excess = Δ dS steps that is assumed to balance the combined entropy increase due to both changing the degree of melting and temperature, respectively, . This equation can be approximately solved 39,49 for the increase of the degree of melt depletion, ΔF, on each marker as: Tm F P Latent heating is given by L = T 3 ΔS melt , where ΔS melt , the entropy change due to melting, is assumed to be constant and independent of the degree of melting 49 . Melting is assumed to be irreversible at time scales greater than individual time steps, and the degree of melting therefore does not decrease if ΔT < 0. This assumption is based on yet another simplifying assumption that melts are rapidly removed from the source 50 , even though the melting parameterization is based on batch melting experiments. The calculated melt productivity can therefore only be regarded as a rough approximation.
Following Phipps Morgan & Morgan 49 , the temperature change due to latent heat consumption during melting is updated in accord with the change of T m due to the melt increment ΔF, and the final temperature corrected for both melting and solid state phase changes is therefore (if ΔF > 0): If ΔF > 0, the final temperature is simply T 3 .

Rheological model. Strain rate is related to the velocity field by
and is assumed to be a mutual contribution from brittle, elastic and viscous deformation, where both diffusion and dislocation creep 24,38 are permitted: The plastic multiplier, χ is nonzero only if the second invariant of the deviatoric stress tensor, σ σ = II kk , reaches the pressure-dependent Mohr-Coulomb failure limit: Where C is cohesion. Viscosity depends on composition, pressure and temperature: . The continental crust follows a 'plagioclase' power-law 51 . The mafic lower crust initially present in the orogen at the centre of the model and within the mantle transition zone is assumed to follow an eclogite power-law 52 at pressures less than 17.5 GPa. At higher pressures, a garnetite 53 flow law is assumed.
For the pyrolite upper mantle, flow laws 54 for both disclocation and diffusion creep are employed. For the pyrolite lower mantle, a simple linear flow law constrained by slab sinking rates 55 is employed (assuming that . A density threshold of 4300 kg m 3 is used to discriminate between upper and lower mantle rheology for pyrolite. The procedure for calculating density is described in the above section. Mechanical parameters are given in Table S1. The presence of interstitial melt reduces the effective viscosity of mantle rocks 56 . However, the change the composition of the mantle source upon melt extraction tends to increase viscosity 56,57 . This strengthening effect is assumed to be of a larger magnitude than that of melt weakening, as we assume that melt is continuously extracted from the mantle source 47,58 . Following Phipps Morgan 57 , we assume two components of strengthening due to melt depletion (at a given pressure and temperature): (1) A 30-fold viscosity increase during the first 5% of melting where most water extraction is expected to occur, and (2) A third of an order of magnitude increase for each 15% increment of melt extraction due to the generally higher viscosity of refractory mantle with higher Mg-content. These effects are approximated by assuming that the depletion strengthening factor, d(F), is: (20 )) exp(8 02 ) (9) This is consistent with the parameterisation of Phipps Morgan 57 , and approximately consistent with the viscosity-depletion increase employed by Afonso et al. 59 even though the latter only assume the effect of dehydration of viscosity.
Melt depletion is also known to decrease the density of the mantle residue 57,60 . This effect is approximated by correcting the density, ρ (Gibbs) , calculated from Gibb's free energy minimisation: Following the Boussinesq approximation 51 , density changes due to changing pressure and temperature conditions are assumed to affect only the momentum equation (2), but are ignored in the continuity equation (3). Similarly, mass exchange between solid and melt phase and the relative advection between the two phases (compaction) is not accounted for in the continuity equation 47 .
Thermal conductivity and radiogenic heat productivity. For continental crustal material, constant values of thermal conductivity are assumed (see Table S1). For pyrolite and for the MORB material assumed in the orogenic lower crust and in the mantle transition zone, a pressure and temperature-dependent parameterisation 61 is assumed: Here q is an exponent assumed to be 0.406 for both pyrolite and MORB. Similarly, for both compositions, a constant pressure gradient, = .
− a G Pa 0 04 1 is assumed 62 . The standard conductivity value, k 298 , is however assumed to differ between pyrolite and MORB compositions, such that for pyrolite at all conditions and MORB at pressures less than 17.GPa, it is = .
⋅ k 4 08 J m K 298 . At higher pressure, MORB is assumed behave similarly to a garnet rich assemblage that has a conductivity 63 Radiogenic heat production is assumed to occur in the crust at a rate of = . ⋅ − H 1 0 10 r W m 6 3 . For MORB material, heat production is assumed to be negligible due to both small amounts and low productivities. Due to the great volumes of pyrolite, a nonzero heat production of ρ = . ⋅ ⋅ − H 1 0 10 r W kg 12 is assumed. This heat production corresponds to that reported for an average MORB source in Table 9 of ref. 32 and references therein.
Initial model state and boundary conditions. The model domain is a 2000 km by 2000 km square box (Fig. S1). A free surface is approximated using a weak 'sticky air' layer 64 in the upper 80 km of the domain. The sticky air is visco-plastic with a viscosity of 10 20 Pa s and a finite strength of 0.1 MPa (similar to the numerical precision of the modelling method). This renders the surface effectively stress-free even during rapid surface movements and dynamically limits the effective viscosity only when needed 39 . Below the sticky air is a layer of continental crust with a thickness of 40-50 km. The crustal thickness is 50 km in a 600 km wide region in the lateral centre of the model domain. Elsewhere, the crust is 40 km. Below, the thickened crust is a 20 km layer of MORB-like material. For reference, a model without an initial lower crustal MORB layer is calculated (Extended data Fig. 3). The rest of the model domain consists of mantle with a pyrolite composition with two exceptions: (1) A lithosphere/asthenosphere boundary is assumed at a depth of 150 km. Mantle within the lithosphere is possibly more depleted, stronger and more buoyant than the asthenospheric mantle below 57,59,60 . To approximate this, it is assumed that lithospheric mantle is depleted by F = 30%. This has the effect that the mantle lithosphere has ~1.8% less density than pyrolite at similar conditions (equation 10), and ~330 times higher viscosity (equation 9). Consistently with previously published models of long-term stability of continental lithosphere [65][66][67] , the mantle lithosphere remains stable against small-scale convective erosion during the first 300 Myr years of model evolution prior to rifting and retains an approximately constant thermal structure (Fig. 1). (2) At the mantle transition zone, a layer of MORB-material 25 is imposed. The thermodynamic model implies that this material is neutrally buoyant relative to pyrolite at this depth 26 . This is because MORB-compositions are garnet rich and therefore denser than pyrolite above the 660 km discontinuity (where pyrolite is predominantly rich in ringwoodite), but less dense below the discontinuity where pyrolite is a perovskite-rich assemblage. Further below (ca. 100 km, depending on temperature), the MORB composition is again denser than pyrolite, due to garnet to perovskite transformation. Rather than imposing the garnetite layer at a fixed depth, we allow for it to form at a level of neutral buoyancy by initializing this material as a number of circular bodies that are capable of freely moving due to their buoyancy. They are initially placed at a depth of 660 km with 200 km lateral distance between their centres and have a radius of 50 km (Fig. S1). Within the first Myr of model evolution, they sink to a level of neutral buoyancy and spread out to a single layer due their finite viscosity.
Initial thermal structure and thermal boundary conditions. The temperature structure of the sub-lithospheric mantle is assumed to be initially adiabatic and is calculated as the contour (isentrope) of the pyrolite entropy map that intersects the surface pressure (P = 1 Pa), at a temperature of 1325 °C. At depths greater than 660 km, temperatures are elevated relative to this adiabat by 200 °C (Fig. S1), assuming that mantle convection is regionally in a layered state 21,22 with a thermal boundary layer separating the upper and lower mantle (see also discussion in main paper). The 150 km lithosphere is assumed to be in a 2D conductive thermal steady state with the adiabatic mantle temperature as a lower boundary condition and 0 °C as an upper boundary. Vertical boundaries at the side of the model are everywhere isolating with zero lateral heat flow. The lower boundary is a constant temperature condition that equals the initial temperature at that depth.
Mechanical boundary conditions. During the first 300 Myr model time, the model has free slip at all 4 bounding sides and closed sides. This means that velocity perpendicular to a side is 0, and that traction parallel to the side is also zero. At 300 Myr, rifting of the lithosphere is kinematically forced by imposing plate separation, v sep , at a rate of v sep = 1 cm/yr. This is employed by imposing outwards perpendicular velocities at both left and right boundaries of v sep 1 2 at depths from 0 to 240 km. To obey the volume-conserving continuity equation (3) globally, an inwards perpendicular flow of 0.4 cm/yr is imposed at the same boundaries at depths from 240 to 540 km. The material flowing into the model domain during the rift stage is pyrolite with F = 0% and a temperature that corresponds to the initial adiabat.
Modelled melt productivity. The intensity of melt generation is calculated following an approach 68 where an effective igneous crustal thickness is found by calculating the total mass of melt (per meter perpendicular to the 2D model cross section) produced during a given time step by integrating the change of depletion, ΔF, in a given time step, Δt, over the whole modelling domain: , where ρ source is an average density of the melt source region. If this mass of melt would be extracted to form igneous crust with average density, ρ crust , its volume (i.e. area because the model is 2D) would be If the total amount of plate separation during a time step, would be accommodated by the formation of a single dike with a width of v sep ∆t, and if the total volume of crust produced would be emplaced within that dike as a single column, the effective igneous crustal thickness, h c , can be defined as the corresponding height of that column: This measure of productivity corresponds to the left y-axis of Fig. 3. Observationally constrained North Atlantic productivity estimates used in this study are not directly reported as magmatic crustal thickness values (e.g. in meters), but instead reflect the cross-sectional area of melt generated within certain (rift-perpendicular) profiles pr. time and are reported as e.g. km 3 /km/Ma 1 . To compare such data with modelling results, modelled productivities are also calculated as For the modelled productivities in Fig. 3, this rate corresponds to the right y-axis, and happens to be proportional to magmatic crustal thickness (equation 13), because the imposed boundary plate separation velocity, v sep , does not change at times where melts are produced.
The ratio of melt source density and melt density is assumed to be = The modelled melt productivity is shown in Fig. 3 where it is decomposed into melt originating from a source that was initially (at the model onset time) part of the upper or lower mantle, respectively. This is specifically archived by splitting the integral Δ ∬ F dxdz appearing in equations (13 and 14) into two terms that each either integrate over the domain comprised of upper mantle material or the domain comprised by lower mantle material, respectively. Fig. 3. The North Atlantic melt productivity curves shown in Fig. 3 Table S2 where productivities for Greenland-Faroe, Aegir/Kolbeinsey and Reykjanes ridge segments are given in units of volume pr. time. In the same table, characteristic lengths of each of these segments are given as 300, 500 and 1100 km, respectively. For comparison with model predictions and the other productivity curves, we normalise the volumetric productivities by these lengths to present productivities given in terms of volume pr. length pr. time. Furthermore, Mjelde & Faleide 4 subtracted a 'normal' oceanic melt productivity of 5.5 km from their volumetric productivity estimates by assuming the time-dependent full-spreading rates given in their Table S1. We have re-added this background production, as our modelled melt productivities are not only a result of excess mantle upwelling and delamination, but also encompass melt generated from passive upwelling due to plate separation. Fig. 3

Code Availability
Numerical modelling code is available upon request from the authors.