Melt-induced buoyancy may explain the elevated rift-rapid sag paradox during breakup of continental plates

The division of the earth’s surface into continents and oceans is a consequence of plate tectonics but a geological paradox exists at continent-ocean boundaries. Continental plate is thicker and lighter than oceanic plate, floating higher on the mantle asthenosphere, but it can rift apart by thinning and heating to form new oceans. In theory, continental plate subsides in proportion to the amount it is thinned and subsequently by the rate it cools down. However, seismic and borehole data from continental margins like the Atlantic show that the upper surface of many plates remains close to sea-level during rifting, inconsistent with its thickness, and subsides after breakup more rapidly than cooling predicts. Here we use numerical models to investigate the origin and nature of this puzzling behaviour with data from the Kwanza Basin, offshore Angola. We explore an idea where the continental plate is made increasingly buoyant during rifting by melt produced and trapped in the asthenosphere. Using finite element simulation, we demonstrate that partially molten asthenosphere combined with other mantle processes can counteract the subsidence effect of thinning plate, keeping it elevated by 2-3 km until breakup. Rapid subsidence occurs after breakup when melt is lost to the embryonic ocean ridge.


Information on the ages of rifting
The best estimate of the onset of rifting is 131 Ma based on the oldest syn-rift sediments which are Hauterivian in age 16 . The precursor Paraná volcanics are precisely constrained by 40 Ar/ 39 Ar dating to 135 Ma +/-1 My 31 and uniform layering within them shows that they were erupted before extension began. Information on the timing of the end of rifting (123 Ma) is reported in the main article and in Refs. 2,5,8,16,18,30 .

Seismic interpretation
The nature of the seismic units in the central South Atlantic region has previously been described by Ref. 5 . The syn-rift interval is characterized by tilted and fanning reflection geometries. Strong parallel reflections at the base of the syn-rift interval correlate with the top of basalt (the major Paraná eruptive event, the precursor to rifting in the central South Atlantic region 8,31 ). The top of the syn-rift interval is marked by an angular unconformity onto which flatlying sag strata progressively overstep (Fig. S3). The top of the sag is a flat horizon above which the salt has a clear but chaotic seismic expression resulting from salt tectonic movement. The top of the salt is marked by a strong seismic reflection and is overlain by Albian carbonates characterized by growth, fault and fold geometries. These structures were formed during gravity spreading as the salt tilted towards the ocean as a result of enhanced subsidence in that direction 18 .

Are anomalous subsidence patterns during breakup related to mantle plumes?
The aim of our work is to find a generic solution for the common problem of elevated breakup and subsequent collapse at rifted margins 6,11 . We acknowledge that the thermal and dynamic effect of a mantle plume can have a semi-regional effect [32][33][34] . However, we feel that this is not a universal explanation for the subsidence issue associated with McKenzie's model 1 during the breakup of continental plates. Nonetheless, one of us (D.G.Q., Ref. 5 ) has previously promoted this idea for a region of the central South Atlantic lying south of the current study area (Fig.S3), where syn-rift strata are particularly thin. It is therefore feasible that anomalous uplift caused by mantle plumes may enhance buoyancy caused by melt and phase changes in the mantle but probably not along every rift margin.

Information on melt produced during rifting
In the accompanying article, we explore the idea that melt is retained in the asthenosphere during rifting and lost rapidly after breakup. That melt exists in the mantle beneath modern continental and oceanic rifts is recognized from an increase in electrical conductivity, a reduction in seismic velocities and a decrease in gravity  . The geophysical evidence for melt is also supported by theoretical models, experiments and geochemical data [61][62][63][64][65][66][67][68][69][70][71][72][73] . Unfortunately, the amount of melt is difficult to constrain because many properties of the asthenosphere such as interconnectivity of melt-filled pores and permeability are poorly constrained with large uncertainties in physical characteristics, including the shape and alignment of melt inclusions 40,43,[74][75][76][77] and whether the peridotite is in textural equilibrium [78][79][80] .
As geophysical data and geodynamic models grow, two different hypotheses have emerged on the amount of melt present beneath modern rifts, with continental and oceanic rifts often treated similarly. 1) A common view is that dihedral angles of melt formed in olivine-rich peridotite are significantly less than 60 o 79,81 and the peridotite is texturally equilibriated 78 meaning that melt flows through the matrix at very low porosities. Consequently, the mantle will retain less than c.2% melt 35,62,70,[83][84][85][86] , migrating out of the source region faster than it is generated 67,87 . This view is supported by various geophysical interpretations of <1%-3% melt in regions of upwelling asthenosphere 43 , more commonly in oceanic domains 48,50,[54][55][56]59,88 .
2) Alternative experimental work and theoretical models indicate that lherzolite can actually retain 5%-c.20% melt 79,80,89,90 , consistent with other geophysical interpretations of c. 5-16%  melt in zones 30-100 km thick and 30->150 km wide in the upper mantle 44,46,47,91 , including  continental regions 52 . For example, Ref. 92 , using magnetotellurics, and Ref. 77 , using P-S wave velocities and anisotropy, have independently produced similar estimates of 6-15% melt over an interval 30-60 km thick in the upper mantle beneath the Dabbahu rift segment in Afar, Ethiopia, equivalent to a melt column of 3-10 km. Assuming these interpretations are valid, there are five potential reasons for high melt retention. a) Where peridotite is neither mono-mineralic nor equi-granular, melt does not become interconnected until porosities of somewhere between 5-20% 79,80,90 . b) Textural equilibrium is difficult to establish in flowing mantle 77 . c) Interfacial tension presents a significant barrier to compaction-and buoyancy-driven flow at porosities of less than c.20%, even where grain boundaries are wet 89,93 . d) Viscous resistance to flow means migration can be slower than generation 86,94 . e) Melt flowing between crystals or carried upwards in circulating asthenosphere is trapped by a viscous lid of mantle lithosphere which represents a permeability barrier 39,56,88,95,96 , ultimately resulting in crystal fractionation [96][97][98] .
One of the key points is that melting experiments on real peridotite show there is a drastic reduction in the interconnectivity of pores when minerals other than olivine are present 80 . As much as 21-29% basaltic melt can be retained in lherzolite containing less than c.80% olivine 90 up to the point of disaggregation 89 . A typical fertile lherzolite comprises 60% olivine, 20% orthopyroxene, 10% clinopyroxene and 10% garnet or spinel or plagioclase 99 , so theoretically all the aluminium-rich phases, all the clinopyroxene and some orthopyroxene have to melt before the pores become fully connected. Thus, it may be that melt is retained in asthenosphere composed of heterogeneous lherzolite but will tend to flow out when it is mostly olivine, provided interfacial tension is overcome 89 . However, in the continental domain a more significant permeability barrier exists at the top of the asthenosphere in the form of the continental plate 95,96,98 . Magma can only ascend through the overlying lithospheric mantle in discrete, fracture-assisted intrusions rather than by porous flow 9,43,74,84 . These considerations offer a way of aligning models of low versus high melt retention beneath rifted margins.
 Prior to breakup, when rifting reaches a peak, large amounts of melt have been generated and a significant proportion is temporarily retained in the upwelling asthenosphere, initially because porous flow is inhibited in multi-mineralic, non-equilibriated lherzolite and thereafter because the overlying lithosphere forms a closed lid. The zone of meltimpregnated asthenosphere will tend to drain downwards by capillary action, a wetting profile driven by grain-grain tension.
 After breakup, the mantle down to c.90 km below the inactive rifted margins is largely devoid of garnet-spinel-plagioclase and clinopyroxene, these minerals having already been lost to melt. The depleted asthenosphere consists predominantly of olivine-rich, texturallyequilibriated harzburgite-dunite with a dihedral angle significantly less than 60 o meaning most of the remaining melt is interconnected at low porosities (<1-3%) and migrates relatively quickly to the base of the lithosphere. Provided the melt does not solidify, it will then flow up-dip towards the edge of the continental plate to be incorporated in thick oceanic crust created at the new spreading axis. As the depleted asthenosphere cools, it is gradually accreted to the base of the continental plate to form new lithosphere together with melt which has solidified before reaching the oceanic domain.
With these processes in mind, we suggest it is possible that asthenosphere retains 5-16% melt during the latter stages of continental rifting prior to breakup 52,77,92 , more than below mid-ocean ridges where there is no mantle lid. In our melt models, we actually use average amounts at the lower end of this range (those prefixed M and F in Figs. 1 and 2 in the accompanying article, Figs.S1 and S2 and Movies S1-S3).
It is also worth noting that the presence of melt beneath an extending plate is likely to cause periodic intrusions of sills or underplates into the lithosphere 97, , an effect which may partly compensate for thinning of the plate, allowing certain rifted margins to become surprisingly wide (e.g. Fig.8 in ref. 5 ). Underplates are quite well documented beneath central South Atlantic rifted margins with interpreted thicknesses up to 6 km 14,25 .
Based on the arguments above, we test the effect of relatively high melt retention in the asthenosphere. The idea is that beneath continental rifts up to breakup i) significant porosity is required for melt to become interconnected in the asthenosphere 76,79,80,89,90 and ii) the base of the continental lithosphere acts as a lid 39,57,88,95,96,98 .
Once it becomes interconnected, melt migrates upwards through the asthenosphere by porous flow, controlled by compaction and buoyancy, viscosity, permeability and surface tensions 62,67,89,93,94,99 . The upward migrating melt will be impeded by continental mantle lithosphere and the zone of melt-impregnated asthenosphere will then increase in thickness downwards by capillary action, with wetting driven by grain-grain tension to give a smoothly decreasing melt profile 89,93 , as used in our models ( Fig. 1c in main article).
We choose to limit the average amount of melt in the asthenosphere to 8%, starting with 0% at 90 km depth up to a maximum of 16% at 20 km close to the point of breakup where the lithosphere is highly attenuated (stretching factor of >7.5, Figs. S4c and S4d). This is equivalent to a maximum column of 5.6 km melt. However, over most of the Kwanza Basin, where stretching factors are <5 (Fig. S4c), a column less than half this amount has developed in our models by the end of rifting. For comparison, more melt is involved in the construction of oceanic crust but here it ascends relatively quickly to a mid-ocean ridge, unimpeded by lithosphere. A total column of 7.5 km melt is needed for average oceanic crust 73,82,97 but as much as 12-24 km of melt may be involved during the formation of volcanic rifted margins 9,87,96 .
We do not argue for unusual mantle conditions of temperature, only that a significant proportion of melt cannot migrate out of the asthenosphere until after breakup. In our models, the asthenosphere wells up passively while the lithosphere extends, a general assumption in most models of continental rifts 13,83 . In reality, the presence of melt in the asthenosphere may enhance convection due to buoyancy and a possible reduction in viscosity 60,61,83,99,100 although Ref. 101 suggests that viscosity may initially increase as the asthenosphere becomes drier with melting. If active upwelling or a mantle plume are associated with rifting, then higher degrees of melt are expected due to replenishment of fertile peridotite 96 or higher temperatures 64,82 which will lead to more volcanism at breakup 87 .

Information on loss of melt after breakup
During rifting, significant columns of melt may be trapped beneath continental lithosphere until it breaks, after which it is likely to drain rapidly to the new volcanic ridge 87,94,102 . This may explain the release of huge volumes of magma over short periods of time forming thick crust on volcanic margins 82,96,97,102 as well as outer volcanic highs on non-volcanic margins 5,23,103 . In the North Atlantic Igneous Province, Ref. 104 suggests that a c.7 km column of melt per unit area was released to the crust in the 5 My up to breakup followed by 6-17 km soon after breakup, which we suggest was partly due to lateral migration of melt stored in the asthenosphere (cf. Ref. 105 ).
After breakup, the asthenosphere beneath the rifted margins ceases to flow and textural equilibrium is re-established. Melt production also stops and porosity decreases causing the rate of melt migration to decline. The density of the mantle will increase as melt is lost, causing the rifted margins to subside. Based on empirical observations of post-breakup volcanism and thicknesses of early post-rift sag 5 , we expect melt to be lost at an exponential rate with only small percentages left after 10 My, although the lithosphere continues to thicken and subside by cooling.

Basin reconstruction approach
TecMod is modeling software for lithosphere-scale basin tectonic, isostatic, stratigraphic, thermal modeling which accounts for the feedbacks between sedimentary, crustal, and mantle processes 106 . It is designed to produce automated 2D basin models using present day stratigraphy as input, from which the structural and thermal evolution of a rifted basin is automatically reconstructed. This is achieved through the coupling of a lithosphere scale forward model with an inverse algorithm for automated parameter optimization. This kind of algorithm allows the feedbacks and interactions between shallow sedimentary and deep lithosphere/asthenosphere processes to be investigated, which are difficult to predict without iterative modeling.
Given the number of physical processes involved and the range of parameters, the evolution of sedimentary basins cannot be uniquely reconstructed. The strategy employed in TecMod is scenario analysis performed on the (user-specified) geologically most likely end member cases, which yields a possible range in target parameters (mainly temperature, heatflow and related aspects though time and space). As demonstrated by Ref. 106 , the actual fitting process of crustal and mantle stretching factors, water depth, sediment fill, thermal conditioning and isostatic adjustment can be fully automated by casting the problem in terms of a constrained optimization problem 107,108 . The goal function to be minimized is the misfit between observed and modeled stratigraphic data where additional constraints can be imposed (e.g. finite rifting durations, maximum stretching factors, weighting between different events). The inversion consists of an iterative search for the optimal set of crustal and mantle stretching factors and palaeo-water depth values, which minimize the chosen goal function 109,110 . The forward run requires an initial guess for stretching factors and water depth which are successively improved by iteration with the inverse algorithm according to the scheme described by Ref. 106 .
The 2D model is locally based on pure shear kinematics and resolves simultaneously differential thinning, flexural isostasy, heat transfer, sedimentation, erosion and compaction on a Lagrangian finite-element mesh. The general modeling strategy is that the structural and thermal solutions are in equilibrium after every simulation time step. This is achieved by resolving simultaneously for deep lithosphere and shallow sedimentary basin processes. During every time step, the entire computational mesh is deformed with the mantle and crustal stretching factors, new sediment packages are deposited and become part of the computational domain. Flexural isostasy is ensured by further movement of the mesh in the vertical direction.
Heat transfer occurs by advection and diffusion: where t is time and vx and vz are the velocities in the horizontal x-and vertical z-direction, respectively. In the employed Lagrangian reference frame, the numerical grid is deformed with the pure shear velocities so that the advection terms in equation (1) are automatically accounted for. In the case of differential thinning, re-meshing of the mantle ensures geometric consistency between the crust/sediments and the mantle portions of the modeling domain. ρ, ci, T, ki and Ai are density, specific heat, temperature, thermal conductivity and radiogenic heat production rate, respectively. The index i illustrates that the property is dependent on the local rock type. The thermal conductivity of the sediments is an effective material property that results from geometric averaging of the local matrix conductivity ( km) and pore-fluid conductivity (kw) 111,112 : where φ is porosity. The temperature dependence of thermal conductivity is accounted for by assuming the following relations 113,114 : with c parameters defined in Table S1. Crustal and mantle densities evolve during basin formation as functions of thinning, melt production, magmatic intrusion, heat transfer and sediment and water loading as well as phase transitions 115,116 . The resulting stresses on the lithosphere are projected onto a plate with an assumed elastic thickness and are compensated by isostasy [117][118][119] . The equation describing the vertical deflection of the top of the crust, w, due to flexural isostasy is: where x, D, ρm, ρin, g, and S are the horizontal coordinate, the flexural rigidity, the average density of the lithospheric mantle, the average density of the basin infill (compacted sediments), the acceleration due to gravity and the kinematic deflection due to necking. The differential load, q, represents differences in load relative to the load in the equilibrium state. The flexural rigidity, D, is a cubic function of the effective elastic thickness (Te), which can be kept either constant or determined by the depth of a specified isotherm 120 or the thickness of the crust. The kinematic deflection S, i.e. the surface depression that results from stretching and not from isostatic compensation, is calculated using the crustal thinning factors and the depth of necking 118 . Both ends of the elastic beam are either matched to an analytical solution of a beam clamped at infinity or the modeling domain is extended laterally to minimize boundary effects. Although elastic thickness is believed to be close to zero during late rifting as the plate is thin and hot 12 , flexural isostasy was nonetheless tested (e.g. Model ELT10, Table 1 in main article). The results were similar to those assuming elastic thickness is 0 km and the effect of flexure is therefore regarded as insignificant during breakup.
Sediments are incorporated in the thermal solver to resolve blanketing effects [121][122][123][124] . Sediment deposition is controlled by the time dependent water depth which comes out of the inversion (see below). For each time step, a local sedimentation rate is iteratively computed so that at the end of a time step, i.e. after sedimentation, compaction, isostatic compensation and flexure, the basin is filled to the pre-computed water depth. The deposited sediments are compacted using empirical compaction laws, specifically those published by Refs. 125,126 that have the form: where φ is porosity, c is a compaction length scale and z the burial depth.
The boundary conditions for the thermal solver are user-specified fixed temperatures at the base and top of the numerical domain and zero horizontal heat flow at the sides utilizing a finite element method-based code. This allows the use of a body-fitting Lagrangian mesh strategy and ensures that all material contrasts are accurately resolved. The model and its individual components have been tested versus a number of analytical solutions including tectonic subsidence, thermal subsidence and steady state heat flow, with radioactive heat production rates decaying exponentially with depth.

Modifications and parameterization of melting and melt migration.
We have introduced new features to TecMod in order to test our main hypothesis that melt retention is an important source of buoyancy during and shortly after breakup of continental plates. This has required i) tracking of the lithosphere-asthenosphere boundary (LAB) and accounting for possible compositional changes across it, ii) resolving melt generation, iii) melt migration, and iv) melt emplacement.
Tracking of the LAB has been implemented by assigning distinct rock properties to the oceanic and continental lithospheres and asthenospheres (Table S1). The employed Lagrangian mesh allows accurate tracking of the LAB during active syn-rift and passive post-rift deformation.
Simulations of decompression melting and melt migration continue to be numerically challenging. Models explicitly designed to resolve these processes typically use published solidus functions to solve for melt generation and the equations of compaction driven fluid flow to track the produced melts during their ascent towards the surface. We have deliberately not taken this approach. Models and observations on melt generation, retention and migration are still difficult to reconcile with geophysically inferred melt fractions in various rift settings as they are often significantly higher than predicted by models. This and the considerable technical challenges associated with implementing melt migration into basin modeling software like TecMod has led us to derive a parameterization of melting and melt migration which allows us to investigate different scenarios that are consistent with geophysically inferred melt fractions.
We have parameterized decompression mantle melting by assuming that the degree of melting scales with the heights to which the LAB rises. The onset of melting is assumed to occur at 90 km depth, consistent with published solidus functions 65,69,94,98 . The concentration of melt is then assumed to increase linearly up to a specified depth of 20 km, where a maximum degree of melting of 16% is reached (average melt content of 8 %), which is consistent with observations at active rift zones and melt fractions calculated for adiabatic mantle upwelling 72 . Melt extraction is parameterized by assuming that the continental lithosphere is a permeability barrier for rising melts. As long as it is intact, melts are retained within the underlying asthenosphere. After breakup, loss of melt occurs exponentially at a prescribed rate of, for example, 15% every 1 My, meaning that more than half has been lost after 5 My and less than 2% melt remains in the asthenosphere after 10 My.
In addition, we have implemented the option to intrude igneous underplates. This allows us to explore scenarios in which a gabbroic sill of thickness equivalent to a specified fraction of the melt retained in the asthenosphere (e.g. 25%) is emplaced below the crust during every time step. For the Kwanza model reported in the accompanying article, this results in an underplate 5 km thick at the point of breakup tapering to 0 km towards the inboard edge of the rifted margin.
Finally, we use the models of Ref. 116 (specifically R123 peridotite) to test the effect that mineral phase changes in the mantle (both lithosphere and asthenosphere) have on the isostatic equilibrium of the plate.

Plate reconstructions
Recent plate reconstructions for the South Atlantic (Refs. 29,30 and Movie S4) use tectonostratigraphic ages that fit with the comprehensive well and seismic data acquired by the oil industry in Brazil and Angola 5,8,16 . The first oceanic crust in the central South Atlantic area is interpreted to have formed around 123 Ma (early Aptian), correlating with the breakup unconformity at the base of the pre-salt sag strata and corresponding to the age of the outer volcanic high which marks the inner edge of oceanic crust in Angola 5 . The sag strata  and the overlying salt  were deposited inboard and thinly on top of the high, which is interpreted to represent the embryonic ridge 18 . The favoured plate reconstructions show a relatively simple breakup (e.g. Movie S4). These reconstructions are significantly different to the models of Refs. 24,28 where large amounts of trans-continental shearing are required in order to fit plates which remain attached for another 10 My, until the end of salt deposition. This shear movement is speculative, we believe resulting from an erroneous breakup age. Figure S1. Depth of the syn-rift sediment surface at time of breakup, 123 Ma, for representative models in the accompanying article (Models Base, MA and Fdd15o, Fig.1a) plus additional models where the stratigraphic match is less than perfect.   10 as well as others noted in the same reference. Figure S4. Cross-sections from TecMod model Fdd15o, supplemental to those shown in Fig.1