Lithospheric mantle buoyancy: the role of tectonic convergence and mantle composition

Plate subduction and delamination, two key processes driving plate tectonics, are thought to be controlled by the buoyancy of the lithospheric mantle relative to the underlying asthenosphere. Most mantle delamination models consider a lithospheric density higher than the asthenosphere to ensure negative buoyancy (slab pull). However, mineral physics show that the continental lithospheric mantle density is lighter than the asthenosphere, and that only its pressure-temperature-composition dependence makes it become denser and unstable when sinking adiabatically. Here, we explore the controls on buoyancy using a 2D thermal-diffusive model of plate convergence, considering five chemical compositions and tectonothermal ages, namely Archon (>2.5 Ga), Proton (2.5–1.0 Ga), Tecton (<1.0 Ga), and two oceanic lithospheric plates of 30 Ma and 120 Ma. While the advection of colder rock in oceanic-like plates always results in negative buoyancy, Protons and Tectons exhibit an ability to slowly flip from negative to positive buoyancy at low convergence rates: they first favour the sinking due to advection and then become more buoyant because they are thinner and heat up faster during subduction. In contrast, the lighter density of cratons overprints this effect and hinders delamination or subduction, regardless of the convergence rate. This may explain why Archons are more stable during the Wilson cycle.

www.nature.com/scientificreports www.nature.com/scientificreports/ mobile belts (Tectons, <1 Ga) is closer to the composition of the Primitive Upper Mantle (PUM). The SCLM beneath Meso-and Paleo-proterozoic shields and mobile belts (Protons), with tectonothermal ages of 2.5-1.0 Ga, is typically intermediate in composition. These compositional variations affect the bulk density of the SCLM and the greater the degree of depletion, the lower the density 27 . In contrast, the composition of the oceanic lithospheric mantle corresponds to that of PUM after melt extraction at mid-ocean ridges (MOR), being relatively homogeneous except beneath large oceanic plateaus.
Petrological and geochemical studies show that at identical P-T conditions as at the LAB, the density of the SCLM is lower than that of the PUM 27 (see Supplementary Fig. S1), which is at odds with the aforementioned density distribution adopted in most geodynamic and static models 20 . In fact, the density of the sinking SCLM increases with pressure and decreases with temperature, thus P and T having competing effects on the depth-dependence of density. Whether mantle delamination or subduction are promoted by a negative buoyancy forced by plate convergence depends on the gradients of density relative to both parameters (Table 1).
Here we explore the idea that, whenever the lithosphere is incipiently forced to sink into the asthenosphere, it can become positively or negatively buoyant depending on its composition and on how the pressure-temperature evolution, imposed by plate convergence, affects its density evolution. This may explain why the older regions of the continental lithosphere (Archons) become tectonically more stable, self-prolonging their life span and favoring the Wilson cycle. To this purpose, we calculate the buoyancy of the sinking mantle in a kinematic model that accounts for a thermodynamically-consistent density dependence on temperature, pressure, and chemical composition.

Method
We developed a 2D kinematic approach to model plate convergence and subduction of the lithospheric mantle. The temperature distribution through time in absence of heat sources is given by where T is temperature (°C), t is time (s), α is thermal diffusivity (m 2 s −1 ), and → v is velocity (plate convergence rate, ms −1 ). The first term on the right side of Eq. (1) is related to thermal diffusion and the second term is related to thermal advection and Eq. (1) is solved with an explicit finite-difference scheme on a rectangular grid with a resolution of 5 × 5 km in a domain of 1500 × 600 km. Figure 2 shows the model setup and the related boundary  www.nature.com/scientificreports www.nature.com/scientificreports/ conditions. As we focus on the density changes of the sinking lithospheric mantle, the model excludes the crustal layer.
The initial density distribution is calculated by considering a depth-dependent density according to where the pressure ∂ρ ∂ ( ) P and temperature ∂ρ ∂ ( ) T derivatives for each lithosphere type are calculated by computing the stable mineral assemblages from major oxide compositions using Perple_X 28,29 and averaging the obtained values of the derivatives over the lithospheric mantle thickness ( Table 1). The considered lithospheric mantle types correspond to Arc_3 from the Slave craton, Pr_1 garnet-average, and Tc_1 garnet-average 27 , with their respective thicknesses to calculate representative density contrast at the LAB for Archean, Proterozoic, and Phanerozoic lithospheres, respectively ( Table 1). The Moho temperature in continental lithospheres varies between 550 °C and 750 °C (see Supplementary Fig. S1), and in our kinematic model we assume a mean value of 650 °C. For oceanic lithosphere we assume a mean Moho temperature of 300 °C (see Supplementary Fig. S1).
The down-going lithospheric mantle is submitted to a prescribed velocity field, which is directed downward with an angle of 30°, parallel to the subduction plane (Fig. 2). Only the temperature of the lithospheric mantle is recalculated using Eq. (1). The temperature in the asthenosphere is fixed to the initial adiabatic gradient (0.5 °C km −1 ). In this way, we ensure that the asthenosphere does not cool down due to the heat transfer to the lithosphere, mimicking the effects of sublithospheric convection (for simplicity, it is not explicitly implemented in our model). Density variations are recalculated for each moving particle according to Eq. (2).
The buoyancy force at each time-step is calculated as the integral of the density changes relative to the initial stage over the entire subducting slab, F b = ∫gΔρ dx dy, where g is gravity and Δρ is the density difference between the slab and the asthenosphere. The buoyancy force, measured in N m −1 , results from two components: diffusive (F d ) and advective (F a ), such that the total F b = F d + F a is directed downwards when F b < 0 (corresponding with slab pull).

Results and Discussion
As lithospheric convergence initiates, the downward advection of cold lithospheric material generates a negative temperature anomaly that extends along the lithospheric slab, from the axial plane (Fig. 2) to a growing depth below the undisturbed LAB, depending on the convergence velocity. Simultaneously, thermal diffusion is responsible for the reduction of that negative temperature anomaly by transferring heat from the asthenosphere to the subducting slab. Therefore, heat advection and thermal diffusion have opposite effects on temperature anomaly and buoyancy. The prevalence of each depends on the convergence rate, and the composition/age of the lithospheric mantle and its thickness. Figure 3a illustrates the contribution of advection and diffusion on the total buoyancy through time for an average Proton-type lithospheric mantle with a thickness of 110 km and a convergence rate of 40 mm yr −1 . In this case (Fig. 3a inset), the buoyancy results negative from the beginning of the convergence up to 4.0 Myr of evolution, when it reaches a minimum value of −1.1 TN m −1 (downwards directed force) increasing to neutral buoyancy at 9 Myr. From here on, the buoyancy increases to a maximum of +2 TN m −1 at 17 Myr and decreases again. These time variations of the total buoyancy are related to the temperature and density variations occurring along the down-going slab (Fig. 3b,c). www.nature.com/scientificreports www.nature.com/scientificreports/ Effect of convergence rate and mantle composition. The initial buoyancy of the lithospheric mantle is determined by the density contrast across the LAB Δρ = ρ − ρ ( ) LAB a sth(LAB) lith(LAB) , provided that P-and Tpartial derivatives do not differ much for different compositions. Regardless of the convergence rate, an average Archon-like composition (160 km LM thickness) always shows positive buoyancy (F b > 0) due to the large density contrast with the sublithospheric mantle at given P-T conditions with Δρ LAB = +68 kg m −3 (Fig. 4a, Table 1). For an average Proton-like composition (Δρ LAB = +39 kg m −3 ) with 110 km LM thickness, buoyancy initially has a negative trend (F b < 0) for all convergence rates, after which F b trends change and tending towards positive buoyancy after a time period spanning between 3 Myr for v = 80 mm yr −1 , and 40 Myr for v = 1 mm yr −1 (Fig. 4b). A younger average Tecton-like composition (Δρ LAB = +19 kg m −3 ) and 80 km LM thickness shows trends with negative buoyancy of higher absolute values than those of Proton's, F b < 0 always (Fig. 4c). In the case of oceanic settings, both young (30 Myr) and old (120 Myr) lithospheres have monotonously decreasing buoyancies (F b < 0) regardless of the convergence velocity (Fig. 4d), due to their intrinsically lower Moho temperature and lower density contrast across the LAB (Δρ LAB = +17 kg m −3 ; Table 1).
The changing trend of the total buoyancy observed in Proton-and Tecton-type lithospheres is related to the low convergence rate, which allows for the subducting slab to thermally re-equilibrate with the surrounding asthenosphere, favouring thermal diffusion against advection and causing the slab to become more buoyant. On the contrary, fast convergence rates prevent the down-going slab from thermal re-equilibration because advection prevails on thermal diffusion. F b becomes more negative (downwards slab pull) for higher convergence rates. The initiation of delamination or slab retreating should begin during this negative buoyancy stage, provided that the magnitude of the force and its duration suffice to trigger the process. The minimum (most negative) force for Protons occurs after 4 Myr (v = 40 mm yr −1 ) and 40 Myr (v = 1 mm yr −1 ) (Fig. 4b).
We explore the time needed to attain a certain value for F b that triggers delamination. We adopt a reference value of −3 TN m −1 which is a representative value for the plate tectonics driving forces [30][31][32] . Figure 5 displays the time needed to reach this value as a function of the density contrast across the LAB for two lithospheric thicknesses, 80 km and 160 km. The thicker lithospheric mantle takes less time to reach the reference pulling force (F b ) and over a wider range of density contrast (Δρ LAB ) than the thinner one. A faster convergence systematically leads to a wider compositional range for the triggering of delamination.
Finally, we isolate the effect of the lithospheric mantle thickness on the total buoyancy force, we vary this parameter while keeping a constant average Tecton-type composition with a convergence rate of 40 mm yr −1 . The results in Supplementary Fig. S2 show that a thicker lithospheric mantle results in a more negative buoyancy, reaching the −3 TN m −1 first. As the thickness of the lithospheric mantle decreases the tendency for the positive buoyancy increases. The temperature and density distributions are similar for the three cases (Tecton-type), therefore, the diffusion of the slab occurs with a similar rate, leading to little difference in the buoyancy forces (see Supplementary Fig. S2). The advective components vary more since this component is due to the difference in how much material is being advected down i.e. thicker lithospheric mantle has a more negatively buoyant advective component. www.nature.com/scientificreports www.nature.com/scientificreports/ tectonic relevance. The ability of the lithospheric mantle to delaminate or retreat is governed by its composition determining the density contrast relative to the underlying asthenosphere, and by the convergence rate. Our study shows that the oceanic lithosphere acquires negative buoyancy even at very low convergence rates and predicts ever-increasing slab pull during oceanic plate subduction, proportionally to the accumulated shortening and to the age of the plate. This is consistent with the observation that flat oceanic subduction is rare, requiring additional kinematic and rheological conditions to occur 33 . It is also consistent with the strong load (4 to 10 TN m −1 ) inferred from the flexural bending of the Pacific plate near the Tonga and Kermadec trenches, aged at 105 Ma 30 .
In contrast, the average Archon-type lithosphere always retains its positive buoyancy, whereas intermediate mantle compositions, corresponding to average Proton-type, can attain negative buoyancy depending on the convergence rate and the elapsed time. For convergence rates <70 mm yr −1 , the maximum amplitude of negative buoyancy is attained after 80-250 km shortening (Fig. 4). This is similar to results from models of delamination of  www.nature.com/scientificreports www.nature.com/scientificreports/ a lithospheric root during orogeny 16 , where it was concluded that lithospheric roots (equivalent to our 'subducted slab' portion) of at least 100-170 km length are needed to generate sufficient negative buoyancy for delamination and detachment to proceed.
The presented results are based on average chemical compositions representative of five lithospheric mantle types, and based on the assumption that convergence rates and compositions are constant through time. Mantle compositions and densities actually vary over a wide range even within given age and P-T conditions. The positive buoyancy of Archon's is not warranted for lherzolitic mantle compositions, as in Kaapvaal (South Africa) and Almklovdalen (western Norway) cratons, showing densities 30 kg m −3 larger than average Archon-type and close to average garnet Proton-type 27 . Similarly, Proton's with harzburgite compositions can have densities 32 kg m −3 lower than the average Proton-type and will never get negative buoyancy. Moreover, metasomatic processes can change the major-element mantle composition converting harzburgite/dunite (lighter end-member) to lherzolite (heavier end-member) and back, modifying its buoyancy 34 . Convergence rates may also vary through time, especially for slow speeds (v ≤ 20 mm yr −1 ) since they involve long lasting processes  to attain significant shortening. These buoyancy variations, caused either by lateral compositional heterogeneities, metasomatic and differentiation processes, far-field tectonic stresses, or a combination of these, can contribute to isostatic destabilization of cratons by dripping and detachment of heavier mantle portions, as can be the case of western Gondwana craton 35 .
An outstanding outcome of the model is the predicted change from negative to positive buoyancy through time for continental subduction. The initial density contrast across the LAB determines how soon diffusion overcomes advection such that intermediate density contrasts 20 kg m −3 ≤ Δρ LAB ≤ 50 kg m −3 and convergence velocities ≤50 mm yr −1 eventually lead to a change from negative to positive buoyancy (Fig. 4). This may promote a rising of the subducting continental lithosphere and a subsequent flattening below the overriding plate producing lithosphere underthrusting as proposed in the India-Eurasia collision region ( Fig. 1 and Fig. 6). According to seismic tomography and potential field modelling 36,37 , underthrusting beneath the western Himalaya-Tibet extends 200 km beyond the suture, whereas to the east, the convergence is accommodated by a steep subducting slab. This two-fold behaviour is compatible with a 160-km-thick lithospheric mantle of Proton-type and subducting at a high velocity of 45 mm yr −1 (Fig. 5). Clearly, our model disregards other relevant mechanisms such as the weight of oceanic lithospheric slabs subducted before collision or the role of subduction of the continental crust, and hence quantitative comparisons with real scenarios must be taken with caution 19,20,24,25,38,39 . Figure 5 also shows three geological scenarios (Alboran 40 , Zagros 41 , and Tibet 41 , located in Fig. 1) that can be quantitatively compared to our study. The Alboran Basin is a back-arc setting characterised by a thinner lithospheric mantle 42 (80 km) undergoing Cenozoic mantle subduction, and possibly late Miocene delamination 43 , whereas Zagros and Tibet are in the thicker lithospheric mantle group 37 (red lines from Fig. 5, 160 km). The three regions are clustered on the left side of the plot (low Δρ LAB ). Our model suggests that in order to develop negative buoyancy under higher Δρ LAB as in cratons, the convergence rate should be higher than 80 mm yr −1 or remain forced for tens of million years, making it a rare phenomenon (Fig. 6).
In summary, the results show that incorporating realistic mineralogy-based densities to geodynamic models rises up an unforeseen control on the development of negative buoyancy. Whereas assigning a constant or temperature-dependent higher density to the lithosphere always results in slab pull, accounting for the effects of composition and pressure reveals that the plate convergence velocity is key to determine the development of negative buoyancy, delamination, and subduction.
The model designed here provides a methodological framework for understanding the stability of the lithosphere during the convergence of tectonic plates, and suggests a simple thermodynamic explanation for the long-term preservation of older continental regions (cratons) in the Wilson cycle. As continents aggregated during the early Earth evolution, their average buoyancy relative to the asthenosphere increased, making them less prone to subduct or delaminate and hence more stable. However, even Archean and Proterozoic lithospheric plates retain a chance to become recycled into the mantle if they are forced to sink by a fast-enough plate convergence, depending on the geographical configuration of tectonic plates. Figure 6. Cartoon summarizing the effect of the convergence rate. (a) shows the case with a slow convergence rate which allows for great thermal diffusion, which in turn, reducing the slab's density and the slab becomes more buoyant. (b) shows the case with a faster convergence rate, which allows the slab to keep cooler and denser, leading to a higher density and an increase in negative buoyancy.