Introduction

In contrast to oceanic lithosphere which gets recycled into the mantle at subduction zones, continental lithosphere appears to be quite old, in particular in cratons where rocks more than a billion years are preserved. Such old lithosphere would be expected to be cold and dense. However, if, based on expected density structure, isostatic topography is computed and compared with actual topography, it is found that continents have higher topography than expected. This has led to the isopycnic hypothesis of Jordan (1978, 1988)—that the negative thermal buoyancy of cold continental lithosphere is at least largely, if not completely, compensated by positive chemical buoyancy. More recent work has shown that this compensation is not complete, and there is still a remaining negative buoyancy (Forte and Perry 2000; Wang et al. 2022). Liu et al. (2025) provide a detailed discussion on this topic. Also, there are sub-lithospheric contributions to topography due to the convecting mantle—so-called dynamic topography, which have to be taken into account.

To address the question of total density anomalies in the lithosphere and at what depth these occur we take here an approach similar to Wang et al. (2023) but with different modeling assumptions and covering a wider range of models: We plot the difference between residual and dynamic topography against lithosphere thickness. The slope of this line at different lithosphere thickness values should then tell something about lithosphere density anomalies at those depths: Residual topography is simply computed by subtracting crustal isostatic topography, using a certain crustal model, from observed topography. It therefore contains the effect of density anomalies within the lithosphere. Dynamic topography, on the other hand, is computed here with a dynamic flow calculation based on density anomalies that are inferred from seismic tomography, using a thermal scaling, below the lithosphere only, i.e., density anomalies are set to zero within the continental lithosphere.

At least in theory, the difference, residual minus dynamic topography, should then be the result of density anomalies within the lithosphere only. It also should not matter whether any age-dependent topography in the oceans has been subtracted, as long as it is done exactly the same way for both dynamic and residual topography. Since lithosphere thickness, as well as dynamic and residual topography are all in units of kilometers, the slope is dimensionless and can be expressed in units of per cent. Since we compute (or convert to) all dynamic and residual topography “below air” (with negligible density), it can be interpreted as per cent density anomaly, but one needs to be clear about at what depth and relative to what.

Seismic anomalies are relative to a global mean zero. Hence, if we set density anomalies to zero within the lithosphere, for the computation of dynamic topography, their global mean is no longer zero. So, relative to the new global mean, the density anomalies in continental lithosphere are not zero. One can iteratively correct for that, and we will present results for both cases where this correction has been done, and where it has not. It has also been proposed that mid-ocean ridges, with essentially adiabatic temperatures all the way to the surface are a suitable reference: Relative to this much hotter reference, the continents will obviously be much more negatively buoyant than relative to the global mean chosen here, and the fraction of thermal anomalies that are compensated by chemical anomalies will be much lower in this case. We will also present results where continental lithosphere densities are set to sub-ridge values.

Dynamic topography is computed in a way such that the global mean is zero by definition (as is for residual topography), even if the density mean at lithospheric depths is not zero. Whether or not this should be corrected for does not matter for the purpose of this paper, where we are only concerned about the slope of residual minus dynamic topography versus lithosphere thickness, and any constant shift does not matter.

Methods

To determine to which density anomalies the slope corresponds to, we consider an isostatic balance between topography \(h\) with mantle surface density \(\rho_0\) (as crustal isostasy has been removed) and lithospheric density anomalies \(\delta \rho(z,t_l)\) (i.e., we assume those density anomalies depend on depth and lithosphere thickness): \[\begin{equation} \rho_0 h = \int\limits_0^{t_l} \delta \rho(z,t_l) dz. \label{eq:isostatic1} \end{equation}\] The lithospheric density anomalies consist of thermal and compositional anomalies. If we further assume that the chemical anomalies \(\rho_c\) only depend on depth (not on lithosphere thickness) and the continental lithosphere geotherm follows a complementary error function, equation (\ref{eq:isostatic1}) becomes \[\begin{equation} h = \int\limits_0^{t_l} \delta \rho_c(z)/\rho_0 dz + \delta\rho_0/\rho_0\int\limits_0^{t_l} \mbox{erfc}(z/t_l) dz \label{eq:isostatic2} \end{equation}\] Here \(\delta\rho_0\) is the thermal density anomaly that corresponds to the difference between mantle potential temperature and actual surface temperature, i.e., the reference temperature for thermal density anomaly in equation (\ref{eq:isostatic2}) is the adiabatic temperature profile of mid-ocean ridges. The second integral is simply the product of lithosphere thickness and the integral of the complementary error function from zero to 1, which is \(C=1-\mbox{erf}(1)+(1-e^{-1})/\sqrt\pi=0.514\). According to McKenzie et al. (2005) cratonic mantle temperature instead follows nearly a straight line. In this case, \(C\) would be one half and results would remain essentially the same. The derivative of topography with respect to lithosphere thickness, which corresponds to the slope, then becomes \[\begin{equation} \frac{dh}{dt_l}=\frac{\delta\rho_c(t_l)}{\rho_0}+C\frac{\delta\rho_0}{\rho_0} \label{eq:isostatic3} \end{equation}\] Hence the first term is the chemical density anomaly at the base of the lithosphere, but the second term is not the thermal density anomaly at the base of the lithosphere, but the average thermal density anomaly over the entire lithosphere, which is approximately half its value at the surface. With a thermal expansivity \(3.3\times 10^{-5}\)/K, and a temperature contrast across the lithosphere of \(1325\) K, this means that a purely thermal lithosphere would correspond to a slope of \(-2.25\%\). This is relative to an adiabatic mantle, i.e., if we take the global ridge system as reference. The corresponding slope is less with the global mean density or seismic velocity as reference, as this mean already corresponds to a depth-dependent temperature difference to the adiabat which would have to be subtracted. The interpretation of equation (\ref{[eq:isostatic3}) depends on the assumption that continents with different lithospheric thicknesses share the same chemical density anomaly at the same depth. More generally, if this is not the case, a steeper slope for thicker lithosphere just means that thicker lithosphere is on average less buoyant, but our approach cannot distinguish at what depth.

For lithosphere thickness models, we follow Steinberger and Becker (2018), i.e., we test 5 different tomography models: (1) sl2013 (Schaeffer and Lebedev 2013), (2) semum2 (French et al. 2013), (3) s40rts (Ritsema et al. 2011), (4) savani (Auer et al. 2014), (5) gypsum (Simmons et al. 2010).

Dynamic topography is computed from sub-lithospheric density anomalies derived from the same models. However, since sl2013 and semum2 do not cover the whole mantle, they are merged with Grand10—the 2010 model update of Grand (2002)—in the deeper mantle, with a smooth transition between 150 km and 250 km. All tomography models are interpolated on the same radial layers.

Dynamic topography computation follows Steinberger (2016), with some differences: In most regions we convert relative seismic velocity anomalies to relative density anomalies with a factor that is somewhat depth-dependent and around 0.22 in the upper 400 km, assuming the same conversion coefficients to density for both slow and fast seismic anomalies so that density anomalies follow the same spatial distribution as that of seismic anomalies. However, in a first step we set anomalies within the crust to zero, in line with removing crustal isostasy in residual topography, i.e., we aim at removing effects of the crust as much as possible. In a second step, density anomalies within the continental lithosphere are also modified, using three options: In the first option, they are also set to zero (i.e., the global mean). However, after this replacement the mean is no longer zero (except in our uppermost layer at depth 7 km, which is all within the crust and where hence density anomalies have been set to zero globally), and the lithosphere anomalies no longer corresponds to this new mean. Hence, in our second approach, we replace them by a depth-dependent different value that corresponds to the new mean. This value is found by an iterative trial and error procedure. We first compute global average density at each model depth after replacing densities in the continental lithosphere by zero. Then we replace it by something near that average and re-compute the average at each depth, and repeat that computation until the procedure converges. In our third approach, we replace anomalies in the continental lithosphere by the approximate mean value below mid-ocean ridges. We compute that based on the Müller et al. (2008) age grid version 3.6, available online at https://www.earthbyte.org/webdav/ftp/earthbyte/agegrid/2008/Grids/ as the mean value below all sea floor with age less than or equal to 5 Ma.

Also, in all cases the same radial viscosity model, corresponding to the reference case of Steinberger (2016), is used. A schematic illustration of the calculation procedure is given in Figure 1. Apart from the crust and continental lithosphere, we keep density anomalies from tomography everywhere and not just up to a specific depth (e.g., 250 km) as has sometimes been done before. In this way, dynamic topography can match some smaller-scale features of residual topography, which would not be the case otherwise (Steinberger 2016). Also note that tomography is relative to global average, i.e., including an (average) thermal boundary layer. Relative to an adiabat, the lithosphere is always cold and thermally dense, contrary to the appearance in tomography. However, whether or not a global average thermal boundary layer is included or not does not affect the dynamic topography computation where we always assign zero mean by definition.

Following Steinberger (2016) we also removed the topography signal related to sea floor age, and an ocean-continent contrast by assigning continents an appropriate “equivalent” sea floor age, because the same procedure is also applied to residual topography. For the purpose of our paper, where we are only concerned with continents, this almost does not matter. It only has a slight effect because, due to the spherical harmonic expansion we use, the seafloor age signal slightly “leaks” into the continents. Since the global mean remains at zero, it also shifts continental elevations by a constant value, but since we are only concerned about the slopes, this does not matter either.

Residual topography model 1 was adopted from Steinberger (2016), model 2 comes from a combination of Holdt et al. (2022) in the oceans and Stephenson et al. (2024) on continents. Oceanic data are divided by 3330/(3330–1030)=1.447826, hence converted from “below water” to “below air” to correspond to dynamic topography computations. Also, the way this model combination is constructed corresponds to a different “equivalent” sea floor age than in our dynamic topography computations. Again, since here where we are only concerned with continents, this construction procedure almost does not matter and this is essentially the model of Stephenson et al. (2024), again with slight changes due to spherical harmonic expansion and a possible constant shift in elevation. Model 3 was computed in the same way as in Steinberger (2016) but with a modification of Crust 1.0 (Laske et al. 2013) crustal thickness (T.W. Becker, pers. comm., further developed from Faccenna (2020)).

All models are evaluated on the same latitude-longitude grid with 1 degree spacing in latitude. Longitudinally the grid spacing is also about one degree of arc, i.e., the points are 1/cos(latitude) degrees of longitude apart. In this way, each grid point represents about the same surface area which simplifies the best-fitting line computations, but, more importantly, leads to an improved visual representation, where areas closer to the poles are not over-represented. Each point is also characterized as either 1-oceanic lithosphere or 2-weak plate boundaries or 3-cratons or 4-other continental lithosphere (Cui et al. 2024). Here we only consider cratons and other continental lithosphere, but we include also those plate boundary regions that are on continents (using the -200 m isobath to discriminate, and explicitly excluding Iceland) in the latter. This gives 6966 points on cratons, and 15522 points for the entire continental lithosphere.

Combining five dynamic topography models with three residual topography models gives 15 different combinations for their difference. Among these, we choose the combination of savani dynamic topography, and residual topography model 3 as reference case, as this is the best-fitting combination of Cui et al. (2024). Plotting the difference against lithosphere thickness gives a cloud of points through which we fit either one or two straight lines. In order to not give outliers too much weight, we use the sum of the absolute differences between points and line as misfit which we aim to minimize, using Numerical Recipes (Press et al. 1986) routines amoeba.f and amebsa.f. The one-line fit is determined with the downhill simplex method amoeba.f. To increase chances for achieving a good result for the two-line fit, we initiate the routine with what gives visually a good fit. We then apply amoeba.f which will converge to a local and not necessarily the global minimum. To further improve the results and if possible converge to the global minimum, we subsequently apply the simulated annealing method amebsa.f. We will use the two-line fit if it corresponds to at least 1% improvement compared to the one-line fit, which we will use otherwise.

To reduce spread from the cloud of points, we also compute averages of topographies or their differences and lithosphere thicknesses for points with similar lithosphere thickness. For this, we first sort the points according to lithosphere thickness, using the Numerical Recipes routine hpsort.f. We then average over each 777 consecutive points (slightly less for the final set), yielding 20 averaged points for the entire continents and nine averaged points for the cratons.

Results

In order to be able to assess what difference it makes to set density anomalies in the continental lithosphere to either zero, or to a value such that it is equal to the global mean (after this replacement) for each depth, we first determine iteratively that value. Figure 2 left shows results for our five tomography models (or combinations thereof), with seismic anomalies converted to density anomalies. Given that seismic anomalies in the lithosphere are mostly positive, replacing them by zero would shift the global mean to negative values. If they are to remain equal to the global mean, they need to be (even more) negative, and this is what is shown. At crustal depths, offsets can be both positive and negative, as they are also influenced by negative seismic velocity anomalies within the crust. Figure 2 right shows the average density anomalies below ridges, which are also negative but larger magnitude. If seismic velocity anomalies were only due to temperature anomalies, this curve should approximately correspond to the difference between the adiabatic temperature at a given depth and the global average temperature at that same depth, multiplied with thermal expansivity, i.e., it should continuously increase towards the surface, reaching a large value of more than 4% there. However, the actual profiles only increase up until a certain depth, which also differs among models, and then decreases again. We can see at least two reasons for the discrepancy between what we would expect for thermal anomalies only and what we actually see. The first one is vertical smearing, which also gives rise to differences among models. The second one is that seismic velocities are also influenced by chemical differences—most importantly, velocities are lower in the crust, but there are also the effects of chemical differences in the sub-crustal lithosphere. However, we don’t expect that chemical differences play a major role beneath ridges, which makes them an appropriate reference column, if we want to determine densities of continental lithosphere relative to an adiabatic mantle, and as long as the profiles are not affected by vertical smearing.

In Figure 3 we use our reference case as an example to illustrate our procedure. The top right panel shows the lithosphere thickness on continents determined in Steinberger and Becker (2018) for tomography model savani (Auer et al. 2014). It shows cratons typically having the thickest lithosphere up to around 250 km thick. Another region of very thick lithosphere is around Tibet. Much of the remaining continents have lithosphere thickness in excess of 100 km, but there are also continental regions with less than 50 km thick lithosphere, such as in western North America, East Asia, and a corridor stretching north from Ethiopia to eastern Turkey. Lithosphere thickness on continents for all models in Figure 6 show broadly similar patterns. See Steinberger and Becker (2018) for the global pattern, including oceans. The top center panel of Figure 3 shows that many regions with thin lithosphere are also characterized by high residual topography, including East Asia, western North America and the Afar region. An obvious reason for this relation could be that in regions of thin lithosphere, the hot asthenosphere comes closer to the surface, giving rise to high topography. In other continental regions there is no immediately obvious relation between the two maps; for example there are cratons with both high and low residual topography. However, if we plot residual topography pointwise against lithosphere thickness (top left panel) we can see a pattern: For lithosphere thickness above around 150 km, thicker lithosphere tends to be associated with lower residual topography, which could indicate that lithosphere below a depth of about 150 km is negatively buoyant—something we will further investigate in the rest of the paper. However, there is also quite a spread of about 2–3 km for residual topography with the same lithosphere thickness. When comparing individual cratons, indicated by different colors, we can for example see that African cratons tend to have higher residual topography for the same lithosphere thickness than other cratons—something that can be readily attributed to Africa sitting on a topographic swell above a Large Low Shear Velocity Province (LLSVP). Brown dots indicate that cratons, especially thinner ones, on average have higher residual topography than other continents, which is also due to African cratons. The other models of residual topography, and plots of residual topography versus lithosphere thickness for all 15 combinations are shown in Figure 9, often showing a similar pattern—sometimes clearer, sometimes less clear. Notably, residual topography model 2, largely based on Stephenson et al. (2024) tends to show less variation than the other two, which are based on Crust 1.0 (Laske et al. 2013). In case their crustal thickness is unmodified (model 1) resulting residual topography is strongly negative in some continental regions, in particular Eurasia. The negative topography is less pronounced in our model 3, where crustal thicknesses have been modified (courtesy of Thorsten Becker), except in East Antarctica, where it is now strongly negative instead of positive.

In contrast, if we plot dynamic topography in our reference case in Figure 3, we see a different pattern, which also depends on what density anomalies we assign to the continental lithosphere: If we assign zero densities (i.e., the global mean) there does not appear much of a relation between dynamic topography and lithosphere thickness at all, except that, again, very high topography appears to be associated with very thin lithosphere, as discussed above. Also, African Cratons, and cratons on average, especially thinner ones, tend to again be associated with higher topography, which we also already discussed. In contrast, if continental lithosphere is assigned the same density as mid-ocean ridges, the computed dynamic topography becomes higher, and it also tends to increase with lithosphere thickness. Results for all models shown in Figures 6 and 8, whereby dynamic topography is only plotted against lithosphere thickness derived from the same tomography model, show overall similar patterns. Furthermore, comparison with Figure 7 shows that results are not much affected by modifying densities in the continental lithosphere such that they remain equal to the global mean.

The difference of residual and dynamic topography, which is supposed to be caused by density anomalies in the continental lithosphere again depends on whether we use the global mean (middle right panel of Figure 3 for our reference case) or mid-ocean ridges (lower right panel) as baseline relative to which those anomalies are defined. In the former case, the areal extent of regions with positive and negative difference is similar, whereby regions of thick lithosphere and cratons in particular are more often associated with a negative difference, again suggesting a negatively buoyant lower lithosphere. In the latter case, there are more regions with negative difference, and regions of thick lithosphere and cratons are often associated with a very strongly negative difference, indicating that, relative to an adiabatic mantle, many regions of continental lithosphere, in particular thick lithosphere, are negatively buoyant. Result for all 5 dynamic topographies, three residual topographies and three different “baselines” for continental lithosphere densities (i.e., a total of 45 combinations) are shown in Figures 10, 11 and 12, and summarized in Tables [tab:table1] and [tab:table2].

As already for residual topography, we find that results become clearer if we plot topography difference against lithosphere thickness. In the case with zero (i.e., global average) densities in the continental lithosphere (Figure 4) in the reference case (combination of savani and residual topography 3) the difference appears more or less independent of lithosphere thickness up until around 100–150 km thickness, but above that, and in particular for cratons, it tends to becomes more negative with thickness. This pattern was also found for residual topography, and can therefore be attributed to positive density anomalies in the lithosphere—in particular below 100–150 km depth—which affect residual but not dynamic topography. The large spread of results of around 2 km (\(\pm 1\) km from the best-fit line) is probably due to uncertainties in both dynamic and residual topography models. They become also apparent when comparing the results for different dynamic and residual topography models. The spread is smallest with residual topography 2, due to this model having the lowest amplitude. Despite the spread, we find the overall pattern of topography difference being either independent of, or even slightly increasing with lithosphere thickness, if the thickness is smaller than a certain value, and becoming more strongly negative with increasing thickness for thicknesses larger than that. The corresponding break in slope of a best-fitting line tends to occur around thickness 150 km. Figure 13 shows that results are not much affected if the continental lithosphere density is adjusted such that it remains equal to the global mean after this adjustment.

However, Figure 5 shows that results are strongly affected, if instead continental lithosphere density is set equal to mid-ocean ridge density: In this case, the slopes tend to be more negative and differences become strongly negative for thick lithosphere and cratons in particular, showing that relative to an adiabatic mantle, the lithosphere has excess density, i.e., is negatively buoyant overall. Also, although it is less apparent than in Figure 4, there still tends to be a break in slope at around 150 km in some cases, i.e., the slope of the best-fitting line tends to be even more strongly negative for higher thickness.

Discussion

Our main results are based on plots of residual minus dynamic topography versus lithosphere thickness, but all of these quantities have considerable uncertainties. To at least partly account for uncertainties, we use various models: We use five different lithosphere thickness models, but all of them are derived from seismic tomography, using the same method. There are other methods as well, based on heat flow (Artemieva 2006), elastic thickness (Audet and Bürgmann 2011) and receiver functions (Rychert et al. 2010). A comparison of lithosphere thickness models is given by Steinberger and Becker (2018). It is not necessarily expected that these different models for lithosphere thickness agree, as material may lose its elastic strength at a different (lower) temperature as its rheological strength, and receiver functions could indicate a boundary within the lithosphere rather than the base of the lithosphere (Selway et al. 2015; Rader et al. 2015). In particular Yuan and Romanowicz (2010) propose a boundary between a top chemical layer and a bottom thermal layer at a depth of about 150 km in regions of thick lithosphere, which agrees well with the break in slope we find around that depth. Regarding determining total lithosphere thickness, our tomography-based method should work well as long as seismic velocities below (not necessarily within) the lithosphere are related to temperature, and viscosity is a function of temperature.

Uncertainties in residual topography are caused by uncertainties in crustal density and thickness. For example, if the assumed crustal density is wrong by only 1%, a simple calculation of isostasy yields several hundred meters difference in residual topography. These uncertainties illustrate what is probably the main reason that our plots of topography difference versus lithosphere thickness are rather clouds of points, with a spread of about \(\pm\) 1 km around the best-fitting lines we determine. Nevertheless, both visually and computationally this “break in slope” of the cloud of points can be found for most of the model computations investigated. That is why we regard our results as robust.

Left of the break in slope, which occurs around 150 km, the slope tends to be close to zero, if lithosphere density anomalies are set to the global mean. Hence depth-integrated continental lithospheric density anomaly is close to the global mean for lithosphere of less than around 150 km thickness. Since this global mean includes both young and old seafloor, we expect that this is close to mid-aged seafloor, however, a more definite statement would require further analysis. Since mid-age ocean lithosphere is less than 150 km thick, this means continental lithosphere is chemically buoyant. This is consistent with Artemieva (2009) who finds that velocity anomalies of non-thermal origin contribute significantly to seismic velocity anomalies in the lithosphere, and discusses the causes of these non-thermal anomalies.

Another problem with our approach is that, although it should work in theory, there are well-known discrepancies between dynamic and residual topography models that are probably unrelated to the lithosphere. In particular, dynamic topography tends to be larger than residual topography, with the largest discrepancy in amplitude at spherical harmonic degree two (Steinberger et al. 2019). Richards et al. (2023) suggested a solution for this discrepancy due to chemical anomalies in the lowermost mantle that imply that one cannot simply convert from seismic to density anomalies—something which we do not consider here. It can, however, be hoped that, since the predominant wavelength of dynamic topography from the lower(most) mantle is much larger than craton size, that this discrepancy mainly causes a spread of the data points in the plots of residual minus dynamic topography versus lithosphere thickness, and does not so much affect the slope.

Our results are also in broad agreement with compositional stratification inferred by Steinberger and Becker (2018). In their Figure 6, they plot the apparent temperature (inferred from tomography) versus depth for different ranges of lithosphere thickness. Deviations of this curve from what is expected (i.e., temperatures getting lower and temperature gradient increasing for shallower depth) may indicate compositional effects: For lithosphere of more than 200 km thickness, these deviations become prominent above about 150 km, for thinner lithosphere the deviations occur above a somewhat shallower depth. Similarly, Lekic and Romanowicz (2011) find from a cluster analysis of upper mantle tomography that shear wave velocity versus depth profiles of stable platforms and shields consistently exhibit a mid-lithospheric low velocity zone (LVZ) between 80 and 130 km depth, while the asthenosphere is found at depths greater than 250 km in both regions. A layering is also indicated by anisotropic structure with mostly \(V_{sh} >V_{sv}\) at depth 70 km overlying \(V_{sh} < V_{sv}\) at depth 150 km especially in cratons, as shown by Wang et al. (2023) based on data by Ho et al. (2016).

Liu et al. (2025) argue that cratonic roots are notably denser than the ambient mantle. This result comes from using mid-oceanic ridges, i.e. with adiabatic temperatures all the way to the surface, as reference. This contrast to the isopycnic hypothesis of Jordan (1978, 1988) may be due to using a different reference column: For example, in our previous work (Steinberger 2016) we use the global average, and relative to that, we found that continental lithosphere is denser, although a lot less than it was purely thermal. Similarly, Forte and Perry (2000) found that chemical buoyancy compensates for 80% of (but not all of) thermal buoyancy. The percentage of compensation will be a lot less if mid-ocean ridges instead of global mean is used as reference. Here we address this apparent controversy by computing and comparing results with both reference columns, and a third one where densities in the continental lithosphere are adjusted such that they stay equal to the global mean after adjustment.

Our results imply that positive total density anomalies are prevalent in the lower, and hence hotter and less viscous parts of thick lithosphere. This may, especially in combination with other causes lead to lithospheric dripping of parts with higher density (Hua et al. 2025). Such delaminations happening through Earth history should lead to events of uplift of the overlying lithosphere, which should be detectable in the geologic record (Wang et al. 2023; Liu et al. 2025).

Conclusions

Here we explored a method for determining continental lithosphere densities from the difference between computed dynamic topography and observation-based residual topography. If only density anomalies below the lithosphere are included in the computation of dynamic topography, this difference should indicate density anomalies within the lithosphere. Our results depend on which density structure is chosen for reference. If we choose the global mean, we find close to zero difference for continental lithosphere thickness less than about 150 km, meaning that above that depth, any thermal anomalies are nearly compensated by chemical anomalies. For more than 150 km thickness, there is a difference that increases with thickness, implying that below that depth, density anomalies are either purely thermal or at least negative chemical density anomalies are a lot less pronounced. If, on the other hand, we choose the presumably nearly adiabatic temperatures below mid-ocean ridges for reference we find differences increasing with lithosphere thickness in the entire thickness range, but the increase becomes stronger below approximately 150 km depth. Hence, relative to adiabatic temperatures, lithosphere is negatively buoyant at all depth, but at least above 150 km, this is partially compensated by positive chemical buoyancy. Our results have implications for lithosphere stability, in particular they may be compatible with instability and dripping of a lower lithosphere layer that is denser than the asthenosphere at the same depth but hotter and less viscous than the lithosphere above.

Author contributions

B.S. conceived and wrote the paper. Both authors conducted computations and contributed to discussions.

Data availability

All data required to produce the figures are available at https://doi.org/10.5281/zenodo.17805343.

Competing interests

The authors declare that they have no competing interests.

Supplementary Materials

Supplementary materials contain figures that are supplementary in the sense that they illustrate intermediate steps in the computation and results for all model combinations, whereas in the main text only exemplary figures conveying the main conclusions are shown: In particular, here maps of lithosphere thickness and dynamic topography, and plots of lithosphere thickness vs. dynamic topography are shown for all combinations of five tomography models considered and three assumptions of lithosphere density anomalies. Residual topography and plots of lithosphere thickness vs. residual topography are shown for all possible combinations. Likewise, maps of residual minus dynamic topography are shown for all 45 combinations, and plots of topography difference vs. lithosphere thickness for the 15 combinations not yet shown in the main paper.

References